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