| /* |
| * 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); |
| } |
| |
| } |
