View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math3.dfp;
19  
20  import org.apache.commons.math3.Field;
21  import org.apache.commons.math3.FieldElement;
22  
23  /** Field for Decimal floating point instances.
24   * @since 2.2
25   */
26  public class DfpField implements Field<Dfp> {
27  
28      /** Enumerate for rounding modes. */
29      public enum RoundingMode {
30  
31          /** Rounds toward zero (truncation). */
32          ROUND_DOWN,
33  
34          /** Rounds away from zero if discarded digit is non-zero. */
35          ROUND_UP,
36  
37          /** Rounds towards nearest unless both are equidistant in which case it rounds away from zero. */
38          ROUND_HALF_UP,
39  
40          /** Rounds towards nearest unless both are equidistant in which case it rounds toward zero. */
41          ROUND_HALF_DOWN,
42  
43          /** Rounds towards nearest unless both are equidistant in which case it rounds toward the even neighbor.
44           * This is the default as  specified by IEEE 854-1987
45           */
46          ROUND_HALF_EVEN,
47  
48          /** Rounds towards nearest unless both are equidistant in which case it rounds toward the odd neighbor.  */
49          ROUND_HALF_ODD,
50  
51          /** Rounds towards positive infinity. */
52          ROUND_CEIL,
53  
54          /** Rounds towards negative infinity. */
55          ROUND_FLOOR;
56  
57      }
58  
59      /** IEEE 854-1987 flag for invalid operation. */
60      public static final int FLAG_INVALID   =  1;
61  
62      /** IEEE 854-1987 flag for division by zero. */
63      public static final int FLAG_DIV_ZERO  =  2;
64  
65      /** IEEE 854-1987 flag for overflow. */
66      public static final int FLAG_OVERFLOW  =  4;
67  
68      /** IEEE 854-1987 flag for underflow. */
69      public static final int FLAG_UNDERFLOW =  8;
70  
71      /** IEEE 854-1987 flag for inexact result. */
72      public static final int FLAG_INEXACT   = 16;
73  
74      /** High precision string representation of &radic;2. */
75      private static String sqr2String;
76  
77      // Note: the static strings are set up (once) by the ctor and @GuardedBy("DfpField.class")
78  
79      /** High precision string representation of &radic;2 / 2. */
80      private static String sqr2ReciprocalString;
81  
82      /** High precision string representation of &radic;3. */
83      private static String sqr3String;
84  
85      /** High precision string representation of &radic;3 / 3. */
86      private static String sqr3ReciprocalString;
87  
88      /** High precision string representation of &pi;. */
89      private static String piString;
90  
91      /** High precision string representation of e. */
92      private static String eString;
93  
94      /** High precision string representation of ln(2). */
95      private static String ln2String;
96  
97      /** High precision string representation of ln(5). */
98      private static String ln5String;
99  
100     /** High precision string representation of ln(10). */
101     private static String ln10String;
102 
103     /** The number of radix digits.
104      * Note these depend on the radix which is 10000 digits,
105      * so each one is equivalent to 4 decimal digits.
106      */
107     private final int radixDigits;
108 
109     /** A {@link Dfp} with value 0. */
110     private final Dfp zero;
111 
112     /** A {@link Dfp} with value 1. */
113     private final Dfp one;
114 
115     /** A {@link Dfp} with value 2. */
116     private final Dfp two;
117 
118     /** A {@link Dfp} with value &radic;2. */
119     private final Dfp sqr2;
120 
121     /** A two elements {@link Dfp} array with value &radic;2 split in two pieces. */
122     private final Dfp[] sqr2Split;
123 
124     /** A {@link Dfp} with value &radic;2 / 2. */
125     private final Dfp sqr2Reciprocal;
126 
127     /** A {@link Dfp} with value &radic;3. */
128     private final Dfp sqr3;
129 
130     /** A {@link Dfp} with value &radic;3 / 3. */
131     private final Dfp sqr3Reciprocal;
132 
133     /** A {@link Dfp} with value &pi;. */
134     private final Dfp pi;
135 
136     /** A two elements {@link Dfp} array with value &pi; split in two pieces. */
137     private final Dfp[] piSplit;
138 
139     /** A {@link Dfp} with value e. */
140     private final Dfp e;
141 
142     /** A two elements {@link Dfp} array with value e split in two pieces. */
143     private final Dfp[] eSplit;
144 
145     /** A {@link Dfp} with value ln(2). */
146     private final Dfp ln2;
147 
148     /** A two elements {@link Dfp} array with value ln(2) split in two pieces. */
149     private final Dfp[] ln2Split;
150 
151     /** A {@link Dfp} with value ln(5). */
152     private final Dfp ln5;
153 
154     /** A two elements {@link Dfp} array with value ln(5) split in two pieces. */
155     private final Dfp[] ln5Split;
156 
157     /** A {@link Dfp} with value ln(10). */
158     private final Dfp ln10;
159 
160     /** Current rounding mode. */
161     private RoundingMode rMode;
162 
163     /** IEEE 854-1987 signals. */
164     private int ieeeFlags;
165 
166     /** Create a factory for the specified number of radix digits.
167      * <p>
168      * Note that since the {@link Dfp} class uses 10000 as its radix, each radix
169      * digit is equivalent to 4 decimal digits. This implies that asking for
170      * 13, 14, 15 or 16 decimal digits will really lead to a 4 radix 10000 digits in
171      * all cases.
172      * </p>
173      * @param decimalDigits minimal number of decimal digits.
174      */
175     public DfpField(final int decimalDigits) {
176         this(decimalDigits, true);
177     }
178 
179     /** Create a factory for the specified number of radix digits.
180      * <p>
181      * Note that since the {@link Dfp} class uses 10000 as its radix, each radix
182      * digit is equivalent to 4 decimal digits. This implies that asking for
183      * 13, 14, 15 or 16 decimal digits will really lead to a 4 radix 10000 digits in
184      * all cases.
185      * </p>
186      * @param decimalDigits minimal number of decimal digits
187      * @param computeConstants if true, the transcendental constants for the given precision
188      * must be computed (setting this flag to false is RESERVED for the internal recursive call)
189      */
190     private DfpField(final int decimalDigits, final boolean computeConstants) {
191 
192         this.radixDigits = (decimalDigits < 13) ? 4 : (decimalDigits + 3) / 4;
193         this.rMode       = RoundingMode.ROUND_HALF_EVEN;
194         this.ieeeFlags   = 0;
195         this.zero        = new Dfp(this, 0);
196         this.one         = new Dfp(this, 1);
197         this.two         = new Dfp(this, 2);
198 
199         if (computeConstants) {
200             // set up transcendental constants
201             synchronized (DfpField.class) {
202 
203                 // as a heuristic to circumvent Table-Maker's Dilemma, we set the string
204                 // representation of the constants to be at least 3 times larger than the
205                 // number of decimal digits, also as an attempt to really compute these
206                 // constants only once, we set a minimum number of digits
207                 computeStringConstants((decimalDigits < 67) ? 200 : (3 * decimalDigits));
208 
209                 // set up the constants at current field accuracy
210                 sqr2           = new Dfp(this, sqr2String);
211                 sqr2Split      = split(sqr2String);
212                 sqr2Reciprocal = new Dfp(this, sqr2ReciprocalString);
213                 sqr3           = new Dfp(this, sqr3String);
214                 sqr3Reciprocal = new Dfp(this, sqr3ReciprocalString);
215                 pi             = new Dfp(this, piString);
216                 piSplit        = split(piString);
217                 e              = new Dfp(this, eString);
218                 eSplit         = split(eString);
219                 ln2            = new Dfp(this, ln2String);
220                 ln2Split       = split(ln2String);
221                 ln5            = new Dfp(this, ln5String);
222                 ln5Split       = split(ln5String);
223                 ln10           = new Dfp(this, ln10String);
224 
225             }
226         } else {
227             // dummy settings for unused constants
228             sqr2           = null;
229             sqr2Split      = null;
230             sqr2Reciprocal = null;
231             sqr3           = null;
232             sqr3Reciprocal = null;
233             pi             = null;
234             piSplit        = null;
235             e              = null;
236             eSplit         = null;
237             ln2            = null;
238             ln2Split       = null;
239             ln5            = null;
240             ln5Split       = null;
241             ln10           = null;
242         }
243 
244     }
245 
246     /** Get the number of radix digits of the {@link Dfp} instances built by this factory.
247      * @return number of radix digits
248      */
249     public int getRadixDigits() {
250         return radixDigits;
251     }
252 
253     /** Set the rounding mode.
254      *  If not set, the default value is {@link RoundingMode#ROUND_HALF_EVEN}.
255      * @param mode desired rounding mode
256      * Note that the rounding mode is common to all {@link Dfp} instances
257      * belonging to the current {@link DfpField} in the system and will
258      * affect all future calculations.
259      */
260     public void setRoundingMode(final RoundingMode mode) {
261         rMode = mode;
262     }
263 
264     /** Get the current rounding mode.
265      * @return current rounding mode
266      */
267     public RoundingMode getRoundingMode() {
268         return rMode;
269     }
270 
271     /** Get the IEEE 854 status flags.
272      * @return IEEE 854 status flags
273      * @see #clearIEEEFlags()
274      * @see #setIEEEFlags(int)
275      * @see #setIEEEFlagsBits(int)
276      * @see #FLAG_INVALID
277      * @see #FLAG_DIV_ZERO
278      * @see #FLAG_OVERFLOW
279      * @see #FLAG_UNDERFLOW
280      * @see #FLAG_INEXACT
281      */
282     public int getIEEEFlags() {
283         return ieeeFlags;
284     }
285 
286     /** Clears the IEEE 854 status flags.
287      * @see #getIEEEFlags()
288      * @see #setIEEEFlags(int)
289      * @see #setIEEEFlagsBits(int)
290      * @see #FLAG_INVALID
291      * @see #FLAG_DIV_ZERO
292      * @see #FLAG_OVERFLOW
293      * @see #FLAG_UNDERFLOW
294      * @see #FLAG_INEXACT
295      */
296     public void clearIEEEFlags() {
297         ieeeFlags = 0;
298     }
299 
300     /** Sets the IEEE 854 status flags.
301      * @param flags desired value for the flags
302      * @see #getIEEEFlags()
303      * @see #clearIEEEFlags()
304      * @see #setIEEEFlagsBits(int)
305      * @see #FLAG_INVALID
306      * @see #FLAG_DIV_ZERO
307      * @see #FLAG_OVERFLOW
308      * @see #FLAG_UNDERFLOW
309      * @see #FLAG_INEXACT
310      */
311     public void setIEEEFlags(final int flags) {
312         ieeeFlags = flags & (FLAG_INVALID | FLAG_DIV_ZERO | FLAG_OVERFLOW | FLAG_UNDERFLOW | FLAG_INEXACT);
313     }
314 
315     /** Sets some bits in the IEEE 854 status flags, without changing the already set bits.
316      * <p>
317      * Calling this method is equivalent to call {@code setIEEEFlags(getIEEEFlags() | bits)}
318      * </p>
319      * @param bits bits to set
320      * @see #getIEEEFlags()
321      * @see #clearIEEEFlags()
322      * @see #setIEEEFlags(int)
323      * @see #FLAG_INVALID
324      * @see #FLAG_DIV_ZERO
325      * @see #FLAG_OVERFLOW
326      * @see #FLAG_UNDERFLOW
327      * @see #FLAG_INEXACT
328      */
329     public void setIEEEFlagsBits(final int bits) {
330         ieeeFlags |= bits & (FLAG_INVALID | FLAG_DIV_ZERO | FLAG_OVERFLOW | FLAG_UNDERFLOW | FLAG_INEXACT);
331     }
332 
333     /** Makes a {@link Dfp} with a value of 0.
334      * @return a new {@link Dfp} with a value of 0
335      */
336     public Dfp newDfp() {
337         return new Dfp(this);
338     }
339 
340     /** Create an instance from a byte value.
341      * @param x value to convert to an instance
342      * @return a new {@link Dfp} with the same value as x
343      */
344     public Dfp newDfp(final byte x) {
345         return new Dfp(this, x);
346     }
347 
348     /** Create an instance from an int value.
349      * @param x value to convert to an instance
350      * @return a new {@link Dfp} with the same value as x
351      */
352     public Dfp newDfp(final int x) {
353         return new Dfp(this, x);
354     }
355 
356     /** Create an instance from a long value.
357      * @param x value to convert to an instance
358      * @return a new {@link Dfp} with the same value as x
359      */
360     public Dfp newDfp(final long x) {
361         return new Dfp(this, x);
362     }
363 
364     /** Create an instance from a double value.
365      * @param x value to convert to an instance
366      * @return a new {@link Dfp} with the same value as x
367      */
368     public Dfp newDfp(final double x) {
369         return new Dfp(this, x);
370     }
371 
372     /** Copy constructor.
373      * @param d instance to copy
374      * @return a new {@link Dfp} with the same value as d
375      */
376     public Dfp newDfp(Dfp d) {
377         return new Dfp(d);
378     }
379 
380     /** Create a {@link Dfp} given a String representation.
381      * @param s string representation of the instance
382      * @return a new {@link Dfp} parsed from specified string
383      */
384     public Dfp newDfp(final String s) {
385         return new Dfp(this, s);
386     }
387 
388     /** Creates a {@link Dfp} with a non-finite value.
389      * @param sign sign of the Dfp to create
390      * @param nans code of the value, must be one of {@link Dfp#INFINITE},
391      * {@link Dfp#SNAN},  {@link Dfp#QNAN}
392      * @return a new {@link Dfp} with a non-finite value
393      */
394     public Dfp newDfp(final byte sign, final byte nans) {
395         return new Dfp(this, sign, nans);
396     }
397 
398     /** Get the constant 0.
399      * @return a {@link Dfp} with value 0
400      */
401     public Dfp getZero() {
402         return zero;
403     }
404 
405     /** Get the constant 1.
406      * @return a {@link Dfp} with value 1
407      */
408     public Dfp getOne() {
409         return one;
410     }
411 
412     /** {@inheritDoc} */
413     public Class<? extends FieldElement<Dfp>> getRuntimeClass() {
414         return Dfp.class;
415     }
416 
417     /** Get the constant 2.
418      * @return a {@link Dfp} with value 2
419      */
420     public Dfp getTwo() {
421         return two;
422     }
423 
424     /** Get the constant &radic;2.
425      * @return a {@link Dfp} with value &radic;2
426      */
427     public Dfp getSqr2() {
428         return sqr2;
429     }
430 
431     /** Get the constant &radic;2 split in two pieces.
432      * @return a {@link Dfp} with value &radic;2 split in two pieces
433      */
434     public Dfp[] getSqr2Split() {
435         return sqr2Split.clone();
436     }
437 
438     /** Get the constant &radic;2 / 2.
439      * @return a {@link Dfp} with value &radic;2 / 2
440      */
441     public Dfp getSqr2Reciprocal() {
442         return sqr2Reciprocal;
443     }
444 
445     /** Get the constant &radic;3.
446      * @return a {@link Dfp} with value &radic;3
447      */
448     public Dfp getSqr3() {
449         return sqr3;
450     }
451 
452     /** Get the constant &radic;3 / 3.
453      * @return a {@link Dfp} with value &radic;3 / 3
454      */
455     public Dfp getSqr3Reciprocal() {
456         return sqr3Reciprocal;
457     }
458 
459     /** Get the constant &pi;.
460      * @return a {@link Dfp} with value &pi;
461      */
462     public Dfp getPi() {
463         return pi;
464     }
465 
466     /** Get the constant &pi; split in two pieces.
467      * @return a {@link Dfp} with value &pi; split in two pieces
468      */
469     public Dfp[] getPiSplit() {
470         return piSplit.clone();
471     }
472 
473     /** Get the constant e.
474      * @return a {@link Dfp} with value e
475      */
476     public Dfp getE() {
477         return e;
478     }
479 
480     /** Get the constant e split in two pieces.
481      * @return a {@link Dfp} with value e split in two pieces
482      */
483     public Dfp[] getESplit() {
484         return eSplit.clone();
485     }
486 
487     /** Get the constant ln(2).
488      * @return a {@link Dfp} with value ln(2)
489      */
490     public Dfp getLn2() {
491         return ln2;
492     }
493 
494     /** Get the constant ln(2) split in two pieces.
495      * @return a {@link Dfp} with value ln(2) split in two pieces
496      */
497     public Dfp[] getLn2Split() {
498         return ln2Split.clone();
499     }
500 
501     /** Get the constant ln(5).
502      * @return a {@link Dfp} with value ln(5)
503      */
504     public Dfp getLn5() {
505         return ln5;
506     }
507 
508     /** Get the constant ln(5) split in two pieces.
509      * @return a {@link Dfp} with value ln(5) split in two pieces
510      */
511     public Dfp[] getLn5Split() {
512         return ln5Split.clone();
513     }
514 
515     /** Get the constant ln(10).
516      * @return a {@link Dfp} with value ln(10)
517      */
518     public Dfp getLn10() {
519         return ln10;
520     }
521 
522     /** Breaks a string representation up into two {@link Dfp}'s.
523      * The split is such that the sum of them is equivalent to the input string,
524      * but has higher precision than using a single Dfp.
525      * @param a string representation of the number to split
526      * @return an array of two {@link Dfp Dfp} instances which sum equals a
527      */
528     private Dfp[] split(final String a) {
529       Dfp result[] = new Dfp[2];
530       boolean leading = true;
531       int sp = 0;
532       int sig = 0;
533 
534       char[] buf = new char[a.length()];
535 
536       for (int i = 0; i < buf.length; i++) {
537         buf[i] = a.charAt(i);
538 
539         if (buf[i] >= '1' && buf[i] <= '9') {
540             leading = false;
541         }
542 
543         if (buf[i] == '.') {
544           sig += (400 - sig) % 4;
545           leading = false;
546         }
547 
548         if (sig == (radixDigits / 2) * 4) {
549           sp = i;
550           break;
551         }
552 
553         if (buf[i] >= '0' && buf[i] <= '9' && !leading) {
554             sig ++;
555         }
556       }
557 
558       result[0] = new Dfp(this, new String(buf, 0, sp));
559 
560       for (int i = 0; i < buf.length; i++) {
561         buf[i] = a.charAt(i);
562         if (buf[i] >= '0' && buf[i] <= '9' && i < sp) {
563             buf[i] = '0';
564         }
565       }
566 
567       result[1] = new Dfp(this, new String(buf));
568 
569       return result;
570 
571     }
572 
573     /** Recompute the high precision string constants.
574      * @param highPrecisionDecimalDigits precision at which the string constants mus be computed
575      */
576     private static void computeStringConstants(final int highPrecisionDecimalDigits) {
577         if (sqr2String == null || sqr2String.length() < highPrecisionDecimalDigits - 3) {
578 
579             // recompute the string representation of the transcendental constants
580             final DfpField highPrecisionField = new DfpField(highPrecisionDecimalDigits, false);
581             final Dfp highPrecisionOne        = new Dfp(highPrecisionField, 1);
582             final Dfp highPrecisionTwo        = new Dfp(highPrecisionField, 2);
583             final Dfp highPrecisionThree      = new Dfp(highPrecisionField, 3);
584 
585             final Dfp highPrecisionSqr2 = highPrecisionTwo.sqrt();
586             sqr2String           = highPrecisionSqr2.toString();
587             sqr2ReciprocalString = highPrecisionOne.divide(highPrecisionSqr2).toString();
588 
589             final Dfp highPrecisionSqr3 = highPrecisionThree.sqrt();
590             sqr3String           = highPrecisionSqr3.toString();
591             sqr3ReciprocalString = highPrecisionOne.divide(highPrecisionSqr3).toString();
592 
593             piString   = computePi(highPrecisionOne, highPrecisionTwo, highPrecisionThree).toString();
594             eString    = computeExp(highPrecisionOne, highPrecisionOne).toString();
595             ln2String  = computeLn(highPrecisionTwo, highPrecisionOne, highPrecisionTwo).toString();
596             ln5String  = computeLn(new Dfp(highPrecisionField, 5),  highPrecisionOne, highPrecisionTwo).toString();
597             ln10String = computeLn(new Dfp(highPrecisionField, 10), highPrecisionOne, highPrecisionTwo).toString();
598 
599         }
600     }
601 
602     /** Compute &pi; using Jonathan and Peter Borwein quartic formula.
603      * @param one constant with value 1 at desired precision
604      * @param two constant with value 2 at desired precision
605      * @param three constant with value 3 at desired precision
606      * @return &pi;
607      */
608     private static Dfp computePi(final Dfp one, final Dfp two, final Dfp three) {
609 
610         Dfp sqrt2   = two.sqrt();
611         Dfp yk      = sqrt2.subtract(one);
612         Dfp four    = two.add(two);
613         Dfp two2kp3 = two;
614         Dfp ak      = two.multiply(three.subtract(two.multiply(sqrt2)));
615 
616         // The formula converges quartically. This means the number of correct
617         // digits is multiplied by 4 at each iteration! Five iterations are
618         // sufficient for about 160 digits, eight iterations give about
619         // 10000 digits (this has been checked) and 20 iterations more than
620         // 160 billions of digits (this has NOT been checked).
621         // So the limit here is considered sufficient for most purposes ...
622         for (int i = 1; i < 20; i++) {
623             final Dfp ykM1 = yk;
624 
625             final Dfp y2         = yk.multiply(yk);
626             final Dfp oneMinusY4 = one.subtract(y2.multiply(y2));
627             final Dfp s          = oneMinusY4.sqrt().sqrt();
628             yk = one.subtract(s).divide(one.add(s));
629 
630             two2kp3 = two2kp3.multiply(four);
631 
632             final Dfp p = one.add(yk);
633             final Dfp p2 = p.multiply(p);
634             ak = ak.multiply(p2.multiply(p2)).subtract(two2kp3.multiply(yk).multiply(one.add(yk).add(yk.multiply(yk))));
635 
636             if (yk.equals(ykM1)) {
637                 break;
638             }
639         }
640 
641         return one.divide(ak);
642 
643     }
644 
645     /** Compute exp(a).
646      * @param a number for which we want the exponential
647      * @param one constant with value 1 at desired precision
648      * @return exp(a)
649      */
650     public static Dfp computeExp(final Dfp a, final Dfp one) {
651 
652         Dfp y  = new Dfp(one);
653         Dfp py = new Dfp(one);
654         Dfp f  = new Dfp(one);
655         Dfp fi = new Dfp(one);
656         Dfp x  = new Dfp(one);
657 
658         for (int i = 0; i < 10000; i++) {
659             x = x.multiply(a);
660             y = y.add(x.divide(f));
661             fi = fi.add(one);
662             f = f.multiply(fi);
663             if (y.equals(py)) {
664                 break;
665             }
666             py = new Dfp(y);
667         }
668 
669         return y;
670 
671     }
672 
673 
674     /** Compute ln(a).
675      *
676      *  Let f(x) = ln(x),
677      *
678      *  We know that f'(x) = 1/x, thus from Taylor's theorem we have:
679      *
680      *           -----          n+1         n
681      *  f(x) =   \           (-1)    (x - 1)
682      *           /          ----------------    for 1 <= n <= infinity
683      *           -----             n
684      *
685      *  or
686      *                       2        3       4
687      *                   (x-1)   (x-1)    (x-1)
688      *  ln(x) =  (x-1) - ----- + ------ - ------ + ...
689      *                     2       3        4
690      *
691      *  alternatively,
692      *
693      *                  2    3   4
694      *                 x    x   x
695      *  ln(x+1) =  x - -  + - - - + ...
696      *                 2    3   4
697      *
698      *  This series can be used to compute ln(x), but it converges too slowly.
699      *
700      *  If we substitute -x for x above, we get
701      *
702      *                   2    3    4
703      *                  x    x    x
704      *  ln(1-x) =  -x - -  - -  - - + ...
705      *                  2    3    4
706      *
707      *  Note that all terms are now negative.  Because the even powered ones
708      *  absorbed the sign.  Now, subtract the series above from the previous
709      *  one to get ln(x+1) - ln(1-x).  Note the even terms cancel out leaving
710      *  only the odd ones
711      *
712      *                             3     5      7
713      *                           2x    2x     2x
714      *  ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ...
715      *                            3     5      7
716      *
717      *  By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have:
718      *
719      *                                3        5        7
720      *      x+1           /          x        x        x          \
721      *  ln ----- =   2 *  |  x  +   ----  +  ----  +  ---- + ...  |
722      *      x-1           \          3        5        7          /
723      *
724      *  But now we want to find ln(a), so we need to find the value of x
725      *  such that a = (x+1)/(x-1).   This is easily solved to find that
726      *  x = (a-1)/(a+1).
727      * @param a number for which we want the exponential
728      * @param one constant with value 1 at desired precision
729      * @param two constant with value 2 at desired precision
730      * @return ln(a)
731      */
732 
733     public static Dfp computeLn(final Dfp a, final Dfp one, final Dfp two) {
734 
735         int den = 1;
736         Dfp x = a.add(new Dfp(a.getField(), -1)).divide(a.add(one));
737 
738         Dfp y = new Dfp(x);
739         Dfp num = new Dfp(x);
740         Dfp py = new Dfp(y);
741         for (int i = 0; i < 10000; i++) {
742             num = num.multiply(x);
743             num = num.multiply(x);
744             den += 2;
745             Dfp t = num.divide(den);
746             y = y.add(t);
747             if (y.equals(py)) {
748                 break;
749             }
750             py = new Dfp(y);
751         }
752 
753         return y.multiply(two);
754 
755     }
756 
757 }