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

package org.apache.commons.math3.dfp;

Mathematical routines for use with Dfp. The constants are defined in DfpField
Since:2.2
/** Mathematical routines for use with {@link Dfp}. * The constants are defined in {@link DfpField} * @since 2.2 */
public class DfpMath {
Name for traps triggered by pow.
/** Name for traps triggered by pow. */
private static final String POW_TRAP = "pow";
Private Constructor.
/** * Private Constructor. */
private DfpMath() { }
Breaks a string representation up into two dfp's.

The two dfp are such that the sum of them is equivalent to the input string, but has higher precision than using a single dfp. This is useful for improving accuracy of exponentiation and critical multiplies.

Params:
  • field – field to which the Dfp must belong
  • a – string representation to split
Returns:an array of two Dfp which sum is a
/** Breaks a string representation up into two dfp's. * <p>The two dfp are such that the sum of them is equivalent * to the input string, but has higher precision than using a * single dfp. This is useful for improving accuracy of * exponentiation and critical multiplies. * @param field field to which the Dfp must belong * @param a string representation to split * @return an array of two {@link Dfp} which sum is a */
protected static Dfp[] split(final DfpField field, final String a) { Dfp result[] = new Dfp[2]; char[] buf; boolean leading = true; int sp = 0; int sig = 0; buf = new char[a.length()]; for (int i = 0; i < buf.length; i++) { buf[i] = a.charAt(i); if (buf[i] >= '1' && buf[i] <= '9') { leading = false; } if (buf[i] == '.') { sig += (400 - sig) % 4; leading = false; } if (sig == (field.getRadixDigits() / 2) * 4) { sp = i; break; } if (buf[i] >= '0' && buf[i] <= '9' && !leading) { sig ++; } } result[0] = field.newDfp(new String(buf, 0, sp)); for (int i = 0; i < buf.length; i++) { buf[i] = a.charAt(i); if (buf[i] >= '0' && buf[i] <= '9' && i < sp) { buf[i] = '0'; } } result[1] = field.newDfp(new String(buf)); return result; }
Splits a Dfp into 2 Dfp's such that their sum is equal to the input Dfp.
Params:
  • a – number to split
Returns:two elements array containing the split number
/** Splits a {@link Dfp} into 2 {@link Dfp}'s such that their sum is equal to the input {@link Dfp}. * @param a number to split * @return two elements array containing the split number */
protected static Dfp[] split(final Dfp a) { final Dfp[] result = new Dfp[2]; final Dfp shift = a.multiply(a.power10K(a.getRadixDigits() / 2)); result[0] = a.add(shift).subtract(shift); result[1] = a.subtract(result[0]); return result; }
Multiply two numbers that are split in to two pieces that are meant to be added together. Use binomial multiplication so ab = a0 b0 + a0 b1 + a1 b0 + a1 b1 Store the first term in result0, the rest in result1 @param a first factor of the multiplication, in split form @param b second factor of the multiplication, in split form @return a × b, in split form
/** Multiply two numbers that are split in to two pieces that are * meant to be added together. * Use binomial multiplication so ab = a0 b0 + a0 b1 + a1 b0 + a1 b1 * Store the first term in result0, the rest in result1 * @param a first factor of the multiplication, in split form * @param b second factor of the multiplication, in split form * @return a &times; b, in split form */
protected static Dfp[] splitMult(final Dfp[] a, final Dfp[] b) { final Dfp[] result = new Dfp[2]; result[1] = a[0].getZero(); result[0] = a[0].multiply(b[0]); /* If result[0] is infinite or zero, don't compute result[1]. * Attempting to do so may produce NaNs. */ if (result[0].classify() == Dfp.INFINITE || result[0].equals(result[1])) { return result; } result[1] = a[0].multiply(b[1]).add(a[1].multiply(b[0])).add(a[1].multiply(b[1])); return result; }
Divide two numbers that are split in to two pieces that are meant to be added together. Inverse of split multiply above: (a+b) / (c+d) = (a/c) + ( (bc-ad)/(c**2+cd) ) @param a dividend, in split form @param b divisor, in split form @return a / b, in split form
/** Divide two numbers that are split in to two pieces that are meant to be added together. * Inverse of split multiply above: * (a+b) / (c+d) = (a/c) + ( (bc-ad)/(c**2+cd) ) * @param a dividend, in split form * @param b divisor, in split form * @return a / b, in split form */
protected static Dfp[] splitDiv(final Dfp[] a, final Dfp[] b) { final Dfp[] result; result = new Dfp[2]; result[0] = a[0].divide(b[0]); result[1] = a[1].multiply(b[0]).subtract(a[0].multiply(b[1])); result[1] = result[1].divide(b[0].multiply(b[0]).add(b[0].multiply(b[1]))); return result; }
Raise a split base to the a power.
Params:
  • base – number to raise
  • a – power
Returns:basea
/** Raise a split base to the a power. * @param base number to raise * @param a power * @return base<sup>a</sup> */
protected static Dfp splitPow(final Dfp[] base, int a) { boolean invert = false; Dfp[] r = new Dfp[2]; Dfp[] result = new Dfp[2]; result[0] = base[0].getOne(); result[1] = base[0].getZero(); if (a == 0) { // Special case a = 0 return result[0].add(result[1]); } if (a < 0) { // If a is less than zero invert = true; a = -a; } // Exponentiate by successive squaring do { r[0] = new Dfp(base[0]); r[1] = new Dfp(base[1]); int trial = 1; int prevtrial; while (true) { prevtrial = trial; trial *= 2; if (trial > a) { break; } r = splitMult(r, r); } trial = prevtrial; a -= trial; result = splitMult(result, r); } while (a >= 1); result[0] = result[0].add(result[1]); if (invert) { result[0] = base[0].getOne().divide(result[0]); } return result[0]; }
Raises base to the power a by successive squaring.
Params:
  • base – number to raise
  • a – power
Returns:basea
/** Raises base to the power a by successive squaring. * @param base number to raise * @param a power * @return base<sup>a</sup> */
public static Dfp pow(Dfp base, int a) { boolean invert = false; Dfp result = base.getOne(); if (a == 0) { // Special case return result; } if (a < 0) { invert = true; a = -a; } // Exponentiate by successive squaring do { Dfp r = new Dfp(base); Dfp prevr; int trial = 1; int prevtrial; do { prevr = new Dfp(r); prevtrial = trial; r = r.multiply(r); trial *= 2; } while (a>trial); r = prevr; trial = prevtrial; a -= trial; result = result.multiply(r); } while (a >= 1); if (invert) { result = base.getOne().divide(result); } return base.newInstance(result); }
Computes e to the given power. a is broken into two parts, such that a = n+m where n is an integer. We use pow() to compute en and a Taylor series to compute em. We return e*n × em
Params:
  • a – power at which e should be raised
Returns:ea
/** Computes e to the given power. * a is broken into two parts, such that a = n+m where n is an integer. * We use pow() to compute e<sup>n</sup> and a Taylor series to compute * e<sup>m</sup>. We return e*<sup>n</sup> &times; e<sup>m</sup> * @param a power at which e should be raised * @return e<sup>a</sup> */
public static Dfp exp(final Dfp a) { final Dfp inta = a.rint(); final Dfp fraca = a.subtract(inta); final int ia = inta.intValue(); if (ia > 2147483646) { // return +Infinity return a.newInstance((byte)1, Dfp.INFINITE); } if (ia < -2147483646) { // return 0; return a.newInstance(); } final Dfp einta = splitPow(a.getField().getESplit(), ia); final Dfp efraca = expInternal(fraca); return einta.multiply(efraca); }
Computes e to the given power. Where -1 < a < 1. Use the classic Taylor series. 1 + x**2/2! + x**3/3! + x**4/4! ...
Params:
  • a – power at which e should be raised
Returns:ea
/** Computes e to the given power. * Where -1 < a < 1. Use the classic Taylor series. 1 + x**2/2! + x**3/3! + x**4/4! ... * @param a power at which e should be raised * @return e<sup>a</sup> */
protected static Dfp expInternal(final Dfp a) { Dfp y = a.getOne(); Dfp x = a.getOne(); Dfp fact = a.getOne(); Dfp py = new Dfp(y); for (int i = 1; i < 90; i++) { x = x.multiply(a); fact = fact.divide(i); y = y.add(x.multiply(fact)); if (y.equals(py)) { break; } py = new Dfp(y); } return y; }
Returns the natural logarithm of a. a is first split into three parts such that a = (10000^h)(2^j)k. ln(a) is computed by ln(a) = ln(5)*h + ln(2)*(h+j) + ln(k) k is in the range 2/3 < k <4/3 and is passed on to a series expansion.
Params:
  • a – number from which logarithm is requested
Returns:log(a)
/** Returns the natural logarithm of a. * a is first split into three parts such that a = (10000^h)(2^j)k. * ln(a) is computed by ln(a) = ln(5)*h + ln(2)*(h+j) + ln(k) * k is in the range 2/3 < k <4/3 and is passed on to a series expansion. * @param a number from which logarithm is requested * @return log(a) */
public static Dfp log(Dfp a) { int lr; Dfp x; int ix; int p2 = 0; // Check the arguments somewhat here if (a.equals(a.getZero()) || a.lessThan(a.getZero()) || a.isNaN()) { // negative, zero or NaN a.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); return a.dotrap(DfpField.FLAG_INVALID, "ln", a, a.newInstance((byte)1, Dfp.QNAN)); } if (a.classify() == Dfp.INFINITE) { return a; } x = new Dfp(a); lr = x.log10K(); x = x.divide(pow(a.newInstance(10000), lr)); /* This puts x in the range 0-10000 */ ix = x.floor().intValue(); while (ix > 2) { ix >>= 1; p2++; } Dfp[] spx = split(x); Dfp[] spy = new Dfp[2]; spy[0] = pow(a.getTwo(), p2); // use spy[0] temporarily as a divisor spx[0] = spx[0].divide(spy[0]); spx[1] = spx[1].divide(spy[0]); spy[0] = a.newInstance("1.33333"); // Use spy[0] for comparison while (spx[0].add(spx[1]).greaterThan(spy[0])) { spx[0] = spx[0].divide(2); spx[1] = spx[1].divide(2); p2++; } // X is now in the range of 2/3 < x < 4/3 Dfp[] spz = logInternal(spx); spx[0] = a.newInstance(new StringBuilder().append(p2+4*lr).toString()); spx[1] = a.getZero(); spy = splitMult(a.getField().getLn2Split(), spx); spz[0] = spz[0].add(spy[0]); spz[1] = spz[1].add(spy[1]); spx[0] = a.newInstance(new StringBuilder().append(4*lr).toString()); spx[1] = a.getZero(); spy = splitMult(a.getField().getLn5Split(), spx); spz[0] = spz[0].add(spy[0]); spz[1] = spz[1].add(spy[1]); return a.newInstance(spz[0].add(spz[1])); }
Computes the natural log of a number between 0 and 2. Let f(x) = ln(x), We know that f'(x) = 1/x, thus from Taylor's theorum we have: ----- n+1 n f(x) = \ (-1) (x - 1) / ---------------- for 1 <= n <= infinity ----- n or 2 3 4 (x-1) (x-1) (x-1) ln(x) = (x-1) - ----- + ------ - ------ + ... 2 3 4 alternatively, 2 3 4 x x x ln(x+1) = x - - + - - - + ... 2 3 4 This series can be used to compute ln(x), but it converges too slowly. If we substitute -x for x above, we get 2 3 4 x x x ln(1-x) = -x - - - - - - + ... 2 3 4 Note that all terms are now negative. Because the even powered ones absorbed the sign. Now, subtract the series above from the previous one to get ln(x+1) - ln(1-x). Note the even terms cancel out leaving only the odd ones 3 5 7 2x 2x 2x ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ... 3 5 7 By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have: 3 5 7 x+1 / x x x \ ln ----- = 2 * | x + ---- + ---- + ---- + ... | x-1 \ 3 5 7 / But now we want to find ln(a), so we need to find the value of x such that a = (x+1)/(x-1). This is easily solved to find that x = (a-1)/(a+1).
Params:
  • a – number from which logarithm is requested, in split form
Returns:log(a)
/** Computes the natural log of a number between 0 and 2. * Let f(x) = ln(x), * * We know that f'(x) = 1/x, thus from Taylor's theorum we have: * * ----- n+1 n * f(x) = \ (-1) (x - 1) * / ---------------- for 1 <= n <= infinity * ----- n * * or * 2 3 4 * (x-1) (x-1) (x-1) * ln(x) = (x-1) - ----- + ------ - ------ + ... * 2 3 4 * * alternatively, * * 2 3 4 * x x x * ln(x+1) = x - - + - - - + ... * 2 3 4 * * This series can be used to compute ln(x), but it converges too slowly. * * If we substitute -x for x above, we get * * 2 3 4 * x x x * ln(1-x) = -x - - - - - - + ... * 2 3 4 * * Note that all terms are now negative. Because the even powered ones * absorbed the sign. Now, subtract the series above from the previous * one to get ln(x+1) - ln(1-x). Note the even terms cancel out leaving * only the odd ones * * 3 5 7 * 2x 2x 2x * ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ... * 3 5 7 * * By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have: * * 3 5 7 * x+1 / x x x \ * ln ----- = 2 * | x + ---- + ---- + ---- + ... | * x-1 \ 3 5 7 / * * But now we want to find ln(a), so we need to find the value of x * such that a = (x+1)/(x-1). This is easily solved to find that * x = (a-1)/(a+1). * @param a number from which logarithm is requested, in split form * @return log(a) */
protected static Dfp[] logInternal(final Dfp a[]) { /* Now we want to compute x = (a-1)/(a+1) but this is prone to * loss of precision. So instead, compute x = (a/4 - 1/4) / (a/4 + 1/4) */ Dfp t = a[0].divide(4).add(a[1].divide(4)); Dfp x = t.add(a[0].newInstance("-0.25")).divide(t.add(a[0].newInstance("0.25"))); Dfp y = new Dfp(x); Dfp num = new Dfp(x); Dfp py = new Dfp(y); int den = 1; for (int i = 0; i < 10000; i++) { num = num.multiply(x); num = num.multiply(x); den += 2; t = num.divide(den); y = y.add(t); if (y.equals(py)) { break; } py = new Dfp(y); } y = y.multiply(a[0].getTwo()); return split(y); }
Computes x to the y power.

Uses the following method:

  1. Set u = rint(y), v = y-u
  2. Compute a = v * ln(x)
  3. Compute b = rint( a/ln(2) )
  4. Compute c = a - b*ln(2)
  5. xy = xu * 2b * ec
if |y| > 1e8, then we compute by exp(y*ln(x))

Special Cases

  • if y is 0.0 or -0.0 then result is 1.0
  • if y is 1.0 then result is x
  • if y is NaN then result is NaN
  • if x is NaN and y is not zero then result is NaN
  • if |x| > 1.0 and y is +Infinity then result is +Infinity
  • if |x| < 1.0 and y is -Infinity then result is +Infinity
  • if |x| > 1.0 and y is -Infinity then result is +0
  • if |x| < 1.0 and y is +Infinity then result is +0
  • if |x| = 1.0 and y is +/-Infinity then result is NaN
  • if x = +0 and y > 0 then result is +0
  • if x = +Inf and y < 0 then result is +0
  • if x = +0 and y < 0 then result is +Inf
  • if x = +Inf and y > 0 then result is +Inf
  • if x = -0 and y > 0, finite, not odd integer then result is +0
  • if x = -0 and y < 0, finite, and odd integer then result is -Inf
  • if x = -Inf and y > 0, finite, and odd integer then result is -Inf
  • if x = -0 and y < 0, not finite odd integer then result is +Inf
  • if x = -Inf and y > 0, not finite odd integer then result is +Inf
  • if x < 0 and y > 0, finite, and odd integer then result is -(|x|y)
  • if x < 0 and y > 0, finite, and not integer then result is NaN
@param x base to be raised @param y power to which base should be raised @return xy
/** Computes x to the y power.<p> * * Uses the following method:<p> * * <ol> * <li> Set u = rint(y), v = y-u * <li> Compute a = v * ln(x) * <li> Compute b = rint( a/ln(2) ) * <li> Compute c = a - b*ln(2) * <li> x<sup>y</sup> = x<sup>u</sup> * 2<sup>b</sup> * e<sup>c</sup> * </ol> * if |y| > 1e8, then we compute by exp(y*ln(x)) <p> * * <b>Special Cases</b><p> * <ul> * <li> if y is 0.0 or -0.0 then result is 1.0 * <li> if y is 1.0 then result is x * <li> if y is NaN then result is NaN * <li> if x is NaN and y is not zero then result is NaN * <li> if |x| > 1.0 and y is +Infinity then result is +Infinity * <li> if |x| < 1.0 and y is -Infinity then result is +Infinity * <li> if |x| > 1.0 and y is -Infinity then result is +0 * <li> if |x| < 1.0 and y is +Infinity then result is +0 * <li> if |x| = 1.0 and y is +/-Infinity then result is NaN * <li> if x = +0 and y > 0 then result is +0 * <li> if x = +Inf and y < 0 then result is +0 * <li> if x = +0 and y < 0 then result is +Inf * <li> if x = +Inf and y > 0 then result is +Inf * <li> if x = -0 and y > 0, finite, not odd integer then result is +0 * <li> if x = -0 and y < 0, finite, and odd integer then result is -Inf * <li> if x = -Inf and y > 0, finite, and odd integer then result is -Inf * <li> if x = -0 and y < 0, not finite odd integer then result is +Inf * <li> if x = -Inf and y > 0, not finite odd integer then result is +Inf * <li> if x < 0 and y > 0, finite, and odd integer then result is -(|x|<sup>y</sup>) * <li> if x < 0 and y > 0, finite, and not integer then result is NaN * </ul> * @param x base to be raised * @param y power to which base should be raised * @return x<sup>y</sup> */
public static Dfp pow(Dfp x, final Dfp y) { // make sure we don't mix number with different precision if (x.getField().getRadixDigits() != y.getField().getRadixDigits()) { x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); final Dfp result = x.newInstance(x.getZero()); result.nans = Dfp.QNAN; return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, result); } final Dfp zero = x.getZero(); final Dfp one = x.getOne(); final Dfp two = x.getTwo(); boolean invert = false; int ui; /* Check for special cases */ if (y.equals(zero)) { return x.newInstance(one); } if (y.equals(one)) { if (x.isNaN()) { // Test for NaNs x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x); } return x; } if (x.isNaN() || y.isNaN()) { // Test for NaNs x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN)); } // X == 0 if (x.equals(zero)) { if (Dfp.copysign(one, x).greaterThan(zero)) { // X == +0 if (y.greaterThan(zero)) { return x.newInstance(zero); } else { return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE)); } } else { // X == -0 if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) { // If y is odd integer if (y.greaterThan(zero)) { return x.newInstance(zero.negate()); } else { return x.newInstance(x.newInstance((byte)-1, Dfp.INFINITE)); } } else { // Y is not odd integer if (y.greaterThan(zero)) { return x.newInstance(zero); } else { return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE)); } } } } if (x.lessThan(zero)) { // Make x positive, but keep track of it x = x.negate(); invert = true; } if (x.greaterThan(one) && y.classify() == Dfp.INFINITE) { if (y.greaterThan(zero)) { return y; } else { return x.newInstance(zero); } } if (x.lessThan(one) && y.classify() == Dfp.INFINITE) { if (y.greaterThan(zero)) { return x.newInstance(zero); } else { return x.newInstance(Dfp.copysign(y, one)); } } if (x.equals(one) && y.classify() == Dfp.INFINITE) { x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN)); } if (x.classify() == Dfp.INFINITE) { // x = +/- inf if (invert) { // negative infinity if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) { // If y is odd integer if (y.greaterThan(zero)) { return x.newInstance(x.newInstance((byte)-1, Dfp.INFINITE)); } else { return x.newInstance(zero.negate()); } } else { // Y is not odd integer if (y.greaterThan(zero)) { return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE)); } else { return x.newInstance(zero); } } } else { // positive infinity if (y.greaterThan(zero)) { return x; } else { return x.newInstance(zero); } } } if (invert && !y.rint().equals(y)) { x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN)); } // End special cases Dfp r; if (y.lessThan(x.newInstance(100000000)) && y.greaterThan(x.newInstance(-100000000))) { final Dfp u = y.rint(); ui = u.intValue(); final Dfp v = y.subtract(u); if (v.unequal(zero)) { final Dfp a = v.multiply(log(x)); final Dfp b = a.divide(x.getField().getLn2()).rint(); final Dfp c = a.subtract(b.multiply(x.getField().getLn2())); r = splitPow(split(x), ui); r = r.multiply(pow(two, b.intValue())); r = r.multiply(exp(c)); } else { r = splitPow(split(x), ui); } } else { // very large exponent. |y| > 1e8 r = exp(log(x).multiply(y)); } if (invert && y.rint().equals(y) && !y.remainder(two).equals(zero)) { // if y is odd integer r = r.negate(); } return x.newInstance(r); }
Computes sin(a) Used when 0 < a < pi/4. Uses the classic Taylor series. x - x**3/3! + x**5/5! ...
Params:
  • a – number from which sine is desired, in split form
Returns:sin(a)
/** Computes sin(a) Used when 0 < a < pi/4. * Uses the classic Taylor series. x - x**3/3! + x**5/5! ... * @param a number from which sine is desired, in split form * @return sin(a) */
protected static Dfp sinInternal(Dfp a[]) { Dfp c = a[0].add(a[1]); Dfp y = c; c = c.multiply(c); Dfp x = y; Dfp fact = a[0].getOne(); Dfp py = new Dfp(y); for (int i = 3; i < 90; i += 2) { x = x.multiply(c); x = x.negate(); fact = fact.divide((i-1)*i); // 1 over fact y = y.add(x.multiply(fact)); if (y.equals(py)) { break; } py = new Dfp(y); } return y; }
Computes cos(a) Used when 0 < a < pi/4. Uses the classic Taylor series for cosine. 1 - x**2/2! + x**4/4! ...
Params:
  • a – number from which cosine is desired, in split form
Returns:cos(a)
/** Computes cos(a) Used when 0 < a < pi/4. * Uses the classic Taylor series for cosine. 1 - x**2/2! + x**4/4! ... * @param a number from which cosine is desired, in split form * @return cos(a) */
protected static Dfp cosInternal(Dfp a[]) { final Dfp one = a[0].getOne(); Dfp x = one; Dfp y = one; Dfp c = a[0].add(a[1]); c = c.multiply(c); Dfp fact = one; Dfp py = new Dfp(y); for (int i = 2; i < 90; i += 2) { x = x.multiply(c); x = x.negate(); fact = fact.divide((i - 1) * i); // 1 over fact y = y.add(x.multiply(fact)); if (y.equals(py)) { break; } py = new Dfp(y); } return y; }
computes the sine of the argument.
Params:
  • a – number from which sine is desired
Returns:sin(a)
/** computes the sine of the argument. * @param a number from which sine is desired * @return sin(a) */
public static Dfp sin(final Dfp a) { final Dfp pi = a.getField().getPi(); final Dfp zero = a.getField().getZero(); boolean neg = false; /* First reduce the argument to the range of +/- PI */ Dfp x = a.remainder(pi.multiply(2)); /* if x < 0 then apply identity sin(-x) = -sin(x) */ /* This puts x in the range 0 < x < PI */ if (x.lessThan(zero)) { x = x.negate(); neg = true; } /* Since sine(x) = sine(pi - x) we can reduce the range to * 0 < x < pi/2 */ if (x.greaterThan(pi.divide(2))) { x = pi.subtract(x); } Dfp y; if (x.lessThan(pi.divide(4))) { y = sinInternal(split(x)); } else { final Dfp c[] = new Dfp[2]; final Dfp[] piSplit = a.getField().getPiSplit(); c[0] = piSplit[0].divide(2).subtract(x); c[1] = piSplit[1].divide(2); y = cosInternal(c); } if (neg) { y = y.negate(); } return a.newInstance(y); }
computes the cosine of the argument.
Params:
  • a – number from which cosine is desired
Returns:cos(a)
/** computes the cosine of the argument. * @param a number from which cosine is desired * @return cos(a) */
public static Dfp cos(Dfp a) { final Dfp pi = a.getField().getPi(); final Dfp zero = a.getField().getZero(); boolean neg = false; /* First reduce the argument to the range of +/- PI */ Dfp x = a.remainder(pi.multiply(2)); /* if x < 0 then apply identity cos(-x) = cos(x) */ /* This puts x in the range 0 < x < PI */ if (x.lessThan(zero)) { x = x.negate(); } /* Since cos(x) = -cos(pi - x) we can reduce the range to * 0 < x < pi/2 */ if (x.greaterThan(pi.divide(2))) { x = pi.subtract(x); neg = true; } Dfp y; if (x.lessThan(pi.divide(4))) { Dfp c[] = new Dfp[2]; c[0] = x; c[1] = zero; y = cosInternal(c); } else { final Dfp c[] = new Dfp[2]; final Dfp[] piSplit = a.getField().getPiSplit(); c[0] = piSplit[0].divide(2).subtract(x); c[1] = piSplit[1].divide(2); y = sinInternal(c); } if (neg) { y = y.negate(); } return a.newInstance(y); }
computes the tangent of the argument.
Params:
  • a – number from which tangent is desired
Returns:tan(a)
/** computes the tangent of the argument. * @param a number from which tangent is desired * @return tan(a) */
public static Dfp tan(final Dfp a) { return sin(a).divide(cos(a)); }
computes the arc-tangent of the argument.
Params:
  • a – number from which arc-tangent is desired
Returns:atan(a)
/** computes the arc-tangent of the argument. * @param a number from which arc-tangent is desired * @return atan(a) */
protected static Dfp atanInternal(final Dfp a) { Dfp y = new Dfp(a); Dfp x = new Dfp(y); Dfp py = new Dfp(y); for (int i = 3; i < 90; i += 2) { x = x.multiply(a); x = x.multiply(a); x = x.negate(); y = y.add(x.divide(i)); if (y.equals(py)) { break; } py = new Dfp(y); } return y; }
computes the arc tangent of the argument Uses the typical taylor series but may reduce arguments using the following identity tan(x+y) = (tan(x) + tan(y)) / (1 - tan(x)*tan(y)) since tan(PI/8) = sqrt(2)-1, atan(x) = atan( (x - sqrt(2) + 1) / (1+x*sqrt(2) - x) + PI/8.0
Params:
  • a – number from which arc-tangent is desired
Returns:atan(a)
/** computes the arc tangent of the argument * * Uses the typical taylor series * * but may reduce arguments using the following identity * tan(x+y) = (tan(x) + tan(y)) / (1 - tan(x)*tan(y)) * * since tan(PI/8) = sqrt(2)-1, * * atan(x) = atan( (x - sqrt(2) + 1) / (1+x*sqrt(2) - x) + PI/8.0 * @param a number from which arc-tangent is desired * @return atan(a) */
public static Dfp atan(final Dfp a) { final Dfp zero = a.getField().getZero(); final Dfp one = a.getField().getOne(); final Dfp[] sqr2Split = a.getField().getSqr2Split(); final Dfp[] piSplit = a.getField().getPiSplit(); boolean recp = false; boolean neg = false; boolean sub = false; final Dfp ty = sqr2Split[0].subtract(one).add(sqr2Split[1]); Dfp x = new Dfp(a); if (x.lessThan(zero)) { neg = true; x = x.negate(); } if (x.greaterThan(one)) { recp = true; x = one.divide(x); } if (x.greaterThan(ty)) { Dfp sty[] = new Dfp[2]; sub = true; sty[0] = sqr2Split[0].subtract(one); sty[1] = sqr2Split[1]; Dfp[] xs = split(x); Dfp[] ds = splitMult(xs, sty); ds[0] = ds[0].add(one); xs[0] = xs[0].subtract(sty[0]); xs[1] = xs[1].subtract(sty[1]); xs = splitDiv(xs, ds); x = xs[0].add(xs[1]); //x = x.subtract(ty).divide(dfp.one.add(x.multiply(ty))); } Dfp y = atanInternal(x); if (sub) { y = y.add(piSplit[0].divide(8)).add(piSplit[1].divide(8)); } if (recp) { y = piSplit[0].divide(2).subtract(y).add(piSplit[1].divide(2)); } if (neg) { y = y.negate(); } return a.newInstance(y); }
computes the arc-sine of the argument.
Params:
  • a – number from which arc-sine is desired
Returns:asin(a)
/** computes the arc-sine of the argument. * @param a number from which arc-sine is desired * @return asin(a) */
public static Dfp asin(final Dfp a) { return atan(a.divide(a.getOne().subtract(a.multiply(a)).sqrt())); }
computes the arc-cosine of the argument.
Params:
  • a – number from which arc-cosine is desired
Returns:acos(a)
/** computes the arc-cosine of the argument. * @param a number from which arc-cosine is desired * @return acos(a) */
public static Dfp acos(Dfp a) { Dfp result; boolean negative = false; if (a.lessThan(a.getZero())) { negative = true; } a = Dfp.copysign(a, a.getOne()); // absolute value result = atan(a.getOne().subtract(a.multiply(a)).sqrt().divide(a)); if (negative) { result = a.getField().getPi().subtract(result); } return a.newInstance(result); } }