    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    }