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
018package org.apache.commons.math4.dfp;
019
020import org.apache.commons.math4.Field;
021import org.apache.commons.math4.FieldElement;
022
023/** Field for Decimal floating point instances.
024 * @since 2.2
025 */
026public class DfpField implements Field<Dfp> {
027
028    /** Enumerate for rounding modes. */
029    public enum RoundingMode {
030
031        /** Rounds toward zero (truncation). */
032        ROUND_DOWN,
033
034        /** Rounds away from zero if discarded digit is non-zero. */
035        ROUND_UP,
036
037        /** Rounds towards nearest unless both are equidistant in which case it rounds away from zero. */
038        ROUND_HALF_UP,
039
040        /** Rounds towards nearest unless both are equidistant in which case it rounds toward zero. */
041        ROUND_HALF_DOWN,
042
043        /** Rounds towards nearest unless both are equidistant in which case it rounds toward the even neighbor.
044         * This is the default as  specified by IEEE 854-1987
045         */
046        ROUND_HALF_EVEN,
047
048        /** Rounds towards nearest unless both are equidistant in which case it rounds toward the odd neighbor.  */
049        ROUND_HALF_ODD,
050
051        /** Rounds towards positive infinity. */
052        ROUND_CEIL,
053
054        /** Rounds towards negative infinity. */
055        ROUND_FLOOR;
056
057    }
058
059    /** IEEE 854-1987 flag for invalid operation. */
060    public static final int FLAG_INVALID   =  1;
061
062    /** IEEE 854-1987 flag for division by zero. */
063    public static final int FLAG_DIV_ZERO  =  2;
064
065    /** IEEE 854-1987 flag for overflow. */
066    public static final int FLAG_OVERFLOW  =  4;
067
068    /** IEEE 854-1987 flag for underflow. */
069    public static final int FLAG_UNDERFLOW =  8;
070
071    /** IEEE 854-1987 flag for inexact result. */
072    public static final int FLAG_INEXACT   = 16;
073
074    /** High precision string representation of &radic;2. */
075    private static String sqr2String;
076
077    // Note: the static strings are set up (once) by the ctor and @GuardedBy("DfpField.class")
078
079    /** High precision string representation of &radic;2 / 2. */
080    private static String sqr2ReciprocalString;
081
082    /** High precision string representation of &radic;3. */
083    private static String sqr3String;
084
085    /** High precision string representation of &radic;3 / 3. */
086    private static String sqr3ReciprocalString;
087
088    /** High precision string representation of &pi;. */
089    private static String piString;
090
091    /** High precision string representation of e. */
092    private static String eString;
093
094    /** High precision string representation of ln(2). */
095    private static String ln2String;
096
097    /** High precision string representation of ln(5). */
098    private static String ln5String;
099
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    @Override
402    public Dfp getZero() {
403        return zero;
404    }
405
406    /** Get the constant 1.
407     * @return a {@link Dfp} with value 1
408     */
409    @Override
410    public Dfp getOne() {
411        return one;
412    }
413
414    /** {@inheritDoc} */
415    @Override
416    public Class<? extends FieldElement<Dfp>> getRuntimeClass() {
417        return Dfp.class;
418    }
419
420    /** Get the constant 2.
421     * @return a {@link Dfp} with value 2
422     */
423    public Dfp getTwo() {
424        return two;
425    }
426
427    /** Get the constant &radic;2.
428     * @return a {@link Dfp} with value &radic;2
429     */
430    public Dfp getSqr2() {
431        return sqr2;
432    }
433
434    /** Get the constant &radic;2 split in two pieces.
435     * @return a {@link Dfp} with value &radic;2 split in two pieces
436     */
437    public Dfp[] getSqr2Split() {
438        return sqr2Split.clone();
439    }
440
441    /** Get the constant &radic;2 / 2.
442     * @return a {@link Dfp} with value &radic;2 / 2
443     */
444    public Dfp getSqr2Reciprocal() {
445        return sqr2Reciprocal;
446    }
447
448    /** Get the constant &radic;3.
449     * @return a {@link Dfp} with value &radic;3
450     */
451    public Dfp getSqr3() {
452        return sqr3;
453    }
454
455    /** Get the constant &radic;3 / 3.
456     * @return a {@link Dfp} with value &radic;3 / 3
457     */
458    public Dfp getSqr3Reciprocal() {
459        return sqr3Reciprocal;
460    }
461
462    /** Get the constant &pi;.
463     * @return a {@link Dfp} with value &pi;
464     */
465    public Dfp getPi() {
466        return pi;
467    }
468
469    /** Get the constant &pi; split in two pieces.
470     * @return a {@link Dfp} with value &pi; split in two pieces
471     */
472    public Dfp[] getPiSplit() {
473        return piSplit.clone();
474    }
475
476    /** Get the constant e.
477     * @return a {@link Dfp} with value e
478     */
479    public Dfp getE() {
480        return e;
481    }
482
483    /** Get the constant e split in two pieces.
484     * @return a {@link Dfp} with value e split in two pieces
485     */
486    public Dfp[] getESplit() {
487        return eSplit.clone();
488    }
489
490    /** Get the constant ln(2).
491     * @return a {@link Dfp} with value ln(2)
492     */
493    public Dfp getLn2() {
494        return ln2;
495    }
496
497    /** Get the constant ln(2) split in two pieces.
498     * @return a {@link Dfp} with value ln(2) split in two pieces
499     */
500    public Dfp[] getLn2Split() {
501        return ln2Split.clone();
502    }
503
504    /** Get the constant ln(5).
505     * @return a {@link Dfp} with value ln(5)
506     */
507    public Dfp getLn5() {
508        return ln5;
509    }
510
511    /** Get the constant ln(5) split in two pieces.
512     * @return a {@link Dfp} with value ln(5) split in two pieces
513     */
514    public Dfp[] getLn5Split() {
515        return ln5Split.clone();
516    }
517
518    /** Get the constant ln(10).
519     * @return a {@link Dfp} with value ln(10)
520     */
521    public Dfp getLn10() {
522        return ln10;
523    }
524
525    /** Breaks a string representation up into two {@link Dfp}'s.
526     * The split is such that the sum of them is equivalent to the input string,
527     * but has higher precision than using a single Dfp.
528     * @param a string representation of the number to split
529     * @return an array of two {@link Dfp Dfp} instances which sum equals a
530     */
531    private Dfp[] split(final String a) {
532      Dfp result[] = new Dfp[2];
533      boolean leading = true;
534      int sp = 0;
535      int sig = 0;
536
537      char[] buf = new char[a.length()];
538
539      for (int i = 0; i < buf.length; i++) {
540        buf[i] = a.charAt(i);
541
542        if (buf[i] >= '1' && buf[i] <= '9') {
543            leading = false;
544        }
545
546        if (buf[i] == '.') {
547          sig += (400 - sig) % 4;
548          leading = false;
549        }
550
551        if (sig == (radixDigits / 2) * 4) {
552          sp = i;
553          break;
554        }
555
556        if (buf[i] >= '0' && buf[i] <= '9' && !leading) {
557            sig ++;
558        }
559      }
560
561      result[0] = new Dfp(this, new String(buf, 0, sp));
562
563      for (int i = 0; i < buf.length; i++) {
564        buf[i] = a.charAt(i);
565        if (buf[i] >= '0' && buf[i] <= '9' && i < sp) {
566            buf[i] = '0';
567        }
568      }
569
570      result[1] = new Dfp(this, new String(buf));
571
572      return result;
573
574    }
575
576    /** Recompute the high precision string constants.
577     * @param highPrecisionDecimalDigits precision at which the string constants mus be computed
578     */
579    private static void computeStringConstants(final int highPrecisionDecimalDigits) {
580        if (sqr2String == null || sqr2String.length() < highPrecisionDecimalDigits - 3) {
581
582            // recompute the string representation of the transcendental constants
583            final DfpField highPrecisionField = new DfpField(highPrecisionDecimalDigits, false);
584            final Dfp highPrecisionOne        = new Dfp(highPrecisionField, 1);
585            final Dfp highPrecisionTwo        = new Dfp(highPrecisionField, 2);
586            final Dfp highPrecisionThree      = new Dfp(highPrecisionField, 3);
587
588            final Dfp highPrecisionSqr2 = highPrecisionTwo.sqrt();
589            sqr2String           = highPrecisionSqr2.toString();
590            sqr2ReciprocalString = highPrecisionOne.divide(highPrecisionSqr2).toString();
591
592            final Dfp highPrecisionSqr3 = highPrecisionThree.sqrt();
593            sqr3String           = highPrecisionSqr3.toString();
594            sqr3ReciprocalString = highPrecisionOne.divide(highPrecisionSqr3).toString();
595
596            piString   = computePi(highPrecisionOne, highPrecisionTwo, highPrecisionThree).toString();
597            eString    = computeExp(highPrecisionOne, highPrecisionOne).toString();
598            ln2String  = computeLn(highPrecisionTwo, highPrecisionOne, highPrecisionTwo).toString();
599            ln5String  = computeLn(new Dfp(highPrecisionField, 5),  highPrecisionOne, highPrecisionTwo).toString();
600            ln10String = computeLn(new Dfp(highPrecisionField, 10), highPrecisionOne, highPrecisionTwo).toString();
601
602        }
603    }
604
605    /** Compute &pi; using Jonathan and Peter Borwein quartic formula.
606     * @param one constant with value 1 at desired precision
607     * @param two constant with value 2 at desired precision
608     * @param three constant with value 3 at desired precision
609     * @return &pi;
610     */
611    private static Dfp computePi(final Dfp one, final Dfp two, final Dfp three) {
612
613        Dfp sqrt2   = two.sqrt();
614        Dfp yk      = sqrt2.subtract(one);
615        Dfp four    = two.add(two);
616        Dfp two2kp3 = two;
617        Dfp ak      = two.multiply(three.subtract(two.multiply(sqrt2)));
618
619        // The formula converges quartically. This means the number of correct
620        // digits is multiplied by 4 at each iteration! Five iterations are
621        // sufficient for about 160 digits, eight iterations give about
622        // 10000 digits (this has been checked) and 20 iterations more than
623        // 160 billions of digits (this has NOT been checked).
624        // So the limit here is considered sufficient for most purposes ...
625        for (int i = 1; i < 20; i++) {
626            final Dfp ykM1 = yk;
627
628            final Dfp y2         = yk.multiply(yk);
629            final Dfp oneMinusY4 = one.subtract(y2.multiply(y2));
630            final Dfp s          = oneMinusY4.sqrt().sqrt();
631            yk = one.subtract(s).divide(one.add(s));
632
633            two2kp3 = two2kp3.multiply(four);
634
635            final Dfp p = one.add(yk);
636            final Dfp p2 = p.multiply(p);
637            ak = ak.multiply(p2.multiply(p2)).subtract(two2kp3.multiply(yk).multiply(one.add(yk).add(yk.multiply(yk))));
638
639            if (yk.equals(ykM1)) {
640                break;
641            }
642        }
643
644        return one.divide(ak);
645
646    }
647
648    /** Compute exp(a).
649     * @param a number for which we want the exponential
650     * @param one constant with value 1 at desired precision
651     * @return exp(a)
652     */
653    public static Dfp computeExp(final Dfp a, final Dfp one) {
654
655        Dfp y  = new Dfp(one);
656        Dfp py = new Dfp(one);
657        Dfp f  = new Dfp(one);
658        Dfp fi = new Dfp(one);
659        Dfp x  = new Dfp(one);
660
661        for (int i = 0; i < 10000; i++) {
662            x = x.multiply(a);
663            y = y.add(x.divide(f));
664            fi = fi.add(one);
665            f = f.multiply(fi);
666            if (y.equals(py)) {
667                break;
668            }
669            py = new Dfp(y);
670        }
671
672        return y;
673
674    }
675
676
677    /** Compute ln(a).
678     *
679     *  <pre>{@code
680     *  Let f(x) = ln(x),
681     *
682     *  We know that f'(x) = 1/x, thus from Taylor's theorem we have:
683     *
684     *           -----          n+1         n
685     *  f(x) =   \           (-1)    (x - 1)
686     *           /          ----------------    for 1 <= n <= infinity
687     *           -----             n
688     *
689     *  or
690     *                       2        3       4
691     *                   (x-1)   (x-1)    (x-1)
692     *  ln(x) =  (x-1) - ----- + ------ - ------ + ...
693     *                     2       3        4
694     *
695     *  alternatively,
696     *
697     *                  2    3   4
698     *                 x    x   x
699     *  ln(x+1) =  x - -  + - - - + ...
700     *                 2    3   4
701     *
702     *  This series can be used to compute ln(x), but it converges too slowly.
703     *
704     *  If we substitute -x for x above, we get
705     *
706     *                   2    3    4
707     *                  x    x    x
708     *  ln(1-x) =  -x - -  - -  - - + ...
709     *                  2    3    4
710     *
711     *  Note that all terms are now negative.  Because the even powered ones
712     *  absorbed the sign.  Now, subtract the series above from the previous
713     *  one to get ln(x+1) - ln(1-x).  Note the even terms cancel out leaving
714     *  only the odd ones
715     *
716     *                             3     5      7
717     *                           2x    2x     2x
718     *  ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ...
719     *                            3     5      7
720     *
721     *  By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have:
722     *
723     *                                3        5        7
724     *      x+1           /          x        x        x          \
725     *  ln ----- =   2 *  |  x  +   ----  +  ----  +  ---- + ...  |
726     *      x-1           \          3        5        7          /
727     *
728     *  But now we want to find ln(a), so we need to find the value of x
729     *  such that a = (x+1)/(x-1).   This is easily solved to find that
730     *  x = (a-1)/(a+1).
731     * }</pre>
732     * @param a number for which we want the exponential
733     * @param one constant with value 1 at desired precision
734     * @param two constant with value 2 at desired precision
735     * @return ln(a)
736     */
737
738    public static Dfp computeLn(final Dfp a, final Dfp one, final Dfp two) {
739
740        int den = 1;
741        Dfp x = a.add(new Dfp(a.getField(), -1)).divide(a.add(one));
742
743        Dfp y = new Dfp(x);
744        Dfp num = new Dfp(x);
745        Dfp py = new Dfp(y);
746        for (int i = 0; i < 10000; i++) {
747            num = num.multiply(x);
748            num = num.multiply(x);
749            den += 2;
750            Dfp t = num.divide(den);
751            y = y.add(t);
752            if (y.equals(py)) {
753                break;
754            }
755            py = new Dfp(y);
756        }
757
758        return y.multiply(two);
759
760    }
761
762}