Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

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
This commit is contained in:
Sebastian Bazley 2011-09-11 20:43:27 +00:00
parent 46b35abc2c
commit 806fa1cdd2
2 changed files with 625 additions and 588 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

603
src/main/java/org/apache/commons/math/util/FastMath.java
View File

@ -16,8 +16,6 @@
*/ */
package org.apache.commons.math.util; package org.apache.commons.math.util;
import org.apache.commons.math.exception.DimensionMismatchException;
/** /**
* Faster, more accurate, portable alternative to {@link Math} and * Faster, more accurate, portable alternative to {@link Math} and
* {@link StrictMath} for large scale computation. * {@link StrictMath} for large scale computation.
@ -121,13 +119,13 @@ public class FastMath {
// Populate expIntTable // Populate expIntTable
for (int i = 0; i < FastMath.EXP_INT_TABLE_MAX_INDEX; i++) { 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_A[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[0];
EXP_INT_TABLE_B[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[1]; EXP_INT_TABLE_B[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[1];
if (i != 0) { if (i != 0) {
// Negative integer powers // 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_A[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[0];
EXP_INT_TABLE_B[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[1]; EXP_INT_TABLE_B[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[1];
} }
@ -3167,7 +3165,7 @@ public class FastMath {
// Populate expFracTable // Populate expFracTable
for (int i = 0; i < EXP_FRAC_TABLE_A.length; i++) { 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_A[i] = tmp[0];
EXP_FRAC_TABLE_B[i] = tmp[1]; 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 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 // 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 // Populate lnMant table
for (int i = 0; i < LN_MANT.length; i++) { for (int i = 0; i < LN_MANT.length; i++) {
final double d = Double.longBitsToDouble( (((long) i) << 42) | 0x3ff0000000000000L ); final double d = Double.longBitsToDouble( (((long) i) << 42) | 0x3ff0000000000000L );
LN_MANT[i] = slowLog(d); LN_MANT[i] = FastMathCalc.slowLog(d);
} }
} else { } else {
LN_MANT = new double[][] { LN_MANT = new double[][] {
@ -6309,26 +6282,6 @@ public class FastMath {
/** log(2) (low bits). */ /** log(2) (low bits). */
private static final double LN_2_B = 1.17304635250823482e-7; 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. */ /** Coefficients for log, when input 0.99 < x < 1.01. */
private static final double LN_QUICK_COEF[][] = { private static final double LN_QUICK_COEF[][] = {
{1.0, 5.669184079525E-24}, {1.0, 5.669184079525E-24},
@ -6539,50 +6492,17 @@ public class FastMath {
// } // }
public static void main(String[] a){ public static void main(String[] a){
printarray("EXP_INT_TABLE_A", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_A); FastMathCalc.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); FastMathCalc.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); FastMathCalc.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); FastMathCalc.printarray("EXP_FRAC_TABLE_B", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_B);
printarray("LN_MANT",LN_MANT_LEN, lnMant.LN_MANT); FastMathCalc.printarray("LN_MANT",LN_MANT_LEN, lnMant.LN_MANT);
printarray("SINE_TABLE_A", SINE_TABLE_LEN, SINE_TABLE_A); FastMathCalc.printarray("SINE_TABLE_A", SINE_TABLE_LEN, SINE_TABLE_A);
printarray("SINE_TABLE_B", SINE_TABLE_LEN, SINE_TABLE_B); FastMathCalc.printarray("SINE_TABLE_B", SINE_TABLE_LEN, SINE_TABLE_B);
printarray("COSINE_TABLE_A", SINE_TABLE_LEN, COSINE_TABLE_A); FastMathCalc.printarray("COSINE_TABLE_A", SINE_TABLE_LEN, COSINE_TABLE_A);
printarray("COSINE_TABLE_B", SINE_TABLE_LEN, COSINE_TABLE_B); FastMathCalc.printarray("COSINE_TABLE_B", SINE_TABLE_LEN, COSINE_TABLE_B);
printarray("TANGENT_TABLE_A", SINE_TABLE_LEN, TANGENT_TABLE_A); FastMathCalc.printarray("TANGENT_TABLE_A", SINE_TABLE_LEN, TANGENT_TABLE_A);
printarray("TANGENT_TABLE_B", SINE_TABLE_LEN, TANGENT_TABLE_B); FastMathCalc.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);
}
} }
/** /**
@ -7333,251 +7253,6 @@ public class FastMath {
return ya + yb; 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. * Natural logarithm.
* *
@ -8048,254 +7723,6 @@ public class FastMath {
return result; 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. * Computes sin(x) - x, where |x| < 1/16.

610
src/main/java/org/apache/commons/math/util/FastMathCalc.java Normal file
View File

@ -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);
}
}
}