DfpMath.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. package org.apache.commons.math4.legacy.core.dfp;

  20. /** Mathematical routines for use with {@link Dfp}.
  21.  * The constants are defined in {@link DfpField}
  22.  * @since 2.2
  23.  */
  24. public final class DfpMath {

  26.     /** Name for traps triggered by pow. */
  27.     private static final String POW_TRAP = "pow";

  29.     /**
  30.      * Private Constructor.
  31.      */
  32.     private DfpMath() {
  33.     }

  35.     /** Breaks a string representation up into two dfp's.
  36.      * <p>The two dfp are such that the sum of them is equivalent
  37.      * to the input string, but has higher precision than using a
  38.      * single dfp. This is useful for improving accuracy of
  39.      * exponentiation and critical multiplies.
  40.      * @param field field to which the Dfp must belong
  41.      * @param a string representation to split
  42.      * @return an array of two {@link Dfp} which sum is a
  43.      */
  44.     protected static Dfp[] split(final DfpField field, final String a) {
  45.         Dfp[] result = new Dfp[2];
  46.         char[] buf;
  47.         boolean leading = true;
  48.         int sp = 0;
  49.         int sig = 0;

  51.         buf = new char[a.length()];

  53.         for (int i = 0; i < buf.length; i++) {
  54.             buf[i] = a.charAt(i);

  56.             if (buf[i] >= '1' && buf[i] <= '9') {
  57.                 leading = false;
  58.             }

  60.             if (buf[i] == '.') {
  61.                 sig += (400 - sig) % 4;
  62.                 leading = false;
  63.             }

  65.             if (sig == (field.getRadixDigits() / 2) * 4) {
  66.                 sp = i;
  67.                 break;
  68.             }

  70.             if (buf[i] >= '0' && buf[i] <= '9' && !leading) {
  71.                 sig++;
  72.             }
  73.         }

  75.         result[0] = field.newDfp(String.valueOf(buf, 0, sp));

  77.         for (int i = 0; i < buf.length; i++) {
  78.             buf[i] = a.charAt(i);
  79.             if (buf[i] >= '0' && buf[i] <= '9' && i < sp) {
  80.                 buf[i] = '0';
  81.             }
  82.         }

  84.         result[1] = field.newDfp(String.valueOf(buf));

  86.         return result;
  87.     }

  89.     /** Splits a {@link Dfp} into 2 {@link Dfp}'s such that their sum is equal to the input {@link Dfp}.
  90.      * @param a number to split
  91.      * @return two elements array containing the split number
  92.      */
  93.     protected static Dfp[] split(final Dfp a) {
  94.         final Dfp[] result = new Dfp[2];
  95.         final Dfp shift = a.multiply(a.power10K(a.getRadixDigits() / 2));
  96.         result[0] = a.add(shift).subtract(shift);
  97.         result[1] = a.subtract(result[0]);
  98.         return result;
  99.     }

  101.     /** Multiply two numbers that are split in to two pieces that are
  102.      *  meant to be added together.
  103.      *  Use binomial multiplication so ab = a0 b0 + a0 b1 + a1 b0 + a1 b1
  104.      *  Store the first term in result0, the rest in result1
  105.      *  @param a first factor of the multiplication, in split form
  106.      *  @param b second factor of the multiplication, in split form
  107.      *  @return a &times; b, in split form
  108.      */
  109.     protected static Dfp[] splitMult(final Dfp[] a, final Dfp[] b) {
  110.         final Dfp[] result = new Dfp[2];

  112.         result[1] = a[0].getZero();
  113.         result[0] = a[0].multiply(b[0]);

  115.         /* If result[0] is infinite or zero, don't compute result[1].
  116.          * Attempting to do so may produce NaNs.
  117.          */

  119.         if (result[0].classify() == Dfp.INFINITE || result[0].equals(result[1])) {
  120.             return result;
  121.         }

  123.         result[1] = a[0].multiply(b[1]).add(a[1].multiply(b[0])).add(a[1].multiply(b[1]));

  125.         return result;
  126.     }

  128.     /** Divide two numbers that are split in to two pieces that are meant to be added together.
  129.      * Inverse of split multiply above:
  130.      *  (a+b) / (c+d) = (a/c) + ( (bc-ad)/(c**2+cd) )
  131.      *  @param a dividend, in split form
  132.      *  @param b divisor, in split form
  133.      *  @return a / b, in split form
  134.      */
  135.     protected static Dfp[] splitDiv(final Dfp[] a, final Dfp[] b) {
  136.         final Dfp[] result;

  138.         result = new Dfp[2];

  140.         result[0] = a[0].divide(b[0]);
  141.         result[1] = a[1].multiply(b[0]).subtract(a[0].multiply(b[1]));
  142.         result[1] = result[1].divide(b[0].multiply(b[0]).add(b[0].multiply(b[1])));

  144.         return result;
  145.     }

  147.     /** Raise a split base to the a power.
  148.      * @param base number to raise
  149.      * @param a power
  150.      * @return base<sup>a</sup>
  151.      */
  152.     protected static Dfp splitPow(final Dfp[] base, int a) {
  153.         boolean invert = false;

  155.         Dfp[] r = new Dfp[2];

  157.         Dfp[] result = new Dfp[2];
  158.         result[0] = base[0].getOne();
  159.         result[1] = base[0].getZero();

  161.         if (a == 0) {
  162.             // Special case a = 0
  163.             return result[0].add(result[1]);
  164.         }

  166.         if (a < 0) {
  167.             // If a is less than zero
  168.             invert = true;
  169.             a = -a;
  170.         }

  172.         // Exponentiate by successive squaring
  173.         do {
  174.             r[0] = new Dfp(base[0]);
  175.             r[1] = new Dfp(base[1]);
  176.             int trial = 1;

  178.             int prevtrial;
  179.             while (true) {
  180.                 prevtrial = trial;
  181.                 trial *= 2;
  182.                 if (trial > a) {
  183.                     break;
  184.                 }
  185.                 r = splitMult(r, r);
  186.             }

  188.             trial = prevtrial;

  190.             a -= trial;
  191.             result = splitMult(result, r);
  192.         } while (a >= 1);

  194.         result[0] = result[0].add(result[1]);

  196.         if (invert) {
  197.             result[0] = base[0].getOne().divide(result[0]);
  198.         }

  200.         return result[0];
  201.     }

  203.     /** Raises base to the power a by successive squaring.
  204.      * @param base number to raise
  205.      * @param a power
  206.      * @return base<sup>a</sup>
  207.      */
  208.     public static Dfp pow(Dfp base, int a) {
  209.         boolean invert = false;

  211.         Dfp result = base.getOne();

  213.         if (a == 0) {
  214.             // Special case
  215.             return result;
  216.         }

  218.         if (a < 0) {
  219.             invert = true;
  220.             a = -a;
  221.         }

  223.         // Exponentiate by successive squaring
  224.         do {
  225.             Dfp r = new Dfp(base);
  226.             Dfp prevr;
  227.             int trial = 1;
  228.             int prevtrial;

  230.             do {
  231.                 prevr = new Dfp(r);
  232.                 prevtrial = trial;
  233.                 r = r.multiply(r);
  234.                 trial *= 2;
  235.             } while (a > trial);

  237.             r = prevr;
  238.             trial = prevtrial;

  240.             a -= trial;
  241.             result = result.multiply(r);
  242.         } while (a >= 1);

  244.         if (invert) {
  245.             result = base.getOne().divide(result);
  246.         }

  248.         return base.newInstance(result);
  249.     }

  251.     /** Computes e to the given power.
  252.      * a is broken into two parts, such that a = n+m  where n is an integer.
  253.      * We use pow() to compute e<sup>n</sup> and a Taylor series to compute
  254.      * e<sup>m</sup>.  We return e*<sup>n</sup> &times; e<sup>m</sup>
  255.      * @param a power at which e should be raised
  256.      * @return e<sup>a</sup>
  257.      */
  258.     public static Dfp exp(final Dfp a) {

  260.         final Dfp inta = a.rint();
  261.         final Dfp fraca = a.subtract(inta);

  263.         final int ia = inta.intValue();
  264.         if (ia > 2147483646) {
  265.             // return +Infinity
  266.             return a.newInstance((byte)1, Dfp.INFINITE);
  267.         }

  269.         if (ia < -2147483646) {
  270.             // return 0;
  271.             return a.newInstance();
  272.         }

  274.         final Dfp einta = splitPow(a.getField().getESplit(), ia);
  275.         final Dfp efraca = expInternal(fraca);

  277.         return einta.multiply(efraca);
  278.     }

  280.     /** Computes e to the given power.
  281.      * Where {@code -1 < a < 1}.  Use the classic Taylor series.
  282.      * {@code 1 + x**2/2! + x**3/3! + x**4/4!  ... }
  283.      * @param a power at which e should be raised
  284.      * @return e<sup>a</sup>
  285.      */
  286.     protected static Dfp expInternal(final Dfp a) {
  287.         Dfp y = a.getOne();
  288.         Dfp x = a.getOne();
  289.         Dfp fact = a.getOne();
  290.         Dfp py = new Dfp(y);

  292.         for (int i = 1; i < 90; i++) {
  293.             x = x.multiply(a);
  294.             fact = fact.divide(i);
  295.             y = y.add(x.multiply(fact));
  296.             if (y.equals(py)) {
  297.                 break;
  298.             }
  299.             py = new Dfp(y);
  300.         }

  302.         return y;
  303.     }

  305.     /** Returns the natural logarithm of a.
  306.      * a is first split into three parts such that {@code a = (10000^h)(2^j)k}.
  307.      * ln(a) is computed by {@code ln(a) = ln(5)*4h + ln(2)*(4h+j) + ln(k)}.
  308.      * k is in the range {@code 2/3 < k <4/3} and is passed on to a series expansion.
  309.      * @param a number from which logarithm is requested
  310.      * @return log(a)
  311.      */
  312.     public static Dfp log(Dfp a) {
  313.         int lr;
  314.         Dfp x;
  315.         int ix;
  316.         int p2 = 0;

  318.         // Check the arguments somewhat here
  319.         if (a.equals(a.getZero()) || a.lessThan(a.getZero()) || a.isNaN()) {
  320.             // negative, zero or NaN
  321.             a.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
  322.             return a.dotrap(DfpField.FLAG_INVALID, "ln", a, a.newInstance((byte)1, Dfp.QNAN));
  323.         }

  325.         if (a.classify() == Dfp.INFINITE) {
  326.             return a;
  327.         }

  329.         x = new Dfp(a);
  330.         lr = x.log10K();

  332.         x = x.divide(pow(a.newInstance(10000), lr));  /* This puts x in the range 0-10000 */
  333.         ix = x.floor().intValue();

  335.         while (ix > 2) {
  336.             ix >>= 1;
  337.             p2++;
  338.         }


  341.         Dfp[] spx = split(x);
  342.         Dfp[] spy = new Dfp[2];
  343.         spy[0] = pow(a.getTwo(), p2);          // use spy[0] temporarily as a divisor
  344.         spx[0] = spx[0].divide(spy[0]);
  345.         spx[1] = spx[1].divide(spy[0]);

  347.         spy[0] = a.newInstance("1.33333");    // Use spy[0] for comparison
  348.         while (spx[0].add(spx[1]).greaterThan(spy[0])) {
  349.             spx[0] = spx[0].divide(2);
  350.             spx[1] = spx[1].divide(2);
  351.             p2++;
  352.         }

  354.         // X is now in the range of 2/3 < x < 4/3
  355.         Dfp[] spz = logInternal(spx);

  357.         spx[0] = a.newInstance(new StringBuilder().append(p2 + 4 * lr).toString());
  358.         spx[1] = a.getZero();
  359.         spy = splitMult(a.getField().getLn2Split(), spx);

  361.         spz[0] = spz[0].add(spy[0]);
  362.         spz[1] = spz[1].add(spy[1]);

  364.         spx[0] = a.newInstance(new StringBuilder().append(4 * lr).toString());
  365.         spx[1] = a.getZero();
  366.         spy = splitMult(a.getField().getLn5Split(), spx);

  368.         spz[0] = spz[0].add(spy[0]);
  369.         spz[1] = spz[1].add(spy[1]);

  371.         return a.newInstance(spz[0].add(spz[1]));
  372.     }

  374.     /** Computes the natural log of a number between 0 and 2.
  375.      * <pre>{@code
  376.      *  Let f(x) = ln(x),
  377.      *
  378.      *  We know that f'(x) = 1/x, thus from Taylor's theorum we have:
  379.      *
  380.      *           -----          n+1         n
  381.      *  f(x) =   \           (-1)    (x - 1)
  382.      *           /          ----------------    for 1 <= n <= infinity
  383.      *           -----             n
  384.      *
  385.      *  or
  386.      *                       2        3       4
  387.      *                   (x-1)   (x-1)    (x-1)
  388.      *  ln(x) =  (x-1) - ----- + ------ - ------ + ...
  389.      *                     2       3        4
  390.      *
  391.      *  alternatively,
  392.      *
  393.      *                  2    3   4
  394.      *                 x    x   x
  395.      *  ln(x+1) =  x - -  + - - - + ...
  396.      *                 2    3   4
  397.      *
  398.      *  This series can be used to compute ln(x), but it converges too slowly.
  399.      *
  400.      *  If we substitute -x for x above, we get
  401.      *
  402.      *                   2    3    4
  403.      *                  x    x    x
  404.      *  ln(1-x) =  -x - -  - -  - - + ...
  405.      *                  2    3    4
  406.      *
  407.      *  Note that all terms are now negative.  Because the even powered ones
  408.      *  absorbed the sign.  Now, subtract the series above from the previous
  409.      *  one to get ln(x+1) - ln(1-x).  Note the even terms cancel out leaving
  410.      *  only the odd ones
  411.      *
  412.      *                             3     5      7
  413.      *                           2x    2x     2x
  414.      *  ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ...
  415.      *                            3     5      7
  416.      *
  417.      *  By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have:
  418.      *
  419.      *                                3        5        7
  420.      *      x+1           /          x        x        x          \
  421.      *  ln ----- =   2 *  |  x  +   ----  +  ----  +  ---- + ...  |
  422.      *      x-1           \          3        5        7          /
  423.      *
  424.      *  But now we want to find ln(a), so we need to find the value of x
  425.      *  such that a = (x+1)/(x-1).   This is easily solved to find that
  426.      *  x = (a-1)/(a+1).
  427.      * }</pre>
  428.      * @param a number from which logarithm is requested, in split form
  429.      * @return log(a)
  430.      */
  431.     protected static Dfp[] logInternal(final Dfp[] a) {

  433.         /* Now we want to compute x = (a-1)/(a+1) but this is prone to
  434.          * loss of precision.  So instead, compute x = (a/4 - 1/4) / (a/4 + 1/4)
  435.          */
  436.         Dfp t = a[0].divide(4).add(a[1].divide(4));
  437.         Dfp x = t.add(a[0].newInstance("-0.25")).divide(t.add(a[0].newInstance("0.25")));

  439.         Dfp y = new Dfp(x);
  440.         Dfp num = new Dfp(x);
  441.         Dfp py = new Dfp(y);
  442.         int den = 1;
  443.         for (int i = 0; i < 10000; i++) {
  444.             num = num.multiply(x);
  445.             num = num.multiply(x);
  446.             den += 2;
  447.             t = num.divide(den);
  448.             y = y.add(t);
  449.             if (y.equals(py)) {
  450.                 break;
  451.             }
  452.             py = new Dfp(y);
  453.         }

  455.         y = y.multiply(a[0].getTwo());

  457.         return split(y);
  458.     }

  460.     /** Computes x to the y power.<p>
  461.      *
  462.      *  Uses the following method:
  463.      *
  464.      *  <ol>
  465.      *  <li> Set u = rint(y), v = y-u
  466.      *  <li> Compute a = v * ln(x)
  467.      *  <li> Compute b = rint( a/ln(2) )
  468.      *  <li> Compute c = a - b*ln(2)
  469.      *  <li> x<sup>y</sup> = x<sup>u</sup>  *   2<sup>b</sup> * e<sup>c</sup>
  470.      *  </ol>
  471.      *  if {@code |y| > 1e8}, then we compute by {@code exp(y*ln(x))}<p>
  472.      *
  473.      *  <b>Special Cases</b>
  474.      *  <ul>
  475.      *  <li>  if y is 0.0 or -0.0 then result is 1.0
  476.      *  <li>  if y is 1.0 then result is x
  477.      *  <li>  if y is NaN then result is NaN
  478.      *  <li>  if x is NaN and y is not zero then result is NaN
  479.      *  <li>  if {@code |x| > 1.0} and y is +Infinity then result is +Infinity
  480.      *  <li>  if {@code |x| < 1.0} and y is -Infinity then result is +Infinity
  481.      *  <li>  if {@code |x| > 1.0} and y is -Infinity then result is +0
  482.      *  <li>  if {@code |x| < 1.0} and y is +Infinity then result is +0
  483.      *  <li>  if {@code |x| = 1.0} and y is +/-Infinity then result is NaN
  484.      *  <li>  if {@code x = +0} and {@code y > 0} then result is +0
  485.      *  <li>  if {@code x = +Inf} and {@code y < 0} then result is +0
  486.      *  <li>  if {@code x = +0} and {@code y < 0} then result is +Inf
  487.      *  <li>  if {@code x = +Inf} and {@code y > 0} then result is +Inf
  488.      *  <li>  if {@code x = -0} and {@code y > 0}, finite, not odd integer then result is +0
  489.      *  <li>  if {@code x = -0} and {@code y < 0}, finite, and odd integer then result is -Inf
  490.      *  <li>  if {@code x = -Inf} and {@code y > 0}, finite, and odd integer then result is -Inf
  491.      *  <li>  if {@code x = -0} and {@code y < 0}, not finite odd integer then result is +Inf
  492.      *  <li>  if {@code x = -Inf} and {@code y > 0}, not finite odd integer then result is +Inf
  493.      *  <li>  if {@code x < 0} and {@code y > 0}, finite, and odd integer then result is -(|x|<sup>y</sup>)
  494.      *  <li>  if {@code x < 0} and {@code y > 0}, finite, and not integer then result is NaN
  495.      *  </ul>
  496.      *  @param x base to be raised
  497.      *  @param y power to which base should be raised
  498.      *  @return x<sup>y</sup>
  499.      */
  500.     public static Dfp pow(Dfp x, final Dfp y) {

  502.         // make sure we don't mix number with different precision
  503.         if (x.getField().getRadixDigits() != y.getField().getRadixDigits()) {
  504.             x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
  505.             final Dfp result = x.newInstance(x.getZero());
  506.             result.nans = Dfp.QNAN;
  507.             return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, result);
  508.         }

  510.         final Dfp zero = x.getZero();
  511.         final Dfp one  = x.getOne();
  512.         final Dfp two  = x.getTwo();
  513.         boolean invert = false;
  514.         int ui;

  516.         /* Check for special cases */
  517.         if (y.equals(zero)) {
  518.             return x.newInstance(one);
  519.         }

  521.         if (y.equals(one)) {
  522.             if (x.isNaN()) {
  523.                 // Test for NaNs
  524.                 x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
  525.                 return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x);
  526.             }
  527.             return x;
  528.         }

  530.         if (x.isNaN() || y.isNaN()) {
  531.             // Test for NaNs
  532.             x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
  533.             return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
  534.         }

  536.         // X == 0
  537.         if (x.equals(zero)) {
  538.             if (Dfp.copySign(one, x).greaterThan(zero)) {
  539.                 // X == +0
  540.                 if (y.greaterThan(zero)) {
  541.                     return x.newInstance(zero);
  542.                 } else {
  543.                     return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
  544.                 }
  545.             } else {
  546.                 // X == -0
  547.                 if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
  548.                     // If y is odd integer
  549.                     if (y.greaterThan(zero)) {
  550.                         return x.newInstance(zero.negate());
  551.                     } else {
  552.                         return x.newInstance(x.newInstance((byte)-1, Dfp.INFINITE));
  553.                     }
  554.                 } else {
  555.                     // Y is not odd integer
  556.                     if (y.greaterThan(zero)) {
  557.                         return x.newInstance(zero);
  558.                     } else {
  559.                         return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
  560.                     }
  561.                 }
  562.             }
  563.         }

  565.         if (x.lessThan(zero)) {
  566.             // Make x positive, but keep track of it
  567.             x = x.negate();
  568.             invert = true;
  569.         }

  571.         if (x.greaterThan(one) && y.classify() == Dfp.INFINITE) {
  572.             if (y.greaterThan(zero)) {
  573.                 return y;
  574.             } else {
  575.                 return x.newInstance(zero);
  576.             }
  577.         }

  579.         if (x.lessThan(one) && y.classify() == Dfp.INFINITE) {
  580.             if (y.greaterThan(zero)) {
  581.                 return x.newInstance(zero);
  582.             } else {
  583.                 return x.newInstance(Dfp.copySign(y, one));
  584.             }
  585.         }

  587.         if (x.equals(one) && y.classify() == Dfp.INFINITE) {
  588.             x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
  589.             return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
  590.         }

  592.         if (x.classify() == Dfp.INFINITE) {
  593.             // x = +/- inf
  594.             if (invert) {
  595.                 // negative infinity
  596.                 if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
  597.                     // If y is odd integer
  598.                     if (y.greaterThan(zero)) {
  599.                         return x.newInstance(x.newInstance((byte)-1, Dfp.INFINITE));
  600.                     } else {
  601.                         return x.newInstance(zero.negate());
  602.                     }
  603.                 } else {
  604.                     // Y is not odd integer
  605.                     if (y.greaterThan(zero)) {
  606.                         return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
  607.                     } else {
  608.                         return x.newInstance(zero);
  609.                     }
  610.                 }
  611.             } else {
  612.                 // positive infinity
  613.                 if (y.greaterThan(zero)) {
  614.                     return x;
  615.                 } else {
  616.                     return x.newInstance(zero);
  617.                 }
  618.             }
  619.         }

  621.         if (invert && !y.rint().equals(y)) {
  622.             x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
  623.             return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
  624.         }

  626.         // End special cases

  628.         Dfp r;
  629.         if (y.lessThan(x.newInstance(100000000)) && y.greaterThan(x.newInstance(-100000000))) {
  630.             final Dfp u = y.rint();
  631.             ui = u.intValue();

  633.             final Dfp v = y.subtract(u);

  635.             if (v.unequal(zero)) {
  636.                 final Dfp a = v.multiply(log(x));
  637.                 final Dfp b = a.divide(x.getField().getLn2()).rint();

  639.                 final Dfp c = a.subtract(b.multiply(x.getField().getLn2()));
  640.                 r = splitPow(split(x), ui);
  641.                 r = r.multiply(pow(two, b.intValue()));
  642.                 r = r.multiply(exp(c));
  643.             } else {
  644.                 r = splitPow(split(x), ui);
  645.             }
  646.         } else {
  647.             // very large exponent.  |y| > 1e8
  648.             r = exp(log(x).multiply(y));
  649.         }

  651.         if (invert && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
  652.             // if y is odd integer
  653.             r = r.negate();
  654.         }

  656.         return x.newInstance(r);
  657.     }

  659.     /** Computes sin(a)  Used when {@code {@code 0 < a < pi/4}}.
  660.      * Uses the classic Taylor series.  {@code x - x**3/3! + x**5/5!  ... }
  661.      * @param a number from which sine is desired, in split form
  662.      * @return sin(a)
  663.      */
  664.     protected static Dfp sinInternal(Dfp[] a) {

  666.         Dfp c = a[0].add(a[1]);
  667.         Dfp y = c;
  668.         c = c.multiply(c);
  669.         Dfp x = y;
  670.         Dfp fact = a[0].getOne();
  671.         Dfp py = new Dfp(y);

  673.         for (int i = 3; i < 90; i += 2) {
  674.             x = x.multiply(c);
  675.             x = x.negate();

  677.             fact = fact.divide((i - 1) * i); // 1 over fact
  678.             y = y.add(x.multiply(fact));
  679.             if (y.equals(py)) {
  680.                 break;
  681.             }
  682.             py = new Dfp(y);
  683.         }

  685.         return y;
  686.     }

  688.     /** Computes cos(a)  Used when {@code 0 < a < pi/4}.
  689.      * Uses the classic Taylor series for cosine.  1 - x**2/2! + x**4/4!  ...
  690.      * @param a number from which cosine is desired, in split form
  691.      * @return cos(a)
  692.      */
  693.     protected static Dfp cosInternal(Dfp[] a) {
  694.         final Dfp one = a[0].getOne();


  697.         Dfp x = one;
  698.         Dfp y = one;
  699.         Dfp c = a[0].add(a[1]);
  700.         c = c.multiply(c);

  702.         Dfp fact = one;
  703.         Dfp py = new Dfp(y);

  705.         for (int i = 2; i < 90; i += 2) {
  706.             x = x.multiply(c);
  707.             x = x.negate();

  709.             fact = fact.divide((i - 1) * i);  // 1 over fact

  711.             y = y.add(x.multiply(fact));
  712.             if (y.equals(py)) {
  713.                 break;
  714.             }
  715.             py = new Dfp(y);
  716.         }

  718.         return y;
  719.     }

  721.     /** computes the sine of the argument.
  722.      * @param a number from which sine is desired
  723.      * @return sin(a)
  724.      */
  725.     public static Dfp sin(final Dfp a) {
  726.         final Dfp pi = a.getField().getPi();
  727.         final Dfp zero = a.getField().getZero();
  728.         boolean neg = false;

  730.         /* First reduce the argument to the range of +/- PI */
  731.         Dfp x = a.remainder(pi.multiply(2));

  733.         /* if x < 0 then apply identity sin(-x) = -sin(x) */
  734.         /* This puts x in the range 0 < x < PI            */
  735.         if (x.lessThan(zero)) {
  736.             x = x.negate();
  737.             neg = true;
  738.         }

  740.         /* Since sine(x) = sine(pi - x) we can reduce the range to
  741.          * 0 < x < pi/2
  742.          */

  744.         if (x.greaterThan(pi.divide(2))) {
  745.             x = pi.subtract(x);
  746.         }

  748.         Dfp y;
  749.         if (x.lessThan(pi.divide(4))) {
  750.             y = sinInternal(split(x));
  751.         } else {
  752.             final Dfp[] c = new Dfp[2];
  753.             final Dfp[] piSplit = a.getField().getPiSplit();
  754.             c[0] = piSplit[0].divide(2).subtract(x);
  755.             c[1] = piSplit[1].divide(2);
  756.             y = cosInternal(c);
  757.         }

  759.         if (neg) {
  760.             y = y.negate();
  761.         }

  763.         return a.newInstance(y);
  764.     }

  766.     /** computes the cosine of the argument.
  767.      * @param a number from which cosine is desired
  768.      * @return cos(a)
  769.      */
  770.     public static Dfp cos(Dfp a) {
  771.         final Dfp pi = a.getField().getPi();
  772.         final Dfp zero = a.getField().getZero();
  773.         boolean neg = false;

  775.         /* First reduce the argument to the range of +/- PI */
  776.         Dfp x = a.remainder(pi.multiply(2));

  778.         /* if x < 0 then apply identity cos(-x) = cos(x) */
  779.         /* This puts x in the range 0 < x < PI           */
  780.         if (x.lessThan(zero)) {
  781.             x = x.negate();
  782.         }

  784.         /* Since cos(x) = -cos(pi - x) we can reduce the range to
  785.          * 0 < x < pi/2
  786.          */

  788.         if (x.greaterThan(pi.divide(2))) {
  789.             x = pi.subtract(x);
  790.             neg = true;
  791.         }

  793.         Dfp y;
  794.         if (x.lessThan(pi.divide(4))) {
  795.             Dfp[] c = new Dfp[2];
  796.             c[0] = x;
  797.             c[1] = zero;

  799.             y = cosInternal(c);
  800.         } else {
  801.             final Dfp[] c = new Dfp[2];
  802.             final Dfp[] piSplit = a.getField().getPiSplit();
  803.             c[0] = piSplit[0].divide(2).subtract(x);
  804.             c[1] = piSplit[1].divide(2);
  805.             y = sinInternal(c);
  806.         }

  808.         if (neg) {
  809.             y = y.negate();
  810.         }

  812.         return a.newInstance(y);
  813.     }

  815.     /** computes the tangent of the argument.
  816.      * @param a number from which tangent is desired
  817.      * @return tan(a)
  818.      */
  819.     public static Dfp tan(final Dfp a) {
  820.         return sin(a).divide(cos(a));
  821.     }

  823.     /** computes the arc-tangent of the argument.
  824.      * @param a number from which arc-tangent is desired
  825.      * @return atan(a)
  826.      */
  827.     protected static Dfp atanInternal(final Dfp a) {

  829.         Dfp y = new Dfp(a);
  830.         Dfp x = new Dfp(y);
  831.         Dfp py = new Dfp(y);

  833.         for (int i = 3; i < 90; i += 2) {
  834.             x = x.multiply(a);
  835.             x = x.multiply(a);
  836.             x = x.negate();
  837.             y = y.add(x.divide(i));
  838.             if (y.equals(py)) {
  839.                 break;
  840.             }
  841.             py = new Dfp(y);
  842.         }

  844.         return y;
  845.     }

  847.     /** Computes the arc tangent of the argument.
  848.      *
  849.      *  Uses the typical taylor series
  850.      *
  851.      *  but may reduce arguments using the following identity
  852.      * tan(x+y) = (tan(x) + tan(y)) / (1 - tan(x)*tan(y))
  853.      *
  854.      * since tan(PI/8) = sqrt(2)-1,
  855.      *
  856.      * atan(x) = atan( (x - sqrt(2) + 1) / (1+x*sqrt(2) - x) + PI/8.0
  857.      * @param a number from which arc-tangent is desired
  858.      * @return atan(a)
  859.      */
  860.     public static Dfp atan(final Dfp a) {
  861.         final Dfp   zero      = a.getField().getZero();
  862.         final Dfp   one       = a.getField().getOne();
  863.         final Dfp[] sqr2Split = a.getField().getSqr2Split();
  864.         final Dfp[] piSplit   = a.getField().getPiSplit();
  865.         boolean recp = false;
  866.         boolean neg = false;
  867.         boolean sub = false;

  869.         final Dfp ty = sqr2Split[0].subtract(one).add(sqr2Split[1]);

  871.         Dfp x = new Dfp(a);
  872.         if (x.lessThan(zero)) {
  873.             neg = true;
  874.             x = x.negate();
  875.         }

  877.         if (x.greaterThan(one)) {
  878.             recp = true;
  879.             x = one.divide(x);
  880.         }

  882.         if (x.greaterThan(ty)) {
  883.             Dfp[] sty = new Dfp[2];
  884.             sub = true;

  886.             sty[0] = sqr2Split[0].subtract(one);
  887.             sty[1] = sqr2Split[1];

  889.             Dfp[] xs = split(x);

  891.             Dfp[] ds = splitMult(xs, sty);
  892.             ds[0] = ds[0].add(one);

  894.             xs[0] = xs[0].subtract(sty[0]);
  895.             xs[1] = xs[1].subtract(sty[1]);

  897.             xs = splitDiv(xs, ds);
  898.             x = xs[0].add(xs[1]);

  900.             //x = x.subtract(ty).divide(dfp.one.add(x.multiply(ty)));
  901.         }

  903.         Dfp y = atanInternal(x);

  905.         if (sub) {
  906.             y = y.add(piSplit[0].divide(8)).add(piSplit[1].divide(8));
  907.         }

  909.         if (recp) {
  910.             y = piSplit[0].divide(2).subtract(y).add(piSplit[1].divide(2));
  911.         }

  913.         if (neg) {
  914.             y = y.negate();
  915.         }

  917.         return a.newInstance(y);
  918.     }

  920.     /** computes the arc-sine of the argument.
  921.      * @param a number from which arc-sine is desired
  922.      * @return asin(a)
  923.      */
  924.     public static Dfp asin(final Dfp a) {
  925.         return atan(a.divide(a.getOne().subtract(a.multiply(a)).sqrt()));
  926.     }

  928.     /** computes the arc-cosine of the argument.
  929.      * @param a number from which arc-cosine is desired
  930.      * @return acos(a)
  931.      */
  932.     public static Dfp acos(Dfp a) {
  933.         Dfp result;
  934.         boolean negative = a.lessThan(a.getZero());

  936.         a = Dfp.copySign(a, a.getOne());  // absolute value

  938.         result = atan(a.getOne().subtract(a.multiply(a)).sqrt().divide(a));

  940.         if (negative) {
  941.             result = a.getField().getPi().subtract(result);
  942.         }

  944.         return a.newInstance(result);
  945.     }
  946. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT