FastMathCalc.java
/*
* 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
*
* https://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.
*/
/*
* This is not the original file distributed by the Apache Software Foundation
* It has been modified by the Hipparchus project
*/
package org.hipparchus.util;
import java.io.PrintStream;
/**
* Class used to compute the classical functions tables.
*/
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 = {
+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},
};
/** Table start declaration. */
private static final String TABLE_START_DECL = " {";
/** Table end declaration. */
private static final String TABLE_END_DECL = " };";
/**
* Private Constructor.
*/
private FastMathCalc() {
}
/** Build the sine and cosine tables.
* @param SINE_TABLE_A table of the most significant part of the sines
* @param SINE_TABLE_B table of the least significant part of the sines
* @param COSINE_TABLE_A table of the most significant part of the cosines
* @param COSINE_TABLE_B table of the most significant part of the cosines
* @param SINE_TABLE_LEN length of the tables
* @param TANGENT_TABLE_A table of the most significant part of the tangents
* @param TANGENT_TABLE_B table of the most significant part of the tangents
*/
@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 *= 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] -= tmp - result[0] - zs[0];
result[0] = tmp;
tmp = result[0] + zs[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] -= tmp - result[0] - zs[0];
result[0] = tmp;
tmp = result[0] + zs[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] -= tmp - result[0] - zs[0];
result[0] = tmp;
tmp = result[0] + zs[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;
}
/**
* Print an array.
* @param out text output stream where output should be printed
* @param name array name
* @param expectedLen expected length of the array
* @param array2d array data
*/
static void printarray(PrintStream out, String name, int expectedLen, double[][] array2d) {
out.println(name);
MathUtils.checkDimension(expectedLen, array2d.length);
out.println(TABLE_START_DECL + " ");
int i = 0;
for(double[] array : array2d) { // "double array[]" causes PMD parsing error
out.print(" {");
for(double d : array) { // assume inner array has very few entries
out.printf("%-25.25s", format(d)); // multiple entries per line
}
out.println("}, // " + i++);
}
out.println(TABLE_END_DECL);
}
/**
* Print an array.
* @param out text output stream where output should be printed
* @param name array name
* @param expectedLen expected length of the array
* @param array array data
*/
static void printarray(PrintStream out, String name, int expectedLen, double[] array) {
out.println(name + "=");
MathUtils.checkDimension(expectedLen, array.length);
out.println(TABLE_START_DECL);
for(double d : array){
out.printf(" %s%n", format(d)); // one entry per line
}
out.println(TABLE_END_DECL);
}
/** Format a double.
* @param d double number to format
* @return formatted number
*/
static String format(double d) {
if (d != d) {
return "Double.NaN,";
} else {
return ((d >= 0) ? "+" : "") + d + "d,";
}
}
}