math.JafamaFastMath Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of finwhale Show documentation
Show all versions of finwhale Show documentation
Statistical distributions library (in statu nascendi)
/*
* Copyright 2012-2013 Jeff Hain
*
* Licensed 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.
*/
/*
* =============================================================================
* Notice of fdlibm package this program is partially derived from:
*
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* =============================================================================
*/
package math;
/**
* Class providing math treatments that:
* - are meant to be faster than java.lang.Math class equivalents,
* - are still somehow accurate and robust (handling of NaN and such),
* - do not (or not directly) generate objects at run time (no "new").
*
* Use of look-up tables: around 1 MiB total, and initialized on class load.
*
* Depending on JVM, or JVM options, these treatments can actually be slower
* than Math ones.
* In particular, they can be slower if not optimized by the JIT, which you
* can see with -Xint JVM option.
* Another cause of slowness can be cache-misses on look-up tables.
* Also, look-up tables initialization, done on class load, typically
* takes multiple hundreds of milliseconds (and is about twice slower
* in J6 than in J5, and in J7 than in J6, possibly due to intrinsifications
* preventing optimizations such as use of hardware sqrt, and Math delegating
* to StrictMath with JIT optimizations not yet up during class load).
*
* These treatments are not strictfp: if you want identical results
* across various architectures, you must roll your own implementation
* by adding strictfp and using StrictMath instead of Math.
*
* Methods with same signature than Math ones, are meant to return
* "good" approximations on all range.
*
* --- words, words, words ---
*
* "0x42BE0000 percents of the folks out there
* are completely clueless about floating-point."
*
* The difference between precision and accuracy:
* "3.177777777777777 is a precise (16 digits)
* but inaccurate (only correct up to the second digit)
* approximation of PI=3.141592653589793(etc.)."
*/
final class JafamaFastMath {
/*
* For trigonometric functions, use of look-up tables and Taylor-Lagrange formula
* with 4 derivatives (more take longer to compute and don't add much accuracy,
* less require larger tables (which use more memory, take more time to initialize,
* and are slower to access (at least on the machine they were developed on))).
*
* For angles reduction of cos/sin/tan functions:
* - for small values, instead of reducing angles, and then computing the best index
* for look-up tables, we compute this index right away, and use it for reduction,
* - for large values, treatments derived from fdlibm package are used, as done in
* java.lang.Math. They are faster but still "slow", so if you work with
* large numbers and need speed over accuracy for them, you might want to use
* normalizeXXXFast treatments before your function, or modify cos/sin/tan
* so that they call the fast normalization treatments instead of the accurate ones.
* NB: If an angle is huge (like PI*1e20), in double precision format its last digits
* are zeros, which most likely is not the case for the intended value, and doing
* an accurate reduction on a very inaccurate value is most likely pointless.
* But it gives some sort of coherence that could be needed in some cases.
*
* Multiplication on double appears to be about as fast (or not much slower) than call
* to [], and regrouping some doubles in a private class, to use
* index only once, does not seem to speed things up, so:
* - for uniformly tabulated values, to retrieve the parameter corresponding to
* an index, we recompute it rather than using an array to store it,
* - for cos/sin, we recompute derivatives divided by (multiplied by inverse of)
* factorial each time, rather than storing them in arrays.
*
* Lengths of look-up tables are usually of the form 2^n+1, for their values to be
* of the form ( * k/2^n, k in 0 .. 2^n), so that particular values
* (PI/2, etc.) are "exactly" computed, as well as for other reasons.
*
* Most math treatments I could find on the web, including "fast" ones,
* usually take care of special cases (NaN, etc.) at the beginning, and
* then deal with the general case, which adds a useless overhead for the
* general (and common) case. In this class, special cases are only dealt
* with when needed, and if the general case does not already handle them.
*/
//--------------------------------------------------------------------------
// CONFIGURATION
//--------------------------------------------------------------------------
/**
* Using two pow tab can just make things barely faster,
* but could relatively hurt in case of cache-misses,
* especially for methods that otherwise wouldn't rely
* on any tab, so we don't use it.
*/
private static final boolean USE_TWO_POW_TAB = false;
//--------------------------------------------------------------------------
// GENERAL CONSTANTS
//--------------------------------------------------------------------------
/**
* High approximation of PI, which is further from PI
* than the low approximation Math.PI:
* PI ~= 3.14159265358979323846...
* Math.PI ~= 3.141592653589793
* JafamaFastMath.PI_SUP ~= 3.1415926535897936
*/
public static final double PI_SUP = Double.longBitsToDouble(Double.doubleToRawLongBits(Math.PI)+1);
private static final double ONE_DIV_F2 = 1/2.0;
private static final double ONE_DIV_F3 = 1/6.0;
private static final double ONE_DIV_F4 = 1/24.0;
private static final double TWO_POW_24 = twoPow(24);
private static final double TWO_POW_N24 = twoPow(-24);
private static final double TWO_POW_66 = twoPow(66);
private static final double TWO_POW_450 = twoPow(450);
private static final double TWO_POW_N450 = twoPow(-450);
private static final double TWO_POW_750 = twoPow(750);
private static final double TWO_POW_N750 = twoPow(-750);
// Not storing float/double mantissa size in constants,
// for 23 and 52 are shorter to read and more
// bitwise-explicit than some constant's name.
private static final int MIN_DOUBLE_EXPONENT = -1074;
private static final int MAX_DOUBLE_EXPONENT = 1023;
private static final double LOG_2 = StrictMath.log(2.0);
//--------------------------------------------------------------------------
// CONSTANTS FOR NORMALIZATIONS
//--------------------------------------------------------------------------
/*
* Table of constants for 1/(2*PI), 282 Hex digits (enough for normalizing doubles).
* 1/(2*PI) approximation = sum of ONE_OVER_TWOPI_TAB[i]*2^(-24*(i+1)).
*/
private static final double ONE_OVER_TWOPI_TAB[] = {
0x28BE60, 0xDB9391, 0x054A7F, 0x09D5F4, 0x7D4D37, 0x7036D8,
0xA5664F, 0x10E410, 0x7F9458, 0xEAF7AE, 0xF1586D, 0xC91B8E,
0x909374, 0xB80192, 0x4BBA82, 0x746487, 0x3F877A, 0xC72C4A,
0x69CFBA, 0x208D7D, 0x4BAED1, 0x213A67, 0x1C09AD, 0x17DF90,
0x4E6475, 0x8E60D4, 0xCE7D27, 0x2117E2, 0xEF7E4A, 0x0EC7FE,
0x25FFF7, 0x816603, 0xFBCBC4, 0x62D682, 0x9B47DB, 0x4D9FB3,
0xC9F2C2, 0x6DD3D1, 0x8FD9A7, 0x97FA8B, 0x5D49EE, 0xB1FAF9,
0x7C5ECF, 0x41CE7D, 0xE294A4, 0xBA9AFE, 0xD7EC47};
/*
* Constants for 2*PI. Only the 23 most significant bits of each mantissa are used.
* 2*PI approximation = sum of TWOPI_TAB.
*/
private static final double TWOPI_TAB0 = Double.longBitsToDouble(0x401921FB40000000L);
private static final double TWOPI_TAB1 = Double.longBitsToDouble(0x3E94442D00000000L);
private static final double TWOPI_TAB2 = Double.longBitsToDouble(0x3D18469880000000L);
private static final double TWOPI_TAB3 = Double.longBitsToDouble(0x3B98CC5160000000L);
private static final double TWOPI_TAB4 = Double.longBitsToDouble(0x3A101B8380000000L);
private static final double INVPIO2 = Double.longBitsToDouble(0x3FE45F306DC9C883L); // 6.36619772367581382433e-01 53 bits of 2/pi
private static final double PIO2_HI = Double.longBitsToDouble(0x3FF921FB54400000L); // 1.57079632673412561417e+00 first 33 bits of pi/2
private static final double PIO2_LO = Double.longBitsToDouble(0x3DD0B4611A626331L); // 6.07710050650619224932e-11 pi/2 - PIO2_HI
private static final double INVTWOPI = INVPIO2/4;
private static final double TWOPI_HI = 4*PIO2_HI;
private static final double TWOPI_LO = 4*PIO2_LO;
// fdlibm uses 2^19*PI/2 here, but we normalize with % 2*PI instead of % PI/2,
// and we can bear some more error.
private static final double NORMALIZE_ANGLE_MAX_MEDIUM_DOUBLE = StrictMath.pow(2,20)*(2*Math.PI);
//--------------------------------------------------------------------------
// CONSTANTS AND TABLES FOR SIN AND COS
//--------------------------------------------------------------------------
private static final int SIN_COS_TABS_SIZE = (1<>9) / SIN_COS_INDEXER) * 0.99;
//--------------------------------------------------------------------------
// CONSTANTS AND TABLES FOR ACOS, ASIN
//--------------------------------------------------------------------------
// We use the following formula:
// 1) acos(x) = PI/2 - asin(x)
// 2) asin(-x) = -asin(x)
// ---> we only have to compute asin(x) on [0,1].
// For values not close to +-1, we use look-up tables;
// for values near +-1, we use code derived from fdlibm.
// Supposed to be >= sin(77.2deg), as fdlibm code is supposed to work with values > 0.975,
// but seems to work well enough as long as value >= sin(25deg).
private static final double ASIN_MAX_VALUE_FOR_TABS = StrictMath.sin(Math.toRadians(73.0));
private static final int ASIN_TABS_SIZE = (1< we only have to compute atan(x) on [0,+infinity[.
// For values corresponding to angles not close to +-PI/2, we use look-up tables;
// for values corresponding to angles near +-PI/2, we use code derived from fdlibm.
// Supposed to be >= tan(67.7deg), as fdlibm code is supposed to work with values > 2.4375.
private static final double ATAN_MAX_VALUE_FOR_TABS = StrictMath.tan(Math.toRadians(74.0));
private static final int ATAN_TABS_SIZE = (1< SIN_COS_MAX_VALUE_FOR_INT_MODULO) {
// Faster than using normalizeZeroTwoPi.
angle = remainderTwoPi(angle);
if (angle < 0.0) {
angle += 2*Math.PI;
}
}
int index = (int)(angle * SIN_COS_INDEXER + 0.5);
double delta = (angle - index * SIN_COS_DELTA_HI) - index * SIN_COS_DELTA_LO;
index &= (SIN_COS_TABS_SIZE-2); // index % (SIN_COS_TABS_SIZE-1)
double indexSin = sinTab[index];
double indexCos = cosTab[index];
double result = indexSin + delta * (indexCos + delta * (-indexSin * ONE_DIV_F2 + delta * (-indexCos * ONE_DIV_F3 + delta * indexSin * ONE_DIV_F4)));
return negateResult ? -result : result;
}
/**
* @param angle Angle in radians.
* @return Angle cosine.
*/
static double cos(double angle) {
angle = Math.abs(angle);
if (angle > SIN_COS_MAX_VALUE_FOR_INT_MODULO) {
// Faster than using normalizeZeroTwoPi.
angle = remainderTwoPi(angle);
if (angle < 0.0) {
angle += 2*Math.PI;
}
}
// index: possibly outside tables range.
int index = (int)(angle * SIN_COS_INDEXER + 0.5);
double delta = (angle - index * SIN_COS_DELTA_HI) - index * SIN_COS_DELTA_LO;
// Making sure index is within tables range.
// Last value of each table is the same than first, so we ignore it (tabs size minus one) for modulo.
index &= (SIN_COS_TABS_SIZE-2); // index % (SIN_COS_TABS_SIZE-1)
double indexCos = cosTab[index];
double indexSin = sinTab[index];
return indexCos + delta * (-indexSin + delta * (-indexCos * ONE_DIV_F2 + delta * (indexSin * ONE_DIV_F3 + delta * indexCos * ONE_DIV_F4)));
}
/**
* @param value Value in [-1,1].
* @return Value arcsine, in radians, in [-PI/2,PI/2].
*/
static double asin(double value) {
boolean negateResult;
if (value < 0.0) {
value = -value;
negateResult = true;
} else {
negateResult = false;
}
if (value <= ASIN_MAX_VALUE_FOR_TABS) {
int index = (int)(value * ASIN_INDEXER + 0.5);
double delta = value - index * ASIN_DELTA;
double result = asinTab[index] + delta * (asinDer1DivF1Tab[index] + delta * (asinDer2DivF2Tab[index] + delta * (asinDer3DivF3Tab[index] + delta * asinDer4DivF4Tab[index])));
return negateResult ? -result : result;
} else { // value > ASIN_MAX_VALUE_FOR_TABS, or value is NaN
// This part is derived from fdlibm.
if (value < 1.0) {
double t = (1.0 - value)*0.5;
double p = t*(ASIN_PS0+t*(ASIN_PS1+t*(ASIN_PS2+t*(ASIN_PS3+t*(ASIN_PS4+t*ASIN_PS5)))));
double q = 1.0+t*(ASIN_QS1+t*(ASIN_QS2+t*(ASIN_QS3+t*ASIN_QS4)));
double s = Math.sqrt(t);
double z = s+s*(p/q);
double result = ASIN_PIO2_HI-((z+z)-ASIN_PIO2_LO);
return negateResult ? -result : result;
} else { // value >= 1.0, or value is NaN
if (value == 1.0) {
return negateResult ? -Math.PI/2 : Math.PI/2;
} else {
return Double.NaN;
}
}
}
}
/**
* @param value Value in [-1,1].
* @return Value arccosine, in radians, in [0,PI].
*/
static double acos(double value) {
return Math.PI/2 - JafamaFastMath.asin(value);
}
/**
* @param value A double value.
* @return Value arctangent, in radians, in [-PI/2,PI/2].
*/
static double atan(double value) {
boolean negateResult;
if (value < 0.0) {
value = -value;
negateResult = true;
} else {
negateResult = false;
}
if (value == 1.0) {
// We want "exact" result for 1.0.
return negateResult ? -Math.PI/4 : Math.PI/4;
} else if (value <= ATAN_MAX_VALUE_FOR_TABS) {
int index = (int)(value * ATAN_INDEXER + 0.5);
double delta = value - index * ATAN_DELTA;
double result = atanTab[index] + delta * (atanDer1DivF1Tab[index] + delta * (atanDer2DivF2Tab[index] + delta * (atanDer3DivF3Tab[index] + delta * atanDer4DivF4Tab[index])));
return negateResult ? -result : result;
} else { // value > ATAN_MAX_VALUE_FOR_TABS, or value is NaN
// This part is derived from fdlibm.
if (value < TWO_POW_66) {
double x = -1/value;
double x2 = x*x;
double x4 = x2*x2;
double s1 = x2*(ATAN_AT0+x4*(ATAN_AT2+x4*(ATAN_AT4+x4*(ATAN_AT6+x4*(ATAN_AT8+x4*ATAN_AT10)))));
double s2 = x4*(ATAN_AT1+x4*(ATAN_AT3+x4*(ATAN_AT5+x4*(ATAN_AT7+x4*ATAN_AT9))));
double result = ATAN_HI3-((x*(s1+s2)-ATAN_LO3)-x);
return negateResult ? -result : result;
} else { // value >= 2^66, or value is NaN
if (value != value) {
return Double.NaN;
} else {
return negateResult ? -Math.PI/2 : Math.PI/2;
}
}
}
}
/**
* For special values for which multiple conventions could be adopted,
* behaves like Math.atan2(double,double).
*
* @param y Coordinate on y axis.
* @param x Coordinate on x axis.
* @return Angle from x axis positive side to (x,y) position, in radians, in [-PI,PI].
* Angle measure is positive when going from x axis to y axis (positive sides).
*/
static double atan2(double y, double x) {
/*
* Using sub-methods, to make method lighter for general case,
* and to avoid JIT-optimization crash on NaN.
*/
if (x > 0.0) {
if (y == 0.0) {
// +-0.0
return y;
}
if (x == Double.POSITIVE_INFINITY) {
return atan2_pinf_yyy(y);
} else {
return JafamaFastMath.atan(y/x);
}
} else if (x < 0.0) {
if (y == 0.0) {
return JafamaFastMath.signFromBit(y) * Math.PI;
}
if (x == Double.NEGATIVE_INFINITY) {
return atan2_ninf_yyy(y);
} else if (y > 0.0) {
return Math.PI/2 + JafamaFastMath.atan(-x/y);
} else if (y < 0.0) {
return -Math.PI/2 - JafamaFastMath.atan(x/y);
} else {
return Double.NaN;
}
} else {
return atan2_zeroOrNaN_yyy(x, y);
}
}
/*
* logarithms
*/
/**
* Much more accurate than log(1+value),
* for arguments (and results) close to zero.
*
* @param value A double value.
* @return Logarithm (base e) of (1+value).
*/
static double log1p(double value) {
if (value > -1.0) {
if (value == Double.POSITIVE_INFINITY) {
return Double.POSITIVE_INFINITY;
}
// ln'(x) = 1/x
// so
// log(x+epsilon) ~= log(x) + epsilon/x
//
// Let u be 1+value rounded:
// 1+value = u+epsilon
//
// log(1+value)
// = log(u+epsilon)
// ~= log(u) + epsilon/value
// We compute log(u) as done in log(double), and then add the corrective term.
double valuePlusOne = 1.0+value;
if (valuePlusOne == 1.0) {
return value;
} else if (Math.abs(value) < 0.15) {
double z = value/(value+2.0);
double z2 = z*z;
return z*(2+z2*((2.0/3)+z2*((2.0/5)+z2*((2.0/7)+z2*((2.0/9)+z2*((2.0/11)))))));
}
int valuePlusOneBitsHi = (int)(Double.doubleToRawLongBits(valuePlusOne)>>32) & 0x7FFFFFFF;
int valuePlusOneExp = (valuePlusOneBitsHi>>20)-MAX_DOUBLE_EXPONENT;
// Getting the first LOG_BITS bits of the mantissa.
int xIndex = ((valuePlusOneBitsHi<<12)>>>(32-LOG_BITS));
// 1.mantissa/1.mantissaApprox - 1
double z = (valuePlusOne * twoPowNormalOrSubnormal(-valuePlusOneExp)) * logXInvTab[xIndex] - 1;
z *= (1-z*((1.0/2)-z*(1.0/3)));
// Adding epsilon/valuePlusOne to z,
// with
// epsilon = value - (valuePlusOne-1)
// (valuePlusOne + epsilon ~= 1+value (not rounded))
return valuePlusOneExp * LOG_2 + logXLogTab[xIndex] + (z + (value - (valuePlusOne-1))/valuePlusOne);
} else if (value == -1.0) {
return Double.NEGATIVE_INFINITY;
} else { // value < -1.0, or value is NaN
return Double.NaN;
}
}
/**
* Returns sqrt(x^2+y^2) without intermediate overflow or underflow.
*/
static double hypot(double x, double y) {
x = Math.abs(x);
y = Math.abs(y);
if (y < x) {
double a = x;
x = y;
y = a;
} else if (!(y >= x)) { // Testing if we have some NaN.
if ((x == Double.POSITIVE_INFINITY) || (y == Double.POSITIVE_INFINITY)) {
return Double.POSITIVE_INFINITY;
} else {
return Double.NaN;
}
}
if (y-x == y) { // x too small to substract from y
return y;
} else {
double factor;
if (x > TWO_POW_450) { // 2^450 < x < y
x *= TWO_POW_N750;
y *= TWO_POW_N750;
factor = TWO_POW_750;
} else if (y < TWO_POW_N450) { // x < y < 2^-450
x *= TWO_POW_750;
y *= TWO_POW_750;
factor = TWO_POW_N750;
} else {
factor = 1.0;
}
return factor * Math.sqrt(x*x+y*y);
}
}
//--------------------------------------------------------------------------
// PRIVATE TREATMENTS
//--------------------------------------------------------------------------
/**
* JafamaFastMath is non-instantiable.
*/
private JafamaFastMath() {
}
/**
* Use look-up tables size power through this method,
* to make sure is it small in case java.lang.Math
* is directly used.
*/
private static int getTabSizePower(int tabSizePower) {
return tabSizePower;
}
/**
* @param value A double value.
* @return -1 if sign bit if 1, 1 if sign bit if 0.
*/
private static long signFromBit(double value) {
// Returning a long, to avoid useless cast into int.
return ((Double.doubleToRawLongBits(value)>>62)|1);
}
/**
* Returns the exact result, provided it's in double range,
* i.e. if power is in [-1074,1023].
*
* @param power A power.
* @return 2^power.
*/
private static double twoPow(int power) {
if (power <= -MAX_DOUBLE_EXPONENT) { // Not normal.
if (power >= MIN_DOUBLE_EXPONENT) { // Subnormal.
return Double.longBitsToDouble(0x0008000000000000L>>(-(power+MAX_DOUBLE_EXPONENT)));
} else { // Underflow.
return 0.0;
}
} else if (power > MAX_DOUBLE_EXPONENT) { // Overflow.
return Double.POSITIVE_INFINITY;
} else { // Normal.
return Double.longBitsToDouble(((long)(power+MAX_DOUBLE_EXPONENT))<<52);
}
}
/**
* @param power Must be in normal values range.
*/
private static double twoPowNormal(int power) {
if (USE_TWO_POW_TAB) {
return twoPowTab[power-MIN_DOUBLE_EXPONENT];
} else {
return Double.longBitsToDouble(((long)(power+MAX_DOUBLE_EXPONENT))<<52);
}
}
/**
* @param power Must be in normal or subnormal values range.
*/
private static double twoPowNormalOrSubnormal(int power) {
if (USE_TWO_POW_TAB) {
return twoPowTab[power-MIN_DOUBLE_EXPONENT];
} else {
if (power <= -MAX_DOUBLE_EXPONENT) { // Not normal.
return Double.longBitsToDouble(0x0008000000000000L>>(-(power+MAX_DOUBLE_EXPONENT)));
} else { // Normal.
return Double.longBitsToDouble(((long)(power+MAX_DOUBLE_EXPONENT))<<52);
}
}
}
private static double atan2_pinf_yyy(double y) {
if (y == Double.POSITIVE_INFINITY) {
return Math.PI/4;
} else if (y == Double.NEGATIVE_INFINITY) {
return -Math.PI/4;
} else if (y > 0.0) {
return 0.0;
} else if (y < 0.0) {
return -0.0;
} else {
return Double.NaN;
}
}
private static double atan2_ninf_yyy(double y) {
if (y == Double.POSITIVE_INFINITY) {
return 3*Math.PI/4;
} else if (y == Double.NEGATIVE_INFINITY) {
return -3*Math.PI/4;
} else if (y > 0.0) {
return Math.PI;
} else if (y < 0.0) {
return -Math.PI;
} else {
return Double.NaN;
}
}
private static double atan2_zeroOrNaN_yyy(double x, double y) {
if (x == 0.0) {
if (y == 0.0) {
if (JafamaFastMath.signFromBit(x) < 0) {
// x is -0.0
return JafamaFastMath.signFromBit(y) * Math.PI;
} else {
// +-0.0
return y;
}
}
if (y > 0.0) {
return Math.PI/2;
} else if (y < 0.0) {
return -Math.PI/2;
} else {
return Double.NaN;
}
} else {
return Double.NaN;
}
}
/**
* Remainder using an accurate definition of PI.
* Derived from a fdlibm treatment called __ieee754_rem_pio2.
*
* This method can return values slightly (like one ULP or so) outside [-Math.PI,Math.PI] range.
*
* @param angle Angle in radians.
* @return Remainder of (angle % (2*PI)), which is in [-PI,PI] range.
*/
private static double remainderTwoPi(double angle) {
boolean negateResult;
if (angle < 0.0) {
negateResult = true;
angle = -angle;
} else {
negateResult = false;
}
if (angle <= NORMALIZE_ANGLE_MAX_MEDIUM_DOUBLE) {
double fn = (int)(angle*INVTWOPI+0.5);
double result = (angle - fn*TWOPI_HI) - fn*TWOPI_LO;
return negateResult ? -result : result;
} else if (angle < Double.POSITIVE_INFINITY) {
// Reworking exponent to have a value < 2^24.
long lx = Double.doubleToRawLongBits(angle);
long exp = ((lx>>52)&0x7FF) - 1046;
double z = Double.longBitsToDouble(lx - (exp<<52));
double x0 = ((int)z);
z = (z-x0)*TWO_POW_24;
double x1 = ((int)z);
double x2 = (z-x1)*TWO_POW_24;
double result = subRemainderTwoPi(x0, x1, x2, (int)exp, (x2 == 0) ? 2 : 3);
return negateResult ? -result : result;
} else { // angle is +infinity or NaN
return Double.NaN;
}
}
/**
* Remainder using an accurate definition of PI.
* Derived from a fdlibm treatment called __kernel_rem_pio2.
*
* @param x0 Most significant part of the value, as an integer < 2^24, in double precision format. Must be >= 0.
* @param x1 Following significant part of the value, as an integer < 2^24, in double precision format.
* @param x2 Least significant part of the value, as an integer < 2^24, in double precision format.
* @param e0 Exponent of x0 (value is (2^e0)*(x0+(2^-24)*(x1+(2^-24)*x2))). Must be >= -20.
* @param nx Number of significant parts to take into account. Must be 2 or 3.
* @return Remainder of (value % (2*PI)), which is in [-PI,PI] range.
*/
private static double subRemainderTwoPi(double x0, double x1, double x2, int e0, int nx) {
int ih;
double z,fw;
double f0,f1,f2,f3,f4,f5,f6 = 0.0,f7;
double q0,q1,q2,q3,q4,q5;
int iq0,iq1,iq2,iq3,iq4;
final int jx = nx - 1; // jx in [1,2] (nx in [2,3])
// Could use a table to avoid division, but the gain isn't worth it most likely...
final int jv = (e0-3)/24; // We do not handle the case (e0-3 < -23).
int q = e0-((jv<<4)+(jv<<3))-24; // e0-24*(jv+1)
final int j = jv + 4;
if (jx == 1) {
f5 = (j >= 0) ? ONE_OVER_TWOPI_TAB[j]: 0.0;
f4 = (j >= 1) ? ONE_OVER_TWOPI_TAB[j-1]: 0.0;
f3 = (j >= 2) ? ONE_OVER_TWOPI_TAB[j-2]: 0.0;
f2 = (j >= 3) ? ONE_OVER_TWOPI_TAB[j-3]: 0.0;
f1 = (j >= 4) ? ONE_OVER_TWOPI_TAB[j-4]: 0.0;
f0 = (j >= 5) ? ONE_OVER_TWOPI_TAB[j-5]: 0.0;
q0 = x0*f1 + x1*f0;
q1 = x0*f2 + x1*f1;
q2 = x0*f3 + x1*f2;
q3 = x0*f4 + x1*f3;
q4 = x0*f5 + x1*f4;
} else { // jx == 2
f6 = (j >= 0) ? ONE_OVER_TWOPI_TAB[j]: 0.0;
f5 = (j >= 1) ? ONE_OVER_TWOPI_TAB[j-1]: 0.0;
f4 = (j >= 2) ? ONE_OVER_TWOPI_TAB[j-2]: 0.0;
f3 = (j >= 3) ? ONE_OVER_TWOPI_TAB[j-3]: 0.0;
f2 = (j >= 4) ? ONE_OVER_TWOPI_TAB[j-4]: 0.0;
f1 = (j >= 5) ? ONE_OVER_TWOPI_TAB[j-5]: 0.0;
f0 = (j >= 6) ? ONE_OVER_TWOPI_TAB[j-6]: 0.0;
q0 = x0*f2 + x1*f1 + x2*f0;
q1 = x0*f3 + x1*f2 + x2*f1;
q2 = x0*f4 + x1*f3 + x2*f2;
q3 = x0*f5 + x1*f4 + x2*f3;
q4 = x0*f6 + x1*f5 + x2*f4;
}
z = q4;
fw = ((int)(TWO_POW_N24*z));
iq0 = (int)(z-TWO_POW_24*fw);
z = q3+fw;
fw = ((int)(TWO_POW_N24*z));
iq1 = (int)(z-TWO_POW_24*fw);
z = q2+fw;
fw = ((int)(TWO_POW_N24*z));
iq2 = (int)(z-TWO_POW_24*fw);
z = q1+fw;
fw = ((int)(TWO_POW_N24*z));
iq3 = (int)(z-TWO_POW_24*fw);
z = q0+fw;
// Here, q is in [-25,2] range or so,
// so in normal exponents range.
double twoPowQ = twoPowNormal(q);
z = (z*twoPowQ) % 8.0;
z -= ((int)z);
if (q > 0) {
iq3 &= 0xFFFFFF>>q;
ih = iq3>>(23-q);
} else if (q == 0) {
ih = iq3>>23;
} else if (z >= 0.5) {
ih = 2;
} else {
ih = 0;
}
if (ih > 0) {
int carry;
if (iq0 != 0) {
carry = 1;
iq0 = 0x1000000 - iq0;
iq1 = 0x0FFFFFF - iq1;
iq2 = 0x0FFFFFF - iq2;
iq3 = 0x0FFFFFF - iq3;
} else {
if (iq1 != 0) {
carry = 1;
iq1 = 0x1000000 - iq1;
iq2 = 0x0FFFFFF - iq2;
iq3 = 0x0FFFFFF - iq3;
} else {
if (iq2 != 0) {
carry = 1;
iq2 = 0x1000000 - iq2;
iq3 = 0x0FFFFFF - iq3;
} else {
if (iq3 != 0) {
carry = 1;
iq3 = 0x1000000 - iq3;
} else {
carry = 0;
}
}
}
}
if (q > 0) {
switch (q) {
case 1:
iq3 &= 0x7FFFFF;
break;
case 2:
iq3 &= 0x3FFFFF;
break;
}
}
if (ih == 2) {
z = 1.0 - z;
if (carry != 0) {
z -= twoPowQ;
}
}
}
if (z == 0.0) {
if (jx == 1) {
f6 = ONE_OVER_TWOPI_TAB[jv+5];
q5 = x0*f6 + x1*f5;
} else { // jx == 2
f7 = ONE_OVER_TWOPI_TAB[jv+5];
q5 = x0*f7 + x1*f6 + x2*f5;
}
z = q5;
fw = ((int)(TWO_POW_N24*z));
iq0 = (int)(z-TWO_POW_24*fw);
z = q4+fw;
fw = ((int)(TWO_POW_N24*z));
iq1 = (int)(z-TWO_POW_24*fw);
z = q3+fw;
fw = ((int)(TWO_POW_N24*z));
iq2 = (int)(z-TWO_POW_24*fw);
z = q2+fw;
fw = ((int)(TWO_POW_N24*z));
iq3 = (int)(z-TWO_POW_24*fw);
z = q1+fw;
fw = ((int)(TWO_POW_N24*z));
iq4 = (int)(z-TWO_POW_24*fw);
z = q0+fw;
z = (z*twoPowQ) % 8.0;
z -= ((int)z);
if (q > 0) {
// some parentheses for Eclipse formatter's weaknesses with bits shifts
iq4 &= (0xFFFFFF>>q);
ih = (iq4>>(23-q));
} else if (q == 0) {
ih = iq4>>23;
} else if (z >= 0.5) {
ih = 2;
} else {
ih = 0;
}
if (ih > 0) {
if (iq0 != 0) {
iq0 = 0x1000000 - iq0;
iq1 = 0x0FFFFFF - iq1;
iq2 = 0x0FFFFFF - iq2;
iq3 = 0x0FFFFFF - iq3;
iq4 = 0x0FFFFFF - iq4;
} else {
if (iq1 != 0) {
iq1 = 0x1000000 - iq1;
iq2 = 0x0FFFFFF - iq2;
iq3 = 0x0FFFFFF - iq3;
iq4 = 0x0FFFFFF - iq4;
} else {
if (iq2 != 0) {
iq2 = 0x1000000 - iq2;
iq3 = 0x0FFFFFF - iq3;
iq4 = 0x0FFFFFF - iq4;
} else {
if (iq3 != 0) {
iq3 = 0x1000000 - iq3;
iq4 = 0x0FFFFFF - iq4;
} else {
if (iq4 != 0) {
iq4 = 0x1000000 - iq4;
}
}
}
}
}
if (q > 0) {
switch (q) {
case 1:
iq4 &= 0x7FFFFF;
break;
case 2:
iq4 &= 0x3FFFFF;
break;
}
}
}
fw = twoPowQ * TWO_POW_N24; // q -= 24, so initializing fw with ((2^q)*(2^-24)=2^(q-24))
} else {
iq4 = (int)(z/twoPowQ);
fw = twoPowQ;
}
q4 = fw*iq4;
fw *= TWO_POW_N24;
q3 = fw*iq3;
fw *= TWO_POW_N24;
q2 = fw*iq2;
fw *= TWO_POW_N24;
q1 = fw*iq1;
fw *= TWO_POW_N24;
q0 = fw*iq0;
fw *= TWO_POW_N24;
fw = TWOPI_TAB0*q4;
fw += TWOPI_TAB0*q3 + TWOPI_TAB1*q4;
fw += TWOPI_TAB0*q2 + TWOPI_TAB1*q3 + TWOPI_TAB2*q4;
fw += TWOPI_TAB0*q1 + TWOPI_TAB1*q2 + TWOPI_TAB2*q3 + TWOPI_TAB3*q4;
fw += TWOPI_TAB0*q0 + TWOPI_TAB1*q1 + TWOPI_TAB2*q2 + TWOPI_TAB3*q3 + TWOPI_TAB4*q4;
return (ih == 0) ? fw : -fw;
}
//--------------------------------------------------------------------------
// STATIC INITIALIZATIONS
//--------------------------------------------------------------------------
static {
init();
}
/**
* Initializes look-up tables.
*
* Using redefined pure Java treatments in this method, instead of Math
* or StrictMath ones (even asin(double)), can make this class load much
* slower, because class loading is likely not to be optimized.
*
* Could use strictfp and StrictMath here, to always end up with the same tables,
* but this class doesn't guarantee identical results across various architectures,
* so to make it simple (less keywords/dependencies), we don't use strictfp,
* and use Math - which could also help if used StrictMath methods are slow.
*/
private static void init() {
/*
* sin and cos
*/
final int SIN_COS_PI_INDEX = (SIN_COS_TABS_SIZE-1)/2;
final int SIN_COS_PI_MUL_2_INDEX = 2*SIN_COS_PI_INDEX;
final int SIN_COS_PI_MUL_0_5_INDEX = SIN_COS_PI_INDEX/2;
final int SIN_COS_PI_MUL_1_5_INDEX = 3*SIN_COS_PI_INDEX/2;
for (int i=0;i