commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From s...@apache.org
Subject svn commit: r1169527 - in /commons/proper/math/trunk/src/main/java/org/apache/commons/math/util: FastMath.java FastMathCalc.java
Date Sun, 11 Sep 2011 20:43:27 GMT
Author: sebb
Date: Sun Sep 11 20:43:27 2011
New Revision: 1169527

URL: http://svn.apache.org/viewvc?rev=1169527&view=rev
Log:
MATH-650 move static setup methods to helper class

Added:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMathCalc.java
  (with props)
Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMath.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMath.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMath.java?rev=1169527&r1=1169526&r2=1169527&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMath.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMath.java Sun
Sep 11 20:43:27 2011
@@ -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);
     }
 
     /**
@@ -7334,251 +7254,6 @@ public class FastMath {
     }
 
     /**
-     *  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.
      *
      * @param x   a double
@@ -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.

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMathCalc.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMathCalc.java?rev=1169527&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMathCalc.java
(added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMathCalc.java
Sun Sep 11 20:43:27 2011
@@ -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);
+        }
+    }
+
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMathCalc.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMathCalc.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision



Mime
View raw message