package jafama;
/*
* =============================================================================
* Copyright (C) 2009 oma
*
* This program is free software: you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public License
* as published by the Free Software Foundation, either version 3 of
* the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see .
* =============================================================================
* 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.
* =============================================================================
*/
/**
* Class providing math treatments that:
* - are meant to be faster than those of java.lang.Math class (depending on
* JVM or JVM options, they might be slower),
* - are still somehow accurate and robust (handling of NaN and such),
* - do not (or not directly) generate objects at run time (no "new").
*
* Other than optimized treatments, a valuable feature of this class is the
* presence of angles normalization methods, derived from those used in
* java.lang.Math (for which, sadly, no API is provided, letting everyone
* with the terrible responsibility to write their own ones).
*
* Non-redefined methods of java.lang.Math class are also available,
* for easy replacement.
*
* Use of look-up tables: around 1 Mo total, and initialized on class load.
*
* - Methods with same signature than Math ones, are meant to return
* "good" approximations on all range.
* - Methods terminating with "Fast" are meant to return "good" approximation
* on a reduced range only.
* - Methods terminating with "Quick" are meant to be quick, but do not
* return a good approximation, and might only work on a reduced 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.)."
*/
public strictfp final class FastMath {
/*
* You might want to modify this code some, for example to use redefined sqrt if you
* don't have a hardware one, or to use fast angle normalization treatments by default.
*
* 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 unaccurate 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.
*/
//--------------------------------------------------------------------------
// STATIC CONFIGURATION
//--------------------------------------------------------------------------
// Set it to true if Math.sqrt is slow (i.e. if it's not a hardware sqrt).
private static final boolean USE_REDEFINED_SQRT = false;
// Set it to true if Math.log is slow (i.e. if it's not a hardware log).
private static final boolean USE_REDEFINED_LOG = false;
// Set it to true if FastMath.sqrt(double) is slow (more tables, but less calls to FastMath.sqrt(double)).
private static final boolean USE_POWTABS_FOR_ASIN = false;
//--------------------------------------------------------------------------
// GENERAL CONSTANTS
//--------------------------------------------------------------------------
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 = Double.longBitsToDouble(0x4170000000000000L);
private static final double TWO_POW_N24 = Double.longBitsToDouble(0x3E70000000000000L);
private static final double TWO_POW_26 = Double.longBitsToDouble(0x4190000000000000L);
private static final double TWO_POW_N26 = Double.longBitsToDouble(0x3E50000000000000L);
// First double value (from zero) such as (value+-1/value == value).
private static final double TWO_POW_27 = Double.longBitsToDouble(0x41A0000000000000L);
private static final double TWO_POW_N27 = Double.longBitsToDouble(0x3E40000000000000L);
private static final double TWO_POW_N28 = Double.longBitsToDouble(0x3E30000000000000L);
private static final double TWO_POW_52 = Double.longBitsToDouble(0x4330000000000000L);
private static final double TWO_POW_N54 = Double.longBitsToDouble(0x3C90000000000000L);
private static final double TWO_POW_N55 = Double.longBitsToDouble(0x3C80000000000000L);
private static final double TWO_POW_66 = Double.longBitsToDouble(0x4410000000000000L);
private static final double TWO_POW_450 = Double.longBitsToDouble(0x5C10000000000000L);
private static final double TWO_POW_N450 = Double.longBitsToDouble(0x23D0000000000000L);
private static final double TWO_POW_750 = Double.longBitsToDouble(0x6ED0000000000000L);
private static final double TWO_POW_N750 = Double.longBitsToDouble(0x1110000000000000L);
// Smallest double normal value.
private static final double MIN_DOUBLE_NORMAL = Double.longBitsToDouble(0x0010000000000000L); // 2.2250738585072014E-308
private static final int MIN_DOUBLE_EXPONENT = -1074;
private static final int MAX_DOUBLE_EXPONENT = 1023;
private static final double LOG_2 = Math.log(2);
private static final double LOG_TWO_POW_27 = Math.log(TWO_POW_27);
private static final double LOG_DOUBLE_MAX_VALUE = Math.log(Double.MAX_VALUE);
private static final double DOUBLE_BEFORE_60 = Math.nextAfter(60.0, 0.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 = Math.pow(2,20)*(2*Math.PI);
/**
* 2*Math.PI, normalized into [-PI,PI].
*/
private static final double TWO_MATH_PI_IN_MINUS_PI_PI = -2.449293598153844E-16;
//--------------------------------------------------------------------------
// CONSTANTS AND TABLES FOR COS, SIN
//--------------------------------------------------------------------------
private static final int SIN_COS_TABS_SIZE = (1<<11) + 1;
private static final double SIN_COS_DELTA_HI = TWOPI_HI/(SIN_COS_TABS_SIZE-1);
private static final double SIN_COS_DELTA_LO = TWOPI_LO/(SIN_COS_TABS_SIZE-1);
private static final double SIN_COS_INDEXER = 1/(SIN_COS_DELTA_HI+SIN_COS_DELTA_LO);
private static final double[] sinTab = new double[SIN_COS_TABS_SIZE];
private static final double[] cosTab = new double[SIN_COS_TABS_SIZE];
// Max abs value for fast modulo, above which we use regular angle normalization.
// This value must be < (Integer.MAX_VALUE / SIN_COS_INDEXER), to stay in range of int type.
// The higher it is, the higher the error, but also the faster it is for lower values.
// If you set it to ((Integer.MAX_VALUE / SIN_COS_INDEXER) * 0.99), worse accuracy on double range is about 1e-10.
private static final double SIN_COS_MAX_VALUE_FOR_INT_MODULO = ((Integer.MAX_VALUE>>9) / SIN_COS_INDEXER) * 0.99;
//--------------------------------------------------------------------------
// CONSTANTS AND TABLES FOR TAN
//--------------------------------------------------------------------------
// We use the following formula:
// 1) tan(-x) = -tan(x)
// 2) tan(x) = 1/tan(PI/2-x)
// ---> we only have to compute tan(x) on [0,A] with PI/4<=A= 45deg, and supposed to be >= 51.4deg, as fdlibm code is not
// supposed to work with values inferior to that (51.4deg is about
// (PI/2-Double.longBitsToDouble(0x3FE5942800000000L))).
private static final double TAN_MAX_VALUE_FOR_TABS = Math.toRadians(77.0);
private static final int TAN_TABS_SIZE = (int)((TAN_MAX_VALUE_FOR_TABS/(Math.PI/2)) * (TAN_VIRTUAL_TABS_SIZE-1)) + 1;
private static final double TAN_DELTA_HI = PIO2_HI/(TAN_VIRTUAL_TABS_SIZE-1);
private static final double TAN_DELTA_LO = PIO2_LO/(TAN_VIRTUAL_TABS_SIZE-1);
private static final double TAN_INDEXER = 1/(TAN_DELTA_HI+TAN_DELTA_LO);
private static final double[] tanTab = new double[TAN_TABS_SIZE];
private static final double[] tanDer1DivF1Tab = new double[TAN_TABS_SIZE];
private static final double[] tanDer2DivF2Tab = new double[TAN_TABS_SIZE];
private static final double[] tanDer3DivF3Tab = new double[TAN_TABS_SIZE];
private static final double[] tanDer4DivF4Tab = new double[TAN_TABS_SIZE];
// Max abs value for fast modulo, above which we use regular angle normalization.
// This value must be < (Integer.MAX_VALUE / SIN_COS_INDEXER), to stay in range of int type.
// The higher it is, the higher the error, but also the faster it is for lower values.
private static final double TAN_MAX_VALUE_FOR_INT_MODULO = (((Integer.MAX_VALUE>>9) / TAN_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 is >= sin(25deg).
private static final double ASIN_MAX_VALUE_FOR_TABS = Math.sin(Math.toRadians(73.0));
private static final int ASIN_TABS_SIZE = (1<<13) + 1;
private static final double ASIN_DELTA = ASIN_MAX_VALUE_FOR_TABS/(ASIN_TABS_SIZE - 1);
private static final double ASIN_INDEXER = 1/ASIN_DELTA;
private static final double[] asinTab = new double[ASIN_TABS_SIZE];
private static final double[] asinDer1DivF1Tab = new double[ASIN_TABS_SIZE];
private static final double[] asinDer2DivF2Tab = new double[ASIN_TABS_SIZE];
private static final double[] asinDer3DivF3Tab = new double[ASIN_TABS_SIZE];
private static final double[] asinDer4DivF4Tab = new double[ASIN_TABS_SIZE];
private static final double ASIN_MAX_VALUE_FOR_POWTABS = Math.sin(Math.toRadians(88.6));
private static final int ASIN_POWTABS_POWER = 84;
private static final double ASIN_POWTABS_ONE_DIV_MAX_VALUE = 1/ASIN_MAX_VALUE_FOR_POWTABS;
private static final int ASIN_POWTABS_SIZE = USE_POWTABS_FOR_ASIN ? (1<<12) + 1 : 0;
private static final int ASIN_POWTABS_SIZE_MINUS_ONE = ASIN_POWTABS_SIZE - 1;
private static final double[] asinParamPowTab = new double[ASIN_POWTABS_SIZE];
private static final double[] asinPowTab = new double[ASIN_POWTABS_SIZE];
private static final double[] asinDer1DivF1PowTab = new double[ASIN_POWTABS_SIZE];
private static final double[] asinDer2DivF2PowTab = new double[ASIN_POWTABS_SIZE];
private static final double[] asinDer3DivF3PowTab = new double[ASIN_POWTABS_SIZE];
private static final double[] asinDer4DivF4PowTab = new double[ASIN_POWTABS_SIZE];
private static final double ASIN_PIO2_HI = Double.longBitsToDouble(0x3FF921FB54442D18L); // 1.57079632679489655800e+00
private static final double ASIN_PIO2_LO = Double.longBitsToDouble(0x3C91A62633145C07L); // 6.12323399573676603587e-17
private static final double ASIN_PS0 = Double.longBitsToDouble(0x3fc5555555555555L); // 1.66666666666666657415e-01
private static final double ASIN_PS1 = Double.longBitsToDouble(0xbfd4d61203eb6f7dL); // -3.25565818622400915405e-01
private static final double ASIN_PS2 = Double.longBitsToDouble(0x3fc9c1550e884455L); // 2.01212532134862925881e-01
private static final double ASIN_PS3 = Double.longBitsToDouble(0xbfa48228b5688f3bL); // -4.00555345006794114027e-02
private static final double ASIN_PS4 = Double.longBitsToDouble(0x3f49efe07501b288L); // 7.91534994289814532176e-04
private static final double ASIN_PS5 = Double.longBitsToDouble(0x3f023de10dfdf709L); // 3.47933107596021167570e-05
private static final double ASIN_QS1 = Double.longBitsToDouble(0xc0033a271c8a2d4bL); // -2.40339491173441421878e+00
private static final double ASIN_QS2 = Double.longBitsToDouble(0x40002ae59c598ac8L); // 2.02094576023350569471e+00
private static final double ASIN_QS3 = Double.longBitsToDouble(0xbfe6066c1b8d0159L); // -6.88283971605453293030e-01
private static final double ASIN_QS4 = Double.longBitsToDouble(0x3fb3b8c5b12e9282L); // 7.70381505559019352791e-02
//--------------------------------------------------------------------------
// CONSTANTS AND TABLES FOR ATAN
//--------------------------------------------------------------------------
// We use the formula atan(-x) = -atan(x)
// ---> 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 = Math.tan(Math.toRadians(74.0));
private static final int ATAN_TABS_SIZE = (1<<12) + 1;
private static final double ATAN_DELTA = ATAN_MAX_VALUE_FOR_TABS/(ATAN_TABS_SIZE - 1);
private static final double ATAN_INDEXER = 1/ATAN_DELTA;
private static final double[] atanTab = new double[ATAN_TABS_SIZE];
private static final double[] atanDer1DivF1Tab = new double[ATAN_TABS_SIZE];
private static final double[] atanDer2DivF2Tab = new double[ATAN_TABS_SIZE];
private static final double[] atanDer3DivF3Tab = new double[ATAN_TABS_SIZE];
private static final double[] atanDer4DivF4Tab = new double[ATAN_TABS_SIZE];
private static final double ATAN_HI3 = Double.longBitsToDouble(0x3ff921fb54442d18L); // 1.57079632679489655800e+00 atan(inf)hi
private static final double ATAN_LO3 = Double.longBitsToDouble(0x3c91a62633145c07L); // 6.12323399573676603587e-17 atan(inf)lo
private static final double ATAN_AT0 = Double.longBitsToDouble(0x3fd555555555550dL); // 3.33333333333329318027e-01
private static final double ATAN_AT1 = Double.longBitsToDouble(0xbfc999999998ebc4L); // -1.99999999998764832476e-01
private static final double ATAN_AT2 = Double.longBitsToDouble(0x3fc24924920083ffL); // 1.42857142725034663711e-01
private static final double ATAN_AT3 = Double.longBitsToDouble(0xbfbc71c6fe231671L); // -1.11111104054623557880e-01
private static final double ATAN_AT4 = Double.longBitsToDouble(0x3fb745cdc54c206eL); // 9.09088713343650656196e-02
private static final double ATAN_AT5 = Double.longBitsToDouble(0xbfb3b0f2af749a6dL); // -7.69187620504482999495e-02
private static final double ATAN_AT6 = Double.longBitsToDouble(0x3fb10d66a0d03d51L); // 6.66107313738753120669e-02
private static final double ATAN_AT7 = Double.longBitsToDouble(0xbfadde2d52defd9aL); // -5.83357013379057348645e-02
private static final double ATAN_AT8 = Double.longBitsToDouble(0x3fa97b4b24760debL); // 4.97687799461593236017e-02
private static final double ATAN_AT9 = Double.longBitsToDouble(0xbfa2b4442c6a6c2fL); // -3.65315727442169155270e-02
private static final double ATAN_AT10 = Double.longBitsToDouble(0x3f90ad3ae322da11L); // 1.62858201153657823623e-02
//--------------------------------------------------------------------------
// TABLE FOR POWERS OF TWO
//--------------------------------------------------------------------------
private static final double[] twoPowTab = new double[(MAX_DOUBLE_EXPONENT-MIN_DOUBLE_EXPONENT)+1];
//--------------------------------------------------------------------------
// CONSTANTS AND TABLES FOR EXP AND EXPM1
//--------------------------------------------------------------------------
private static final double EXP_OVERFLOW_LIMIT = Double.longBitsToDouble(0x40862E42FEFA39EFL); // 7.09782712893383973096e+02
private static final double EXP_UNDERFLOW_LIMIT = Double.longBitsToDouble(0xC0874910D52D3051L); // -7.45133219101941108420e+02
private static final double EXP_MIN_INT_LIMIT = -705;
private static final int EXP_LO_DISTANCE_TO_ZERO_POT = 0;
private static final int EXP_LO_DISTANCE_TO_ZERO = (1< 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)));
}
/**
* Quick cosine, with accuracy of about 1.5e-3 (PI/)
* for |angle| < 6588397.0 (Integer.MAX_VALUE * (2*PI/)),
* and no accuracy at all for larger values.
*
* @param angle Angle in radians.
* @return Angle cosine.
*/
public static double cosQuick(double angle) {
return cosTab[((int)(FastMath.abs(angle) * SIN_COS_INDEXER + 0.5)) & (SIN_COS_TABS_SIZE-2)];
}
/**
* @param angle Angle in radians.
* @return Angle sine.
*/
public static double sin(double angle) {
boolean negateResult;
if (angle < 0.0) {
angle = -angle;
negateResult = true;
} else {
negateResult = false;
}
if (angle > 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;
}
/**
* Quick sine, with accuracy of about 1.5e-3 (PI/)
* for |angle| < 6588397.0 (Integer.MAX_VALUE * (2*PI/)),
* and no accuracy at all for larger values.
*
* @param angle Angle in radians.
* @return Angle sine.
*/
public static double sinQuick(double angle) {
return cosTab[((int)(FastMath.abs(angle-Math.PI/2) * SIN_COS_INDEXER + 0.5)) & (SIN_COS_TABS_SIZE-2)];
}
/**
* Computes sine and cosine together, at the cost of... a dependency of this class with DoubleWrapper.
*
* @param angle Angle in radians.
* @param sine Angle sine.
* @param cosine Angle cosine.
*/
public static void sinAndCos(double angle, DoubleWrapper sine, DoubleWrapper cosine) {
// Using the same algorithm than sin(double) method, and computing also cosine at the end.
boolean negateResult;
if (angle < 0.0) {
angle = -angle;
negateResult = true;
} else {
negateResult = false;
}
if (angle > 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)));
sine.value = (negateResult) ? -result : result;
cosine.value = indexCos + delta * (-indexSin + delta * (-indexCos * ONE_DIV_F2 + delta * (indexSin * ONE_DIV_F3 + delta * indexCos * ONE_DIV_F4)));
}
/**
* @param angle Angle in radians.
* @return Angle tangent.
*/
public static double tan(double angle) {
if (FastMath.abs(angle) > TAN_MAX_VALUE_FOR_INT_MODULO) {
// Faster than using normalizeMinusHalfPiHalfPi.
angle = remainderTwoPi(angle);
if (angle < -Math.PI/2) {
angle += Math.PI;
} else if (angle > Math.PI/2) {
angle -= Math.PI;
}
}
boolean negateResult;
if (angle < 0.0) {
angle = -angle;
negateResult = true;
} else {
negateResult = false;
}
int index = (int)(angle * TAN_INDEXER + 0.5);
double delta = (angle - index * TAN_DELTA_HI) - index * TAN_DELTA_LO;
// index modulo PI, i.e. 2*(virtual tab size minus one).
index &= (2*(TAN_VIRTUAL_TABS_SIZE-1)-1); // index % (2*(TAN_VIRTUAL_TABS_SIZE-1))
// Here, index is in [0,2*(TAN_VIRTUAL_TABS_SIZE-1)-1], i.e. indicates an angle in [0,PI[.
if (index > (TAN_VIRTUAL_TABS_SIZE-1)) {
index = (2*(TAN_VIRTUAL_TABS_SIZE-1)) - index;
delta = -delta;
negateResult = !negateResult;
}
double result;
if (index < TAN_TABS_SIZE) {
result = tanTab[index] + delta * (tanDer1DivF1Tab[index] + delta * (tanDer2DivF2Tab[index] + delta * (tanDer3DivF3Tab[index] + delta * tanDer4DivF4Tab[index])));
} else { // angle in ]TAN_MAX_VALUE_FOR_TABS,TAN_MAX_VALUE_FOR_INT_MODULO], or angle is NaN
// Using tan(angle) == 1/tan(PI/2-angle) formula: changing angle (index and delta), and inverting.
index = (TAN_VIRTUAL_TABS_SIZE-1) - index;
result = 1/(tanTab[index] - delta * (tanDer1DivF1Tab[index] - delta * (tanDer2DivF2Tab[index] - delta * (tanDer3DivF3Tab[index] - delta * tanDer4DivF4Tab[index]))));
}
return (negateResult) ? -result : result;
}
/**
* @param value Value in [-1,1].
* @return Value arccosine, in radians, in [0,PI].
*/
public static double acos(double value) {
return Math.PI/2 - FastMath.asin(value);
}
/**
* If value is not NaN and is outside [-1,1] range, closest value in this range is used.
*
* @param value Value in [-1,1].
* @return Value arccosine, in radians, in [0,PI].
*/
public static double acosInRange(double value) {
if (value < -1) {
return Math.PI;
} else if (value > 1) {
return 0.0;
} else {
return FastMath.acos(value);
}
}
/**
* @param value Value in [-1,1].
* @return Value arcsine, in radians, in [-PI/2,PI/2].
*/
public 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 if (USE_POWTABS_FOR_ASIN && (value <= ASIN_MAX_VALUE_FOR_POWTABS)) {
int index = (int)(FastMath.pow(value * ASIN_POWTABS_ONE_DIV_MAX_VALUE, ASIN_POWTABS_POWER) * ASIN_POWTABS_SIZE_MINUS_ONE + 0.5);
double delta = value - asinParamPowTab[index];
double result = asinPowTab[index] + delta * (asinDer1DivF1PowTab[index] + delta * (asinDer2DivF2PowTab[index] + delta * (asinDer3DivF3PowTab[index] + delta * asinDer4DivF4PowTab[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 = FastMath.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;
}
}
}
}
/**
* If value is not NaN and is outside [-1,1] range, closest value in this range is used.
*
* @param value Value in [-1,1].
* @return Value arcsine, in radians, in [-PI/2,PI/2].
*/
public static double asinInRange(double value) {
if (value < -1) {
return -Math.PI/2;
} else if (value > 1) {
return Math.PI/2;
} else {
return FastMath.asin(value);
}
}
/**
* @param value A double value.
* @return Value arctangent, in radians, in [-PI/2,PI/2].
*/
public 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 (Double.isNaN(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).
*/
public static double atan2(double y, double x) {
if (x > 0.0) {
if (y == 0.0) {
return (1/y == Double.NEGATIVE_INFINITY) ? -0.0 : 0.0;
}
if (x == Double.POSITIVE_INFINITY) {
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;
}
} else {
return FastMath.atan(y/x);
}
} else if (x < 0.0) {
if (y == 0.0) {
return (1/y == Double.NEGATIVE_INFINITY) ? -Math.PI : Math.PI;
}
if (x == Double.NEGATIVE_INFINITY) {
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;
}
} else if (y > 0.0) {
return Math.PI/2 + FastMath.atan(-x/y);
} else if (y < 0.0) {
return -Math.PI/2 - FastMath.atan(x/y);
} else {
return Double.NaN;
}
} else if (x == 0.0) {
if (y == 0.0) {
if (1/x == Double.NEGATIVE_INFINITY) {
return (1/y == Double.NEGATIVE_INFINITY) ? -Math.PI : Math.PI;
} else {
return (1/y == Double.NEGATIVE_INFINITY) ? -0.0 : 0.0;
}
}
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;
}
}
/**
* @param value A double value.
* @return Value hyperbolic cosine.
*/
public static double cosh(double value) {
// cosh(x) = (exp(x)+exp(-x))/2
if (value < 0.0) {
value = -value;
}
if (value < LOG_TWO_POW_27) {
if (value < TWO_POW_N27) {
// cosh(x)
// = (exp(x)+exp(-x))/2
// = ((1+x+x^2/2!+...) + (1-x+x^2/2!-...))/2
// = 1+x^2/2!+x^4/4!+...
// For value of x small in magnitude, the sum of the terms does not add to 1.
return 1;
} else {
double t = FastMath.exp(value);
return 0.5 * (t+1/t);
}
} else if (value < LOG_DOUBLE_MAX_VALUE) {
return 0.5 * FastMath.exp(value);
} else {
double t = FastMath.exp(value*0.5);
return (0.5*t)*t;
}
}
/**
* @param value A double value.
* @return Value hyperbolic sine.
*/
public static double sinh(double value) {
// sinh(x) = (exp(x)-exp(-x))/2
double h;
if (value < 0.0) {
value = -value;
h = -0.5;
} else {
h = 0.5;
}
if (value < 22.0) {
if (value < TWO_POW_N28) {
return (h < 0.0) ? -value : value;
} else {
double t = FastMath.expm1(value);
// Might be more accurate, if value < 1: return h*((t+t)-t*t/(t+1.0)).
return h * (t + t/(t+1.0));
}
} else if (value < LOG_DOUBLE_MAX_VALUE) {
return h * FastMath.exp(value);
} else {
double t = FastMath.exp(value*0.5);
return (h*t)*t;
}
}
/**
* Computes hyperbolic sine and hyperbolic cosine together, at the cost of... a dependency of this class with DoubleWrapper.
*
* @param value A double value.
* @param hsine Value hyperbolic sine.
* @param hcosine Value hyperbolic cosine.
*/
public static void sinhAndCosh(double value, DoubleWrapper hsine, DoubleWrapper hcosine) {
// Mixup of sinh and cosh treatments: if you modify them,
// you might want to also modify this.
double h;
if (value < 0.0) {
value = -value;
h = -0.5;
} else {
h = 0.5;
}
// LOG_TWO_POW_27 = 18.714973875118524
if (value < LOG_TWO_POW_27) { // test from cosh
// sinh
if (value < TWO_POW_N28) {
hsine.value = (h < 0.0) ? -value : value;
} else {
double t = FastMath.expm1(value);
hsine.value = h * (t + t/(t+1.0));
}
// cosh
if (value < TWO_POW_N27) {
hcosine.value = 1;
} else {
double t = FastMath.exp(value);
hcosine.value = 0.5 * (t+1/t);
}
} else if (value < 22.0) { // test from sinh
// Here, value is in [18.714973875118524,22.0[.
double t = FastMath.expm1(value);
hsine.value = h * (t + t/(t+1.0));
hcosine.value = 0.5 * (t+1.0);
} else {
if (value < LOG_DOUBLE_MAX_VALUE) {
hsine.value = h * FastMath.exp(value);
} else {
double t = FastMath.exp(value*0.5);
hsine.value = (h*t)*t;
}
hcosine.value = Math.abs(hsine.value);
}
}
/**
* @param value A double value.
* @return Value hyperbolic tangent.
*/
public static double tanh(double value) {
// tanh(x) = sinh(x)/cosh(x)
// = (exp(x)-exp(-x))/(exp(x)+exp(-x))
// = (exp(2*x)-1)/(exp(2*x)+1)
boolean negateResult;
if (value < 0.0) {
value = -value;
negateResult = true;
} else {
negateResult = false;
}
double z;
if (value < 22.0) {
if (value < TWO_POW_N55) {
return (negateResult) ? -value*(1.0-value) : value*(1.0+value);
} else if (value >= 1) {
z = 1.0-2.0/(FastMath.expm1(value+value)+2.0);
} else {
double t = FastMath.expm1(-(value+value));
z = -t/(t+2.0);
}
} else {
z = (value != value) ? Double.NaN : 1.0;
}
return (negateResult) ? -z : z;
}
/**
* 1e-13ish accuracy (or better) on whole double range.
*
* @param value A double value.
* @param power A power.
* @return value^power.
*/
public static double pow(double value, double power) {
if (power == 0.0) {
return 1.0;
} else if (power == 1.0) {
return value;
}
if (value <= 0.0) {
// powerInfo: 0 if not integer, 1 if even integer, -1 if odd integer
int powerInfo;
if (Math.abs(power) >= (TWO_POW_52*2)) {
// The binary digit just before coma is outside mantissa,
// thus it is always 0: power is an even integer.
powerInfo = 1;
} else {
// If power's magnitude permits, we cast into int instead of into long,
// as it is faster (TODO check it on 64 bits machines...).
if (Math.abs(power) <= (double)Integer.MAX_VALUE) {
int powerAsInt = (int)power;
if (power == (double)powerAsInt) {
powerInfo = ((powerAsInt & 1) == 0) ? 1 : -1;
} else { // power is not an integer (and not NaN, due to test against Integer.MAX_VALUE)
powerInfo = 0;
}
} else {
long powerAsLong = (long)power;
if (power == (double)powerAsLong) {
powerInfo = ((powerAsLong & 1) == 0) ? 1 : -1;
} else { // power is not an integer, or is NaN
if (power != power) {
return Double.NaN;
}
powerInfo = 0;
}
}
}
if (value == 0.0) {
if (power < 0.0) {
return (powerInfo < 0) ? 1/value : Double.POSITIVE_INFINITY;
} else { // power > 0.0 (0 and NaN cases already treated)
return (powerInfo < 0) ? value : 0.0;
}
} else { // value < 0.0
if (value == Double.NEGATIVE_INFINITY) {
if (powerInfo < 0) { // power odd integer
return (power < 0.0) ? -0.0 : Double.NEGATIVE_INFINITY;
} else { // power even integer, or not an integer
return (power < 0.0) ? 0.0 : Double.POSITIVE_INFINITY;
}
} else {
return (powerInfo != 0) ? powerInfo * FastMath.exp(power*FastMath.log(-value)) : Double.NaN;
}
}
} else { // a > 0.0, or a is NaN
return FastMath.exp(power*FastMath.log(value));
}
}
/**
* This treatment is somehow accurate for low values of |power|.
*
* @param value A double value.
* @param power A power.
* @return value^power.
*/
public static double powFast(double value, int power) {
if (power > 5) { // Most common case first.
double oddRemains = 1.0;
do {
// Test if power is odd.
if ((power & 1) != 0) {
oddRemains *= value;
}
value *= value;
power >>= 1; // power = power / 2
} while (power > 5);
// Here, power is in [3,5]: faster to finish outside the loop.
if (power == 3) {
return oddRemains * value * value * value;
} else {
double v2 = value * value;
if (power == 4) {
return oddRemains * v2 * v2;
} else { // power == 5
return oddRemains * v2 * v2 * value;
}
}
} else if (power >= 0) { // power in [0,5]
if (power < 3) { // power in [0,2]
if (power == 2) { // Most common case first.
return value * value;
} else if (power != 0) { // faster than == 1
return value;
} else { // power == 0
return 1.0;
}
} else { // power in [3,5]
if (power == 3) {
return value * value * value;
} else { // power in [4,5]
double v2 = value * value;
if (power == 4) {
return v2 * v2;
} else { // power == 5
return v2 * v2 * value;
}
}
}
} else { // power < 0
// Opposite of Integer.MIN_VALUE does not exist as int.
if (power == Integer.MIN_VALUE) {
return 1.0/(FastMath.powFast(value,-(power+1)) * value);
} else {
return 1.0/FastMath.powFast(value,-power);
}
}
}
/**
* Returns the exact result, provided it's in double range.
*
* @param power A power.
* @return 2^power.
*/
public static double twoPow(int power) {
if (power >= 0) {
if (power <= MAX_DOUBLE_EXPONENT) {
return twoPowTab[power-MIN_DOUBLE_EXPONENT];
} else {
// Overflow.
return Double.POSITIVE_INFINITY;
}
} else {
if (power >= MIN_DOUBLE_EXPONENT) {
return twoPowTab[power-MIN_DOUBLE_EXPONENT];
} else {
// Underflow.
return 0.0;
}
}
/* Version without table.
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 value An int value.
* @return value*value.
*/
public static int pow2(int value) {
return value*value;
}
/**
* @param value A long value.
* @return value*value.
*/
public static long pow2(long value) {
return value*value;
}
/**
* @param value A float value.
* @return value*value.
*/
public static float pow2(float value) {
return value*value;
}
/**
* @param value A double value.
* @return value*value.
*/
public static double pow2(double value) {
return value*value;
}
/**
* @param value An int value.
* @return value*value*value.
*/
public static int pow3(int value) {
return value*value*value;
}
/**
* @param value A long value.
* @return value*value*value.
*/
public static long pow3(long value) {
return value*value*value;
}
/**
* @param value A float value.
* @return value*value*value.
*/
public static float pow3(float value) {
return value*value*value;
}
/**
* @param value A double value.
* @return value*value*value.
*/
public static double pow3(double value) {
return value*value*value;
}
/**
* @param value A double value.
* @return e^value.
*/
public static double exp(double value) {
// exp(x) = exp([x])*exp(y)
// with [x] the integer part of x, and y = x-[x]
// ===>
// We find an approximation of y, called z.
// ===>
// exp(x) = exp([x])*(exp(z)*exp(epsilon))
// ===>
// We have exp([x]) and exp(z) pre-computed in tables, we "just" have to compute exp(epsilon).
//
// We use the same indexing (cast to int) to compute x integer part and the
// table index corresponding to z, to avoit two int casts.
// Also, to optimize index multiplication and division, we use powers of two,
// so that we can do it with bits shifts.
if (value >= 0.0) {
if (value > EXP_OVERFLOW_LIMIT) {
return Double.POSITIVE_INFINITY;
}
int i = (int)(value*EXP_LO_INDEXING);
int valueInt = (i>>EXP_LO_INDEXING_DIV_SHIFT);
i -= (valueInt<= EXP_UNDERFLOW_LIMIT)) { // value < EXP_UNDERFLOW_LIMIT, or value is NaN
return (value < EXP_UNDERFLOW_LIMIT) ? 0.0 : Double.NaN;
}
// TODO JVM bug with -server option: test with values of all magnitudes
// is very slow, if using (int)x instead of -(int)-x or (int)(long)x (which give the same result).
// The guessed cause is that when the same expression is used to define "i" in
// both sides of the above "else", some (desastrous) optimization is done which factorizes
// it above the first "if" statement, making it computed all the time, without the protecting "sub-ifs".
// Since cast from double to int with huge values is extremely slow (another buggy feature...),
// this makes this whole treatment extremely slow for huge values.
// The solution is therefore to modify a bit the expression for the "optimization" not to occur.
int i = -(int)-(value*EXP_LO_INDEXING);
int valueInt = -((-i)>>EXP_LO_INDEXING_DIV_SHIFT);
i -= ((valueInt)<= EXP_MIN_INT_LIMIT) ? tmp : tmp * TWO_POW_N54;
}
}
/**
* Much more accurate than exp(value)-1, for values close to zero.
*
* @param value A double value.
* @return e^value-1.
*/
public static double expm1(double value) {
// If value is far from zero, we use exp(value)-1.
//
// If value is close to zero, we use the following formula:
// exp(value)-1
// = exp(valueApprox)*exp(epsilon)-1
// = exp(valueApprox)*(exp(epsilon)-exp(-valueApprox))
// = exp(valueApprox)*(1+epsilon+epsilon^2/2!+...-exp(-valueApprox))
// = exp(valueApprox)*((1-exp(-valueApprox))+epsilon+epsilon^2/2!+...)
// exp(valueApprox) and exp(-valueApprox) being stored in tables.
if (Math.abs(value) < EXP_LO_DISTANCE_TO_ZERO) {
// Taking int part instead of rounding, which takes too long.
int i = (int)(value*EXP_LO_INDEXING);
double delta = value-i*(1.0/EXP_LO_INDEXING);
return expLoPosTab[i+EXP_LO_TAB_MID_INDEX]*(expLoNegTab[i+EXP_LO_TAB_MID_INDEX]+delta*(1+delta*(1.0/2+delta*(1.0/6+delta*(1.0/24+delta*(1.0/120))))));
} else {
return FastMath.exp(value)-1;
}
}
/**
* @param value A double value.
* @return Value logarithm (base e).
*/
public static double log(double value) {
if (USE_REDEFINED_LOG) {
if (value > 0.0) {
if (value == Double.POSITIVE_INFINITY) {
return Double.POSITIVE_INFINITY;
}
// For normal values not close to 1.0, we use the following formula:
// log(value)
// = log(2^exponent*1.mantissa)
// = log(2^exponent) + log(1.mantissa)
// = exponent * log(2) + log(1.mantissa)
// = exponent * log(2) + log(1.mantissaApprox) + log(1.mantissa/1.mantissaApprox)
// = exponent * log(2) + log(1.mantissaApprox) + log(1+epsilon)
// = exponent * log(2) + log(1.mantissaApprox) + epsilon-epsilon^2/2+epsilon^3/3-epsilon^4/4+...
// with:
// 1.mantissaApprox <= 1.mantissa,
// log(1.mantissaApprox) in table,
// epsilon = (1.mantissa/1.mantissaApprox)-1
//
// To avoid bad relative error for small results,
// values close to 1.0 are treated apart, with the formula:
// log(x) = z*(2+z^2*((2.0/3)+z^2*((2.0/5))+z^2*((2.0/7))+...)))
// with z=(x-1)/(x+1)
double h;
if (value > 0.95) {
if (value < 1.14) {
double z = (value-1.0)/(value+1.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)))))));
}
h = 0.0;
} else if (value < MIN_DOUBLE_NORMAL) {
// Ensuring value is normal.
value *= TWO_POW_52;
// log(x*2^52)
// = log(x)-ln(2^52)
// = log(x)-52*ln(2)
h = -52*LOG_2;
} else {
h = 0.0;
}
int valueBitsHi = (int)(Double.doubleToRawLongBits(value)>>32);
int valueExp = (valueBitsHi>>20)-MAX_DOUBLE_EXPONENT;
// Getting the first LOG_BITS bits of the mantissa.
int xIndex = ((valueBitsHi<<12)>>>(32-LOG_BITS));
// 1.mantissa/1.mantissaApprox - 1
double z = (value * twoPowTab[-valueExp-MIN_DOUBLE_EXPONENT]) * logXInvTab[xIndex] - 1;
z *= (1-z*((1.0/2)-z*((1.0/3))));
return h + valueExp * LOG_2 + (logXLogTab[xIndex] + z);
} else if (value == 0.0) {
return Double.NEGATIVE_INFINITY;
} else { // value < 0.0, or value is NaN
return Double.NaN;
}
} else {
return Math.log(value);
}
}
/**
* Much more accurate than log(1+value), for values close to zero.
*
* @param value A double value.
* @return Logarithm (base e) of (1+value).
*/
public static double log1p(double value) {
if (false) {
// This also works. Simpler but a bit slower.
if (value == Double.POSITIVE_INFINITY) {
return Double.POSITIVE_INFINITY;
}
double valuePlusOne = 1+value;
if (valuePlusOne == 1.0) {
return value;
} else {
return Math.log(valuePlusOne)*(value/(valuePlusOne-1.0));
}
}
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 * twoPowTab[-valuePlusOneExp-MIN_DOUBLE_EXPONENT]) * 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;
}
}
/**
* @param value A double value.
* @return Value square root.
*/
public static double sqrt(double value) {
if (USE_REDEFINED_SQRT) {
// See cbrt for comments, sqrt uses the same ideas.
if (!(value > 0.0)) { // value <= 0.0, or value is NaN
return (value == 0.0) ? value : Double.NaN;
} else if (value == Double.POSITIVE_INFINITY) {
return Double.POSITIVE_INFINITY;
}
double h;
if (value < MIN_DOUBLE_NORMAL) {
value *= TWO_POW_52;
h = 2.0*TWO_POW_N26;
} else {
h = 2.0;
}
int valueBitsHi = (int)(Double.doubleToRawLongBits(value)>>32);
int valueExponentIndex = (valueBitsHi>>20)+(-MAX_DOUBLE_EXPONENT-MIN_DOUBLE_EXPONENT);
int xIndex = ((valueBitsHi<<12)>>>(32-SQRT_LO_BITS));
double result = sqrtXSqrtHiTab[valueExponentIndex] * sqrtXSqrtLoTab[xIndex];
double slope = sqrtSlopeHiTab[valueExponentIndex] * sqrtSlopeLoTab[xIndex];
value *= 0.25;
result += (value - result * result) * slope;
result += (value - result * result) * slope;
return h*(result + (value - result * result) * slope);
} else {
return Math.sqrt(value);
}
}
/**
* @param value A double value.
* @return Value cubic root.
*/
public static double cbrt(double value) {
double h;
if (value < 0.0) {
if (value == Double.NEGATIVE_INFINITY) {
return Double.NEGATIVE_INFINITY;
}
value = -value;
// Making sure value is normal.
if (value < MIN_DOUBLE_NORMAL) {
value *= (TWO_POW_52*TWO_POW_26);
// h = * /
h = -2*TWO_POW_N26;
} else {
h = -2.0;
}
} else {
if (!(value < Double.POSITIVE_INFINITY)) { // value is +infinity, or value is NaN
return value;
}
// Making sure value is normal.
if (value < MIN_DOUBLE_NORMAL) {
if (value == 0.0) {
// cbrt(0.0) = 0.0, cbrt(-0.0) = -0.0
return value;
}
value *= (TWO_POW_52*TWO_POW_26);
h = 2*TWO_POW_N26;
} else {
h = 2.0;
}
}
// Normal value is (2^ * ).
// First member cubic root is computed, and multiplied with an approximation
// of the cubic root of the second member, to end up with a good guess of
// the result before using Newton's method.
// To compute the cubic root approximation, we use the formula "cbrt(value) = cbrt(x) * cbrt(value/x)",
// choosing x as close to value as possible but inferior to it, so that cbrt(value/x) is close to 1
// (we could iterate on this method, using value/x as new value for each iteration,
// but finishing with Newton's method is faster).
// TODO for everywhere: check if slower to work on int like this with 64 bits machines
// (adding a shift and a cast here for that)...
int valueBitsHi = (int)(Double.doubleToRawLongBits(value)>>32);
int valueExponentIndex = (valueBitsHi>>20)+(-MAX_DOUBLE_EXPONENT-MIN_DOUBLE_EXPONENT);
// Getting the first CBRT_LO_BITS bits of the mantissa.
int xIndex = ((valueBitsHi<<12)>>>(32-CBRT_LO_BITS));
double result = cbrtXCbrtHiTab[valueExponentIndex] * cbrtXCbrtLoTab[xIndex];
double slope = cbrtSlopeHiTab[valueExponentIndex] * cbrtSlopeLoTab[xIndex];
// Lowering values to avoid overflows when using Newton's method
// (we will then just have to return twice the result).
// result^3 = value
// (result/2)^3 = value/8
value *= 0.125;
// No need to divide result here, as division is factorized in result computation tables.
// result *= 0.5;
// Newton's method, looking for y = x^(1/p):
// y(n) = y(n-1) + (x-y(n-1)^p) * slope(y(n-1))
// y(n) = y(n-1) + (x-y(n-1)^p) * (1/p)*(x(n-1)^(1/p-1))
// y(n) = y(n-1) + (x-y(n-1)^p) * (1/p)*(x(n-1)^((1-p)/p))
// with x(n-1)=y(n-1)^p, i.e.:
// y(n) = y(n-1) + (x-y(n-1)^p) * (1/p)*(y(n-1)^(1-p))
//
// For p=3:
// y(n) = y(n-1) + (x-y(n-1)^3) * (1/(3*y(n-1)^2))
// To save time, we don't recompute the slope between Newton's method steps,
// as initial slope is good enough for a few iterations.
//
// NB: slope = 1/(3*trueResult*trueResult)
// As we have result = trueResult/2 (to avoid overflows), we have:
// slope = 4/(3*result*result)
// = (4/3)*resultInv*resultInv
// with newResultInv = 1/newResult
// = 1/(oldResult+resultDelta)
// = (oldResultInv)*1/(1+resultDelta/oldResult)
// = (oldResultInv)*1/(1+resultDelta*oldResultInv)
// ~= (oldResultInv)*(1-resultDelta*oldResultInv)
// ===> Successive slopes could be computed without division, if needed,
// by computing resultInv (instead of slope right away) and retrieving
// slopes from it.
result += (value - result * result * result) * slope;
result += (value - result * result * result) * slope;
return h*(result + (value - result * result * result) * slope);
}
/**
* Behaves about (exactly?) like Math.IEEEremainder(double,double).
*
* @param dividend Dividend.
* @param divisor Divisor.
* @return Remainder of dividend/divisor.
*/
public static double remainder(double dividend, double divisor) {
if (Double.isInfinite(divisor)) {
if (Double.isInfinite(dividend)) {
return Double.NaN;
} else {
return dividend;
}
}
double value = dividend % divisor;
if (FastMath.abs(value+value) > FastMath.abs(divisor)) {
return value + ((value > 0.0) ? -FastMath.abs(divisor) : FastMath.abs(divisor));
} else {
return value;
}
}
/**
* @param angle Angle in radians.
* @return The same angle, in radians, but in [-Math.PI,Math.PI].
* @published
*/
public static double normalizeMinusPiPi(double angle) {
// Not modifying values in output range.
if ((angle >= -Math.PI) && (angle <= Math.PI)) {
return angle;
}
double angleMinusPiPiOrSo = remainderTwoPi(angle);
if (angleMinusPiPiOrSo < -Math.PI) {
return -Math.PI;
} else if (angleMinusPiPiOrSo > Math.PI) {
return Math.PI;
} else {
return angleMinusPiPiOrSo;
}
}
/**
* Not accurate for large values.
*
* @param angle Angle in radians.
* @return The same angle, in radians, but in [-Math.PI,Math.PI].
*/
public static double normalizeMinusPiPiFast(double angle) {
// Not modifying values in output range.
if ((angle >= -Math.PI) && (angle <= Math.PI)) {
return angle;
}
double angleMinusPiPiOrSo = remainderTwoPiFast(angle);
if (angleMinusPiPiOrSo < -Math.PI) {
return -Math.PI;
} else if (angleMinusPiPiOrSo > Math.PI) {
return Math.PI;
} else {
return angleMinusPiPiOrSo;
}
}
/**
* @param angle Angle in radians.
* @return The same angle, in radians, but in [0,2*Math.PI].
*/
public static double normalizeZeroTwoPi(double angle) {
// Not modifying values in output range.
if ((angle >= 0.0) && (angle <= 2*Math.PI)) {
return angle;
}
double angleMinusPiPiOrSo = remainderTwoPi(angle);
if (angleMinusPiPiOrSo < 0.0) {
// Not a problem if angle is slightly < -Math.PI,
// since result ends up around PI, which is not near output range borders.
return angleMinusPiPiOrSo + 2*Math.PI;
} else {
// Not a problem if angle is slightly > Math.PI,
// since result ends up around PI, which is not near output range borders.
return angleMinusPiPiOrSo;
}
}
/**
* Not accurate for large values.
*
* @param angle Angle in radians.
* @return The same angle, in radians, but in [0,2*Math.PI].
*/
public static double normalizeZeroTwoPiFast(double angle) {
// Not modifying values in output range.
if ((angle >= 0.0) && (angle <= 2*Math.PI)) {
return angle;
}
double angleMinusPiPiOrSo = remainderTwoPiFast(angle);
if (angleMinusPiPiOrSo < 0.0) {
// Not a problem if angle is slightly < -Math.PI,
// since result ends up around PI, which is not near output range borders.
return angleMinusPiPiOrSo + 2*Math.PI;
} else {
// Not a problem if angle is slightly > Math.PI,
// since result ends up around PI, which is not near output range borders.
return angleMinusPiPiOrSo;
}
}
/**
* @param angle Angle in radians.
* @return Angle value modulo PI, in radians, in [-Math.PI/2,Math.PI/2].
*/
public static double normalizeMinusHalfPiHalfPi(double angle) {
// Not modifying values in output range.
if ((angle >= -Math.PI/2) && (angle <= Math.PI/2)) {
return angle;
}
double angleMinusPiPiOrSo = remainderTwoPi(angle);
if (angleMinusPiPiOrSo < -Math.PI/2) {
// Not a problem if angle is slightly < -Math.PI,
// since result ends up around zero, which is not near output range borders.
return angleMinusPiPiOrSo + Math.PI;
} else if (angleMinusPiPiOrSo > Math.PI/2) {
// Not a problem if angle is slightly > Math.PI,
// since result ends up around zero, which is not near output range borders.
return angleMinusPiPiOrSo - Math.PI;
} else {
return angleMinusPiPiOrSo;
}
}
/**
* Not accurate for large values.
*
* @param angle Angle in radians.
* @return Angle value modulo PI, in radians, in [-Math.PI/2,Math.PI/2].
*/
public static double normalizeMinusHalfPiHalfPiFast(double angle) {
// Not modifying values in output range.
if ((angle >= -Math.PI/2) && (angle <= Math.PI/2)) {
return angle;
}
double angleMinusPiPiOrSo = remainderTwoPiFast(angle);
if (angleMinusPiPiOrSo < -Math.PI/2) {
// Not a problem if angle is slightly < -Math.PI,
// since result ends up around zero, which is not near output range borders.
return angleMinusPiPiOrSo + Math.PI;
} else if (angleMinusPiPiOrSo > Math.PI/2) {
// Not a problem if angle is slightly > Math.PI,
// since result ends up around zero, which is not near output range borders.
return angleMinusPiPiOrSo - Math.PI;
} else {
return angleMinusPiPiOrSo;
}
}
/**
* Returns sqrt(x^2+y^2) without intermediate overflow or underflow.
*/
public static double hypot(double x, double y) {
x = FastMath.abs(x);
y = FastMath.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 * FastMath.sqrt(x*x+y*y);
}
}
/**
* Supposed to behave like Math.ceil(double), for safe interchangeability.
*
* @param value A double value.
* @return Ceiling of value.
*/
public static double ceil(double value) {
return -FastMath.floor(-value);
}
/**
* @param value A float value.
* @return Ceiling of value.
*/
public static float ceil(float value) {
return -FastMath.floor(-value);
}
/**
* Supposed to behave like Math.floor(double), for safe interchangeability.
*
* @param value A double value.
* @return Floor of value.
*/
public static double floor(double value) {
// Faster than to work directly on bits, as done for floor(float).
if (FastMath.abs(value) <= (double)Integer.MAX_VALUE) {
if (value > 0.0) {
return (double)(int)value;
} else if (value < 0.0) {
double anteComaDigits = (double)(int)value;
if (value != anteComaDigits) {
return anteComaDigits - 1.0;
} else {
return anteComaDigits;
}
} else { // value is +-0.0 (not NaN due to test against Integer.MAX_VALUE)
return value;
}
} else if (FastMath.abs(value) < TWO_POW_52) {
// Possible repartitions of the 17 decimal digits around coma:
// xxxxxxxxxx.yyyyyyy
// xxxxxxxxxxx.yyyyyy
// (...)
// xxxxxxxxxxxxxxxx.y
double highPart = ((int)(value * TWO_POW_N26)) * TWO_POW_26;
if (value > 0.0) {
return highPart + (double)((int)(value - highPart));
} else {
double anteComaDigits = highPart + (double)((int)(value - highPart));
if (value != anteComaDigits) {
return anteComaDigits - 1.0;
} else {
return anteComaDigits;
}
}
} else { // abs(value) >= 2^52, or value is NaN
return value;
}
}
/**
* @param value A float value.
* @return Floor of value.
*/
public static float floor(float value) {
int exp = FastMath.getExponent(value);
if (exp < 0) {
if (value < 0.0f) {
return -1.0f;
} else { // value in [0.0f,1.0f[
return 0.0f * value; // 0.0f, or -0.0f if value is -0.0f
}
} else {
if (exp < 24) {
int valueBits = Float.floatToRawIntBits(value);
int anteComaDigits = valueBits & (0xFF800000>>exp);
if ((value < 0.0f) && (anteComaDigits != valueBits)) {
return Float.intBitsToFloat(anteComaDigits) - 1.0f;
} else {
return Float.intBitsToFloat(anteComaDigits);
}
} else {
return value;
}
}
}
/**
* Supposed to behave like Math.round(double), for safe interchangeability.
*
* @param value A double value.
* @return Value rounded to nearest long.
*/
public static long round(double value) {
// Would be more coherent with rint, to call rint(double) instead of
// floor(double), but that would not give same results than Math.round(double).
double roundedValue = FastMath.floor(value+0.5);
if (FastMath.abs(roundedValue) <= (double)Integer.MAX_VALUE) {
// Faster with intermediary cast in int.
return (long)(int)roundedValue;
} else {
return (long)roundedValue;
}
}
/**
* Supposed to behave like Math.round(float), for safe interchangeability.
*
* @param value A double value.
* @return Value rounded to nearest int.
*/
public static int round(float value) {
// "return (int)FastMath.floor((float)(value+0.5));" would be more accurate for values in [8388609.0f,16777216.0f]
// (i.e. [0x800001,0x1000000]), but would not give same results than Math.round(float).
return (int)FastMath.floor(value+0.5f);
}
/**
* @param value A double value.
* @return Value unbiased exponent.
*/
public static int getExponent(double value) {
return (((int)(Double.doubleToRawLongBits(value)>>52))&0x7FF)-MAX_DOUBLE_EXPONENT;
}
/**
* @param value A float value.
* @return Value unbiased exponent.
*/
public static int getExponent(float value) {
return ((Float.floatToRawIntBits(value)>>23)&0xFF)-127;
}
/**
* Gives same result as Math.toDegrees for some particular values
* like Math.PI/2, Math.PI or 2*Math.PI, but is faster (no division).
*
* @param angrad Angle value in radians.
* @return Angle value in degrees.
*/
public static double toDegrees(double angrad) {
return angrad * (180/Math.PI);
}
/**
* Gives same result as Math.toRadians for some particular values
* like 90.0, 180.0 or 360.0, but is faster (no division).
*
* @param angdeg Angle value in degrees.
* @return Angle value in radians.
*/
public static double toRadians(double angdeg) {
return angdeg * (Math.PI/180);
}
/**
* @param sign Sign of the angle: true for positive, false for negative.
* @param degrees Degrees, in [0,180].
* @param minutes Minutes, in [0,59].
* @param seconds Seconds, in [0.0,60.0[.
* @return Angle in radians.
*/
public static double toRadians(boolean sign, int degrees, int minutes, double seconds) {
return FastMath.toRadians(FastMath.toDegrees(sign, degrees, minutes, seconds));
}
/**
* @param sign Sign of the angle: true for positive, false for negative.
* @param degrees Degrees, in [0,180].
* @param minutes Minutes, in [0,59].
* @param seconds Seconds, in [0.0,60.0[.
* @return Angle in degrees.
*/
public static double toDegrees(boolean sign, int degrees, int minutes, double seconds) {
double signFactor = sign ? 1.0 : -1.0;
return signFactor * (degrees + (1.0/60)*(minutes + (1.0/60)*seconds));
}
/**
* @param angrad Angle in radians.
* @param degrees (out) Degrees, in [0,180].
* @param minutes (out) Minutes, in [0,59].
* @param seconds (out) Seconds, in [0.0,60.0[.
* @return True if the resulting angle in [-180deg,180deg] is positive, false if it is negative.
*/
public static boolean toDMS(double angrad, IntWrapper degrees, IntWrapper minutes, DoubleWrapper seconds) {
// Computing longitude DMS.
double tmp = FastMath.toDegrees(FastMath.normalizeMinusPiPi(angrad));
boolean isNeg = (tmp < 0.0);
if (isNeg) {
tmp = -tmp;
}
degrees.value = (int)tmp;
tmp = (tmp-degrees.value)*60.0;
minutes.value = (int)tmp;
seconds.value = Math.min((tmp-minutes.value)*60.0,DOUBLE_BEFORE_60);
return !isNeg;
}
/**
* @param value An int value.
* @return The absolute value, except if value is Integer.MIN_VALUE, for which it returns Integer.MIN_VALUE.
*/
public static int abs(int value) {
return (value^(value>>31))-(value>>31);
}
/**
* @param value A long value.
* @return The closest int value in [Integer.MIN_VALUE,Integer.MAX_VALUE] range.
*/
public static int toInt(long value) {
if (value < (long)Integer.MIN_VALUE) {
return Integer.MIN_VALUE;
} else if (value > (long)Integer.MAX_VALUE) {
return Integer.MAX_VALUE;
} else {
return (int)value;
}
}
/**
* @param value A long value.
* @return If value is in [Integer.MIN_VALUE,Integer.MAX_VALUE] range, this value as int,
* otherwise throws an exception.
*/
public static int toIntSafe(long value) {
if ((value < (long)Integer.MIN_VALUE) || (value > (long)Integer.MAX_VALUE)) {
throw new ArithmeticException("overflow: "+value);
} else {
return (int)value;
}
}
/**
* @param a An int value.
* @param b An int value.
* @return The int value of [Integer.MIN_VALUE,Integer.MAX_VALUE] range which is the closest to mathematical result of a+b.
*/
public static int plusNoModulo(int a, int b) {
return FastMath.toInt(((long)a) + ((long)b));
}
/**
* @param a An int value.
* @param b An int value.
* @return The mathematical result of a+b, if it is in [Integer.MIN_VALUE,Integer.MAX_VALUE] range,
* otherwise throws an exception.
*/
public static int plusNoModuloSafe(int a, int b) {
if ((a^b) < 0) { // test if a and b signs are different
return a + b;
} else {
int sum = a + b;
if ((a^sum) < 0) {
throw new ArithmeticException("overflow: "+a+"+"+b);
} else {
return sum;
}
}
}
/**
* @param a A long value.
* @param b A long value.
* @return The long value of [Long.MIN_VALUE,Long.MAX_VALUE] range which is the closest to mathematical result of a+b.
*/
public static long plusNoModulo(long a, long b) {
// Algorithm tested on int type, for which it's faster with cast to long.
if ((a^b) < 0) { // test if a and b signs are different
return a + b;
} else {
long sum = a + b;
if ((a^sum) < 0) {
return (sum >= 0) ? Long.MIN_VALUE : Long.MAX_VALUE;
} else {
return sum;
}
}
}
/**
* @param a A long value.
* @param b A long value.
* @return The mathematical result of a+b, if it is in [Long.MIN_VALUE,Long.MAX_VALUE] range,
* otherwise throws an exception.
*/
public static long plusNoModuloSafe(long a, long b) {
// Algorithm tested on int type.
if ((a^b) < 0) { // test if a and b signs are different
return a + b;
} else {
long sum = a + b;
if ((a^sum) < 0) {
throw new ArithmeticException("overflow: "+a+"+"+b);
} else {
return sum;
}
}
}
/**
* @param a An int value.
* @param b An int value.
* @return The int value of [Integer.MIN_VALUE,Integer.MAX_VALUE] range which is the closest to mathematical result of a-b.
*/
public static int minusNoModulo(int a, int b) {
return FastMath.toInt(((long)a) - ((long)b));
}
/**
* @param a An int value.
* @param b An int value.
* @return The mathematical result of a-b, if it is in [Integer.MIN_VALUE,Integer.MAX_VALUE] range,
* otherwise throws an exception.
*/
public static int minusNoModuloSafe(int a, int b) {
if ((a^b) >= 0) { // test if a and b signs are identical
return a - b;
} else {
int diff = a - b;
if ((a^diff) < 0) {
throw new ArithmeticException("overflow: "+a+"-"+b);
} else {
return diff;
}
}
}
/**
* @param a A long value.
* @param b A long value.
* @return The long value of [Long.MIN_VALUE,Long.MAX_VALUE] range which is the closest to mathematical result of a-b.
*/
public static long minusNoModulo(long a, long b) {
// Algorithm tested on int type, for which it's faster with cast to long.
if ((a^b) >= 0) { // test if a and b signs are identical
return a - b;
} else {
long diff = a - b;
if ((a^diff) < 0) {
return (diff >= 0) ? Long.MIN_VALUE : Long.MAX_VALUE;
} else {
return diff;
}
}
}
/**
* @param a A long value.
* @param b A long value.
* @return The mathematical result of a-b, if it is in [Long.MIN_VALUE,Long.MAX_VALUE] range,
* otherwise throws an exception.
*/
public static long minusNoModuloSafe(long a, long b) {
// Algorithm tested on int type.
if ((a^b) >= 0) { // test if a and b signs are identical
return a - b;
} else {
long diff = a - b;
if ((a^diff) < 0) {
throw new ArithmeticException("overflow: "+a+"-"+b);
} else {
return diff;
}
}
}
/**
* @param a An int value.
* @param b An int value.
* @return The int value of [Integer.MIN_VALUE,Integer.MAX_VALUE] range which is the closest to mathematical result of a*b.
*/
public static int timesNoModulo(int a, int b) {
if (false) {
/*
* slower
*/
return FastMath.toInt(((long)a) * ((long)b));
}
return (int)(a * (double)b);
}
/**
* @param a An int value.
* @param b An int value.
* @return The mathematical result of a*b, if it is in [Integer.MIN_VALUE,Integer.MAX_VALUE] range,
* otherwise throws an exception.
*/
public static int timesNoModuloSafe(int a, int b) {
if (false) {
/*
* slower
*/
if (b == 0) {
return 0;
}
int product = a * b;
if ((product == Integer.MIN_VALUE) && ((a^b) >= 0)) {
// product negative, but a and b have the same sign:
// that means total would be -Integer.MIN_VALUE,
// which does not exist as int.
throw new ArithmeticException("overflow: "+a+"*"+b);
} else {
if (product / b != a) {
throw new ArithmeticException("overflow: "+a+"*"+b);
} else {
return product;
}
}
}
double product = a * (double)b;
if ((product >= (double)Integer.MIN_VALUE) && (product <= (double)Integer.MAX_VALUE)) {
return (int)product;
} else {
throw new ArithmeticException("overflow: "+a+"*"+b);
}
}
/**
* @param a A long value.
* @param b A long value.
* @return The long value of [Long.MIN_VALUE,Long.MAX_VALUE] range which is the closest to mathematical result of a*b.
*/
public static long timesNoModulo(long a, long b) {
// Algorithm tested on int type, for which it's faster to do another way.
if (b == 0) {
return 0;
}
long product = a * b;
if ((product == Long.MIN_VALUE) && ((a^b) >= 0)) {
// product negative, but a and b have the same sign:
// that means total would be -Long.MIN_VALUE,
// which does not exist as long.
return Long.MAX_VALUE;
} else {
if (product / b != a) {
return ((a^b) >= 0) ? Long.MAX_VALUE : Long.MIN_VALUE;
} else {
return product;
}
}
}
/**
* @param a A long value.
* @param b A long value.
* @return The mathematical result of a*b, if it is in [Long.MIN_VALUE,Long.MAX_VALUE] range,
* otherwise throws an exception.
*/
public static long timesNoModuloSafe(long a, long b) {
// Algorithm tested on int type.
if (b == 0) {
return 0;
}
long product = a * b;
if ((product == Long.MIN_VALUE) && ((a^b) >= 0)) {
// product negative, but a and b have the same sign:
// that means total would be -Long.MIN_VALUE,
// which does not exist as long.
throw new ArithmeticException("overflow: "+a+"*"+b);
} else {
if (product / b != a) {
throw new ArithmeticException("overflow: "+a+"*"+b);
} else {
return product;
}
}
}
/**
* @param minValue An int value.
* @param maxValue An int value.
* @param value An int value.
* @return minValue if value < minValue, maxValue if value > maxValue, value otherwise.
*/
public static int inRange(int minValue, int maxValue, int value) {
if (value < minValue) {
return minValue;
} else if (value > maxValue) {
return maxValue;
} else {
return value;
}
}
/**
* @param minValue A long value.
* @param maxValue A long value.
* @param value A long value.
* @return minValue if value < minValue, maxValue if value > maxValue, value otherwise.
*/
public static long inRange(long minValue, long maxValue, long value) {
if (value < minValue) {
return minValue;
} else if (value > maxValue) {
return maxValue;
} else {
return value;
}
}
/**
* @param minValue A float value.
* @param maxValue A float value.
* @param value A float value.
* @return minValue if value < minValue, maxValue if value > maxValue, value otherwise. * @published
*/
public static float inRange(float minValue, float maxValue, float value) {
if (value < minValue) {
return minValue;
} else if (value > maxValue) {
return maxValue;
} else {
return value;
}
}
/**
* @param minValue A double value.
* @param maxValue A double value.
* @param value A double value.
* @return minValue if value < minValue, maxValue if value > maxValue, value otherwise.
*/
public static double inRange(double minValue, double maxValue, double value) {
if (value < minValue) {
return minValue;
} else if (value > maxValue) {
return maxValue;
} else {
return value;
}
}
/**
* NB: Since 2*Math.PI < 2*PI, a span of 2*Math.PI does not mean full angular range.
* ex.: isInClockwiseDomain(0.0, 2*Math.PI, -1e-20) returns false.
* ---> For full angular range, use a span > 2*Math.PI.
*
* @param startAngRad An angle, in radians.
* @param angSpanRad An angular span, >= 0.0, in radians.
* @param angRad An angle, in radians.
* @return True if angRad is in the clockwise angular domain going from startAngRad, over angSpanRad,
* extremities included, false otherwise.
*/
public static boolean isInClockwiseDomain(double startAngRad, double angSpanRad, double angRad) {
if (FastMath.abs(angRad) < -TWO_MATH_PI_IN_MINUS_PI_PI) {
// special case for angular values of small magnitude
if (angSpanRad < 0.0) {
// empty domain
return false;
} else if (angSpanRad <= 2*Math.PI) { // angSpanRad is in [0.0,2*Math.PI]
startAngRad = FastMath.normalizeMinusPiPi(startAngRad);
double endAngRad = FastMath.normalizeMinusPiPi(startAngRad + angSpanRad);
//
if (startAngRad <= endAngRad) {
return (angRad >= startAngRad) && (angRad <= endAngRad);
} else {
return (angRad >= startAngRad) || (angRad <= endAngRad);
}
} else if (angSpanRad != angSpanRad) { // angSpanRad is NaN
return false;
} else { // angSpanRad > 2*Math.PI
// we know angRad is not NaN, due to a previous test
return true;
}
} else {
// general case
return (FastMath.normalizeZeroTwoPi(angRad - startAngRad) <= angSpanRad);
}
}
public static boolean isNaNOrInfinite(double value) {
// value-value is not equal to 0.0 (and is NaN) <-> value is NaN or +-infinity
return !(value-value == 0.0);
}
/*
*
* Not-redefined java.lang.Math public values and treatments, for quick replacement of Math with FastMath.
*
*/
public static final double E = Math.E;
public static final double PI = Math.PI;
public static double abs(double a) {
return Math.abs(a);
}
public static float abs(float a) {
return Math.abs(a);
}
public static long abs(long a) {
return Math.abs(a);
}
public static double copySign(double magnitude, double sign) {
return Math.copySign(magnitude,sign);
}
public static float copySign(float magnitude, float sign) {
return Math.copySign(magnitude,sign);
}
public static double IEEEremainder(double f1, double f2) {
return Math.IEEEremainder(f1,f2);
}
public static double log10(double a) {
return Math.log10(a);
}
public static double max(double a, double b) {
return Math.max(a,b);
}
public static float max(float a, float b) {
return Math.max(a,b);
}
public static int max(int a, int b) {
return Math.max(a,b);
}
public static long max(long a, long b) {
return Math.max(a,b);
}
public static double min(double a, double b) {
return Math.min(a,b);
}
public static float min(float a, float b) {
return Math.min(a,b);
}
public static int min(int a, int b) {
return Math.min(a,b);
}
public static long min(long a, long b) {
return Math.min(a,b);
}
public static double nextAfter(double start, double direction) {
return Math.nextAfter(start,direction);
}
public static float nextAfter(float start, float direction) {
return Math.nextAfter(start,direction);
}
public static double nextUp(double d) {
return Math.nextUp(d);
}
public static float nextUp(float f) {
return Math.nextUp(f);
}
public static double random() {
return Math.random();
}
public static double rint(double a) {
return Math.rint(a);
}
public static double scalb(double d, int scaleFactor) {
return Math.scalb(d,scaleFactor);
}
public static float scalb(float f, int scaleFactor) {
return Math.scalb(f,scaleFactor);
}
public static double signum(double d) {
return Math.signum(d);
}
public static float signum(float f) {
return Math.signum(f);
}
public static double ulp(double d) {
return Math.ulp(d);
}
public static float ulp(float f) {
return Math.ulp(f);
}
//--------------------------------------------------------------------------
// PRIVATE TREATMENTS
//--------------------------------------------------------------------------
/**
* FastMath is non-instantiable.
*/
private FastMath() {
}
/**
* 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 = (double)(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 = (double)((int)z);
z = (z-x0)*TWO_POW_24;
double x1 = (double)((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;
}
}
/**
* Not accurate for large values.
*
* 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 remainderTwoPiFast(double angle) {
boolean negateResult;
if (angle < 0.0) {
negateResult = true;
angle = -angle;
} else {
negateResult = false;
}
// - We don't bother with values higher than (2*PI*(2^52)),
// since they are spaced by 2*PI or more from each other.
// - For large values, we don't use % because it might be very slow,
// and we split computation in two, because cast from double to int
// with large numbers might be very slow also.
if (angle <= TWO_POW_26*(2*Math.PI)) {
double fn = (double)(int)(angle*INVTWOPI+0.5);
double result = (angle - fn*TWOPI_HI) - fn*TWOPI_LO;
return (negateResult) ? -result : result;
} else if (angle <= TWO_POW_52*(2*Math.PI)) {
// 1) Computing remainder of angle modulo TWO_POW_26*(2*PI).
double fn = (double)(int)(angle*(INVTWOPI/TWO_POW_26)+0.5);
double result = (angle - fn*(TWOPI_HI*TWO_POW_26)) - fn*(TWOPI_LO*TWO_POW_26);
// Here, result is in [-TWO_POW_26*Math.PI,TWO_POW_26*Math.PI].
if (result < 0.0) {
result = -result;
negateResult = !negateResult;
}
// 2) Computing remainder of angle modulo 2*PI.
fn = (double)(int)(result*INVTWOPI+0.5);
result = (result - fn*TWOPI_HI) - fn*TWOPI_LO;
return (negateResult) ? -result : result;
} else if (angle < Double.POSITIVE_INFINITY) {
return 0.0;
} 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 = (double)((int)(TWO_POW_N24*z));
iq0 = (int)(z-TWO_POW_24*fw);
z = q3+fw;
fw = (double)((int)(TWO_POW_N24*z));
iq1 = (int)(z-TWO_POW_24*fw);
z = q2+fw;
fw = (double)((int)(TWO_POW_N24*z));
iq2 = (int)(z-TWO_POW_24*fw);
z = q1+fw;
fw = (double)((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 we can use the table right away.
double twoPowQ = twoPowTab[q-MIN_DOUBLE_EXPONENT];
z = (z*twoPowQ) % 8.0;
z -= (double)((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 = (double)((int)(TWO_POW_N24*z));
iq0 = (int)(z-TWO_POW_24*fw);
z = q4+fw;
fw = (double)((int)(TWO_POW_N24*z));
iq1 = (int)(z-TWO_POW_24*fw);
z = q3+fw;
fw = (double)((int)(TWO_POW_N24*z));
iq2 = (int)(z-TWO_POW_24*fw);
z = q2+fw;
fw = (double)((int)(TWO_POW_N24*z));
iq3 = (int)(z-TWO_POW_24*fw);
z = q1+fw;
fw = (double)((int)(TWO_POW_N24*z));
iq4 = (int)(z-TWO_POW_24*fw);
z = q0+fw;
z = (z*twoPowQ) % 8.0;
z -= (double)((int)z);
if (q > 0) {
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 {
// Here, q is in [-25,-2] range or so, so we could use twoPow's table right away with
// iq4 = (int)(z*twoPowTab[-q-TWO_POW_TAB_MIN_POW]);
// but tests show using division is faster...
iq4 = (int)(z/twoPowQ);
fw = twoPowQ;
}
q4 = fw*(double)iq4;
fw *= TWO_POW_N24;
q3 = fw*(double)iq3;
fw *= TWO_POW_N24;
q2 = fw*(double)iq2;
fw *= TWO_POW_N24;
q1 = fw*(double)iq1;
fw *= TWO_POW_N24;
q0 = fw*(double)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
//--------------------------------------------------------------------------
/**
* Initializes look-up tables.
*
* Using some FastMath methods in there, not to spend
* an hour in it. Must take care not to use methods
* which look-up tables have not yet been initialized,
* or that are not accurate enough.
*/
static {
// 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= EXP_MIN_INT_LIMIT) {
expHiInvTab[i] = Math.exp(-i);
} else {
expHiInvTab[i] = Math.exp(54*LOG_2-i);
}
}
// log
for (int i=0;i>SQRT_LO_BITS));
for (int i=1;i>CBRT_LO_BITS));
for (int i=1;i