diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math/dfp/DfpMath.java')
-rw-r--r-- | src/main/java/org/apache/commons/math/dfp/DfpMath.java | 969 |
1 files changed, 969 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math/dfp/DfpMath.java b/src/main/java/org/apache/commons/math/dfp/DfpMath.java new file mode 100644 index 0000000..56af1d2 --- /dev/null +++ b/src/main/java/org/apache/commons/math/dfp/DfpMath.java @@ -0,0 +1,969 @@ +/* + * 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.math.dfp; + +/** Mathematical routines for use with {@link Dfp}. + * The constants are defined in {@link DfpField} + * @version $Revision: 1042376 $ $Date: 2010-12-05 16:54:55 +0100 (dim. 05 déc. 2010) $ + * @since 2.2 + */ +public class DfpMath { + + /** Name for traps triggered by pow. */ + private static final String POW_TRAP = "pow"; + + /** + * Private Constructor. + */ + private DfpMath() { + } + + /** 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 {@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 + */ + 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 + */ + 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. + * @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 = 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. + * @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 = trial * 2; + } while (a>trial); + + r = prevr; + trial = prevtrial; + + a = 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 e<sup>n</sup> and a Taylor series to compute + * e<sup>m</sup>. We return e*<sup>n</sup> × 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! ... + * @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. + * @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). + * @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 = 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.<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) { + // if y is odd integer + if (y.rint().equals(y) && !y.remainder(two).equals(zero)) { + 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! ... + * @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! ... + * @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. + * @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))) { + Dfp c[] = new Dfp[2]; + c[0] = x; + c[1] = zero; + + //y = sinInternal(c); + 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. + * @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. + * @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. + * @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 + * @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. + * @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. + * @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); + } + +} |