    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    /** Mathematical routines for use with {@link Dfp}.
    021     * The constants are defined in {@link DfpField}
    022     * @version $Revision: 1042376 $ $Date: 2010-12-05 16:54:55 +0100 (dim. 05 d??c. 2010) $
    023     * @since 2.2
    024     */
    025    public class DfpMath {
    026
    027        /** Name for traps triggered by pow. */
    028        private static final String POW_TRAP = "pow";
    029
    030        /**
    031         * Private Constructor.
    032         */
    033        private DfpMath() {
    034        }
    035
    036        /** Breaks a string representation up into two dfp's.
    037         * <p>The two dfp are such that the sum of them is equivalent
    038         * to the input string, but has higher precision than using a
    039         * single dfp. This is useful for improving accuracy of
    040         * exponentiation and critical multiplies.
    041         * @param field field to which the Dfp must belong
    042         * @param a string representation to split
    043         * @return an array of two {@link Dfp} which sum is a
    044         */
    045        protected static Dfp[] split(final DfpField field, final String a) {
    046            Dfp result[] = new Dfp[2];
    047            char[] buf;
    048            boolean leading = true;
    049            int sp = 0;
    050            int sig = 0;
    051
    052            buf = new char[a.length()];
    053
    054            for (int i = 0; i < buf.length; i++) {
    055                buf[i] = a.charAt(i);
    056
    057                if (buf[i] >= '1' && buf[i] <= '9') {
    058                    leading = false;
    059                }
    060
    061                if (buf[i] == '.') {
    062                    sig += (400 - sig) % 4;
    063                    leading = false;
    064                }
    065
    066                if (sig == (field.getRadixDigits() / 2) * 4) {
    067                    sp = i;
    068                    break;
    069                }
    070
    071                if (buf[i] >= '0' && buf[i] <= '9' && !leading) {
    072                    sig ++;
    073                }
    074            }
    075
    076            result[0] = field.newDfp(new String(buf, 0, sp));
    077
    078            for (int i = 0; i < buf.length; i++) {
    079                buf[i] = a.charAt(i);
    080                if (buf[i] >= '0' && buf[i] <= '9' && i < sp) {
    081                    buf[i] = '0';
    082                }
    083            }
    084
    085            result[1] = field.newDfp(new String(buf));
    086
    087            return result;
    088        }
    089
    090        /** Splits a {@link Dfp} into 2 {@link Dfp}'s such that their sum is equal to the input {@link Dfp}.
    091         * @param a number to split
    092         * @return two elements array containing the split number
    093         */
    094        protected static Dfp[] split(final Dfp a) {
    095            final Dfp[] result = new Dfp[2];
    096            final Dfp shift = a.multiply(a.power10K(a.getRadixDigits() / 2));
    097            result[0] = a.add(shift).subtract(shift);
    098            result[1] = a.subtract(result[0]);
    099            return result;
    100        }
    101
    102        /** Multiply two numbers that are split in to two pieces that are
    103         *  meant to be added together.
    104         *  Use binomial multiplication so ab = a0 b0 + a0 b1 + a1 b0 + a1 b1
    105         *  Store the first term in result0, the rest in result1
    106         *  @param a first factor of the multiplication, in split form
    107         *  @param b second factor of the multiplication, in split form
    108         *  @return a &times; b, in split form
    109         */
    110        protected static Dfp[] splitMult(final Dfp[] a, final Dfp[] b) {
    111            final Dfp[] result = new Dfp[2];
    112
    113            result[1] = a[0].getZero();
    114            result[0] = a[0].multiply(b[0]);
    115
    116            /* If result[0] is infinite or zero, don't compute result[1].
    117             * Attempting to do so may produce NaNs.
    118             */
    119
    120            if (result[0].classify() == Dfp.INFINITE || result[0].equals(result[1])) {
    121                return result;
    122            }
    123
    124            result[1] = a[0].multiply(b[1]).add(a[1].multiply(b[0])).add(a[1].multiply(b[1]));
    125
    126            return result;
    127        }
    128
    129        /** Divide two numbers that are split in to two pieces that are meant to be added together.
    130         * Inverse of split multiply above:
    131         *  (a+b) / (c+d) = (a/c) + ( (bc-ad)/(c**2+cd) )
    132         *  @param a dividend, in split form
    133         *  @param b divisor, in split form
    134         *  @return a / b, in split form
    135         */
    136        protected static Dfp[] splitDiv(final Dfp[] a, final Dfp[] b) {
    137            final Dfp[] result;
    138
    139            result = new Dfp[2];
    140
    141            result[0] = a[0].divide(b[0]);
    142            result[1] = a[1].multiply(b[0]).subtract(a[0].multiply(b[1]));
    143            result[1] = result[1].divide(b[0].multiply(b[0]).add(b[0].multiply(b[1])));
    144
    145            return result;
    146        }
    147
    148        /** Raise a split base to the a power.
    149         * @param base number to raise
    150         * @param a power
    151         * @return base<sup>a</sup>
    152         */
    153        protected static Dfp splitPow(final Dfp[] base, int a) {
    154            boolean invert = false;
    155
    156            Dfp[] r = new Dfp[2];
    157
    158            Dfp[] result = new Dfp[2];
    159            result[0] = base[0].getOne();
    160            result[1] = base[0].getZero();
    161
    162            if (a == 0) {
    163                // Special case a = 0
    164                return result[0].add(result[1]);
    165            }
    166
    167            if (a < 0) {
    168                // If a is less than zero
    169                invert = true;
    170                a = -a;
    171            }
    172
    173            // Exponentiate by successive squaring
    174            do {
    175                r[0] = new Dfp(base[0]);
    176                r[1] = new Dfp(base[1]);
    177                int trial = 1;
    178
    179                int prevtrial;
    180                while (true) {
    181                    prevtrial = trial;
    182                    trial = trial * 2;
    183                    if (trial > a) {
    184                        break;
    185                    }
    186                    r = splitMult(r, r);
    187                }
    188
    189                trial = prevtrial;
    190
    191                a -= trial;
    192                result = splitMult(result, r);
    193
    194            } while (a >= 1);
    195
    196            result[0] = result[0].add(result[1]);
    197
    198            if (invert) {
    199                result[0] = base[0].getOne().divide(result[0]);
    200            }
    201
    202            return result[0];
    203
    204        }
    205
    206        /** Raises base to the power a by successive squaring.
    207         * @param base number to raise
    208         * @param a power
    209         * @return base<sup>a</sup>
    210         */
    211        public static Dfp pow(Dfp base, int a)
    212        {
    213            boolean invert = false;
    214
    215            Dfp result = base.getOne();
    216
    217            if (a == 0) {
    218                // Special case
    219                return result;
    220            }
    221
    222            if (a < 0) {
    223                invert = true;
    224                a = -a;
    225            }
    226
    227            // Exponentiate by successive squaring
    228            do {
    229                Dfp r = new Dfp(base);
    230                Dfp prevr;
    231                int trial = 1;
    232                int prevtrial;
    233
    234                do {
    235                    prevr = new Dfp(r);
    236                    prevtrial = trial;
    237                    r = r.multiply(r);
    238                    trial = trial * 2;
    239                } while (a>trial);
    240
    241                r = prevr;
    242                trial = prevtrial;
    243
    244                a = a - trial;
    245                result = result.multiply(r);
    246
    247            } while (a >= 1);
    248
    249            if (invert) {
    250                result = base.getOne().divide(result);
    251            }
    252
    253            return base.newInstance(result);
    254
    255        }
    256
    257        /** Computes e to the given power.
    258         * a is broken into two parts, such that a = n+m  where n is an integer.
    259         * We use pow() to compute e<sup>n</sup> and a Taylor series to compute
    260         * e<sup>m</sup>.  We return e*<sup>n</sup> &times; e<sup>m</sup>
    261         * @param a power at which e should be raised
    262         * @return e<sup>a</sup>
    263         */
    264        public static Dfp exp(final Dfp a) {
    265
    266            final Dfp inta = a.rint();
    267            final Dfp fraca = a.subtract(inta);
    268
    269            final int ia = inta.intValue();
    270            if (ia > 2147483646) {
    271                // return +Infinity
    272                return a.newInstance((byte)1, Dfp.INFINITE);
    273            }
    274
    275            if (ia < -2147483646) {
    276                // return 0;
    277                return a.newInstance();
    278            }
    279
    280            final Dfp einta = splitPow(a.getField().getESplit(), ia);
    281            final Dfp efraca = expInternal(fraca);
    282
    283            return einta.multiply(efraca);
    284        }
    285
    286        /** Computes e to the given power.
    287         * Where -1 < a < 1.  Use the classic Taylor series.  1 + x**2/2! + x**3/3! + x**4/4!  ...
    288         * @param a power at which e should be raised
    289         * @return e<sup>a</sup>
    290         */
    291        protected static Dfp expInternal(final Dfp a) {
    292            Dfp y = a.getOne();
    293            Dfp x = a.getOne();
    294            Dfp fact = a.getOne();
    295            Dfp py = new Dfp(y);
    296
    297            for (int i = 1; i < 90; i++) {
    298                x = x.multiply(a);
    299                fact = fact.divide(i);
    300                y = y.add(x.multiply(fact));
    301                if (y.equals(py)) {
    302                    break;
    303                }
    304                py = new Dfp(y);
    305            }
    306
    307            return y;
    308        }
    309
    310        /** Returns the natural logarithm of a.
    311         * a is first split into three parts such that  a = (10000^h)(2^j)k.
    312         * ln(a) is computed by ln(a) = ln(5)*h + ln(2)*(h+j) + ln(k)
    313         * k is in the range 2/3 < k <4/3 and is passed on to a series expansion.
    314         * @param a number from which logarithm is requested
    315         * @return log(a)
    316         */
    317        public static Dfp log(Dfp a) {
    318            int lr;
    319            Dfp x;
    320            int ix;
    321            int p2 = 0;
    322
    323            // Check the arguments somewhat here
    324            if (a.equals(a.getZero()) || a.lessThan(a.getZero()) || a.isNaN()) {
    325                // negative, zero or NaN
    326                a.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
    327                return a.dotrap(DfpField.FLAG_INVALID, "ln", a, a.newInstance((byte)1, Dfp.QNAN));
    328            }
    329
    330            if (a.classify() == Dfp.INFINITE) {
    331                return a;
    332            }
    333
    334            x = new Dfp(a);
    335            lr = x.log10K();
    336
    337            x = x.divide(pow(a.newInstance(10000), lr));  /* This puts x in the range 0-10000 */
    338            ix = x.floor().intValue();
    339
    340            while (ix > 2) {
    341                ix >>= 1;
    342                p2++;
    343            }
    344
    345
    346            Dfp[] spx = split(x);
    347            Dfp[] spy = new Dfp[2];
    348            spy[0] = pow(a.getTwo(), p2);          // use spy[0] temporarily as a divisor
    349            spx[0] = spx[0].divide(spy[0]);
    350            spx[1] = spx[1].divide(spy[0]);
    351
    352            spy[0] = a.newInstance("1.33333");    // Use spy[0] for comparison
    353            while (spx[0].add(spx[1]).greaterThan(spy[0])) {
    354                spx[0] = spx[0].divide(2);
    355                spx[1] = spx[1].divide(2);
    356                p2++;
    357            }
    358
    359            // X is now in the range of 2/3 < x < 4/3
    360            Dfp[] spz = logInternal(spx);
    361
    362            spx[0] = a.newInstance(new StringBuilder().append(p2+4*lr).toString());
    363            spx[1] = a.getZero();
    364            spy = splitMult(a.getField().getLn2Split(), spx);
    365
    366            spz[0] = spz[0].add(spy[0]);
    367            spz[1] = spz[1].add(spy[1]);
    368
    369            spx[0] = a.newInstance(new StringBuilder().append(4*lr).toString());
    370            spx[1] = a.getZero();
    371            spy = splitMult(a.getField().getLn5Split(), spx);
    372
    373            spz[0] = spz[0].add(spy[0]);
    374            spz[1] = spz[1].add(spy[1]);
    375
    376            return a.newInstance(spz[0].add(spz[1]));
    377
    378        }
    379
    380        /** Computes the natural log of a number between 0 and 2.
    381         *  Let f(x) = ln(x),
    382         *
    383         *  We know that f'(x) = 1/x, thus from Taylor's theorum we have:
    384         *
    385         *           -----          n+1         n
    386         *  f(x) =   \           (-1)    (x - 1)
    387         *           /          ----------------    for 1 <= n <= infinity
    388         *           -----             n
    389         *
    390         *  or
    391         *                       2        3       4
    392         *                   (x-1)   (x-1)    (x-1)
    393         *  ln(x) =  (x-1) - ----- + ------ - ------ + ...
    394         *                     2       3        4
    395         *
    396         *  alternatively,
    397         *
    398         *                  2    3   4
    399         *                 x    x   x
    400         *  ln(x+1) =  x - -  + - - - + ...
    401         *                 2    3   4
    402         *
    403         *  This series can be used to compute ln(x), but it converges too slowly.
    404         *
    405         *  If we substitute -x for x above, we get
    406         *
    407         *                   2    3    4
    408         *                  x    x    x
    409         *  ln(1-x) =  -x - -  - -  - - + ...
    410         *                  2    3    4
    411         *
    412         *  Note that all terms are now negative.  Because the even powered ones
    413         *  absorbed the sign.  Now, subtract the series above from the previous
    414         *  one to get ln(x+1) - ln(1-x).  Note the even terms cancel out leaving
    415         *  only the odd ones
    416         *
    417         *                             3     5      7
    418         *                           2x    2x     2x
    419         *  ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ...
    420         *                            3     5      7
    421         *
    422         *  By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have:
    423         *
    424         *                                3        5        7
    425         *      x+1           /          x        x        x          \
    426         *  ln ----- =   2 *  |  x  +   ----  +  ----  +  ---- + ...  |
    427         *      x-1           \          3        5        7          /
    428         *
    429         *  But now we want to find ln(a), so we need to find the value of x
    430         *  such that a = (x+1)/(x-1).   This is easily solved to find that
    431         *  x = (a-1)/(a+1).
    432         * @param a number from which logarithm is requested, in split form
    433         * @return log(a)
    434         */
    435        protected static Dfp[] logInternal(final Dfp a[]) {
    436
    437            /* Now we want to compute x = (a-1)/(a+1) but this is prone to
    438             * loss of precision.  So instead, compute x = (a/4 - 1/4) / (a/4 + 1/4)
    439             */
    440            Dfp t = a[0].divide(4).add(a[1].divide(4));
    441            Dfp x = t.add(a[0].newInstance("-0.25")).divide(t.add(a[0].newInstance("0.25")));
    442
    443            Dfp y = new Dfp(x);
    444            Dfp num = new Dfp(x);
    445            Dfp py = new Dfp(y);
    446            int den = 1;
    447            for (int i = 0; i < 10000; i++) {
    448                num = num.multiply(x);
    449                num = num.multiply(x);
    450                den = den + 2;
    451                t = num.divide(den);
    452                y = y.add(t);
    453                if (y.equals(py)) {
    454                    break;
    455                }
    456                py = new Dfp(y);
    457            }
    458
    459            y = y.multiply(a[0].getTwo());
    460
    461            return split(y);
    462
    463        }
    464
    465        /** Computes x to the y power.<p>
    466         *
    467         *  Uses the following method:<p>
    468         *
    469         *  <ol>
    470         *  <li> Set u = rint(y), v = y-u
    471         *  <li> Compute a = v * ln(x)
    472         *  <li> Compute b = rint( a/ln(2) )
    473         *  <li> Compute c = a - b*ln(2)
    474         *  <li> x<sup>y</sup> = x<sup>u</sup>  *   2<sup>b</sup> * e<sup>c</sup>
    475         *  </ol>
    476         *  if |y| > 1e8, then we compute by exp(y*ln(x))   <p>
    477         *
    478         *  <b>Special Cases</b><p>
    479         *  <ul>
    480         *  <li>  if y is 0.0 or -0.0 then result is 1.0
    481         *  <li>  if y is 1.0 then result is x
    482         *  <li>  if y is NaN then result is NaN
    483         *  <li>  if x is NaN and y is not zero then result is NaN
    484         *  <li>  if |x| > 1.0 and y is +Infinity then result is +Infinity
    485         *  <li>  if |x| < 1.0 and y is -Infinity then result is +Infinity
    486         *  <li>  if |x| > 1.0 and y is -Infinity then result is +0
    487         *  <li>  if |x| < 1.0 and y is +Infinity then result is +0
    488         *  <li>  if |x| = 1.0 and y is +/-Infinity then result is NaN
    489         *  <li>  if x = +0 and y > 0 then result is +0
    490         *  <li>  if x = +Inf and y < 0 then result is +0
    491         *  <li>  if x = +0 and y < 0 then result is +Inf
    492         *  <li>  if x = +Inf and y > 0 then result is +Inf
    493         *  <li>  if x = -0 and y > 0, finite, not odd integer then result is +0
    494         *  <li>  if x = -0 and y < 0, finite, and odd integer then result is -Inf
    495         *  <li>  if x = -Inf and y > 0, finite, and odd integer then result is -Inf
    496         *  <li>  if x = -0 and y < 0, not finite odd integer then result is +Inf
    497         *  <li>  if x = -Inf and y > 0, not finite odd integer then result is +Inf
    498         *  <li>  if x < 0 and y > 0, finite, and odd integer then result is -(|x|<sup>y</sup>)
    499         *  <li>  if x < 0 and y > 0, finite, and not integer then result is NaN
    500         *  </ul>
    501         *  @param x base to be raised
    502         *  @param y power to which base should be raised
    503         *  @return x<sup>y</sup>
    504         */
    505        public static Dfp pow(Dfp x, final Dfp y) {
    506
    507            // make sure we don't mix number with different precision
    508            if (x.getField().getRadixDigits() != y.getField().getRadixDigits()) {
    509                x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
    510                final Dfp result = x.newInstance(x.getZero());
    511                result.nans = Dfp.QNAN;
    512                return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, result);
    513            }
    514
    515            final Dfp zero = x.getZero();
    516            final Dfp one  = x.getOne();
    517            final Dfp two  = x.getTwo();
    518            boolean invert = false;
    519            int ui;
    520
    521            /* Check for special cases */
    522            if (y.equals(zero)) {
    523                return x.newInstance(one);
    524            }
    525
    526            if (y.equals(one)) {
    527                if (x.isNaN()) {
    528                    // Test for NaNs
    529                    x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
    530                    return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x);
    531                }
    532                return x;
    533            }
    534
    535            if (x.isNaN() || y.isNaN()) {
    536                // Test for NaNs
    537                x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
    538                return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
    539            }
    540
    541            // X == 0
    542            if (x.equals(zero)) {
    543                if (Dfp.copysign(one, x).greaterThan(zero)) {
    544                    // X == +0
    545                    if (y.greaterThan(zero)) {
    546                        return x.newInstance(zero);
    547                    } else {
    548                        return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
    549                    }
    550                } else {
    551                    // X == -0
    552                    if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
    553                        // If y is odd integer
    554                        if (y.greaterThan(zero)) {
    555                            return x.newInstance(zero.negate());
    556                        } else {
    557                            return x.newInstance(x.newInstance((byte)-1, Dfp.INFINITE));
    558                        }
    559                    } else {
    560                        // Y is not odd integer
    561                        if (y.greaterThan(zero)) {
    562                            return x.newInstance(zero);
    563                        } else {
    564                            return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
    565                        }
    566                    }
    567                }
    568            }
    569
    570            if (x.lessThan(zero)) {
    571                // Make x positive, but keep track of it
    572                x = x.negate();
    573                invert = true;
    574            }
    575
    576            if (x.greaterThan(one) && y.classify() == Dfp.INFINITE) {
    577                if (y.greaterThan(zero)) {
    578                    return y;
    579                } else {
    580                    return x.newInstance(zero);
    581                }
    582            }
    583
    584            if (x.lessThan(one) && y.classify() == Dfp.INFINITE) {
    585                if (y.greaterThan(zero)) {
    586                    return x.newInstance(zero);
    587                } else {
    588                    return x.newInstance(Dfp.copysign(y, one));
    589                }
    590            }
    591
    592            if (x.equals(one) && y.classify() == Dfp.INFINITE) {
    593                x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
    594                return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
    595            }
    596
    597            if (x.classify() == Dfp.INFINITE) {
    598                // x = +/- inf
    599                if (invert) {
    600                    // negative infinity
    601                    if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
    602                        // If y is odd integer
    603                        if (y.greaterThan(zero)) {
    604                            return x.newInstance(x.newInstance((byte)-1, Dfp.INFINITE));
    605                        } else {
    606                            return x.newInstance(zero.negate());
    607                        }
    608                    } else {
    609                        // Y is not odd integer
    610                        if (y.greaterThan(zero)) {
    611                            return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
    612                        } else {
    613                            return x.newInstance(zero);
    614                        }
    615                    }
    616                } else {
    617                    // positive infinity
    618                    if (y.greaterThan(zero)) {
    619                        return x;
    620                    } else {
    621                        return x.newInstance(zero);
    622                    }
    623                }
    624            }
    625
    626            if (invert && !y.rint().equals(y)) {
    627                x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
    628                return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
    629            }
    630
    631            // End special cases
    632
    633            Dfp r;
    634            if (y.lessThan(x.newInstance(100000000)) && y.greaterThan(x.newInstance(-100000000))) {
    635                final Dfp u = y.rint();
    636                ui = u.intValue();
    637
    638                final Dfp v = y.subtract(u);
    639
    640                if (v.unequal(zero)) {
    641                    final Dfp a = v.multiply(log(x));
    642                    final Dfp b = a.divide(x.getField().getLn2()).rint();
    643
    644                    final Dfp c = a.subtract(b.multiply(x.getField().getLn2()));
    645                    r = splitPow(split(x), ui);
    646                    r = r.multiply(pow(two, b.intValue()));
    647                    r = r.multiply(exp(c));
    648                } else {
    649                    r = splitPow(split(x), ui);
    650                }
    651            } else {
    652                // very large exponent.  |y| > 1e8
    653                r = exp(log(x).multiply(y));
    654            }
    655
    656            if (invert) {
    657                // if y is odd integer
    658                if (y.rint().equals(y) && !y.remainder(two).equals(zero)) {
    659                    r = r.negate();
    660                }
    661            }
    662
    663            return x.newInstance(r);
    664
    665        }
    666
    667        /** Computes sin(a)  Used when 0 < a < pi/4.
    668         * Uses the classic Taylor series.  x - x**3/3! + x**5/5!  ...
    669         * @param a number from which sine is desired, in split form
    670         * @return sin(a)
    671         */
    672        protected static Dfp sinInternal(Dfp a[]) {
    673
    674            Dfp c = a[0].add(a[1]);
    675            Dfp y = c;
    676            c = c.multiply(c);
    677            Dfp x = y;
    678            Dfp fact = a[0].getOne();
    679            Dfp py = new Dfp(y);
    680
    681            for (int i = 3; i < 90; i += 2) {
    682                x = x.multiply(c);
    683                x = x.negate();
    684
    685                fact = fact.divide((i-1)*i);  // 1 over fact
    686                y = y.add(x.multiply(fact));
    687                if (y.equals(py))
    688                    break;
    689                py = new Dfp(y);
    690            }
    691
    692            return y;
    693
    694        }
    695
    696        /** Computes cos(a)  Used when 0 < a < pi/4.
    697         * Uses the classic Taylor series for cosine.  1 - x**2/2! + x**4/4!  ...
    698         * @param a number from which cosine is desired, in split form
    699         * @return cos(a)
    700         */
    701        protected static Dfp cosInternal(Dfp a[]) {
    702            final Dfp one = a[0].getOne();
    703
    704
    705            Dfp x = one;
    706            Dfp y = one;
    707            Dfp c = a[0].add(a[1]);
    708            c = c.multiply(c);
    709
    710            Dfp fact = one;
    711            Dfp py = new Dfp(y);
    712
    713            for (int i = 2; i < 90; i += 2) {
    714                x = x.multiply(c);
    715                x = x.negate();
    716
    717                fact = fact.divide((i - 1) * i);  // 1 over fact
    718
    719                y = y.add(x.multiply(fact));
    720                if (y.equals(py)) {
    721                    break;
    722                }
    723                py = new Dfp(y);
    724            }
    725
    726            return y;
    727
    728        }
    729
    730        /** computes the sine of the argument.
    731         * @param a number from which sine is desired
    732         * @return sin(a)
    733         */
    734        public static Dfp sin(final Dfp a) {
    735            final Dfp pi = a.getField().getPi();
    736            final Dfp zero = a.getField().getZero();
    737            boolean neg = false;
    738
    739            /* First reduce the argument to the range of +/- PI */
    740            Dfp x = a.remainder(pi.multiply(2));
    741
    742            /* if x < 0 then apply identity sin(-x) = -sin(x) */
    743            /* This puts x in the range 0 < x < PI            */
    744            if (x.lessThan(zero)) {
    745                x = x.negate();
    746                neg = true;
    747            }
    748
    749            /* Since sine(x) = sine(pi - x) we can reduce the range to
    750             * 0 < x < pi/2
    751             */
    752
    753            if (x.greaterThan(pi.divide(2))) {
    754                x = pi.subtract(x);
    755            }
    756
    757            Dfp y;
    758            if (x.lessThan(pi.divide(4))) {
    759                Dfp c[] = new Dfp[2];
    760                c[0] = x;
    761                c[1] = zero;
    762
    763                //y = sinInternal(c);
    764                y = sinInternal(split(x));
    765            } else {
    766                final Dfp c[] = new Dfp[2];
    767                final Dfp[] piSplit = a.getField().getPiSplit();
    768                c[0] = piSplit[0].divide(2).subtract(x);
    769                c[1] = piSplit[1].divide(2);
    770                y = cosInternal(c);
    771            }
    772
    773            if (neg) {
    774                y = y.negate();
    775            }
    776
    777            return a.newInstance(y);
    778
    779        }
    780
    781        /** computes the cosine of the argument.
    782         * @param a number from which cosine is desired
    783         * @return cos(a)
    784         */
    785        public static Dfp cos(Dfp a) {
    786            final Dfp pi = a.getField().getPi();
    787            final Dfp zero = a.getField().getZero();
    788            boolean neg = false;
    789
    790            /* First reduce the argument to the range of +/- PI */
    791            Dfp x = a.remainder(pi.multiply(2));
    792
    793            /* if x < 0 then apply identity cos(-x) = cos(x) */
    794            /* This puts x in the range 0 < x < PI           */
    795            if (x.lessThan(zero)) {
    796                x = x.negate();
    797            }
    798
    799            /* Since cos(x) = -cos(pi - x) we can reduce the range to
    800             * 0 < x < pi/2
    801             */
    802
    803            if (x.greaterThan(pi.divide(2))) {
    804                x = pi.subtract(x);
    805                neg = true;
    806            }
    807
    808            Dfp y;
    809            if (x.lessThan(pi.divide(4))) {
    810                Dfp c[] = new Dfp[2];
    811                c[0] = x;
    812                c[1] = zero;
    813
    814                y = cosInternal(c);
    815            } else {
    816                final Dfp c[] = new Dfp[2];
    817                final Dfp[] piSplit = a.getField().getPiSplit();
    818                c[0] = piSplit[0].divide(2).subtract(x);
    819                c[1] = piSplit[1].divide(2);
    820                y = sinInternal(c);
    821            }
    822
    823            if (neg) {
    824                y = y.negate();
    825            }
    826
    827            return a.newInstance(y);
    828
    829        }
    830
    831        /** computes the tangent of the argument.
    832         * @param a number from which tangent is desired
    833         * @return tan(a)
    834         */
    835        public static Dfp tan(final Dfp a) {
    836            return sin(a).divide(cos(a));
    837        }
    838
    839        /** computes the arc-tangent of the argument.
    840         * @param a number from which arc-tangent is desired
    841         * @return atan(a)
    842         */
    843        protected static Dfp atanInternal(final Dfp a) {
    844
    845            Dfp y = new Dfp(a);
    846            Dfp x = new Dfp(y);
    847            Dfp py = new Dfp(y);
    848
    849            for (int i = 3; i < 90; i += 2) {
    850                x = x.multiply(a);
    851                x = x.multiply(a);
    852                x = x.negate();
    853                y = y.add(x.divide(i));
    854                if (y.equals(py)) {
    855                    break;
    856                }
    857                py = new Dfp(y);
    858            }
    859
    860            return y;
    861
    862        }
    863
    864        /** computes the arc tangent of the argument
    865         *
    866         *  Uses the typical taylor series
    867         *
    868         *  but may reduce arguments using the following identity
    869         * tan(x+y) = (tan(x) + tan(y)) / (1 - tan(x)*tan(y))
    870         *
    871         * since tan(PI/8) = sqrt(2)-1,
    872         *
    873         * atan(x) = atan( (x - sqrt(2) + 1) / (1+x*sqrt(2) - x) + PI/8.0
    874         * @param a number from which arc-tangent is desired
    875         * @return atan(a)
    876         */
    877        public static Dfp atan(final Dfp a) {
    878            final Dfp   zero      = a.getField().getZero();
    879            final Dfp   one       = a.getField().getOne();
    880            final Dfp[] sqr2Split = a.getField().getSqr2Split();
    881            final Dfp[] piSplit   = a.getField().getPiSplit();
    882            boolean recp = false;
    883            boolean neg = false;
    884            boolean sub = false;
    885
    886            final Dfp ty = sqr2Split[0].subtract(one).add(sqr2Split[1]);
    887
    888            Dfp x = new Dfp(a);
    889            if (x.lessThan(zero)) {
    890                neg = true;
    891                x = x.negate();
    892            }
    893
    894            if (x.greaterThan(one)) {
    895                recp = true;
    896                x = one.divide(x);
    897            }
    898
    899            if (x.greaterThan(ty)) {
    900                Dfp sty[] = new Dfp[2];
    901                sub = true;
    902
    903                sty[0] = sqr2Split[0].subtract(one);
    904                sty[1] = sqr2Split[1];
    905
    906                Dfp[] xs = split(x);
    907
    908                Dfp[] ds = splitMult(xs, sty);
    909                ds[0] = ds[0].add(one);
    910
    911                xs[0] = xs[0].subtract(sty[0]);
    912                xs[1] = xs[1].subtract(sty[1]);
    913
    914                xs = splitDiv(xs, ds);
    915                x = xs[0].add(xs[1]);
    916
    917                //x = x.subtract(ty).divide(dfp.one.add(x.multiply(ty)));
    918            }
    919
    920            Dfp y = atanInternal(x);
    921
    922            if (sub) {
    923                y = y.add(piSplit[0].divide(8)).add(piSplit[1].divide(8));
    924            }
    925
    926            if (recp) {
    927                y = piSplit[0].divide(2).subtract(y).add(piSplit[1].divide(2));
    928            }
    929
    930            if (neg) {
    931                y = y.negate();
    932            }
    933
    934            return a.newInstance(y);
    935
    936        }
    937
    938        /** computes the arc-sine of the argument.
    939         * @param a number from which arc-sine is desired
    940         * @return asin(a)
    941         */
    942        public static Dfp asin(final Dfp a) {
    943            return atan(a.divide(a.getOne().subtract(a.multiply(a)).sqrt()));
    944        }
    945
    946        /** computes the arc-cosine of the argument.
    947         * @param a number from which arc-cosine is desired
    948         * @return acos(a)
    949         */
    950        public static Dfp acos(Dfp a) {
    951            Dfp result;
    952            boolean negative = false;
    953
    954            if (a.lessThan(a.getZero())) {
    955                negative = true;
    956            }
    957
    958            a = Dfp.copysign(a, a.getOne());  // absolute value
    959
    960            result = atan(a.getOne().subtract(a.multiply(a)).sqrt().divide(a));
    961
    962            if (negative) {
    963                result = a.getField().getPi().subtract(result);
    964            }
    965
    966            return a.newInstance(result);
    967        }
    968
    969    }