001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math3.dfp;
019    
020    import org.apache.commons.math3.Field;
021    import org.apache.commons.math3.FieldElement;
022    
023    /** Field for Decimal floating point instances.
024     * @version $Id: DfpField.java 1416643 2012-12-03 19:37:14Z tn $
025     * @since 2.2
026     */
027    public class DfpField implements Field<Dfp> {
028    
029        /** Enumerate for rounding modes. */
030        public enum RoundingMode {
031    
032            /** Rounds toward zero (truncation). */
033            ROUND_DOWN,
034    
035            /** Rounds away from zero if discarded digit is non-zero. */
036            ROUND_UP,
037    
038            /** Rounds towards nearest unless both are equidistant in which case it rounds away from zero. */
039            ROUND_HALF_UP,
040    
041            /** Rounds towards nearest unless both are equidistant in which case it rounds toward zero. */
042            ROUND_HALF_DOWN,
043    
044            /** Rounds towards nearest unless both are equidistant in which case it rounds toward the even neighbor.
045             * This is the default as  specified by IEEE 854-1987
046             */
047            ROUND_HALF_EVEN,
048    
049            /** Rounds towards nearest unless both are equidistant in which case it rounds toward the odd neighbor.  */
050            ROUND_HALF_ODD,
051    
052            /** Rounds towards positive infinity. */
053            ROUND_CEIL,
054    
055            /** Rounds towards negative infinity. */
056            ROUND_FLOOR;
057    
058        }
059    
060        /** IEEE 854-1987 flag for invalid operation. */
061        public static final int FLAG_INVALID   =  1;
062    
063        /** IEEE 854-1987 flag for division by zero. */
064        public static final int FLAG_DIV_ZERO  =  2;
065    
066        /** IEEE 854-1987 flag for overflow. */
067        public static final int FLAG_OVERFLOW  =  4;
068    
069        /** IEEE 854-1987 flag for underflow. */
070        public static final int FLAG_UNDERFLOW =  8;
071    
072        /** IEEE 854-1987 flag for inexact result. */
073        public static final int FLAG_INEXACT   = 16;
074    
075        /** High precision string representation of &radic;2. */
076        private static String sqr2String;
077    
078        // Note: the static strings are set up (once) by the ctor and @GuardedBy("DfpField.class")
079    
080        /** High precision string representation of &radic;2 / 2. */
081        private static String sqr2ReciprocalString;
082    
083        /** High precision string representation of &radic;3. */
084        private static String sqr3String;
085    
086        /** High precision string representation of &radic;3 / 3. */
087        private static String sqr3ReciprocalString;
088    
089        /** High precision string representation of &pi;. */
090        private static String piString;
091    
092        /** High precision string representation of e. */
093        private static String eString;
094    
095        /** High precision string representation of ln(2). */
096        private static String ln2String;
097    
098        /** High precision string representation of ln(5). */
099        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 = 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    }