package org.math.plot.utils;

/*
 * Copyright 2012 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.
 * =============================================================================
 */

/**
 * 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.
 * 
 * Properties:
 * 
 * - jodk.fastmath.strict (boolean, default is true):
 *   If true, non-redefined Math methods which results could vary between Math and StrictMath,
 *   delegate to StrictMath, and if false, to Math.
 *   Default is true to ensure consistency across various architectures.
 *   
 * - jodk.fastmath.usejdk (boolean, default is false):
 *   If true, redefined Math methods, as well as their "Fast" or "Quick" terminated counterparts,
 *   delegate to StrictMath or Math, depending on jodk.fastmath.strict property.
 *   
 * - jodk.fastmath.fastlog (boolean, default is true):
 *   If true, using redefined log(double), if false using StrictMath.log(double) or
 *   Math.log(double), depending on jodk.fastmath.strict property.
 *   Default is true because jodk.fastmath.strict is true by default, and StrictMath.log(double)
 *   seems usually slow.
 *   
 * - jodk.fastmath.fastsqrt (boolean, default is false):
 *   If true, using redefined sqrt(double), if false using StrictMath.sqrt(double) or
 *   Math.sqrt(double), depending on jodk.fastmath.strict property.
 *   Default is false because StrictMath.sqrt(double) seems to usually delegate to hardware sqrt.
 * 
 * Unless jodk.fastmath.strict is false and jodk.fastmath.usejdk is true, these treatments
 * are consistent across various architectures, for constants and look-up tables are
 * computed with StrictMath, or exact Math methods.
 * 
 * --- 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 {

    /*
     * 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 <double_array>[<index>], 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 (<a_constant> * 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.
     */
    
public class DoubleWrapper {
    public double value;
    @Override
    public String toString() {
        return Double.toString(this.value);
    }
}

public class IntWrapper {
    public int value;
    @Override
    public String toString() {
        return Integer.toString(this.value);
    }
}
    //--------------------------------------------------------------------------
    // CONFIGURATION
    //--------------------------------------------------------------------------

    private static final boolean STRICT_MATH = false;//LangUtils.getBooleanProperty("jodk.fastmath.strict", true);

    private static final boolean USE_JDK_MATH = false;//LangUtils.getBooleanProperty("jodk.fastmath.usejdk", false);

    /**
     * Used for both log(double) and log10(double).
     */
    private static final boolean USE_REDEFINED_LOG = true;//LangUtils.getBooleanProperty("jodk.fastmath.fastlog", true);

    private static final boolean USE_REDEFINED_SQRT = true;//LangUtils.getBooleanProperty("jodk.fastmath.fastsqrt", 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 = true;

    //--------------------------------------------------------------------------
    // GENERAL CONSTANTS
    //--------------------------------------------------------------------------

    /**
     * High approximation of PI, which is further from PI
     * than the low approximation Math.PI:
     *              PI ~= 3.14159265358979323846...
     *         Math.PI ~= 3.141592653589793
     * FastMath.PI_SUP ~= 3.1415926535897936
     */
    public static final double PI_SUP = Math.nextUp(Math.PI);

    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 int MAX_FLOAT_EXPONENT = 127;

    private static final double LOG_2 = StrictMath.log(2.0);
    private static final double LOG_TWO_POW_27 = StrictMath.log(TWO_POW_27);
    private static final double LOG_DOUBLE_MAX_VALUE = StrictMath.log(Double.MAX_VALUE);

    private static final double INV_LOG_10 = 1.0/StrictMath.log(10.0);

    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<i>.
     */
    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);

    /**
     * 2*Math.PI, normalized into [-PI,PI].
     * Computed using normalizeMinusPiPi(double).
     */
    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<<getTabSizePower(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<PI/2.

    // We use indexing past look-up tables, so that indexing information
    // allows for fast recomputation of angle in [0,PI/2] range.
    private static final int TAN_VIRTUAL_TABS_SIZE = (1<<getTabSizePower(12)) + 1;

    // Must be >= 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 / TAN_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 >= 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<<getTabSizePower(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 = StrictMath.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<<getTabSizePower(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 = StrictMath.tan(Math.toRadians(74.0));

    private static final int ATAN_TABS_SIZE = (1<<getTabSizePower(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 

    //--------------------------------------------------------------------------
    // 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<<EXP_LO_DISTANCE_TO_ZERO_POT);
    private static final int EXP_LO_TAB_SIZE_POT = getTabSizePower(11);
    private static final int EXP_LO_TAB_SIZE = (1<<EXP_LO_TAB_SIZE_POT)+1;
    private static final int EXP_LO_TAB_MID_INDEX = ((EXP_LO_TAB_SIZE-1)/2);
    private static final int EXP_LO_INDEXING = EXP_LO_TAB_MID_INDEX/EXP_LO_DISTANCE_TO_ZERO;
    private static final int EXP_LO_INDEXING_DIV_SHIFT = EXP_LO_TAB_SIZE_POT-1-EXP_LO_DISTANCE_TO_ZERO_POT;
    private static final double[] expHiTab = new double[1+(int)EXP_OVERFLOW_LIMIT];
    private static final double[] expHiInvTab = new double[1-(int)EXP_UNDERFLOW_LIMIT];
    private static final double[] expLoPosTab = new double[EXP_LO_TAB_SIZE];
    private static final double[] expLoNegTab = new double[EXP_LO_TAB_SIZE];

    //--------------------------------------------------------------------------
    // CONSTANTS FOR QUICK EXP
    //--------------------------------------------------------------------------

    private static final double EXP_QUICK_A = TWO_POW_52/LOG_2;
    private static final double EXP_QUICK_B = MAX_DOUBLE_EXPONENT * TWO_POW_52;
    private static final double EXP_QUICK_C = Math.ceil((StrictMath.log(LOG_2+2/Math.E) - LOG_2 - StrictMath.log(LOG_2)) * EXP_QUICK_A);

    //--------------------------------------------------------------------------
    // CONSTANTS AND TABLES FOR LOG AND LOG1P
    //--------------------------------------------------------------------------

    private static final int LOG_BITS = getTabSizePower(12);
    private static final int LOG_TAB_SIZE = (1<<LOG_BITS);
    private static final double[] logXLogTab = new double[LOG_TAB_SIZE];
    private static final double[] logXTab = new double[LOG_TAB_SIZE];
    private static final double[] logXInvTab = new double[LOG_TAB_SIZE];

    //--------------------------------------------------------------------------
    // TABLE FOR POWERS OF TWO
    //--------------------------------------------------------------------------

    private static final double[] twoPowTab = new double[(MAX_DOUBLE_EXPONENT-MIN_DOUBLE_EXPONENT)+1];

    //--------------------------------------------------------------------------
    // CONSTANTS AND TABLES FOR SQRT
    //--------------------------------------------------------------------------

    private static final int SQRT_LO_BITS = getTabSizePower(12);
    private static final int SQRT_LO_TAB_SIZE = (1<<SQRT_LO_BITS);
    private static final double[] sqrtXSqrtHiTab = new double[MAX_DOUBLE_EXPONENT-MIN_DOUBLE_EXPONENT+1];
    private static final double[] sqrtXSqrtLoTab = new double[SQRT_LO_TAB_SIZE];
    private static final double[] sqrtSlopeHiTab = new double[MAX_DOUBLE_EXPONENT-MIN_DOUBLE_EXPONENT+1];
    private static final double[] sqrtSlopeLoTab = new double[SQRT_LO_TAB_SIZE];

    //--------------------------------------------------------------------------
    // CONSTANTS AND TABLES FOR CBRT
    //--------------------------------------------------------------------------

    private static final int CBRT_LO_BITS = getTabSizePower(12);
    private static final int CBRT_LO_TAB_SIZE = (1<<CBRT_LO_BITS);
    // For CBRT_LO_BITS = 12:
    // cbrtXCbrtLoTab[0] = 1.0.
    // cbrtXCbrtLoTab[1] = cbrt(1. 000000000000 1111111111111111111111111111111111111111b)
    // cbrtXCbrtLoTab[2] = cbrt(1. 000000000001 1111111111111111111111111111111111111111b)
    // cbrtXCbrtLoTab[3] = cbrt(1. 000000000010 1111111111111111111111111111111111111111b)
    // cbrtXCbrtLoTab[4] = cbrt(1. 000000000011 1111111111111111111111111111111111111111b)
    // etc.
    private static final double[] cbrtXCbrtHiTab = new double[MAX_DOUBLE_EXPONENT-MIN_DOUBLE_EXPONENT+1];
    private static final double[] cbrtXCbrtLoTab = new double[CBRT_LO_TAB_SIZE];
    private static final double[] cbrtSlopeHiTab = new double[MAX_DOUBLE_EXPONENT-MIN_DOUBLE_EXPONENT+1];
    private static final double[] cbrtSlopeLoTab = new double[CBRT_LO_TAB_SIZE];

    //--------------------------------------------------------------------------
    // PUBLIC TREATMENTS
    //--------------------------------------------------------------------------

    /**
     * @param angle Angle in radians.
     * @return Angle cosine.
     */
    public static double cos(double angle) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.cos(angle) : Math.cos(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)));
    }

    /**
     * Quick cosine, with accuracy of about 1.6e-3 (PI/'look-up tabs size')
     * for |angle| &lt; 6588397.0 (Integer.MAX_VALUE * (2*PI/'look-up tabs size')),
     * and no accuracy at all for larger values.
     * 
     * @param angle Angle in radians.
     * @return Angle cosine.
     */
    public static double cosQuick(double angle) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.cos(angle) : Math.cos(angle);
        }
        return cosTab[((int)(Math.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) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.sin(angle) : Math.sin(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.6e-3 (PI/'look-up tabs size')
     * for |angle| &lt; 6588397.0 (Integer.MAX_VALUE * (2*PI/'look-up tabs size')),
     * and no accuracy at all for larger values.
     * 
     * @param angle Angle in radians.
     * @return Angle sine.
     */
    public static double sinQuick(double angle) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.sin(angle) : Math.sin(angle);
        }
        return cosTab[((int)(Math.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) {
        if (USE_JDK_MATH) {
            sine.value = STRICT_MATH ? StrictMath.sin(angle) : Math.sin(angle);
            cosine.value = STRICT_MATH ? StrictMath.cos(angle) : Math.cos(angle);
            return;
        }
        // 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 (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.tan(angle) : Math.tan(angle);
        }
        if (Math.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) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.acos(value) : Math.acos(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) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.asin(value) : Math.asin(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.powFast(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) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.atan(value) : Math.atan(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 (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.atan2(y,x) : Math.atan2(y,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) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.cosh(value) : Math.cosh(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) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.sinh(value) : Math.sinh(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) {
        if (USE_JDK_MATH) {
            hsine.value = STRICT_MATH ? StrictMath.sinh(value) : Math.sinh(value);
            hcosine.value = STRICT_MATH ? StrictMath.cosh(value) : Math.cosh(value);
            return;
        }
        // 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) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.tanh(value) : Math.tanh(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;
    }

    /**
     * @param value A double value.
     * @return e^value.
     */
    public static double exp(double value) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.exp(value) : Math.exp(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 avoid 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_LO_INDEXING_DIV_SHIFT);
            double delta = (value-valueInt)-i*(1.0/EXP_LO_INDEXING);
            return expHiTab[valueInt] * (expLoPosTab[i+EXP_LO_TAB_MID_INDEX]*(1+delta*(1+delta*(1.0/2+delta*(1.0/6+delta*(1.0/24))))));
        } else { // value < 0.0, or value is NaN
            if (!(value >= 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,
            // 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_LO_INDEXING_DIV_SHIFT);
            double delta = (value-valueInt)-i*(1.0/EXP_LO_INDEXING);
            double tmp = expHiInvTab[-valueInt] * (expLoPosTab[i+EXP_LO_TAB_MID_INDEX]*(1+delta*(1+delta*(1.0/2+delta*(1.0/6+delta*(1.0/24))))));
            // We took care not to compute with subnormal values.
            return (valueInt >= EXP_MIN_INT_LIMIT) ? tmp : tmp * TWO_POW_N54;
        }
    }

    /**
     * Quick exp, with a max relative error of about 3e-2 for |value| &lt; 700.0 or so,
     * and no accuracy at all outside this range.
     * Derived from a note by Nicol N. Schraudolph, IDSIA, 1998.
     * 
     * @param value A double value.
     * @return e^value.
     */
    public static double expQuick(double value) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.exp(value) : Math.exp(value);
        }
        if (false) {
            // Schraudolph's original method.
            return Double.longBitsToDouble((long)(EXP_QUICK_A * value + (EXP_QUICK_B - EXP_QUICK_C)));
        }
        /*
         * Cast of double values, even in long range, into long, is slower than
         * from double to int for values in int range, and then from int to long.
         * For that reason, we only work with integer values in int range (corresponding to the 32 first bits of the long,
         * containing sign, exponent, and highest significant bits of double's mantissa), and cast twice.
         */
        return Double.longBitsToDouble(((long)(int)(EXP_QUICK_A/(1L<<32) * value + (EXP_QUICK_B - EXP_QUICK_C)/(1L<<32)))<<32);
    }

    /**
     * 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 (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.expm1(value) : Math.expm1(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_JDK_MATH || (!USE_REDEFINED_LOG)) {
            return STRICT_MATH ? StrictMath.log(value) : Math.log(value);
        } else {
            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 aside, 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;
            }
        }
    }

    /**
     * Quick log, with a max relative error of about 2.8e-4
     * for values in ]0,+infinity[, and no accuracy at all
     * outside this range.
     */
    public static double logQuick(double value) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.log(value) : Math.log(value);
        }
        /*
         * Inverse of Schraudolph's method for exp, is very inaccurate near 1,
         * and not that fast (even using floats), especially with added if's
         * to deal with values near 1, so we don't use it, and use a simplified
         * version of our log's redefined algorithm.
         */

        // Simplified version of log's redefined algorithm:
        // log(value) ~= exponent * log(2) + log(1.mantissaApprox)

        double h;
        if (value > 0.87) {
            if (value < 1.16) {
                return 2.0 * (value-1.0)/(value+1.0);
            }
            h = 0.0;
        } else if (value < MIN_DOUBLE_NORMAL) {
            value *= TWO_POW_52;
            h = -52*LOG_2;
        } else {
            h = 0.0;
        }

        int valueBitsHi = (int)(Double.doubleToRawLongBits(value)>>32);
        int valueExp = (valueBitsHi>>20)-MAX_DOUBLE_EXPONENT;
        int xIndex = ((valueBitsHi<<12)>>>(32-LOG_BITS));

        return h + valueExp * LOG_2 + logXLogTab[xIndex];
    }

    /**
     * @param value A double value.
     * @return Value logarithm (base 10).
     */
    public static double log10(double value) {
        if (USE_JDK_MATH || (!USE_REDEFINED_LOG)) {
            return STRICT_MATH ? StrictMath.log10(value) : Math.log10(value);
        } else {
            // INV_LOG_10 is < 1, but there is no risk of log(double)
            // overflow (positive or negative) while the end result shouldn't,
            // since log(Double.MIN_VALUE) and log(Double.MAX_VALUE) have
            // magnitudes of just a few hundreds.
            return log(value) * INV_LOG_10;
        }
    }

    /**
     * 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 (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.log1p(value) : Math.log1p(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 FastMath.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 An integer value in [1,Integer.MAX_VALUE].
     * @return The integer part of the logarithm, in base 2, of the specified value,
     *         i.e. a result in [0,30]
     * @throws IllegalArgumentException if the specified value is &lt;= 0.
     */
    public static int log2(int value) {
        return NumbersUtils.log2(value);
    }

    /**
     * @param value An integer value in [1,Long.MAX_VALUE].
     * @return The integer part of the logarithm, in base 2, of the specified value,
     *         i.e. a result in [0,62]
     * @throws IllegalArgumentException if the specified value is &lt;= 0.
     */
    public static int log2(long value) {
        return NumbersUtils.log2(value);
    }

    /**
     * 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 (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.pow(value,power) : Math.pow(value,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 comma 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.
                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 { // value > 0.0, or value is NaN
            return FastMath.exp(power*FastMath.log(value));
        }
    }

    /**
     * Quick pow, with a max relative error of about 3.5e-2
     * for |a^b| &lt; 1e10, of about 0.17 for |a^b| &lt; 1e50,
     * and worse accuracy above.
     * 
     * @param value A double value, in ]0,+infinity[ (strictly positive and finite).
     * @param power A double value.
     * @return value^power.
     */
    public static double powQuick(double value, double power) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.pow(value,power) : Math.pow(value,power);
        }
        return FastMath.exp(power*FastMath.logQuick(value));
    }

    /**
     * This treatment is somehow accurate for low values of |power|,
     * and for |power*getExponent(value)| &lt; 1023 or so (to stay away
     * from double extreme magnitudes (large and small)).
     * 
     * @param value A double value.
     * @param power A power.
     * @return value^power.
     */
    public static double powFast(double value, int power) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.pow(value,power) : Math.pow(value,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) {
                // Integer.MAX_VALUE = -(power+1)
                return 1.0/(FastMath.powFast(value,Integer.MAX_VALUE) * 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) {
        /*
         * Using table, to go faster than NumbersUtils.twoPow(int).
         */
        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;
            }
        }
    }

    /**
     * @param value An int value.
     * @return value*value.
     */
    public static int pow2(int value) {
        return NumbersUtils.pow2(value);
    }

    /**
     * @param value A long value.
     * @return value*value.
     */
    public static long pow2(long value) {
        return NumbersUtils.pow2(value);
    }

    /**
     * @param value A float value.
     * @return value*value.
     */
    public static float pow2(float value) {
        return NumbersUtils.pow2(value);
    }

    /**
     * @param value A double value.
     * @return value*value.
     */
    public static double pow2(double value) {
        return NumbersUtils.pow2(value);
    }

    /**
     * @param value An int value.
     * @return value*value*value.
     */
    public static int pow3(int value) {
        return NumbersUtils.pow3(value);
    }

    /**
     * @param value A long value.
     * @return value*value*value.
     */
    public static long pow3(long value) {
        return NumbersUtils.pow3(value);
    }

    /**
     * @param value A float value.
     * @return value*value*value.
     */
    public static float pow3(float value) {
        return NumbersUtils.pow3(value);
    }

    /**
     * @param value A double value.
     * @return value*value*value.
     */
    public static double pow3(double value) {
        return NumbersUtils.pow3(value);
    }

    /**
     * @param value A double value.
     * @return Value square root.
     */
    public static double sqrt(double value) {
        if (USE_JDK_MATH || (!USE_REDEFINED_SQRT)) {
            return STRICT_MATH ? StrictMath.sqrt(value) : Math.sqrt(value);
        } else {
            // 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*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);
        }
    }

    /**
     * @param value A double value.
     * @return Value cubic root.
     */
    public static double cbrt(double value) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.cbrt(value) : Math.cbrt(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 = <result_sign> * <result_multiplicator_to_avoid_overflow> / <cbrt(value_multiplicator_to_avoid_subnormal)>
                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^<value exponent> * <a value in [1,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 (or Archimedes'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).

        // Shift and cast into an int, which overall is faster than working with a long.
        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);
    }

    /**
     * Returns dividend - divisor * n, where n is the mathematical integer
     * closest to dividend/divisor.
     * If dividend/divisor is equally close to surrounding integers,
     * we choose n to be the integer of smallest magnitude, which makes
     * this treatment differ from Math.IEEEremainder(double,double),
     * where n is chosen to be the even integer.
     * Note that the choice of n is not done considering the double
     * approximation of dividend/divisor, because it could cause
     * result to be outside [-|divisor|/2,|divisor|/2] range.
     * The practical effect is that if multiple results would be possible,
     * we always choose the result that is the closest to (and has the same
     * sign as) the dividend.
     * Ex. :
     * - for (-3.0,2.0), this method returns -1.0,
     *   whereas Math.IEEEremainder returns 1.0.
     * - for (-5.0,2.0), both this method and Math.IEEEremainder return -1.0.
     * 
     * If the remainder is zero, its sign is the same as the sign of the first argument.
     * If either argument is NaN, or the first argument is infinite,
     * or the second argument is positive zero or negative zero,
     * then the result is NaN.
     * If the first argument is finite and the second argument is
     * infinite, then the result is the same as the first argument.
     * 
     * NB:
     * - Modulo operator (%) returns a value in ]-|divisor|,|divisor|[,
     *   which sign is the same as dividend.
     * - As for modulo operator, the sign of the divisor has no effect on the result.
     * 
     * @param dividend Dividend.
     * @param divisor Divisor.
     * @return Remainder of dividend/divisor, i.e. a value in [-|divisor|/2,|divisor|/2].
     */
    public static double remainder(double dividend, double divisor) {
        if (USE_JDK_MATH) {
            // no Math equivalent (differs from IEEEremainder(double,double))
        }
        if (Double.isInfinite(divisor)) {
            if (Double.isInfinite(dividend)) {
                return Double.NaN;
            } else {
                return dividend;
            }
        }
        double value = dividend % divisor;
        if (Math.abs(value+value) > Math.abs(divisor)) {
            return value + ((value > 0.0) ? -Math.abs(divisor) : Math.abs(divisor));
        } else {
            return value;
        }
    }

    /**
     * @param angle Angle in radians.
     * @return The same angle, in radians, but in [-Math.PI,Math.PI].
     */
    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) {
        if (USE_JDK_MATH) {
            return STRICT_MATH ? StrictMath.hypot(x,y) : Math.hypot(x,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 * FastMath.sqrt(x*x+y*y);
        }
    }

    /**
     * @param value A float value.
     * @return Ceiling of value.
     */
    public static float ceil(float value) {
        if (USE_JDK_MATH) {
            // TODO use Math.ceil(float) if exists
            return (float)Math.ceil((double)value);
        }
        return -FastMath.floor(-value);
    }

    /**
     * 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) {
        if (USE_JDK_MATH) {
            return Math.ceil(value);
        }
        return -FastMath.floor(-value);
    }

    /**
     * @param value A float value.
     * @return Floor of value.
     */
    public static float floor(float value) {
        if (USE_JDK_MATH) {
            // TODO use Math.floor(float) if exists
            return (float)Math.floor((double)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 anteCommaDigits = valueBits & (0xFF800000>>exp);
                if ((value < 0.0f) && (anteCommaDigits != valueBits)) {
                    return Float.intBitsToFloat(anteCommaDigits) - 1.0f;
                } else {
                    return Float.intBitsToFloat(anteCommaDigits);
                }
            } else {
                return 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) {
        if (USE_JDK_MATH) {
            return Math.floor(value);
        }
        // Faster than to work directly on bits.
        if (Math.abs(value) <= (double)Integer.MAX_VALUE) {
            if (value > 0.0) {
                return (double)(int)value;
            } else if (value < 0.0) {
                double anteCommaDigits = (double)(int)value;
                if (value != anteCommaDigits) {
                    return anteCommaDigits - 1.0;
                } else {
                    return anteCommaDigits;
                }
            } else { // value is +-0.0 (not NaN due to test against Integer.MAX_VALUE)
                return value;
            }
        } else if (Math.abs(value) < TWO_POW_52) {
            // We split the value in two:
            // high part, which is a mathematical integer,
            // and the rest, for which we can get rid of the
            // post comma digits by casting into an int.
            double highPart = ((int)(value * TWO_POW_N26)) * TWO_POW_26;
            if (value > 0.0) {
                return highPart + (double)((int)(value - highPart));
            } else {
                double anteCommaDigits = highPart + (double)((int)(value - highPart));
                if (value != anteCommaDigits) {
                    return anteCommaDigits - 1.0;
                } else {
                    return anteCommaDigits;
                }
            }
        } else { // abs(value) >= 2^52, or value is NaN
            return value;
        }
    }

    /**
     * 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) {
        if (USE_JDK_MATH) {
            return Math.round(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);
    }

    /**
     * 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) {
        if (USE_JDK_MATH) {
            return Math.round(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 (Math.abs(roundedValue) <= (double)Integer.MAX_VALUE) {
            // Faster with intermediary cast in int.
            return (long)(int)roundedValue;
        } else {
            return (long)roundedValue;
        }
    }

    /**
     * @param value A float value.
     * @return Value unbiased exponent.
     */
    public static int getExponent(float value) {
        if (USE_JDK_MATH) {
            return Math.getExponent(value);
        }
        return ((Float.floatToRawIntBits(value)>>23)&0xFF)-MAX_FLOAT_EXPONENT;
    }

    /**
     * @param value A double value.
     * @return Value unbiased exponent.
     */
    public static int getExponent(double value) {
        if (USE_JDK_MATH) {
            return Math.getExponent(value);
        }
        return (((int)(Double.doubleToRawLongBits(value)>>52))&0x7FF)-MAX_DOUBLE_EXPONENT;
    }

    /**
     * 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) {
        if (USE_JDK_MATH) {
            return Math.toDegrees(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) {
        if (USE_JDK_MATH) {
            return Math.toRadians(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) {
        if (USE_JDK_MATH) {
            return Math.abs(value);
        }
        return NumbersUtils.abs(value);
    }

    /**
     * @param value A long value.
     * @return The specified value as int.
     * @throws ArithmeticException if the specified value is not in [Integer.MIN_VALUE,Integer.MAX_VALUE] range.
     */
    public static int toIntExact(long value) {
        return NumbersUtils.asInt(value);
    }

    /**
     * @param value A long value.
     * @return The closest int value in [Integer.MIN_VALUE,Integer.MAX_VALUE] range.
     */
    public static int toInt(long value) {
        return NumbersUtils.toInt(value);
    }

    /**
     * @param a An int value.
     * @param b An int value.
     * @return The mathematical result of a+b.
     * @throws ArithmeticException if the mathematical result of a+b is not in [Integer.MIN_VALUE,Integer.MAX_VALUE] range.
     */
    public static int addExact(int a, int b) {
        return NumbersUtils.plusExact(a, b);
    }

    /**
     * @param a A long value.
     * @param b A long value.
     * @return The mathematical result of a+b.
     * @throws ArithmeticException if the mathematical result of a+b is not in [Long.MIN_VALUE,Long.MAX_VALUE] range.
     */
    public static long addExact(long a, long b) {
        return NumbersUtils.plusExact(a, b);
    }

    /**
     * @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 addBounded(int a, int b) {
        return NumbersUtils.plusBounded(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 addBounded(long a, long b) {
        return NumbersUtils.plusBounded(a, b);
    }

    /**
     * @param a An int value.
     * @param b An int value.
     * @return The mathematical result of a-b.
     * @throws ArithmeticException if the mathematical result of a-b is not in [Integer.MIN_VALUE,Integer.MAX_VALUE] range.
     */
    public static int subtractExact(int a, int b) {
        return NumbersUtils.minusExact(a, b);
    }

    /**
     * @param a A long value.
     * @param b A long value.
     * @return The mathematical result of a-b.
     * @throws ArithmeticException if the mathematical result of a-b is not in [Long.MIN_VALUE,Long.MAX_VALUE] range.
     */
    public static long subtractExact(long a, long b) {
        return NumbersUtils.minusExact(a, b);
    }

    /**
     * @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 subtractBounded(int a, int b) {
        return NumbersUtils.minusBounded(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 subtractBounded(long a, long b) {
        return NumbersUtils.minusBounded(a, b);
    }

    /**
     * @param a An int value.
     * @param b An int value.
     * @return The mathematical result of a*b.
     * @throws ArithmeticException if the mathematical result of a*b is not in [Integer.MIN_VALUE,Integer.MAX_VALUE] range.
     */
    public static int multiplyExact(int a, int b) {
        return NumbersUtils.timesExact(a, b);
    }

    /**
     * @param a A long value.
     * @param b A long value.
     * @return The mathematical result of a*b.
     * @throws ArithmeticException if the mathematical result of a*b is not in [Long.MIN_VALUE,Long.MAX_VALUE] range.
     */
    public static long multiplyExact(long a, long b) {
        return NumbersUtils.timesExact(a, b);
    }

    /**
     * @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 multiplyBounded(int a, int b) {
        return NumbersUtils.timesBounded(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 multiplyBounded(long a, long b) {
        return NumbersUtils.timesBounded(a, b);
    }

    /**
     * @param minValue An int value.
     * @param maxValue An int value.
     * @param value An int value.
     * @return minValue if value &lt; minValue, maxValue if value &gt; maxValue, value otherwise.
     */
    public static int toRange(int minValue, int maxValue, int value) {
        return NumbersUtils.toRange(minValue, maxValue, value);
    }

    /**
     * @param minValue A long value.
     * @param maxValue A long value.
     * @param value A long value.
     * @return minValue if value &lt; minValue, maxValue if value &gt; maxValue, value otherwise.
     */
    public static long toRange(long minValue, long maxValue, long value) {
        return NumbersUtils.toRange(minValue, maxValue, value);
    }

    /**
     * @param minValue A float value.
     * @param maxValue A float value.
     * @param value A float value.
     * @return minValue if value &lt; minValue, maxValue if value &gt; maxValue, value otherwise.
     */
    public static float toRange(float minValue, float maxValue, float value) {
        return NumbersUtils.toRange(minValue, maxValue, value);
    }

    /**
     * @param minValue A double value.
     * @param maxValue A double value.
     * @param value A double value.
     * @return minValue if value &lt; minValue, maxValue if value &gt; maxValue, value otherwise.
     */
    public static double toRange(double minValue, double maxValue, double value) {
        return NumbersUtils.toRange(minValue, maxValue, value);
    }

    /**
     * NB: Since 2*Math.PI &lt; 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 &gt; 2*Math.PI, like 2*PI_SUP constant of this class.
     * 
     * @param startAngRad An angle, in radians.
     * @param angSpanRad An angular span, &gt;= 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 (Math.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 &gt; 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(float value) {
        return NumbersUtils.isNaNOrInfinite(value);
    }

    public static boolean isNaNOrInfinite(double value) {
        return NumbersUtils.isNaNOrInfinite(value);
    }

    /*
     * 
     * 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 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() {
        // StrictMath and Math use different RNG instances,
        // so their random() methods are not equivalent.
        return STRICT_MATH ? StrictMath.random() : 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() {
    }

    /**
     * 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 USE_JDK_MATH ? Math.min(2, tabSizePower) : tabSizePower;
    }

    /**
     * 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) {
        if (USE_JDK_MATH) {
            double y = STRICT_MATH ? StrictMath.sin(angle) : Math.sin(angle);
            double x = STRICT_MATH ? StrictMath.cos(angle) : Math.cos(angle);
            return STRICT_MATH ? StrictMath.atan2(y,x) : Math.atan2(y,x);
        }
        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) {
        if (USE_JDK_MATH) {
            return remainderTwoPi(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) {
                // 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 {
            // 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.
     * 
     * Might use some FastMath methods in there, not to spend
     * an hour in it, but 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<SIN_COS_TABS_SIZE;i++) {
            // angle: in [0,2*PI].
            double angle = i * SIN_COS_DELTA_HI + i * SIN_COS_DELTA_LO;
            double sinAngle = StrictMath.sin(angle);
            double cosAngle = StrictMath.cos(angle);
            // For indexes corresponding to null cosine or sine, we make sure the value is zero
            // and not an epsilon. This allows for a much better accuracy for results close to zero.
            if (i == SIN_COS_PI_INDEX) {
                sinAngle = 0.0;
            } else if (i == SIN_COS_PI_MUL_2_INDEX) {
                sinAngle = 0.0;
            } else if (i == SIN_COS_PI_MUL_0_5_INDEX) {
                cosAngle = 0.0;
            } else if (i == SIN_COS_PI_MUL_1_5_INDEX) {
                cosAngle = 0.0;
            }
            sinTab[i] = sinAngle;
            cosTab[i] = cosAngle;
        }

        // tan

        for (int i=0;i<TAN_TABS_SIZE;i++) {
            // angle: in [0,TAN_MAX_VALUE_FOR_TABS].
            double angle = i * TAN_DELTA_HI + i * TAN_DELTA_LO;
            tanTab[i] = StrictMath.tan(angle);
            double cosAngle = StrictMath.cos(angle);
            double sinAngle = StrictMath.sin(angle);
            double cosAngleInv = 1/cosAngle;
            double cosAngleInv2 = cosAngleInv*cosAngleInv;
            double cosAngleInv3 = cosAngleInv2*cosAngleInv;
            double cosAngleInv4 = cosAngleInv2*cosAngleInv2;
            double cosAngleInv5 = cosAngleInv3*cosAngleInv2;
            tanDer1DivF1Tab[i] = cosAngleInv2;
            tanDer2DivF2Tab[i] = ((2*sinAngle)*cosAngleInv3) * ONE_DIV_F2;
            tanDer3DivF3Tab[i] = ((2*(1+2*sinAngle*sinAngle))*cosAngleInv4) * ONE_DIV_F3;
            tanDer4DivF4Tab[i] = ((8*sinAngle*(2+sinAngle*sinAngle))*cosAngleInv5) * ONE_DIV_F4;
        }

        // asin

        for (int i=0;i<ASIN_TABS_SIZE;i++) {
            // x: in [0,ASIN_MAX_VALUE_FOR_TABS].
            double x = i * ASIN_DELTA;
            asinTab[i] = StrictMath.asin(x);
            double oneMinusXSqInv = 1.0/(1-x*x);
            double oneMinusXSqInv0_5 = StrictMath.sqrt(oneMinusXSqInv);
            double oneMinusXSqInv1_5 = oneMinusXSqInv0_5*oneMinusXSqInv;
            double oneMinusXSqInv2_5 = oneMinusXSqInv1_5*oneMinusXSqInv;
            double oneMinusXSqInv3_5 = oneMinusXSqInv2_5*oneMinusXSqInv;
            asinDer1DivF1Tab[i] = oneMinusXSqInv0_5;
            asinDer2DivF2Tab[i] = (x*oneMinusXSqInv1_5) * ONE_DIV_F2;
            asinDer3DivF3Tab[i] = ((1+2*x*x)*oneMinusXSqInv2_5) * ONE_DIV_F3;
            asinDer4DivF4Tab[i] = ((5+2*x*(2+x*(5-2*x)))*oneMinusXSqInv3_5) * ONE_DIV_F4;
        }

        if (USE_POWTABS_FOR_ASIN) {
            for (int i=0;i<ASIN_POWTABS_SIZE;i++) {
                // x: in [0,ASIN_MAX_VALUE_FOR_POWTABS].
                double x = StrictMath.pow(i*(1.0/ASIN_POWTABS_SIZE_MINUS_ONE), 1.0/ASIN_POWTABS_POWER) * ASIN_MAX_VALUE_FOR_POWTABS;
                asinParamPowTab[i] = x;
                asinPowTab[i] = StrictMath.asin(x);
                double oneMinusXSqInv = 1.0/(1-x*x);
                double oneMinusXSqInv0_5 = StrictMath.sqrt(oneMinusXSqInv);
                double oneMinusXSqInv1_5 = oneMinusXSqInv0_5*oneMinusXSqInv;
                double oneMinusXSqInv2_5 = oneMinusXSqInv1_5*oneMinusXSqInv;
                double oneMinusXSqInv3_5 = oneMinusXSqInv2_5*oneMinusXSqInv;
                asinDer1DivF1PowTab[i] = oneMinusXSqInv0_5;
                asinDer2DivF2PowTab[i] = (x*oneMinusXSqInv1_5) * ONE_DIV_F2;
                asinDer3DivF3PowTab[i] = ((1+2*x*x)*oneMinusXSqInv2_5) * ONE_DIV_F3;
                asinDer4DivF4PowTab[i] = ((5+2*x*(2+x*(5-2*x)))*oneMinusXSqInv3_5) * ONE_DIV_F4;
            }
        }

        // atan

        for (int i=0;i<ATAN_TABS_SIZE;i++) {
            // x: in [0,ATAN_MAX_VALUE_FOR_TABS].
            double x = i * ATAN_DELTA;
            double onePlusXSqInv = 1.0/(1+x*x);
            double onePlusXSqInv2 = onePlusXSqInv*onePlusXSqInv;
            double onePlusXSqInv3 = onePlusXSqInv2*onePlusXSqInv;
            double onePlusXSqInv4 = onePlusXSqInv2*onePlusXSqInv2;
            atanTab[i] = StrictMath.atan(x);
            atanDer1DivF1Tab[i] = onePlusXSqInv;
            atanDer2DivF2Tab[i] = (-2*x*onePlusXSqInv2) * ONE_DIV_F2;
            atanDer3DivF3Tab[i] = ((-2+6*x*x)*onePlusXSqInv3) * ONE_DIV_F3;
            atanDer4DivF4Tab[i] = ((24*x*(1-x*x))*onePlusXSqInv4) * ONE_DIV_F4;
        }

        // exp

        for (int i=0;i<EXP_LO_TAB_SIZE;i++) {
            // x: in [-EXPM1_DISTANCE_TO_ZERO,EXPM1_DISTANCE_TO_ZERO].
            double x = -EXP_LO_DISTANCE_TO_ZERO + i/(double)EXP_LO_INDEXING;
            // exp(x)
            expLoPosTab[i] = StrictMath.exp(x);
            // 1-exp(-x), accurately computed
            expLoNegTab[i] = -StrictMath.expm1(-x);
        }
        for (int i=0;i<=(int)EXP_OVERFLOW_LIMIT;i++) {
            expHiTab[i] = StrictMath.exp(i);
        }
        for (int i=0;i<=-(int)EXP_UNDERFLOW_LIMIT;i++) {
            // We take care not to compute with subnormal values.
            if ((double)-i >= EXP_MIN_INT_LIMIT) {
                expHiInvTab[i] = StrictMath.exp(-i);
            } else {
                expHiInvTab[i] = StrictMath.exp(54*LOG_2-i);
            }
        }

        // log

        for (int i=0;i<LOG_TAB_SIZE;i++) {
            // Exact to use inverse of tab size, since it is a power of two.
            double x = 1+i*(1.0/LOG_TAB_SIZE);
            logXLogTab[i] = StrictMath.log(x);
            logXTab[i] = x;
            logXInvTab[i] = 1/x;
        }

        // twoPow

        for (int i=MIN_DOUBLE_EXPONENT;i<=MAX_DOUBLE_EXPONENT;i++) {
            twoPowTab[i-MIN_DOUBLE_EXPONENT] = StrictMath.pow(2.0,i);
        }

        // sqrt

        for (int i=MIN_DOUBLE_EXPONENT;i<=MAX_DOUBLE_EXPONENT;i++) {
            double twoPowExpDiv2 = StrictMath.pow(2.0,i*0.5);
            sqrtXSqrtHiTab[i-MIN_DOUBLE_EXPONENT] = twoPowExpDiv2 * 0.5; // Half sqrt, to avoid overflows.
            sqrtSlopeHiTab[i-MIN_DOUBLE_EXPONENT] = 1/twoPowExpDiv2;
        }
        sqrtXSqrtLoTab[0] = 1.0;
        sqrtSlopeLoTab[0] = 1.0;
        final long SQRT_LO_MASK = (0x3FF0000000000000L | (0x000FFFFFFFFFFFFFL>>SQRT_LO_BITS));
        for (int i=1;i<SQRT_LO_TAB_SIZE;i++) {
            long xBits = SQRT_LO_MASK | (((long)(i-1))<<(52-SQRT_LO_BITS));
            double sqrtX = StrictMath.sqrt(Double.longBitsToDouble(xBits));
            sqrtXSqrtLoTab[i] = sqrtX;
            sqrtSlopeLoTab[i] = 1/sqrtX;
        }

        // cbrt

        for (int i=MIN_DOUBLE_EXPONENT;i<=MAX_DOUBLE_EXPONENT;i++) {
            double twoPowExpDiv3 = StrictMath.pow(2.0,i/3.0);
            cbrtXCbrtHiTab[i-MIN_DOUBLE_EXPONENT] = twoPowExpDiv3 * 0.5; // Half cbrt, to avoid overflows.
            double tmp = 1/twoPowExpDiv3;
            cbrtSlopeHiTab[i-MIN_DOUBLE_EXPONENT] = (4.0/3)*tmp*tmp;
        }
        cbrtXCbrtLoTab[0] = 1.0;
        cbrtSlopeLoTab[0] = 1.0;
        final long CBRT_LO_MASK = (0x3FF0000000000000L | (0x000FFFFFFFFFFFFFL>>CBRT_LO_BITS));
        for (int i=1;i<CBRT_LO_TAB_SIZE;i++) {
            long xBits = CBRT_LO_MASK | (((long)(i-1))<<(52-CBRT_LO_BITS));
            double cbrtX = StrictMath.cbrt(Double.longBitsToDouble(xBits));
            cbrtXCbrtLoTab[i] = cbrtX;
            cbrtSlopeLoTab[i] = 1/(cbrtX*cbrtX);
        }
    }
}