From 806fa1cdd2945ef2234666270f6aa2886c9e636f Mon Sep 17 00:00:00 2001 From: Sebastian Bazley Date: Sun, 11 Sep 2011 20:43:27 +0000 Subject: [PATCH] MATH-650 move static setup methods to helper class git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1169527 13f79535-47bb-0310-9956-ffa450edef68 --- .../apache/commons/math/util/FastMath.java | 603 +---------------- .../commons/math/util/FastMathCalc.java | 610 ++++++++++++++++++ 2 files changed, 625 insertions(+), 588 deletions(-) create mode 100644 src/main/java/org/apache/commons/math/util/FastMathCalc.java diff --git a/src/main/java/org/apache/commons/math/util/FastMath.java b/src/main/java/org/apache/commons/math/util/FastMath.java index 0e2e9a0f6..462fca189 100644 --- a/src/main/java/org/apache/commons/math/util/FastMath.java +++ b/src/main/java/org/apache/commons/math/util/FastMath.java @@ -16,8 +16,6 @@ */ package org.apache.commons.math.util; -import org.apache.commons.math.exception.DimensionMismatchException; - /** * Faster, more accurate, portable alternative to {@link Math} and * {@link StrictMath} for large scale computation. @@ -121,13 +119,13 @@ public class FastMath { // Populate expIntTable for (int i = 0; i < FastMath.EXP_INT_TABLE_MAX_INDEX; i++) { - expint(i, tmp); + FastMathCalc.expint(i, tmp); EXP_INT_TABLE_A[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[0]; EXP_INT_TABLE_B[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[1]; if (i != 0) { // Negative integer powers - splitReciprocal(tmp, recip); + FastMathCalc.splitReciprocal(tmp, recip); EXP_INT_TABLE_A[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[0]; EXP_INT_TABLE_B[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[1]; } @@ -3167,7 +3165,7 @@ public class FastMath { // Populate expFracTable for (int i = 0; i < EXP_FRAC_TABLE_A.length; i++) { - slowexp(i/1024.0, tmp); // TWO_POWER_10 + FastMathCalc.slowexp(i/1024.0, tmp); // TWO_POWER_10 EXP_FRAC_TABLE_A[i] = tmp[0]; EXP_FRAC_TABLE_B[i] = tmp[1]; } @@ -5231,31 +5229,6 @@ public class FastMath { } } - /** Factorial table, for Taylor series expansions. 0!, 1!, 2!, ... 19! */ - private static final double FACT[] = new double[] - { - +1.0d, // 0 - +1.0d, // 1 - +2.0d, // 2 - +6.0d, // 3 - +24.0d, // 4 - +120.0d, // 5 - +720.0d, // 6 - +5040.0d, // 7 - +40320.0d, // 8 - +362880.0d, // 9 - +3628800.0d, // 10 - +39916800.0d, // 11 - +479001600.0d, // 12 - +6227020800.0d, // 13 - +87178291200.0d, // 14 - +1307674368000.0d, // 15 - +20922789888000.0d, // 16 - +355687428096000.0d, // 17 - +6402373705728000.0d, // 18 - +121645100408832000.0d, // 19 - }; - private static final int LN_MANT_LEN = 1024; // MAGIC NUMBER // Enclose large data table in nested static class so it's only loaded on first access @@ -5270,7 +5243,7 @@ public class FastMath { // Populate lnMant table for (int i = 0; i < LN_MANT.length; i++) { final double d = Double.longBitsToDouble( (((long) i) << 42) | 0x3ff0000000000000L ); - LN_MANT[i] = slowLog(d); + LN_MANT[i] = FastMathCalc.slowLog(d); } } else { LN_MANT = new double[][] { @@ -6309,26 +6282,6 @@ public class FastMath { /** log(2) (low bits). */ private static final double LN_2_B = 1.17304635250823482e-7; - /** Coefficients for slowLog. */ - private static final double LN_SPLIT_COEF[][] = { - {2.0, 0.0}, - {0.6666666269302368, 3.9736429850260626E-8}, - {0.3999999761581421, 2.3841857910019882E-8}, - {0.2857142686843872, 1.7029898543501842E-8}, - {0.2222222089767456, 1.3245471311735498E-8}, - {0.1818181574344635, 2.4384203044354907E-8}, - {0.1538461446762085, 9.140260083262505E-9}, - {0.13333332538604736, 9.220590270857665E-9}, - {0.11764700710773468, 1.2393345855018391E-8}, - {0.10526403784751892, 8.251545029714408E-9}, - {0.0952233225107193, 1.2675934823758863E-8}, - {0.08713622391223907, 1.1430250008909141E-8}, - {0.07842259109020233, 2.404307984052299E-9}, - {0.08371849358081818, 1.176342548272881E-8}, - {0.030589580535888672, 1.2958646899018938E-9}, - {0.14982303977012634, 1.225743062930824E-8}, - }; - /** Coefficients for log, when input 0.99 < x < 1.01. */ private static final double LN_QUICK_COEF[][] = { {1.0, 5.669184079525E-24}, @@ -6539,50 +6492,17 @@ public class FastMath { // } public static void main(String[] a){ - printarray("EXP_INT_TABLE_A", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_A); - printarray("EXP_INT_TABLE_B", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_B); - printarray("EXP_FRAC_TABLE_A", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_A); - printarray("EXP_FRAC_TABLE_B", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_B); - printarray("LN_MANT",LN_MANT_LEN, lnMant.LN_MANT); - printarray("SINE_TABLE_A", SINE_TABLE_LEN, SINE_TABLE_A); - printarray("SINE_TABLE_B", SINE_TABLE_LEN, SINE_TABLE_B); - printarray("COSINE_TABLE_A", SINE_TABLE_LEN, COSINE_TABLE_A); - printarray("COSINE_TABLE_B", SINE_TABLE_LEN, COSINE_TABLE_B); - printarray("TANGENT_TABLE_A", SINE_TABLE_LEN, TANGENT_TABLE_A); - printarray("TANGENT_TABLE_B", SINE_TABLE_LEN, TANGENT_TABLE_B); - } - - private static void printarray(String string, int expectedLen, double[][] array2d) { - System.out.println(string); - checkLen(expectedLen, array2d.length); - System.out.println(" { "); - int i = 0; - for(double array[] : array2d) { - System.out.print(" {"); - for(double d : array) { // assume inner array has very few entries - String ds = d >= 0 ? "+"+Double.toString(d)+"d," : Double.toString(d)+"d,"; - System.out.printf("%-25.25s",ds); // multiple entries per line - } - System.out.println("}, // "+i++); - } - System.out.println(" };"); - } - - private static void printarray(String string, int expectedLen, double[] array) { - System.out.println(string+"="); - checkLen(expectedLen, array.length); - System.out.println(" {"); - for(double d : array){ - String ds = d!=d ? "Double.NaN," : d >= 0 ? "+"+Double.toString(d)+"d," : Double.toString(d)+"d,"; - System.out.printf(" %s%n",ds); // one entry per line - } - System.out.println(" };"); - } - - private static void checkLen(int expectedLen, int actual) { - if (expectedLen != actual) { - throw new DimensionMismatchException(actual, expectedLen); - } + FastMathCalc.printarray("EXP_INT_TABLE_A", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_A); + FastMathCalc.printarray("EXP_INT_TABLE_B", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_B); + FastMathCalc.printarray("EXP_FRAC_TABLE_A", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_A); + FastMathCalc.printarray("EXP_FRAC_TABLE_B", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_B); + FastMathCalc.printarray("LN_MANT",LN_MANT_LEN, lnMant.LN_MANT); + FastMathCalc.printarray("SINE_TABLE_A", SINE_TABLE_LEN, SINE_TABLE_A); + FastMathCalc.printarray("SINE_TABLE_B", SINE_TABLE_LEN, SINE_TABLE_B); + FastMathCalc.printarray("COSINE_TABLE_A", SINE_TABLE_LEN, COSINE_TABLE_A); + FastMathCalc.printarray("COSINE_TABLE_B", SINE_TABLE_LEN, COSINE_TABLE_B); + FastMathCalc.printarray("TANGENT_TABLE_A", SINE_TABLE_LEN, TANGENT_TABLE_A); + FastMathCalc.printarray("TANGENT_TABLE_B", SINE_TABLE_LEN, TANGENT_TABLE_B); } /** @@ -7333,251 +7253,6 @@ public class FastMath { return ya + yb; } - /** - * For x between 0 and 1, returns exp(x), uses extended precision - * @param x argument of exponential - * @param result placeholder where to place exp(x) split in two terms - * for extra precision (i.e. exp(x) = result[0] + result[1] - * @return exp(x) - */ - private static double slowexp(final double x, final double result[]) { - final double xs[] = new double[2]; - final double ys[] = new double[2]; - final double facts[] = new double[2]; - final double as[] = new double[2]; - split(x, xs); - ys[0] = ys[1] = 0.0; - - for (int i = FACT.length-1; i >= 0; i--) { - splitMult(xs, ys, as); - ys[0] = as[0]; - ys[1] = as[1]; - - split(FACT[i], as); - splitReciprocal(as, facts); - - splitAdd(ys, facts, as); - ys[0] = as[0]; - ys[1] = as[1]; - } - - if (result != null) { - result[0] = ys[0]; - result[1] = ys[1]; - } - - return ys[0] + ys[1]; - } - - /** Compute split[0], split[1] such that their sum is equal to d, - * and split[0] has its 30 least significant bits as zero. - * @param d number to split - * @param split placeholder where to place the result - */ - private static void split(final double d, final double split[]) { - if (d < 8e298 && d > -8e298) { - final double a = d * HEX_40000000; - split[0] = (d + a) - a; - split[1] = d - split[0]; - } else { - final double a = d * 9.31322574615478515625E-10; - split[0] = (d + a - d) * HEX_40000000; - split[1] = d - split[0]; - } - } - - /** Recompute a split. - * @param a input/out array containing the split, changed - * on output - */ - private static void resplit(final double a[]) { - final double c = a[0] + a[1]; - final double d = -(c - a[0] - a[1]); - - if (c < 8e298 && c > -8e298) { // MAGIC NUMBER - double z = c * HEX_40000000; - a[0] = (c + z) - z; - a[1] = c - a[0] + d; - } else { - double z = c * 9.31322574615478515625E-10; - a[0] = (c + z - c) * HEX_40000000; - a[1] = c - a[0] + d; - } - } - - /** Multiply two numbers in split form. - * @param a first term of multiplication - * @param b second term of multiplication - * @param ans placeholder where to put the result - */ - private static void splitMult(double a[], double b[], double ans[]) { - ans[0] = a[0] * b[0]; - ans[1] = a[0] * b[1] + a[1] * b[0] + a[1] * b[1]; - - /* Resplit */ - resplit(ans); - } - - /** Add two numbers in split form. - * @param a first term of addition - * @param b second term of addition - * @param ans placeholder where to put the result - */ - private static void splitAdd(final double a[], final double b[], final double ans[]) { - ans[0] = a[0] + b[0]; - ans[1] = a[1] + b[1]; - - resplit(ans); - } - - /** Compute the reciprocal of in. Use the following algorithm. - * in = c + d. - * want to find x + y such that x+y = 1/(c+d) and x is much - * larger than y and x has several zero bits on the right. - * - * Set b = 1/(2^22), a = 1 - b. Thus (a+b) = 1. - * Use following identity to compute (a+b)/(c+d) - * - * (a+b)/(c+d) = a/c + (bc - ad) / (c^2 + cd) - * set x = a/c and y = (bc - ad) / (c^2 + cd) - * This will be close to the right answer, but there will be - * some rounding in the calculation of X. So by carefully - * computing 1 - (c+d)(x+y) we can compute an error and - * add that back in. This is done carefully so that terms - * of similar size are subtracted first. - * @param in initial number, in split form - * @param result placeholder where to put the result - */ - private static void splitReciprocal(final double in[], final double result[]) { - final double b = 1.0/4194304.0; - final double a = 1.0 - b; - - if (in[0] == 0.0) { - in[0] = in[1]; - in[1] = 0.0; - } - - result[0] = a / in[0]; - result[1] = (b*in[0]-a*in[1]) / (in[0]*in[0] + in[0]*in[1]); - - if (result[1] != result[1]) { // can happen if result[1] is NAN - result[1] = 0.0; - } - - /* Resplit */ - resplit(result); - - for (int i = 0; i < 2; i++) { - /* this may be overkill, probably once is enough */ - double err = 1.0 - result[0] * in[0] - result[0] * in[1] - - result[1] * in[0] - result[1] * in[1]; - /*err = 1.0 - err; */ - err = err * (result[0] + result[1]); - /*printf("err = %16e\n", err); */ - result[1] += err; - } - } - - /** Compute (a[0] + a[1]) * (b[0] + b[1]) in extended precision. - * @param a first term of the multiplication - * @param b second term of the multiplication - * @param result placeholder where to put the result - */ - private static void quadMult(final double a[], final double b[], final double result[]) { - final double xs[] = new double[2]; - final double ys[] = new double[2]; - final double zs[] = new double[2]; - - /* a[0] * b[0] */ - split(a[0], xs); - split(b[0], ys); - splitMult(xs, ys, zs); - - result[0] = zs[0]; - result[1] = zs[1]; - - /* a[0] * b[1] */ - split(b[1], ys); - splitMult(xs, ys, zs); - - double tmp = result[0] + zs[0]; - result[1] = result[1] - (tmp - result[0] - zs[0]); - result[0] = tmp; - tmp = result[0] + zs[1]; - result[1] = result[1] - (tmp - result[0] - zs[1]); - result[0] = tmp; - - /* a[1] * b[0] */ - split(a[1], xs); - split(b[0], ys); - splitMult(xs, ys, zs); - - tmp = result[0] + zs[0]; - result[1] = result[1] - (tmp - result[0] - zs[0]); - result[0] = tmp; - tmp = result[0] + zs[1]; - result[1] = result[1] - (tmp - result[0] - zs[1]); - result[0] = tmp; - - /* a[1] * b[0] */ - split(a[1], xs); - split(b[1], ys); - splitMult(xs, ys, zs); - - tmp = result[0] + zs[0]; - result[1] = result[1] - (tmp - result[0] - zs[0]); - result[0] = tmp; - tmp = result[0] + zs[1]; - result[1] = result[1] - (tmp - result[0] - zs[1]); - result[0] = tmp; - } - - /** Compute exp(p) for a integer p in extended precision. - * @param p integer whose exponential is requested - * @param result placeholder where to put the result in extended precision - * @return exp(p) in standard precision (equal to result[0] + result[1]) - */ - private static double expint(int p, final double result[]) { - //double x = M_E; - final double xs[] = new double[2]; - final double as[] = new double[2]; - final double ys[] = new double[2]; - //split(x, xs); - //xs[1] = (double)(2.7182818284590452353602874713526625L - xs[0]); - //xs[0] = 2.71827697753906250000; - //xs[1] = 4.85091998273542816811e-06; - //xs[0] = Double.longBitsToDouble(0x4005bf0800000000L); - //xs[1] = Double.longBitsToDouble(0x3ed458a2bb4a9b00L); - - /* E */ - xs[0] = 2.718281828459045; - xs[1] = 1.4456468917292502E-16; - - split(1.0, ys); - - while (p > 0) { - if ((p & 1) != 0) { - quadMult(ys, xs, as); - ys[0] = as[0]; ys[1] = as[1]; - } - - quadMult(xs, xs, as); - xs[0] = as[0]; xs[1] = as[1]; - - p >>= 1; - } - - if (result != null) { - result[0] = ys[0]; - result[1] = ys[1]; - - resplit(result); - } - - return ys[0] + ys[1]; - } - - /** * Natural logarithm. * @@ -8048,254 +7723,6 @@ public class FastMath { return result; } - /** xi in the range of [1, 2]. - * 3 5 7 - * x+1 / x x x \ - * ln ----- = 2 * | x + ---- + ---- + ---- + ... | - * 1-x \ 3 5 7 / - * - * So, compute a Remez approximation of the following function - * - * ln ((sqrt(x)+1)/(1-sqrt(x))) / x - * - * This will be an even function with only positive coefficents. - * x is in the range [0 - 1/3]. - * - * Transform xi for input to the above function by setting - * x = (xi-1)/(xi+1). Input to the polynomial is x^2, then - * the result is multiplied by x. - * @param xi number from which log is requested - * @return log(xi) - */ - private static double[] slowLog(double xi) { - double x[] = new double[2]; - double x2[] = new double[2]; - double y[] = new double[2]; - double a[] = new double[2]; - - split(xi, x); - - /* Set X = (x-1)/(x+1) */ - x[0] += 1.0; - resplit(x); - splitReciprocal(x, a); - x[0] -= 2.0; - resplit(x); - splitMult(x, a, y); - x[0] = y[0]; - x[1] = y[1]; - - /* Square X -> X2*/ - splitMult(x, x, x2); - - - //x[0] -= 1.0; - //resplit(x); - - y[0] = LN_SPLIT_COEF[LN_SPLIT_COEF.length-1][0]; - y[1] = LN_SPLIT_COEF[LN_SPLIT_COEF.length-1][1]; - - for (int i = LN_SPLIT_COEF.length-2; i >= 0; i--) { - splitMult(y, x2, a); - y[0] = a[0]; - y[1] = a[1]; - splitAdd(y, LN_SPLIT_COEF[i], a); - y[0] = a[0]; - y[1] = a[1]; - } - - splitMult(y, x, a); - y[0] = a[0]; - y[1] = a[1]; - - return y; - } - - /** - * For x between 0 and pi/4 compute sine using Taylor expansion: - * sin(x) = x - x^3/3! + x^5/5! - x^7/7! ... - * @param x number from which sine is requested - * @param result placeholder where to put the result in extended precision - * (may be null) - * @return sin(x) - */ - private static double slowSin(final double x, final double result[]) { - final double xs[] = new double[2]; - final double ys[] = new double[2]; - final double facts[] = new double[2]; - final double as[] = new double[2]; - split(x, xs); - ys[0] = ys[1] = 0.0; - - for (int i = FACT.length-1; i >= 0; i--) { - splitMult(xs, ys, as); - ys[0] = as[0]; ys[1] = as[1]; - - if ( (i & 1) == 0) { // Ignore even numbers - continue; - } - - split(FACT[i], as); - splitReciprocal(as, facts); - - if ( (i & 2) != 0 ) { // alternate terms are negative - facts[0] = -facts[0]; - facts[1] = -facts[1]; - } - - splitAdd(ys, facts, as); - ys[0] = as[0]; ys[1] = as[1]; - } - - if (result != null) { - result[0] = ys[0]; - result[1] = ys[1]; - } - - return ys[0] + ys[1]; - } - - /** - * For x between 0 and pi/4 compute cosine using Talor series - * cos(x) = 1 - x^2/2! + x^4/4! ... - * @param x number from which cosine is requested - * @param result placeholder where to put the result in extended precision - * (may be null) - * @return cos(x) - */ - private static double slowCos(final double x, final double result[]) { - - final double xs[] = new double[2]; - final double ys[] = new double[2]; - final double facts[] = new double[2]; - final double as[] = new double[2]; - split(x, xs); - ys[0] = ys[1] = 0.0; - - for (int i = FACT.length-1; i >= 0; i--) { - splitMult(xs, ys, as); - ys[0] = as[0]; ys[1] = as[1]; - - if ( (i & 1) != 0) { // skip odd entries - continue; - } - - split(FACT[i], as); - splitReciprocal(as, facts); - - if ( (i & 2) != 0 ) { // alternate terms are negative - facts[0] = -facts[0]; - facts[1] = -facts[1]; - } - - splitAdd(ys, facts, as); - ys[0] = as[0]; ys[1] = as[1]; - } - - if (result != null) { - result[0] = ys[0]; - result[1] = ys[1]; - } - - return ys[0] + ys[1]; - } - - /** Build the sine and cosine tables. - */ - @SuppressWarnings("unused") - private static void buildSinCosTables() { - final double result[] = new double[2]; - - /* Use taylor series for 0 <= x <= 6/8 */ - for (int i = 0; i < 7; i++) { - double x = i / 8.0; - - slowSin(x, result); - SINE_TABLE_A[i] = result[0]; - SINE_TABLE_B[i] = result[1]; - - slowCos(x, result); - COSINE_TABLE_A[i] = result[0]; - COSINE_TABLE_B[i] = result[1]; - } - - /* Use angle addition formula to complete table to 13/8, just beyond pi/2 */ - for (int i = 7; i < SINE_TABLE_LEN; i++) { - double xs[] = new double[2]; - double ys[] = new double[2]; - double as[] = new double[2]; - double bs[] = new double[2]; - double temps[] = new double[2]; - - if ( (i & 1) == 0) { - // Even, use double angle - xs[0] = SINE_TABLE_A[i/2]; - xs[1] = SINE_TABLE_B[i/2]; - ys[0] = COSINE_TABLE_A[i/2]; - ys[1] = COSINE_TABLE_B[i/2]; - - /* compute sine */ - splitMult(xs, ys, result); - SINE_TABLE_A[i] = result[0] * 2.0; - SINE_TABLE_B[i] = result[1] * 2.0; - - /* Compute cosine */ - splitMult(ys, ys, as); - splitMult(xs, xs, temps); - temps[0] = -temps[0]; - temps[1] = -temps[1]; - splitAdd(as, temps, result); - COSINE_TABLE_A[i] = result[0]; - COSINE_TABLE_B[i] = result[1]; - } else { - xs[0] = SINE_TABLE_A[i/2]; - xs[1] = SINE_TABLE_B[i/2]; - ys[0] = COSINE_TABLE_A[i/2]; - ys[1] = COSINE_TABLE_B[i/2]; - as[0] = SINE_TABLE_A[i/2+1]; - as[1] = SINE_TABLE_B[i/2+1]; - bs[0] = COSINE_TABLE_A[i/2+1]; - bs[1] = COSINE_TABLE_B[i/2+1]; - - /* compute sine */ - splitMult(xs, bs, temps); - splitMult(ys, as, result); - splitAdd(result, temps, result); - SINE_TABLE_A[i] = result[0]; - SINE_TABLE_B[i] = result[1]; - - /* Compute cosine */ - splitMult(ys, bs, result); - splitMult(xs, as, temps); - temps[0] = -temps[0]; - temps[1] = -temps[1]; - splitAdd(result, temps, result); - COSINE_TABLE_A[i] = result[0]; - COSINE_TABLE_B[i] = result[1]; - } - } - - /* Compute tangent = sine/cosine */ - for (int i = 0; i < SINE_TABLE_LEN; i++) { - double xs[] = new double[2]; - double ys[] = new double[2]; - double as[] = new double[2]; - - as[0] = COSINE_TABLE_A[i]; - as[1] = COSINE_TABLE_B[i]; - - splitReciprocal(as, ys); - - xs[0] = SINE_TABLE_A[i]; - xs[1] = SINE_TABLE_B[i]; - - splitMult(xs, ys, as); - - TANGENT_TABLE_A[i] = as[0]; - TANGENT_TABLE_B[i] = as[1]; - } - - } /** * Computes sin(x) - x, where |x| < 1/16. diff --git a/src/main/java/org/apache/commons/math/util/FastMathCalc.java b/src/main/java/org/apache/commons/math/util/FastMathCalc.java new file mode 100644 index 000000000..88f8aa342 --- /dev/null +++ b/src/main/java/org/apache/commons/math/util/FastMathCalc.java @@ -0,0 +1,610 @@ +/* + * 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.util; + +import org.apache.commons.math.exception.DimensionMismatchException; + +class FastMathCalc { + + /** + * 0x40000000 - used to split a double into two parts, both with the low order bits cleared. + * Equivalent to 2^30. + */ + private static final long HEX_40000000 = 0x40000000L; // 1073741824L + + /** Factorial table, for Taylor series expansions. 0!, 1!, 2!, ... 19! */ + private static final double FACT[] = new double[] + { + +1.0d, // 0 + +1.0d, // 1 + +2.0d, // 2 + +6.0d, // 3 + +24.0d, // 4 + +120.0d, // 5 + +720.0d, // 6 + +5040.0d, // 7 + +40320.0d, // 8 + +362880.0d, // 9 + +3628800.0d, // 10 + +39916800.0d, // 11 + +479001600.0d, // 12 + +6227020800.0d, // 13 + +87178291200.0d, // 14 + +1307674368000.0d, // 15 + +20922789888000.0d, // 16 + +355687428096000.0d, // 17 + +6402373705728000.0d, // 18 + +121645100408832000.0d, // 19 + }; + + /** Coefficients for slowLog. */ + private static final double LN_SPLIT_COEF[][] = { + {2.0, 0.0}, + {0.6666666269302368, 3.9736429850260626E-8}, + {0.3999999761581421, 2.3841857910019882E-8}, + {0.2857142686843872, 1.7029898543501842E-8}, + {0.2222222089767456, 1.3245471311735498E-8}, + {0.1818181574344635, 2.4384203044354907E-8}, + {0.1538461446762085, 9.140260083262505E-9}, + {0.13333332538604736, 9.220590270857665E-9}, + {0.11764700710773468, 1.2393345855018391E-8}, + {0.10526403784751892, 8.251545029714408E-9}, + {0.0952233225107193, 1.2675934823758863E-8}, + {0.08713622391223907, 1.1430250008909141E-8}, + {0.07842259109020233, 2.404307984052299E-9}, + {0.08371849358081818, 1.176342548272881E-8}, + {0.030589580535888672, 1.2958646899018938E-9}, + {0.14982303977012634, 1.225743062930824E-8}, + }; + + /** Build the sine and cosine tables. + * @param SINE_TABLE_A + * @param SINE_TABLE_B + * @param COSINE_TABLE_A + * @param COSINE_TABLE_B + * @param SINE_TABLE_LEN + * @param TANGENT_TABLE_A + * @param TANGENT_TABLE_B + */ + @SuppressWarnings("unused") + private static void buildSinCosTables(double[] SINE_TABLE_A, double[] SINE_TABLE_B, double[] COSINE_TABLE_A, double[] COSINE_TABLE_B, int SINE_TABLE_LEN, double[] TANGENT_TABLE_A, double[] TANGENT_TABLE_B) { + final double result[] = new double[2]; + + /* Use taylor series for 0 <= x <= 6/8 */ + for (int i = 0; i < 7; i++) { + double x = i / 8.0; + + slowSin(x, result); + SINE_TABLE_A[i] = result[0]; + SINE_TABLE_B[i] = result[1]; + + slowCos(x, result); + COSINE_TABLE_A[i] = result[0]; + COSINE_TABLE_B[i] = result[1]; + } + + /* Use angle addition formula to complete table to 13/8, just beyond pi/2 */ + for (int i = 7; i < SINE_TABLE_LEN; i++) { + double xs[] = new double[2]; + double ys[] = new double[2]; + double as[] = new double[2]; + double bs[] = new double[2]; + double temps[] = new double[2]; + + if ( (i & 1) == 0) { + // Even, use double angle + xs[0] = SINE_TABLE_A[i/2]; + xs[1] = SINE_TABLE_B[i/2]; + ys[0] = COSINE_TABLE_A[i/2]; + ys[1] = COSINE_TABLE_B[i/2]; + + /* compute sine */ + splitMult(xs, ys, result); + SINE_TABLE_A[i] = result[0] * 2.0; + SINE_TABLE_B[i] = result[1] * 2.0; + + /* Compute cosine */ + splitMult(ys, ys, as); + splitMult(xs, xs, temps); + temps[0] = -temps[0]; + temps[1] = -temps[1]; + splitAdd(as, temps, result); + COSINE_TABLE_A[i] = result[0]; + COSINE_TABLE_B[i] = result[1]; + } else { + xs[0] = SINE_TABLE_A[i/2]; + xs[1] = SINE_TABLE_B[i/2]; + ys[0] = COSINE_TABLE_A[i/2]; + ys[1] = COSINE_TABLE_B[i/2]; + as[0] = SINE_TABLE_A[i/2+1]; + as[1] = SINE_TABLE_B[i/2+1]; + bs[0] = COSINE_TABLE_A[i/2+1]; + bs[1] = COSINE_TABLE_B[i/2+1]; + + /* compute sine */ + splitMult(xs, bs, temps); + splitMult(ys, as, result); + splitAdd(result, temps, result); + SINE_TABLE_A[i] = result[0]; + SINE_TABLE_B[i] = result[1]; + + /* Compute cosine */ + splitMult(ys, bs, result); + splitMult(xs, as, temps); + temps[0] = -temps[0]; + temps[1] = -temps[1]; + splitAdd(result, temps, result); + COSINE_TABLE_A[i] = result[0]; + COSINE_TABLE_B[i] = result[1]; + } + } + + /* Compute tangent = sine/cosine */ + for (int i = 0; i < SINE_TABLE_LEN; i++) { + double xs[] = new double[2]; + double ys[] = new double[2]; + double as[] = new double[2]; + + as[0] = COSINE_TABLE_A[i]; + as[1] = COSINE_TABLE_B[i]; + + splitReciprocal(as, ys); + + xs[0] = SINE_TABLE_A[i]; + xs[1] = SINE_TABLE_B[i]; + + splitMult(xs, ys, as); + + TANGENT_TABLE_A[i] = as[0]; + TANGENT_TABLE_B[i] = as[1]; + } + + } + + /** + * For x between 0 and pi/4 compute cosine using Talor series + * cos(x) = 1 - x^2/2! + x^4/4! ... + * @param x number from which cosine is requested + * @param result placeholder where to put the result in extended precision + * (may be null) + * @return cos(x) + */ + static double slowCos(final double x, final double result[]) { + + final double xs[] = new double[2]; + final double ys[] = new double[2]; + final double facts[] = new double[2]; + final double as[] = new double[2]; + split(x, xs); + ys[0] = ys[1] = 0.0; + + for (int i = FACT.length-1; i >= 0; i--) { + splitMult(xs, ys, as); + ys[0] = as[0]; ys[1] = as[1]; + + if ( (i & 1) != 0) { // skip odd entries + continue; + } + + split(FACT[i], as); + splitReciprocal(as, facts); + + if ( (i & 2) != 0 ) { // alternate terms are negative + facts[0] = -facts[0]; + facts[1] = -facts[1]; + } + + splitAdd(ys, facts, as); + ys[0] = as[0]; ys[1] = as[1]; + } + + if (result != null) { + result[0] = ys[0]; + result[1] = ys[1]; + } + + return ys[0] + ys[1]; + } + + /** + * For x between 0 and pi/4 compute sine using Taylor expansion: + * sin(x) = x - x^3/3! + x^5/5! - x^7/7! ... + * @param x number from which sine is requested + * @param result placeholder where to put the result in extended precision + * (may be null) + * @return sin(x) + */ + static double slowSin(final double x, final double result[]) { + final double xs[] = new double[2]; + final double ys[] = new double[2]; + final double facts[] = new double[2]; + final double as[] = new double[2]; + split(x, xs); + ys[0] = ys[1] = 0.0; + + for (int i = FACT.length-1; i >= 0; i--) { + splitMult(xs, ys, as); + ys[0] = as[0]; ys[1] = as[1]; + + if ( (i & 1) == 0) { // Ignore even numbers + continue; + } + + split(FACT[i], as); + splitReciprocal(as, facts); + + if ( (i & 2) != 0 ) { // alternate terms are negative + facts[0] = -facts[0]; + facts[1] = -facts[1]; + } + + splitAdd(ys, facts, as); + ys[0] = as[0]; ys[1] = as[1]; + } + + if (result != null) { + result[0] = ys[0]; + result[1] = ys[1]; + } + + return ys[0] + ys[1]; + } + + + /** + * For x between 0 and 1, returns exp(x), uses extended precision + * @param x argument of exponential + * @param result placeholder where to place exp(x) split in two terms + * for extra precision (i.e. exp(x) = result[0] + result[1] + * @return exp(x) + */ + static double slowexp(final double x, final double result[]) { + final double xs[] = new double[2]; + final double ys[] = new double[2]; + final double facts[] = new double[2]; + final double as[] = new double[2]; + split(x, xs); + ys[0] = ys[1] = 0.0; + + for (int i = FACT.length-1; i >= 0; i--) { + splitMult(xs, ys, as); + ys[0] = as[0]; + ys[1] = as[1]; + + split(FACT[i], as); + splitReciprocal(as, facts); + + splitAdd(ys, facts, as); + ys[0] = as[0]; + ys[1] = as[1]; + } + + if (result != null) { + result[0] = ys[0]; + result[1] = ys[1]; + } + + return ys[0] + ys[1]; + } + + /** Compute split[0], split[1] such that their sum is equal to d, + * and split[0] has its 30 least significant bits as zero. + * @param d number to split + * @param split placeholder where to place the result + */ + private static void split(final double d, final double split[]) { + if (d < 8e298 && d > -8e298) { + final double a = d * HEX_40000000; + split[0] = (d + a) - a; + split[1] = d - split[0]; + } else { + final double a = d * 9.31322574615478515625E-10; + split[0] = (d + a - d) * HEX_40000000; + split[1] = d - split[0]; + } + } + + /** Recompute a split. + * @param a input/out array containing the split, changed + * on output + */ + private static void resplit(final double a[]) { + final double c = a[0] + a[1]; + final double d = -(c - a[0] - a[1]); + + if (c < 8e298 && c > -8e298) { // MAGIC NUMBER + double z = c * HEX_40000000; + a[0] = (c + z) - z; + a[1] = c - a[0] + d; + } else { + double z = c * 9.31322574615478515625E-10; + a[0] = (c + z - c) * HEX_40000000; + a[1] = c - a[0] + d; + } + } + + /** Multiply two numbers in split form. + * @param a first term of multiplication + * @param b second term of multiplication + * @param ans placeholder where to put the result + */ + private static void splitMult(double a[], double b[], double ans[]) { + ans[0] = a[0] * b[0]; + ans[1] = a[0] * b[1] + a[1] * b[0] + a[1] * b[1]; + + /* Resplit */ + resplit(ans); + } + + /** Add two numbers in split form. + * @param a first term of addition + * @param b second term of addition + * @param ans placeholder where to put the result + */ + private static void splitAdd(final double a[], final double b[], final double ans[]) { + ans[0] = a[0] + b[0]; + ans[1] = a[1] + b[1]; + + resplit(ans); + } + + /** Compute the reciprocal of in. Use the following algorithm. + * in = c + d. + * want to find x + y such that x+y = 1/(c+d) and x is much + * larger than y and x has several zero bits on the right. + * + * Set b = 1/(2^22), a = 1 - b. Thus (a+b) = 1. + * Use following identity to compute (a+b)/(c+d) + * + * (a+b)/(c+d) = a/c + (bc - ad) / (c^2 + cd) + * set x = a/c and y = (bc - ad) / (c^2 + cd) + * This will be close to the right answer, but there will be + * some rounding in the calculation of X. So by carefully + * computing 1 - (c+d)(x+y) we can compute an error and + * add that back in. This is done carefully so that terms + * of similar size are subtracted first. + * @param in initial number, in split form + * @param result placeholder where to put the result + */ + static void splitReciprocal(final double in[], final double result[]) { + final double b = 1.0/4194304.0; + final double a = 1.0 - b; + + if (in[0] == 0.0) { + in[0] = in[1]; + in[1] = 0.0; + } + + result[0] = a / in[0]; + result[1] = (b*in[0]-a*in[1]) / (in[0]*in[0] + in[0]*in[1]); + + if (result[1] != result[1]) { // can happen if result[1] is NAN + result[1] = 0.0; + } + + /* Resplit */ + resplit(result); + + for (int i = 0; i < 2; i++) { + /* this may be overkill, probably once is enough */ + double err = 1.0 - result[0] * in[0] - result[0] * in[1] - + result[1] * in[0] - result[1] * in[1]; + /*err = 1.0 - err; */ + err = err * (result[0] + result[1]); + /*printf("err = %16e\n", err); */ + result[1] += err; + } + } + + /** Compute (a[0] + a[1]) * (b[0] + b[1]) in extended precision. + * @param a first term of the multiplication + * @param b second term of the multiplication + * @param result placeholder where to put the result + */ + private static void quadMult(final double a[], final double b[], final double result[]) { + final double xs[] = new double[2]; + final double ys[] = new double[2]; + final double zs[] = new double[2]; + + /* a[0] * b[0] */ + split(a[0], xs); + split(b[0], ys); + splitMult(xs, ys, zs); + + result[0] = zs[0]; + result[1] = zs[1]; + + /* a[0] * b[1] */ + split(b[1], ys); + splitMult(xs, ys, zs); + + double tmp = result[0] + zs[0]; + result[1] = result[1] - (tmp - result[0] - zs[0]); + result[0] = tmp; + tmp = result[0] + zs[1]; + result[1] = result[1] - (tmp - result[0] - zs[1]); + result[0] = tmp; + + /* a[1] * b[0] */ + split(a[1], xs); + split(b[0], ys); + splitMult(xs, ys, zs); + + tmp = result[0] + zs[0]; + result[1] = result[1] - (tmp - result[0] - zs[0]); + result[0] = tmp; + tmp = result[0] + zs[1]; + result[1] = result[1] - (tmp - result[0] - zs[1]); + result[0] = tmp; + + /* a[1] * b[0] */ + split(a[1], xs); + split(b[1], ys); + splitMult(xs, ys, zs); + + tmp = result[0] + zs[0]; + result[1] = result[1] - (tmp - result[0] - zs[0]); + result[0] = tmp; + tmp = result[0] + zs[1]; + result[1] = result[1] - (tmp - result[0] - zs[1]); + result[0] = tmp; + } + + /** Compute exp(p) for a integer p in extended precision. + * @param p integer whose exponential is requested + * @param result placeholder where to put the result in extended precision + * @return exp(p) in standard precision (equal to result[0] + result[1]) + */ + static double expint(int p, final double result[]) { + //double x = M_E; + final double xs[] = new double[2]; + final double as[] = new double[2]; + final double ys[] = new double[2]; + //split(x, xs); + //xs[1] = (double)(2.7182818284590452353602874713526625L - xs[0]); + //xs[0] = 2.71827697753906250000; + //xs[1] = 4.85091998273542816811e-06; + //xs[0] = Double.longBitsToDouble(0x4005bf0800000000L); + //xs[1] = Double.longBitsToDouble(0x3ed458a2bb4a9b00L); + + /* E */ + xs[0] = 2.718281828459045; + xs[1] = 1.4456468917292502E-16; + + split(1.0, ys); + + while (p > 0) { + if ((p & 1) != 0) { + quadMult(ys, xs, as); + ys[0] = as[0]; ys[1] = as[1]; + } + + quadMult(xs, xs, as); + xs[0] = as[0]; xs[1] = as[1]; + + p >>= 1; + } + + if (result != null) { + result[0] = ys[0]; + result[1] = ys[1]; + + resplit(result); + } + + return ys[0] + ys[1]; + } + /** xi in the range of [1, 2]. + * 3 5 7 + * x+1 / x x x \ + * ln ----- = 2 * | x + ---- + ---- + ---- + ... | + * 1-x \ 3 5 7 / + * + * So, compute a Remez approximation of the following function + * + * ln ((sqrt(x)+1)/(1-sqrt(x))) / x + * + * This will be an even function with only positive coefficents. + * x is in the range [0 - 1/3]. + * + * Transform xi for input to the above function by setting + * x = (xi-1)/(xi+1). Input to the polynomial is x^2, then + * the result is multiplied by x. + * @param xi number from which log is requested + * @return log(xi) + */ + static double[] slowLog(double xi) { + double x[] = new double[2]; + double x2[] = new double[2]; + double y[] = new double[2]; + double a[] = new double[2]; + + split(xi, x); + + /* Set X = (x-1)/(x+1) */ + x[0] += 1.0; + resplit(x); + splitReciprocal(x, a); + x[0] -= 2.0; + resplit(x); + splitMult(x, a, y); + x[0] = y[0]; + x[1] = y[1]; + + /* Square X -> X2*/ + splitMult(x, x, x2); + + + //x[0] -= 1.0; + //resplit(x); + + y[0] = LN_SPLIT_COEF[LN_SPLIT_COEF.length-1][0]; + y[1] = LN_SPLIT_COEF[LN_SPLIT_COEF.length-1][1]; + + for (int i = LN_SPLIT_COEF.length-2; i >= 0; i--) { + splitMult(y, x2, a); + y[0] = a[0]; + y[1] = a[1]; + splitAdd(y, LN_SPLIT_COEF[i], a); + y[0] = a[0]; + y[1] = a[1]; + } + + splitMult(y, x, a); + y[0] = a[0]; + y[1] = a[1]; + + return y; + } + + + static void printarray(String string, int expectedLen, double[][] array2d) { + System.out.println(string); + checkLen(expectedLen, array2d.length); + System.out.println(" { "); + int i = 0; + for(double array[] : array2d) { + System.out.print(" {"); + for(double d : array) { // assume inner array has very few entries + String ds = d >= 0 ? "+"+Double.toString(d)+"d," : Double.toString(d)+"d,"; + System.out.printf("%-25.25s",ds); // multiple entries per line + } + System.out.println("}, // "+i++); + } + System.out.println(" };"); + } + + static void printarray(String string, int expectedLen, double[] array) { + System.out.println(string+"="); + checkLen(expectedLen, array.length); + System.out.println(" {"); + for(double d : array){ + String ds = d!=d ? "Double.NaN," : d >= 0 ? "+"+Double.toString(d)+"d," : Double.toString(d)+"d,"; + System.out.printf(" %s%n",ds); // one entry per line + } + System.out.println(" };"); + } + + private static void checkLen(int expectedLen, int actual) { + if (expectedLen != actual) { + throw new DimensionMismatchException(actual, expectedLen); + } + } + +}