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