    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math.util;
    018
    019    /**
    020     * Faster, more accurate, portable alternative to {@link StrictMath}.
    021     * <p>
    022     * Additionally implements the following methods not found in StrictMath:
    023     * <ul>
    024     * <li>{@link #asinh(double)}</li>
    025     * <li>{@link #acosh(double)}</li>
    026     * <li>{@link #atanh(double)}</li>
    027     * </ul>
    028     * The following methods are found in StrictMath since 1.6 only
    029     * <ul>
    030     * <li>{@link #copySign(double, double)}</li>
    031     * <li>{@link #getExponent(double)}</li>
    032     * <li>{@link #nextAfter(double,double)}</li>
    033     * <li>{@link #nextUp(double)}</li>
    034     * <li>{@link #scalb(double, int)}</li>
    035     * <li>{@link #copySign(float, float)}</li>
    036     * <li>{@link #getExponent(float)}</li>
    037     * <li>{@link #nextAfter(float,double)}</li>
    038     * <li>{@link #nextUp(float)}</li>
    039     * <li>{@link #scalb(float, int)}</li>
    040     * </ul>
    041     * @version $Revision: 1074294 $ $Date: 2011-02-24 22:18:59 +0100 (jeu. 24 f??vr. 2011) $
    042     * @since 2.2
    043     */
    044    public class FastMath {
    045
    046        /** Archimede's constant PI, ratio of circle circumference to diameter. */
    047        public static final double PI = 105414357.0 / 33554432.0 + 1.984187159361080883e-9;
    048
    049        /** Napier's constant e, base of the natural logarithm. */
    050        public static final double E = 2850325.0 / 1048576.0 + 8.254840070411028747e-8;
    051
    052        /** Exponential evaluated at integer values,
    053         * exp(x) =  expIntTableA[x + 750] + expIntTableB[x+750].
    054         */
    055        private static final double EXP_INT_TABLE_A[] = new double[1500];
    056
    057        /** Exponential evaluated at integer values,
    058         * exp(x) =  expIntTableA[x + 750] + expIntTableB[x+750]
    059         */
    060        private static final double EXP_INT_TABLE_B[] = new double[1500];
    061
    062        /** Exponential over the range of 0 - 1 in increments of 2^-10
    063         * exp(x/1024) =  expFracTableA[x] + expFracTableB[x].
    064         */
    065        private static final double EXP_FRAC_TABLE_A[] = new double[1025];
    066
    067        /** Exponential over the range of 0 - 1 in increments of 2^-10
    068         * exp(x/1024) =  expFracTableA[x] + expFracTableB[x].
    069         */
    070        private static final double EXP_FRAC_TABLE_B[] = new double[1025];
    071
    072        /** Factorial table, for Taylor series expansions. */
    073        private static final double FACT[] = new double[20];
    074
    075        /** Extended precision logarithm table over the range 1 - 2 in increments of 2^-10. */
    076        private static final double LN_MANT[][] = new double[1024][];
    077
    078        /** log(2) (high bits). */
    079        private static final double LN_2_A = 0.693147063255310059;
    080
    081        /** log(2) (low bits). */
    082        private static final double LN_2_B = 1.17304635250823482e-7;
    083
    084        /** Coefficients for slowLog. */
    085        private static final double LN_SPLIT_COEF[][] = {
    086            {2.0, 0.0},
    087            {0.6666666269302368, 3.9736429850260626E-8},
    088            {0.3999999761581421, 2.3841857910019882E-8},
    089            {0.2857142686843872, 1.7029898543501842E-8},
    090            {0.2222222089767456, 1.3245471311735498E-8},
    091            {0.1818181574344635, 2.4384203044354907E-8},
    092            {0.1538461446762085, 9.140260083262505E-9},
    093            {0.13333332538604736, 9.220590270857665E-9},
    094            {0.11764700710773468, 1.2393345855018391E-8},
    095            {0.10526403784751892, 8.251545029714408E-9},
    096            {0.0952233225107193, 1.2675934823758863E-8},
    097            {0.08713622391223907, 1.1430250008909141E-8},
    098            {0.07842259109020233, 2.404307984052299E-9},
    099            {0.08371849358081818, 1.176342548272881E-8},
    100            {0.030589580535888672, 1.2958646899018938E-9},
    101            {0.14982303977012634, 1.225743062930824E-8},
    102        };
    103
    104        /** Coefficients for log, when input 0.99 < x < 1.01. */
    105        private static final double LN_QUICK_COEF[][] = {
    106            {1.0, 5.669184079525E-24},
    107            {-0.25, -0.25},
    108            {0.3333333134651184, 1.986821492305628E-8},
    109            {-0.25, -6.663542893624021E-14},
    110            {0.19999998807907104, 1.1921056801463227E-8},
    111            {-0.1666666567325592, -7.800414592973399E-9},
    112            {0.1428571343421936, 5.650007086920087E-9},
    113            {-0.12502530217170715, -7.44321345601866E-11},
    114            {0.11113807559013367, 9.219544613762692E-9},
    115        };
    116
    117        /** Coefficients for log in the range of 1.0 < x < 1.0 + 2^-10. */
    118        private static final double LN_HI_PREC_COEF[][] = {
    119            {1.0, -6.032174644509064E-23},
    120            {-0.25, -0.25},
    121            {0.3333333134651184, 1.9868161777724352E-8},
    122            {-0.2499999701976776, -2.957007209750105E-8},
    123            {0.19999954104423523, 1.5830993332061267E-10},
    124            {-0.16624879837036133, -2.6033824355191673E-8}
    125        };
    126
    127        /** Sine table (high bits). */
    128        private static final double SINE_TABLE_A[] = new double[14];
    129
    130        /** Sine table (low bits). */
    131        private static final double SINE_TABLE_B[] = new double[14];
    132
    133        /** Cosine table (high bits). */
    134        private static final double COSINE_TABLE_A[] = new double[14];
    135
    136        /** Cosine table (low bits). */
    137        private static final double COSINE_TABLE_B[] = new double[14];
    138
    139        /** Tangent table, used by atan() (high bits). */
    140        private static final double TANGENT_TABLE_A[] = new double[14];
    141
    142        /** Tangent table, used by atan() (low bits). */
    143        private static final double TANGENT_TABLE_B[] = new double[14];
    144
    145        /** Bits of 1/(2*pi), need for reducePayneHanek(). */
    146        private static final long RECIP_2PI[] = new long[] {
    147            (0x28be60dbL << 32) | 0x9391054aL,
    148            (0x7f09d5f4L << 32) | 0x7d4d3770L,
    149            (0x36d8a566L << 32) | 0x4f10e410L,
    150            (0x7f9458eaL << 32) | 0xf7aef158L,
    151            (0x6dc91b8eL << 32) | 0x909374b8L,
    152            (0x01924bbaL << 32) | 0x82746487L,
    153            (0x3f877ac7L << 32) | 0x2c4a69cfL,
    154            (0xba208d7dL << 32) | 0x4baed121L,
    155            (0x3a671c09L << 32) | 0xad17df90L,
    156            (0x4e64758eL << 32) | 0x60d4ce7dL,
    157            (0x272117e2L << 32) | 0xef7e4a0eL,
    158            (0xc7fe25ffL << 32) | 0xf7816603L,
    159            (0xfbcbc462L << 32) | 0xd6829b47L,
    160            (0xdb4d9fb3L << 32) | 0xc9f2c26dL,
    161            (0xd3d18fd9L << 32) | 0xa797fa8bL,
    162            (0x5d49eeb1L << 32) | 0xfaf97c5eL,
    163            (0xcf41ce7dL << 32) | 0xe294a4baL,
    164             0x9afed7ecL << 32  };
    165
    166        /** Bits of pi/4, need for reducePayneHanek(). */
    167        private static final long PI_O_4_BITS[] = new long[] {
    168            (0xc90fdaa2L << 32) | 0x2168c234L,
    169            (0xc4c6628bL << 32) | 0x80dc1cd1L };
    170
    171        /** Eighths.
    172         * This is used by sinQ, because its faster to do a table lookup than
    173         * a multiply in this time-critical routine
    174         */
    175        private static final double EIGHTHS[] = {0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625};
    176
    177        /** Table of 2^((n+2)/3) */
    178        private static final double CBRTTWO[] = { 0.6299605249474366,
    179                                                0.7937005259840998,
    180                                                1.0,
    181                                                1.2599210498948732,
    182                                                1.5874010519681994 };
    183
    184        /*
    185         *  There are 52 bits in the mantissa of a double.
    186         *  For additional precision, the code splits double numbers into two parts,
    187         *  by clearing the low order 30 bits if possible, and then performs the arithmetic
    188         *  on each half separately.
    189         */
    190
    191        /**
    192         * 0x40000000 - used to split a double into two parts, both with the low order bits cleared.
    193         * Equivalent to 2^30.
    194         */
    195        private static final long HEX_40000000 = 0x40000000L; // 1073741824L
    196
    197        /** Mask used to clear low order 30 bits */
    198        private static final long MASK_30BITS = -1L - (HEX_40000000 -1); // 0xFFFFFFFFC0000000L;
    199
    200        /** 2^52 - double numbers this large must be integral (no fraction) or NaN or Infinite */
    201        private static final double TWO_POWER_52 = 4503599627370496.0;
    202
    203        // Initialize tables
    204        static {
    205            int i;
    206
    207            // Generate an array of factorials
    208            FACT[0] = 1.0;
    209            for (i = 1; i < FACT.length; i++) {
    210                FACT[i] = FACT[i-1] * i;
    211            }
    212
    213            double tmp[] = new double[2];
    214            double recip[] = new double[2];
    215
    216            // Populate expIntTable
    217            for (i = 0; i < 750; i++) {
    218                expint(i, tmp);
    219                EXP_INT_TABLE_A[i+750] = tmp[0];
    220                EXP_INT_TABLE_B[i+750] = tmp[1];
    221
    222                if (i != 0) {
    223                    // Negative integer powers
    224                    splitReciprocal(tmp, recip);
    225                    EXP_INT_TABLE_A[750-i] = recip[0];
    226                    EXP_INT_TABLE_B[750-i] = recip[1];
    227                }
    228            }
    229
    230            // Populate expFracTable
    231            for (i = 0; i < EXP_FRAC_TABLE_A.length; i++) {
    232                slowexp(i/1024.0, tmp);
    233                EXP_FRAC_TABLE_A[i] = tmp[0];
    234                EXP_FRAC_TABLE_B[i] = tmp[1];
    235            }
    236
    237            // Populate lnMant table
    238            for (i = 0; i < LN_MANT.length; i++) {
    239                double d = Double.longBitsToDouble( (((long) i) << 42) | 0x3ff0000000000000L );
    240                LN_MANT[i] = slowLog(d);
    241            }
    242
    243            // Build the sine and cosine tables
    244            buildSinCosTables();
    245        }
    246
    247        /**
    248         * Private Constructor
    249         */
    250        private FastMath() {
    251        }
    252
    253        // Generic helper methods
    254
    255        /**
    256         * Get the high order bits from the mantissa.
    257         * Equivalent to adding and subtracting HEX_40000 but also works for very large numbers
    258         *
    259         * @param d the value to split
    260         * @return the high order part of the mantissa
    261         */
    262        private static double doubleHighPart(double d) {
    263            if (d > -MathUtils.SAFE_MIN && d < MathUtils.SAFE_MIN){
    264                return d; // These are un-normalised - don't try to convert
    265            }
    266            long xl = Double.doubleToLongBits(d);
    267            xl = xl & MASK_30BITS; // Drop low order bits
    268            return Double.longBitsToDouble(xl);
    269        }
    270
    271        /** Compute the square root of a number.
    272         * <p><b>Note:</b> this implementation currently delegates to {@link Math#sqrt}
    273         * @param a number on which evaluation is done
    274         * @return square root of a
    275         */
    276        public static double sqrt(final double a) {
    277            return Math.sqrt(a);
    278        }
    279
    280        /** Compute the hyperbolic cosine of a number.
    281         * @param x number on which evaluation is done
    282         * @return hyperbolic cosine of x
    283         */
    284        public static double cosh(double x) {
    285          if (x != x) {
    286              return x;
    287          }
    288
    289          if (x > 20.0) {
    290              return exp(x)/2.0;
    291          }
    292
    293          if (x < -20) {
    294              return exp(-x)/2.0;
    295          }
    296
    297          double hiPrec[] = new double[2];
    298          if (x < 0.0) {
    299              x = -x;
    300          }
    301          exp(x, 0.0, hiPrec);
    302
    303          double ya = hiPrec[0] + hiPrec[1];
    304          double yb = -(ya - hiPrec[0] - hiPrec[1]);
    305
    306          double temp = ya * HEX_40000000;
    307          double yaa = ya + temp - temp;
    308          double yab = ya - yaa;
    309
    310          // recip = 1/y
    311          double recip = 1.0/ya;
    312          temp = recip * HEX_40000000;
    313          double recipa = recip + temp - temp;
    314          double recipb = recip - recipa;
    315
    316          // Correct for rounding in division
    317          recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
    318          // Account for yb
    319          recipb += -yb * recip * recip;
    320
    321          // y = y + 1/y
    322          temp = ya + recipa;
    323          yb += -(temp - ya - recipa);
    324          ya = temp;
    325          temp = ya + recipb;
    326          yb += -(temp - ya - recipb);
    327          ya = temp;
    328
    329          double result = ya + yb;
    330          result *= 0.5;
    331          return result;
    332        }
    333
    334        /** Compute the hyperbolic sine of a number.
    335         * @param x number on which evaluation is done
    336         * @return hyperbolic sine of x
    337         */
    338        public static double sinh(double x) {
    339          boolean negate = false;
    340          if (x != x) {
    341              return x;
    342          }
    343
    344          if (x > 20.0) {
    345              return exp(x)/2.0;
    346          }
    347
    348          if (x < -20) {
    349              return -exp(-x)/2.0;
    350          }
    351
    352          if (x == 0) {
    353              return x;
    354          }
    355
    356          if (x < 0.0) {
    357              x = -x;
    358              negate = true;
    359          }
    360
    361          double result;
    362
    363          if (x > 0.25) {
    364              double hiPrec[] = new double[2];
    365              exp(x, 0.0, hiPrec);
    366
    367              double ya = hiPrec[0] + hiPrec[1];
    368              double yb = -(ya - hiPrec[0] - hiPrec[1]);
    369
    370              double temp = ya * HEX_40000000;
    371              double yaa = ya + temp - temp;
    372              double yab = ya - yaa;
    373
    374              // recip = 1/y
    375              double recip = 1.0/ya;
    376              temp = recip * HEX_40000000;
    377              double recipa = recip + temp - temp;
    378              double recipb = recip - recipa;
    379
    380              // Correct for rounding in division
    381              recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
    382              // Account for yb
    383              recipb += -yb * recip * recip;
    384
    385              recipa = -recipa;
    386              recipb = -recipb;
    387
    388              // y = y + 1/y
    389              temp = ya + recipa;
    390              yb += -(temp - ya - recipa);
    391              ya = temp;
    392              temp = ya + recipb;
    393              yb += -(temp - ya - recipb);
    394              ya = temp;
    395
    396              result = ya + yb;
    397              result *= 0.5;
    398          }
    399          else {
    400              double hiPrec[] = new double[2];
    401              expm1(x, hiPrec);
    402
    403              double ya = hiPrec[0] + hiPrec[1];
    404              double yb = -(ya - hiPrec[0] - hiPrec[1]);
    405
    406              /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
    407              double denom = 1.0 + ya;
    408              double denomr = 1.0 / denom;
    409              double denomb = -(denom - 1.0 - ya) + yb;
    410              double ratio = ya * denomr;
    411              double temp = ratio * HEX_40000000;
    412              double ra = ratio + temp - temp;
    413              double rb = ratio - ra;
    414
    415              temp = denom * HEX_40000000;
    416              double za = denom + temp - temp;
    417              double zb = denom - za;
    418
    419              rb += (ya - za*ra - za*rb - zb*ra - zb*rb) * denomr;
    420
    421              // Adjust for yb
    422              rb += yb*denomr;                        // numerator
    423              rb += -ya * denomb * denomr * denomr;   // denominator
    424
    425              // y = y - 1/y
    426              temp = ya + ra;
    427              yb += -(temp - ya - ra);
    428              ya = temp;
    429              temp = ya + rb;
    430              yb += -(temp - ya - rb);
    431              ya = temp;
    432
    433              result = ya + yb;
    434              result *= 0.5;
    435          }
    436
    437          if (negate) {
    438              result = -result;
    439          }
    440
    441          return result;
    442        }
    443
    444        /** Compute the hyperbolic tangent of a number.
    445         * @param x number on which evaluation is done
    446         * @return hyperbolic tangent of x
    447         */
    448        public static double tanh(double x) {
    449          boolean negate = false;
    450
    451          if (x != x) {
    452              return x;
    453          }
    454
    455          if (x > 20.0) {
    456              return 1.0;
    457          }
    458
    459          if (x < -20) {
    460              return -1.0;
    461          }
    462
    463          if (x == 0) {
    464              return x;
    465          }
    466
    467          if (x < 0.0) {
    468              x = -x;
    469              negate = true;
    470          }
    471
    472          double result;
    473          if (x >= 0.5) {
    474              double hiPrec[] = new double[2];
    475              // tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
    476              exp(x*2.0, 0.0, hiPrec);
    477
    478              double ya = hiPrec[0] + hiPrec[1];
    479              double yb = -(ya - hiPrec[0] - hiPrec[1]);
    480
    481              /* Numerator */
    482              double na = -1.0 + ya;
    483              double nb = -(na + 1.0 - ya);
    484              double temp = na + yb;
    485              nb += -(temp - na - yb);
    486              na = temp;
    487
    488              /* Denominator */
    489              double da = 1.0 + ya;
    490              double db = -(da - 1.0 - ya);
    491              temp = da + yb;
    492              db += -(temp - da - yb);
    493              da = temp;
    494
    495              temp = da * HEX_40000000;
    496              double daa = da + temp - temp;
    497              double dab = da - daa;
    498
    499              // ratio = na/da
    500              double ratio = na/da;
    501              temp = ratio * HEX_40000000;
    502              double ratioa = ratio + temp - temp;
    503              double ratiob = ratio - ratioa;
    504
    505              // Correct for rounding in division
    506              ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;
    507
    508              // Account for nb
    509              ratiob += nb / da;
    510              // Account for db
    511              ratiob += -db * na / da / da;
    512
    513              result = ratioa + ratiob;
    514          }
    515          else {
    516              double hiPrec[] = new double[2];
    517              // tanh(x) = expm1(2x) / (expm1(2x) + 2)
    518              expm1(x*2.0, hiPrec);
    519
    520              double ya = hiPrec[0] + hiPrec[1];
    521              double yb = -(ya - hiPrec[0] - hiPrec[1]);
    522
    523              /* Numerator */
    524              double na = ya;
    525              double nb = yb;
    526
    527              /* Denominator */
    528              double da = 2.0 + ya;
    529              double db = -(da - 2.0 - ya);
    530              double temp = da + yb;
    531              db += -(temp - da - yb);
    532              da = temp;
    533
    534              temp = da * HEX_40000000;
    535              double daa = da + temp - temp;
    536              double dab = da - daa;
    537
    538              // ratio = na/da
    539              double ratio = na/da;
    540              temp = ratio * HEX_40000000;
    541              double ratioa = ratio + temp - temp;
    542              double ratiob = ratio - ratioa;
    543
    544              // Correct for rounding in division
    545              ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;
    546
    547              // Account for nb
    548              ratiob += nb / da;
    549              // Account for db
    550              ratiob += -db * na / da / da;
    551
    552              result = ratioa + ratiob;
    553          }
    554
    555          if (negate) {
    556              result = -result;
    557          }
    558
    559          return result;
    560        }
    561
    562        /** Compute the inverse hyperbolic cosine of a number.
    563         * @param a number on which evaluation is done
    564         * @return inverse hyperbolic cosine of a
    565         */
    566        public static double acosh(final double a) {
    567            return FastMath.log(a + FastMath.sqrt(a * a - 1));
    568        }
    569
    570        /** Compute the inverse hyperbolic sine of a number.
    571         * @param a number on which evaluation is done
    572         * @return inverse hyperbolic sine of a
    573         */
    574        public static double asinh(double a) {
    575
    576            boolean negative = false;
    577            if (a < 0) {
    578                negative = true;
    579                a = -a;
    580            }
    581
    582            double absAsinh;
    583            if (a > 0.167) {
    584                absAsinh = FastMath.log(FastMath.sqrt(a * a + 1) + a);
    585            } else {
    586                final double a2 = a * a;
    587                if (a > 0.097) {
    588                    absAsinh = a * (1 - a2 * (1 / 3.0 - a2 * (1 / 5.0 - a2 * (1 / 7.0 - a2 * (1 / 9.0 - a2 * (1.0 / 11.0 - a2 * (1.0 / 13.0 - a2 * (1.0 / 15.0 - a2 * (1.0 / 17.0) * 15.0 / 16.0) * 13.0 / 14.0) * 11.0 / 12.0) * 9.0 / 10.0) * 7.0 / 8.0) * 5.0 / 6.0) * 3.0 / 4.0) / 2.0);
    589                } else if (a > 0.036) {
    590                    absAsinh = a * (1 - a2 * (1 / 3.0 - a2 * (1 / 5.0 - a2 * (1 / 7.0 - a2 * (1 / 9.0 - a2 * (1.0 / 11.0 - a2 * (1.0 / 13.0) * 11.0 / 12.0) * 9.0 / 10.0) * 7.0 / 8.0) * 5.0 / 6.0) * 3.0 / 4.0) / 2.0);
    591                } else if (a > 0.0036) {
    592                    absAsinh = a * (1 - a2 * (1 / 3.0 - a2 * (1 / 5.0 - a2 * (1 / 7.0 - a2 * (1 / 9.0) * 7.0 / 8.0) * 5.0 / 6.0) * 3.0 / 4.0) / 2.0);
    593                } else {
    594                    absAsinh = a * (1 - a2 * (1 / 3.0 - a2 * (1 / 5.0) * 3.0 / 4.0) / 2.0);
    595                }
    596            }
    597
    598            return negative ? -absAsinh : absAsinh;
    599
    600        }
    601
    602        /** Compute the inverse hyperbolic tangent of a number.
    603         * @param a number on which evaluation is done
    604         * @return inverse hyperbolic tangent of a
    605         */
    606        public static double atanh(double a) {
    607
    608            boolean negative = false;
    609            if (a < 0) {
    610                negative = true;
    611                a = -a;
    612            }
    613
    614            double absAtanh;
    615            if (a > 0.15) {
    616                absAtanh = 0.5 * FastMath.log((1 + a) / (1 - a));
    617            } else {
    618                final double a2 = a * a;
    619                if (a > 0.087) {
    620                    absAtanh = a * (1 + a2 * (1.0 / 3.0 + a2 * (1.0 / 5.0 + a2 * (1.0 / 7.0 + a2 * (1.0 / 9.0 + a2 * (1.0 / 11.0 + a2 * (1.0 / 13.0 + a2 * (1.0 / 15.0 + a2 * (1.0 / 17.0)))))))));
    621                } else if (a > 0.031) {
    622                    absAtanh = a * (1 + a2 * (1.0 / 3.0 + a2 * (1.0 / 5.0 + a2 * (1.0 / 7.0 + a2 * (1.0 / 9.0 + a2 * (1.0 / 11.0 + a2 * (1.0 / 13.0)))))));
    623                } else if (a > 0.003) {
    624                    absAtanh = a * (1 + a2 * (1.0 / 3.0 + a2 * (1.0 / 5.0 + a2 * (1.0 / 7.0 + a2 * (1.0 / 9.0)))));
    625                } else {
    626                    absAtanh = a * (1 + a2 * (1.0 / 3.0 + a2 * (1.0 / 5.0)));
    627                }
    628            }
    629
    630            return negative ? -absAtanh : absAtanh;
    631
    632        }
    633
    634        /** Compute the signum of a number.
    635         * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise
    636         * @param a number on which evaluation is done
    637         * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
    638         */
    639        public static double signum(final double a) {
    640            return (a < 0.0) ? -1.0 : ((a > 0.0) ? 1.0 : a); // return +0.0/-0.0/NaN depending on a
    641        }
    642
    643        /** Compute the signum of a number.
    644         * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise
    645         * @param a number on which evaluation is done
    646         * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
    647         */
    648        public static float signum(final float a) {
    649            return (a < 0.0f) ? -1.0f : ((a > 0.0f) ? 1.0f : a); // return +0.0/-0.0/NaN depending on a
    650        }
    651
    652        /** Compute next number towards positive infinity.
    653         * @param a number to which neighbor should be computed
    654         * @return neighbor of a towards positive infinity
    655         */
    656        public static double nextUp(final double a) {
    657            return nextAfter(a, Double.POSITIVE_INFINITY);
    658        }
    659
    660        /** Compute next number towards positive infinity.
    661         * @param a number to which neighbor should be computed
    662         * @return neighbor of a towards positive infinity
    663         */
    664        public static float nextUp(final float a) {
    665            return nextAfter(a, Float.POSITIVE_INFINITY);
    666        }
    667
    668        /** Returns a pseudo-random number between 0.0 and 1.0.
    669         * <p><b>Note:</b> this implementation currently delegates to {@link Math#random}
    670         * @return a random number between 0.0 and 1.0
    671         */
    672        public static double random() {
    673            return Math.random();
    674        }
    675
    676        /**
    677         * Exponential function.
    678         *
    679         * Computes exp(x), function result is nearly rounded.   It will be correctly
    680         * rounded to the theoretical value for 99.9% of input values, otherwise it will
    681         * have a 1 UPL error.
    682         *
    683         * Method:
    684         *    Lookup intVal = exp(int(x))
    685         *    Lookup fracVal = exp(int(x-int(x) / 1024.0) * 1024.0 );
    686         *    Compute z as the exponential of the remaining bits by a polynomial minus one
    687         *    exp(x) = intVal * fracVal * (1 + z)
    688         *
    689         * Accuracy:
    690         *    Calculation is done with 63 bits of precision, so result should be correctly
    691         *    rounded for 99.9% of input values, with less than 1 ULP error otherwise.
    692         *
    693         * @param x   a double
    694         * @return double e<sup>x</sup>
    695         */
    696        public static double exp(double x) {
    697            return exp(x, 0.0, null);
    698        }
    699
    700        /**
    701         * Internal helper method for exponential function.
    702         * @param x original argument of the exponential function
    703         * @param extra extra bits of precision on input (To Be Confirmed)
    704         * @param hiPrec extra bits of precision on output (To Be Confirmed)
    705         * @return exp(x)
    706         */
    707        private static double exp(double x, double extra, double[] hiPrec) {
    708            double intPartA;
    709            double intPartB;
    710            int intVal;
    711
    712            /* Lookup exp(floor(x)).
    713             * intPartA will have the upper 22 bits, intPartB will have the lower
    714             * 52 bits.
    715             */
    716            if (x < 0.0) {
    717                intVal = (int) -x;
    718
    719                if (intVal > 746) {
    720                    if (hiPrec != null) {
    721                        hiPrec[0] = 0.0;
    722                        hiPrec[1] = 0.0;
    723                    }
    724                    return 0.0;
    725                }
    726
    727                if (intVal > 709) {
    728                    /* This will produce a subnormal output */
    729                    final double result = exp(x+40.19140625, extra, hiPrec) / 285040095144011776.0;
    730                    if (hiPrec != null) {
    731                        hiPrec[0] /= 285040095144011776.0;
    732                        hiPrec[1] /= 285040095144011776.0;
    733                    }
    734                    return result;
    735                }
    736
    737                if (intVal == 709) {
    738                    /* exp(1.494140625) is nearly a machine number... */
    739                    final double result = exp(x+1.494140625, extra, hiPrec) / 4.455505956692756620;
    740                    if (hiPrec != null) {
    741                        hiPrec[0] /= 4.455505956692756620;
    742                        hiPrec[1] /= 4.455505956692756620;
    743                    }
    744                    return result;
    745                }
    746
    747                intVal++;
    748
    749                intPartA = EXP_INT_TABLE_A[750-intVal];
    750                intPartB = EXP_INT_TABLE_B[750-intVal];
    751
    752                intVal = -intVal;
    753            } else {
    754                intVal = (int) x;
    755
    756                if (intVal > 709) {
    757                    if (hiPrec != null) {
    758                        hiPrec[0] = Double.POSITIVE_INFINITY;
    759                        hiPrec[1] = 0.0;
    760                    }
    761                    return Double.POSITIVE_INFINITY;
    762                }
    763
    764                intPartA = EXP_INT_TABLE_A[750+intVal];
    765                intPartB = EXP_INT_TABLE_B[750+intVal];
    766            }
    767
    768            /* Get the fractional part of x, find the greatest multiple of 2^-10 less than
    769             * x and look up the exp function of it.
    770             * fracPartA will have the upper 22 bits, fracPartB the lower 52 bits.
    771             */
    772            final int intFrac = (int) ((x - intVal) * 1024.0);
    773            final double fracPartA = EXP_FRAC_TABLE_A[intFrac];
    774            final double fracPartB = EXP_FRAC_TABLE_B[intFrac];
    775
    776            /* epsilon is the difference in x from the nearest multiple of 2^-10.  It
    777             * has a value in the range 0 <= epsilon < 2^-10.
    778             * Do the subtraction from x as the last step to avoid possible loss of percison.
    779             */
    780            final double epsilon = x - (intVal + intFrac / 1024.0);
    781
    782            /* Compute z = exp(epsilon) - 1.0 via a minimax polynomial.  z has
    783           full double precision (52 bits).  Since z < 2^-10, we will have
    784           62 bits of precision when combined with the contant 1.  This will be
    785           used in the last addition below to get proper rounding. */
    786
    787            /* Remez generated polynomial.  Converges on the interval [0, 2^-10], error
    788           is less than 0.5 ULP */
    789            double z = 0.04168701738764507;
    790            z = z * epsilon + 0.1666666505023083;
    791            z = z * epsilon + 0.5000000000042687;
    792            z = z * epsilon + 1.0;
    793            z = z * epsilon + -3.940510424527919E-20;
    794
    795            /* Compute (intPartA+intPartB) * (fracPartA+fracPartB) by binomial
    796           expansion.
    797           tempA is exact since intPartA and intPartB only have 22 bits each.
    798           tempB will have 52 bits of precision.
    799             */
    800            double tempA = intPartA * fracPartA;
    801            double tempB = intPartA * fracPartB + intPartB * fracPartA + intPartB * fracPartB;
    802
    803            /* Compute the result.  (1+z)(tempA+tempB).  Order of operations is
    804           important.  For accuracy add by increasing size.  tempA is exact and
    805           much larger than the others.  If there are extra bits specified from the
    806           pow() function, use them. */
    807            final double tempC = tempB + tempA;
    808            final double result;
    809            if (extra != 0.0) {
    810                result = tempC*extra*z + tempC*extra + tempC*z + tempB + tempA;
    811            } else {
    812                result = tempC*z + tempB + tempA;
    813            }
    814
    815            if (hiPrec != null) {
    816                // If requesting high precision
    817                hiPrec[0] = tempA;
    818                hiPrec[1] = tempC*extra*z + tempC*extra + tempC*z + tempB;
    819            }
    820
    821            return result;
    822        }
    823
    824        /** Compute exp(x) - 1
    825         * @param x number to compute shifted exponential
    826         * @return exp(x) - 1
    827         */
    828        public static double expm1(double x) {
    829          return expm1(x, null);
    830        }
    831
    832        /** Internal helper method for expm1
    833         * @param x number to compute shifted exponential
    834         * @param hiPrecOut receive high precision result for -1.0 < x < 1.0
    835         * @return exp(x) - 1
    836         */
    837        private static double expm1(double x, double hiPrecOut[]) {
    838            if (x != x || x == 0.0) { // NaN or zero
    839                return x;
    840            }
    841
    842            if (x <= -1.0 || x >= 1.0) {
    843                // If not between +/- 1.0
    844                //return exp(x) - 1.0;
    845                double hiPrec[] = new double[2];
    846                exp(x, 0.0, hiPrec);
    847                if (x > 0.0) {
    848                    return -1.0 + hiPrec[0] + hiPrec[1];
    849                } else {
    850                    final double ra = -1.0 + hiPrec[0];
    851                    double rb = -(ra + 1.0 - hiPrec[0]);
    852                    rb += hiPrec[1];
    853                    return ra + rb;
    854                }
    855            }
    856
    857            double baseA;
    858            double baseB;
    859            double epsilon;
    860            boolean negative = false;
    861
    862            if (x < 0.0) {
    863                x = -x;
    864                negative = true;
    865            }
    866
    867            {
    868                int intFrac = (int) (x * 1024.0);
    869                double tempA = EXP_FRAC_TABLE_A[intFrac] - 1.0;
    870                double tempB = EXP_FRAC_TABLE_B[intFrac];
    871
    872                double temp = tempA + tempB;
    873                tempB = -(temp - tempA - tempB);
    874                tempA = temp;
    875
    876                temp = tempA * HEX_40000000;
    877                baseA = tempA + temp - temp;
    878                baseB = tempB + (tempA - baseA);
    879
    880                epsilon = x - intFrac/1024.0;
    881            }
    882
    883
    884            /* Compute expm1(epsilon) */
    885            double zb = 0.008336750013465571;
    886            zb = zb * epsilon + 0.041666663879186654;
    887            zb = zb * epsilon + 0.16666666666745392;
    888            zb = zb * epsilon + 0.49999999999999994;
    889            zb = zb * epsilon;
    890            zb = zb * epsilon;
    891
    892            double za = epsilon;
    893            double temp = za + zb;
    894            zb = -(temp - za - zb);
    895            za = temp;
    896
    897            temp = za * HEX_40000000;
    898            temp = za + temp - temp;
    899            zb += za - temp;
    900            za = temp;
    901
    902            /* Combine the parts.   expm1(a+b) = expm1(a) + expm1(b) + expm1(a)*expm1(b) */
    903            double ya = za * baseA;
    904            //double yb = za*baseB + zb*baseA + zb*baseB;
    905            temp = ya + za * baseB;
    906            double yb = -(temp - ya - za * baseB);
    907            ya = temp;
    908
    909            temp = ya + zb * baseA;
    910            yb += -(temp - ya - zb * baseA);
    911            ya = temp;
    912
    913            temp = ya + zb * baseB;
    914            yb += -(temp - ya - zb*baseB);
    915            ya = temp;
    916
    917            //ya = ya + za + baseA;
    918            //yb = yb + zb + baseB;
    919            temp = ya + baseA;
    920            yb += -(temp - baseA - ya);
    921            ya = temp;
    922
    923            temp = ya + za;
    924            //yb += (ya > za) ? -(temp - ya - za) : -(temp - za - ya);
    925            yb += -(temp - ya - za);
    926            ya = temp;
    927
    928            temp = ya + baseB;
    929            //yb += (ya > baseB) ? -(temp - ya - baseB) : -(temp - baseB - ya);
    930            yb += -(temp - ya - baseB);
    931            ya = temp;
    932
    933            temp = ya + zb;
    934            //yb += (ya > zb) ? -(temp - ya - zb) : -(temp - zb - ya);
    935            yb += -(temp - ya - zb);
    936            ya = temp;
    937
    938            if (negative) {
    939                /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
    940                double denom = 1.0 + ya;
    941                double denomr = 1.0 / denom;
    942                double denomb = -(denom - 1.0 - ya) + yb;
    943                double ratio = ya * denomr;
    944                temp = ratio * HEX_40000000;
    945                final double ra = ratio + temp - temp;
    946                double rb = ratio - ra;
    947
    948                temp = denom * HEX_40000000;
    949                za = denom + temp - temp;
    950                zb = denom - za;
    951
    952                rb += (ya - za * ra - za * rb - zb * ra - zb * rb) * denomr;
    953
    954                // f(x) = x/1+x
    955                // Compute f'(x)
    956                // Product rule:  d(uv) = du*v + u*dv
    957                // Chain rule:  d(f(g(x)) = f'(g(x))*f(g'(x))
    958                // d(1/x) = -1/(x*x)
    959                // d(1/1+x) = -1/( (1+x)^2) *  1 =  -1/((1+x)*(1+x))
    960                // d(x/1+x) = -x/((1+x)(1+x)) + 1/1+x = 1 / ((1+x)(1+x))
    961
    962                // Adjust for yb
    963                rb += yb * denomr;                      // numerator
    964                rb += -ya * denomb * denomr * denomr;   // denominator
    965
    966                // negate
    967                ya = -ra;
    968                yb = -rb;
    969            }
    970
    971            if (hiPrecOut != null) {
    972                hiPrecOut[0] = ya;
    973                hiPrecOut[1] = yb;
    974            }
    975
    976            return ya + yb;
    977        }
    978
    979        /**
    980         *  For x between 0 and 1, returns exp(x), uses extended precision
    981         *  @param x argument of exponential
    982         *  @param result placeholder where to place exp(x) split in two terms
    983         *  for extra precision (i.e. exp(x) = result[0] ?? result[1]
    984         *  @return exp(x)
    985         */
    986        private static double slowexp(final double x, final double result[]) {
    987            final double xs[] = new double[2];
    988            final double ys[] = new double[2];
    989            final double facts[] = new double[2];
    990            final double as[] = new double[2];
    991            split(x, xs);
    992            ys[0] = ys[1] = 0.0;
    993
    994            for (int i = 19; i >= 0; i--) {
    995                splitMult(xs, ys, as);
    996                ys[0] = as[0];
    997                ys[1] = as[1];
    998
    999                split(FACT[i], as);
    1000                splitReciprocal(as, facts);
    1001
    1002                splitAdd(ys, facts, as);
    1003                ys[0] = as[0];
    1004                ys[1] = as[1];
    1005            }
    1006
    1007            if (result != null) {
    1008                result[0] = ys[0];
    1009                result[1] = ys[1];
    1010            }
    1011
    1012            return ys[0] + ys[1];
    1013        }
    1014
    1015        /** Compute split[0], split[1] such that their sum is equal to d,
    1016         * and split[0] has its 30 least significant bits as zero.
    1017         * @param d number to split
    1018         * @param split placeholder where to place the result
    1019         */
    1020        private static void split(final double d, final double split[]) {
    1021            if (d < 8e298 && d > -8e298) {
    1022                final double a = d * HEX_40000000;
    1023                split[0] = (d + a) - a;
    1024                split[1] = d - split[0];
    1025            } else {
    1026                final double a = d * 9.31322574615478515625E-10;
    1027                split[0] = (d + a - d) * HEX_40000000;
    1028                split[1] = d - split[0];
    1029            }
    1030        }
    1031
    1032        /** Recompute a split.
    1033         * @param a input/out array containing the split, changed
    1034         * on output
    1035         */
    1036        private static void resplit(final double a[]) {
    1037            final double c = a[0] + a[1];
    1038            final double d = -(c - a[0] - a[1]);
    1039
    1040            if (c < 8e298 && c > -8e298) {
    1041                double z = c * HEX_40000000;
    1042                a[0] = (c + z) - z;
    1043                a[1] = c - a[0] + d;
    1044            } else {
    1045                double z = c * 9.31322574615478515625E-10;
    1046                a[0] = (c + z - c) * HEX_40000000;
    1047                a[1] = c - a[0] + d;
    1048            }
    1049        }
    1050
    1051        /** Multiply two numbers in split form.
    1052         * @param a first term of multiplication
    1053         * @param b second term of multiplication
    1054         * @param ans placeholder where to put the result
    1055         */
    1056        private static void splitMult(double a[], double b[], double ans[]) {
    1057            ans[0] = a[0] * b[0];
    1058            ans[1] = a[0] * b[1] + a[1] * b[0] + a[1] * b[1];
    1059
    1060            /* Resplit */
    1061            resplit(ans);
    1062        }
    1063
    1064        /** Add two numbers in split form.
    1065         * @param a first term of addition
    1066         * @param b second term of addition
    1067         * @param ans placeholder where to put the result
    1068         */
    1069        private static void splitAdd(final double a[], final double b[], final double ans[]) {
    1070            ans[0] = a[0] + b[0];
    1071            ans[1] = a[1] + b[1];
    1072
    1073            resplit(ans);
    1074        }
    1075
    1076        /** Compute the reciprocal of in.  Use the following algorithm.
    1077         *  in = c + d.
    1078         *  want to find x + y such that x+y = 1/(c+d) and x is much
    1079         *  larger than y and x has several zero bits on the right.
    1080         *
    1081         *  Set b = 1/(2^22),  a = 1 - b.  Thus (a+b) = 1.
    1082         *  Use following identity to compute (a+b)/(c+d)
    1083         *
    1084         *  (a+b)/(c+d)  =   a/c   +    (bc - ad) / (c^2 + cd)
    1085         *  set x = a/c  and y = (bc - ad) / (c^2 + cd)
    1086         *  This will be close to the right answer, but there will be
    1087         *  some rounding in the calculation of X.  So by carefully
    1088         *  computing 1 - (c+d)(x+y) we can compute an error and
    1089         *  add that back in.   This is done carefully so that terms
    1090         *  of similar size are subtracted first.
    1091         *  @param in initial number, in split form
    1092         *  @param result placeholder where to put the result
    1093         */
    1094        private static void splitReciprocal(final double in[], final double result[]) {
    1095            final double b = 1.0/4194304.0;
    1096            final double a = 1.0 - b;
    1097
    1098            if (in[0] == 0.0) {
    1099                in[0] = in[1];
    1100                in[1] = 0.0;
    1101            }
    1102
    1103            result[0] = a / in[0];
    1104            result[1] = (b*in[0]-a*in[1]) / (in[0]*in[0] + in[0]*in[1]);
    1105
    1106            if (result[1] != result[1]) { // can happen if result[1] is NAN
    1107                result[1] = 0.0;
    1108            }
    1109
    1110            /* Resplit */
    1111            resplit(result);
    1112
    1113            for (int i = 0; i < 2; i++) {
    1114                /* this may be overkill, probably once is enough */
    1115                double err = 1.0 - result[0] * in[0] - result[0] * in[1] -
    1116                result[1] * in[0] - result[1] * in[1];
    1117                /*err = 1.0 - err; */
    1118                err = err * (result[0] + result[1]);
    1119                /*printf("err = %16e\n", err); */
    1120                result[1] += err;
    1121            }
    1122        }
    1123
    1124        /** Compute (a[0] + a[1]) * (b[0] + b[1]) in extended precision.
    1125         * @param a first term of the multiplication
    1126         * @param b second term of the multiplication
    1127         * @param result placeholder where to put the result
    1128         */
    1129        private static void quadMult(final double a[], final double b[], final double result[]) {
    1130            final double xs[] = new double[2];
    1131            final double ys[] = new double[2];
    1132            final double zs[] = new double[2];
    1133
    1134            /* a[0] * b[0] */
    1135            split(a[0], xs);
    1136            split(b[0], ys);
    1137            splitMult(xs, ys, zs);
    1138
    1139            result[0] = zs[0];
    1140            result[1] = zs[1];
    1141
    1142            /* a[0] * b[1] */
    1143            split(b[1], ys);
    1144            splitMult(xs, ys, zs);
    1145
    1146            double tmp = result[0] + zs[0];
    1147            result[1] = result[1] - (tmp - result[0] - zs[0]);
    1148            result[0] = tmp;
    1149            tmp = result[0] + zs[1];
    1150            result[1] = result[1] - (tmp - result[0] - zs[1]);
    1151            result[0] = tmp;
    1152
    1153            /* a[1] * b[0] */
    1154            split(a[1], xs);
    1155            split(b[0], ys);
    1156            splitMult(xs, ys, zs);
    1157
    1158            tmp = result[0] + zs[0];
    1159            result[1] = result[1] - (tmp - result[0] - zs[0]);
    1160            result[0] = tmp;
    1161            tmp = result[0] + zs[1];
    1162            result[1] = result[1] - (tmp - result[0] - zs[1]);
    1163            result[0] = tmp;
    1164
    1165            /* a[1] * b[0] */
    1166            split(a[1], xs);
    1167            split(b[1], ys);
    1168            splitMult(xs, ys, zs);
    1169
    1170            tmp = result[0] + zs[0];
    1171            result[1] = result[1] - (tmp - result[0] - zs[0]);
    1172            result[0] = tmp;
    1173            tmp = result[0] + zs[1];
    1174            result[1] = result[1] - (tmp - result[0] - zs[1]);
    1175            result[0] = tmp;
    1176        }
    1177
    1178        /** Compute exp(p) for a integer p in extended precision.
    1179         * @param p integer whose exponential is requested
    1180         * @param result placeholder where to put the result in extended precision
    1181         * @return exp(p) in standard precision (equal to result[0] + result[1])
    1182         */
    1183        private static double expint(int p, final double result[]) {
    1184            //double x = M_E;
    1185            final double xs[] = new double[2];
    1186            final double as[] = new double[2];
    1187            final double ys[] = new double[2];
    1188            //split(x, xs);
    1189            //xs[1] = (double)(2.7182818284590452353602874713526625L - xs[0]);
    1190            //xs[0] = 2.71827697753906250000;
    1191            //xs[1] = 4.85091998273542816811e-06;
    1192            //xs[0] = Double.longBitsToDouble(0x4005bf0800000000L);
    1193            //xs[1] = Double.longBitsToDouble(0x3ed458a2bb4a9b00L);
    1194
    1195            /* E */
    1196            xs[0] = 2.718281828459045;
    1197            xs[1] = 1.4456468917292502E-16;
    1198
    1199            split(1.0, ys);
    1200
    1201            while (p > 0) {
    1202                if ((p & 1) != 0) {
    1203                    quadMult(ys, xs, as);
    1204                    ys[0] = as[0]; ys[1] = as[1];
    1205                }
    1206
    1207                quadMult(xs, xs, as);
    1208                xs[0] = as[0]; xs[1] = as[1];
    1209
    1210                p >>= 1;
    1211            }
    1212
    1213            if (result != null) {
    1214                result[0] = ys[0];
    1215                result[1] = ys[1];
    1216
    1217                resplit(result);
    1218            }
    1219
    1220            return ys[0] + ys[1];
    1221        }
    1222
    1223
    1224        /**
    1225         * Natural logarithm.
    1226         *
    1227         * @param x   a double
    1228         * @return log(x)
    1229         */
    1230        public static double log(final double x) {
    1231            return log(x, null);
    1232        }
    1233
    1234        /**
    1235         * Internal helper method for natural logarithm function.
    1236         * @param x original argument of the natural logarithm function
    1237         * @param hiPrec extra bits of precision on output (To Be Confirmed)
    1238         * @return log(x)
    1239         */
    1240        private static double log(final double x, final double[] hiPrec) {
    1241            if (x==0) { // Handle special case of +0/-0
    1242                return Double.NEGATIVE_INFINITY;
    1243            }
    1244            long bits = Double.doubleToLongBits(x);
    1245
    1246            /* Handle special cases of negative input, and NaN */
    1247            if ((bits & 0x8000000000000000L) != 0 || x != x) {
    1248                if (x != 0.0) {
    1249                    if (hiPrec != null) {
    1250                        hiPrec[0] = Double.NaN;
    1251                    }
    1252
    1253                    return Double.NaN;
    1254                }
    1255            }
    1256
    1257            /* Handle special cases of Positive infinity. */
    1258            if (x == Double.POSITIVE_INFINITY) {
    1259                if (hiPrec != null) {
    1260                    hiPrec[0] = Double.POSITIVE_INFINITY;
    1261                }
    1262
    1263                return Double.POSITIVE_INFINITY;
    1264            }
    1265
    1266            /* Extract the exponent */
    1267            int exp = (int)(bits >> 52)-1023;
    1268
    1269            if ((bits & 0x7ff0000000000000L) == 0) {
    1270                // Subnormal!
    1271                if (x == 0) {
    1272                    // Zero
    1273                    if (hiPrec != null) {
    1274                        hiPrec[0] = Double.NEGATIVE_INFINITY;
    1275                    }
    1276
    1277                    return Double.NEGATIVE_INFINITY;
    1278                }
    1279
    1280                /* Normalize the subnormal number. */
    1281                bits <<= 1;
    1282                while ( (bits & 0x0010000000000000L) == 0) {
    1283                    exp--;
    1284                    bits <<= 1;
    1285                }
    1286            }
    1287
    1288
    1289            if (exp == -1 || exp == 0) {
    1290                if (x < 1.01 && x > 0.99 && hiPrec == null) {
    1291                    /* The normal method doesn't work well in the range [0.99, 1.01], so call do a straight
    1292               polynomial expansion in higer precision. */
    1293
    1294                   /* Compute x - 1.0 and split it */
    1295                    double xa = x - 1.0;
    1296                    double xb = xa - x + 1.0;
    1297                    double tmp = xa * HEX_40000000;
    1298                    double aa = xa + tmp - tmp;
    1299                    double ab = xa - aa;
    1300                    xa = aa;
    1301                    xb = ab;
    1302
    1303                    double ya = LN_QUICK_COEF[LN_QUICK_COEF.length-1][0];
    1304                    double yb = LN_QUICK_COEF[LN_QUICK_COEF.length-1][1];
    1305
    1306                    for (int i = LN_QUICK_COEF.length - 2; i >= 0; i--) {
    1307                        /* Multiply a = y * x */
    1308                        aa = ya * xa;
    1309                        ab = ya * xb + yb * xa + yb * xb;
    1310                        /* split, so now y = a */
    1311                        tmp = aa * HEX_40000000;
    1312                        ya = aa + tmp - tmp;
    1313                        yb = aa - ya + ab;
    1314
    1315                        /* Add  a = y + lnQuickCoef */
    1316                        aa = ya + LN_QUICK_COEF[i][0];
    1317                        ab = yb + LN_QUICK_COEF[i][1];
    1318                        /* Split y = a */
    1319                        tmp = aa * HEX_40000000;
    1320                        ya = aa + tmp - tmp;
    1321                        yb = aa - ya + ab;
    1322                    }
    1323
    1324                    /* Multiply a = y * x */
    1325                    aa = ya * xa;
    1326                    ab = ya * xb + yb * xa + yb * xb;
    1327                    /* split, so now y = a */
    1328                    tmp = aa * HEX_40000000;
    1329                    ya = aa + tmp - tmp;
    1330                    yb = aa - ya + ab;
    1331
    1332                    return ya + yb;
    1333                }
    1334            }
    1335
    1336            // lnm is a log of a number in the range of 1.0 - 2.0, so 0 <= lnm < ln(2)
    1337            double lnm[] = LN_MANT[(int)((bits & 0x000ffc0000000000L) >> 42)];
    1338
    1339            /*
    1340        double epsilon = x / Double.longBitsToDouble(bits & 0xfffffc0000000000L);
    1341
    1342        epsilon -= 1.0;
    1343             */
    1344
    1345            // y is the most significant 10 bits of the mantissa
    1346            //double y = Double.longBitsToDouble(bits & 0xfffffc0000000000L);
    1347            //double epsilon = (x - y) / y;
    1348            double epsilon = (bits & 0x3ffffffffffL) / (TWO_POWER_52 + (bits & 0x000ffc0000000000L));
    1349
    1350            double lnza = 0.0;
    1351            double lnzb = 0.0;
    1352
    1353            if (hiPrec != null) {
    1354                /* split epsilon -> x */
    1355                double tmp = epsilon * HEX_40000000;
    1356                double aa = epsilon + tmp - tmp;
    1357                double ab = epsilon - aa;
    1358                double xa = aa;
    1359                double xb = ab;
    1360
    1361                /* Need a more accurate epsilon, so adjust the division. */
    1362                double numer = bits & 0x3ffffffffffL;
    1363                double denom = TWO_POWER_52 + (bits & 0x000ffc0000000000L);
    1364                aa = numer - xa*denom - xb * denom;
    1365                xb += aa / denom;
    1366
    1367                /* Remez polynomial evaluation */
    1368                double ya = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1][0];
    1369                double yb = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1][1];
    1370
    1371                for (int i = LN_HI_PREC_COEF.length - 2; i >= 0; i--) {
    1372                    /* Multiply a = y * x */
    1373                    aa = ya * xa;
    1374                    ab = ya * xb + yb * xa + yb * xb;
    1375                    /* split, so now y = a */
    1376                    tmp = aa * HEX_40000000;
    1377                    ya = aa + tmp - tmp;
    1378                    yb = aa - ya + ab;
    1379
    1380                    /* Add  a = y + lnHiPrecCoef */
    1381                    aa = ya + LN_HI_PREC_COEF[i][0];
    1382                    ab = yb + LN_HI_PREC_COEF[i][1];
    1383                    /* Split y = a */
    1384                    tmp = aa * HEX_40000000;
    1385                    ya = aa + tmp - tmp;
    1386                    yb = aa - ya + ab;
    1387                }
    1388
    1389                /* Multiply a = y * x */
    1390                aa = ya * xa;
    1391                ab = ya * xb + yb * xa + yb * xb;
    1392
    1393                /* split, so now lnz = a */
    1394                /*
    1395          tmp = aa * 1073741824.0;
    1396          lnza = aa + tmp - tmp;
    1397          lnzb = aa - lnza + ab;
    1398                 */
    1399                lnza = aa + ab;
    1400                lnzb = -(lnza - aa - ab);
    1401            } else {
    1402                /* High precision not required.  Eval Remez polynomial
    1403             using standard double precision */
    1404                lnza = -0.16624882440418567;
    1405                lnza = lnza * epsilon + 0.19999954120254515;
    1406                lnza = lnza * epsilon + -0.2499999997677497;
    1407                lnza = lnza * epsilon + 0.3333333333332802;
    1408                lnza = lnza * epsilon + -0.5;
    1409                lnza = lnza * epsilon + 1.0;
    1410                lnza = lnza * epsilon;
    1411            }
    1412
    1413            /* Relative sizes:
    1414             * lnzb     [0, 2.33E-10]
    1415             * lnm[1]   [0, 1.17E-7]
    1416             * ln2B*exp [0, 1.12E-4]
    1417             * lnza      [0, 9.7E-4]
    1418             * lnm[0]   [0, 0.692]
    1419             * ln2A*exp [0, 709]
    1420             */
    1421
    1422            /* Compute the following sum:
    1423             * lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
    1424             */
    1425
    1426            //return lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
    1427            double a = LN_2_A*exp;
    1428            double b = 0.0;
    1429            double c = a+lnm[0];
    1430            double d = -(c-a-lnm[0]);
    1431            a = c;
    1432            b = b + d;
    1433
    1434            c = a + lnza;
    1435            d = -(c - a - lnza);
    1436            a = c;
    1437            b = b + d;
    1438
    1439            c = a + LN_2_B*exp;
    1440            d = -(c - a - LN_2_B*exp);
    1441            a = c;
    1442            b = b + d;
    1443
    1444            c = a + lnm[1];
    1445            d = -(c - a - lnm[1]);
    1446            a = c;
    1447            b = b + d;
    1448
    1449            c = a + lnzb;
    1450            d = -(c - a - lnzb);
    1451            a = c;
    1452            b = b + d;
    1453
    1454            if (hiPrec != null) {
    1455                hiPrec[0] = a;
    1456                hiPrec[1] = b;
    1457            }
    1458
    1459            return a + b;
    1460        }
    1461
    1462        /** Compute log(1 + x).
    1463         * @param x a number
    1464         * @return log(1 + x)
    1465         */
    1466        public static double log1p(final double x) {
    1467            double xpa = 1.0 + x;
    1468            double xpb = -(xpa - 1.0 - x);
    1469
    1470            if (x == -1) {
    1471                return x/0.0;   // -Infinity
    1472            }
    1473
    1474            if (x > 0 && 1/x == 0) { // x = Infinity
    1475                return x;
    1476            }
    1477
    1478            if (x>1e-6 || x<-1e-6) {
    1479                double hiPrec[] = new double[2];
    1480
    1481                final double lores = log(xpa, hiPrec);
    1482                if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN
    1483                    return lores;
    1484                }
    1485
    1486                /* Do a taylor series expansion around xpa */
    1487                /* f(x+y) = f(x) + f'(x)*y + f''(x)/2 y^2 */
    1488                double fx1 = xpb/xpa;
    1489
    1490                double epsilon = 0.5 * fx1 + 1.0;
    1491                epsilon = epsilon * fx1;
    1492
    1493                return epsilon + hiPrec[1] + hiPrec[0];
    1494            }
    1495
    1496            /* Value is small |x| < 1e6, do a Taylor series centered on 1.0 */
    1497            double y = x * 0.333333333333333 - 0.5;
    1498            y = y * x + 1.0;
    1499            y = y * x;
    1500
    1501            return y;
    1502        }
    1503
    1504        /** Compute the base 10 logarithm.
    1505         * @param x a number
    1506         * @return log10(x)
    1507         */
    1508        public static double log10(final double x) {
    1509            final double hiPrec[] = new double[2];
    1510
    1511            final double lores = log(x, hiPrec);
    1512            if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN
    1513                return lores;
    1514            }
    1515
    1516            final double tmp = hiPrec[0] * HEX_40000000;
    1517            final double lna = hiPrec[0] + tmp - tmp;
    1518            final double lnb = hiPrec[0] - lna + hiPrec[1];
    1519
    1520            final double rln10a = 0.4342944622039795;
    1521            final double rln10b = 1.9699272335463627E-8;
    1522
    1523            return rln10b * lnb + rln10b * lna + rln10a * lnb + rln10a * lna;
    1524        }
    1525
    1526        /**
    1527         * Power function.  Compute x^y.
    1528         *
    1529         * @param x   a double
    1530         * @param y   a double
    1531         * @return double
    1532         */
    1533        public static double pow(double x, double y) {
    1534            final double lns[] = new double[2];
    1535
    1536            if (y == 0.0) {
    1537                return 1.0;
    1538            }
    1539
    1540            if (x != x) { // X is NaN
    1541                return x;
    1542            }
    1543
    1544
    1545            if (x == 0) {
    1546                long bits = Double.doubleToLongBits(x);
    1547                if ((bits & 0x8000000000000000L) != 0) {
    1548                    // -zero
    1549                    long yi = (long) y;
    1550
    1551                    if (y < 0 && y == yi && (yi & 1) == 1) {
    1552                        return Double.NEGATIVE_INFINITY;
    1553                    }
    1554
    1555                    if (y < 0 && y == yi && (yi & 1) == 1) {
    1556                        return -0.0;
    1557                    }
    1558
    1559                    if (y > 0 && y == yi && (yi & 1) == 1) {
    1560                        return -0.0;
    1561                    }
    1562                }
    1563
    1564                if (y < 0) {
    1565                    return Double.POSITIVE_INFINITY;
    1566                }
    1567                if (y > 0) {
    1568                    return 0.0;
    1569                }
    1570
    1571                return Double.NaN;
    1572            }
    1573
    1574            if (x == Double.POSITIVE_INFINITY) {
    1575                if (y != y) { // y is NaN
    1576                    return y;
    1577                }
    1578                if (y < 0.0) {
    1579                    return 0.0;
    1580                } else {
    1581                    return Double.POSITIVE_INFINITY;
    1582                }
    1583            }
    1584
    1585            if (y == Double.POSITIVE_INFINITY) {
    1586                if (x * x == 1.0)
    1587                  return Double.NaN;
    1588
    1589                if (x * x > 1.0) {
    1590                    return Double.POSITIVE_INFINITY;
    1591                } else {
    1592                    return 0.0;
    1593                }
    1594            }
    1595
    1596            if (x == Double.NEGATIVE_INFINITY) {
    1597                if (y != y) { // y is NaN
    1598                    return y;
    1599                }
    1600
    1601                if (y < 0) {
    1602                    long yi = (long) y;
    1603                    if (y == yi && (yi & 1) == 1) {
    1604                        return -0.0;
    1605                    }
    1606
    1607                    return 0.0;
    1608                }
    1609
    1610                if (y > 0)  {
    1611                    long yi = (long) y;
    1612                    if (y == yi && (yi & 1) == 1) {
    1613                        return Double.NEGATIVE_INFINITY;
    1614                    }
    1615
    1616                    return Double.POSITIVE_INFINITY;
    1617                }
    1618            }
    1619
    1620            if (y == Double.NEGATIVE_INFINITY) {
    1621
    1622                if (x * x == 1.0) {
    1623                    return Double.NaN;
    1624                }
    1625
    1626                if (x * x < 1.0) {
    1627                    return Double.POSITIVE_INFINITY;
    1628                } else {
    1629                    return 0.0;
    1630                }
    1631            }
    1632
    1633            /* Handle special case x<0 */
    1634            if (x < 0) {
    1635                // y is an even integer in this case
    1636                if (y >= TWO_POWER_52 || y <= -TWO_POWER_52) {
    1637                    return pow(-x, y);
    1638                }
    1639
    1640                if (y == (long) y) {
    1641                    // If y is an integer
    1642                    return ((long)y & 1) == 0 ? pow(-x, y) : -pow(-x, y);
    1643                } else {
    1644                    return Double.NaN;
    1645                }
    1646            }
    1647
    1648            /* Split y into ya and yb such that y = ya+yb */
    1649            double ya;
    1650            double yb;
    1651            if (y < 8e298 && y > -8e298) {
    1652                double tmp1 = y * HEX_40000000;
    1653                ya = y + tmp1 - tmp1;
    1654                yb = y - ya;
    1655            } else {
    1656                double tmp1 = y * 9.31322574615478515625E-10;
    1657                double tmp2 = tmp1 * 9.31322574615478515625E-10;
    1658                ya = (tmp1 + tmp2 - tmp1) * HEX_40000000 * HEX_40000000;
    1659                yb = y - ya;
    1660            }
    1661
    1662            /* Compute ln(x) */
    1663            final double lores = log(x, lns);
    1664            if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN
    1665                return lores;
    1666            }
    1667
    1668            double lna = lns[0];
    1669            double lnb = lns[1];
    1670
    1671            /* resplit lns */
    1672            double tmp1 = lna * HEX_40000000;
    1673            double tmp2 = lna + tmp1 - tmp1;
    1674            lnb += lna - tmp2;
    1675            lna = tmp2;
    1676
    1677            // y*ln(x) = (aa+ab)
    1678            final double aa = lna * ya;
    1679            final double ab = lna * yb + lnb * ya + lnb * yb;
    1680
    1681            lna = aa+ab;
    1682            lnb = -(lna - aa - ab);
    1683
    1684            double z = 1.0 / 120.0;
    1685            z = z * lnb + (1.0 / 24.0);
    1686            z = z * lnb + (1.0 / 6.0);
    1687            z = z * lnb + 0.5;
    1688            z = z * lnb + 1.0;
    1689            z = z * lnb;
    1690
    1691            final double result = exp(lna, z, null);
    1692            //result = result + result * z;
    1693            return result;
    1694        }
    1695
    1696        /** xi in the range of [1, 2].
    1697         *                                3        5        7
    1698         *      x+1           /          x        x        x          \
    1699         *  ln ----- =   2 *  |  x  +   ----  +  ----  +  ---- + ...  |
    1700         *      1-x           \          3        5        7          /
    1701         *
    1702         * So, compute a Remez approximation of the following function
    1703         *
    1704         *  ln ((sqrt(x)+1)/(1-sqrt(x)))  /  x
    1705         *
    1706         * This will be an even function with only positive coefficents.
    1707         * x is in the range [0 - 1/3].
    1708         *
    1709         * Transform xi for input to the above function by setting
    1710         * x = (xi-1)/(xi+1).   Input to the polynomial is x^2, then
    1711         * the result is multiplied by x.
    1712         * @param xi number from which log is requested
    1713         * @return log(xi)
    1714         */
    1715        private static double[] slowLog(double xi) {
    1716            double x[] = new double[2];
    1717            double x2[] = new double[2];
    1718            double y[] = new double[2];
    1719            double a[] = new double[2];
    1720
    1721            split(xi, x);
    1722
    1723            /* Set X = (x-1)/(x+1) */
    1724            x[0] += 1.0;
    1725            resplit(x);
    1726            splitReciprocal(x, a);
    1727            x[0] -= 2.0;
    1728            resplit(x);
    1729            splitMult(x, a, y);
    1730            x[0] = y[0];
    1731            x[1] = y[1];
    1732
    1733            /* Square X -> X2*/
    1734            splitMult(x, x, x2);
    1735
    1736
    1737            //x[0] -= 1.0;
    1738            //resplit(x);
    1739
    1740            y[0] = LN_SPLIT_COEF[LN_SPLIT_COEF.length-1][0];
    1741            y[1] = LN_SPLIT_COEF[LN_SPLIT_COEF.length-1][1];
    1742
    1743            for (int i = LN_SPLIT_COEF.length-2; i >= 0; i--) {
    1744                splitMult(y, x2, a);
    1745                y[0] = a[0];
    1746                y[1] = a[1];
    1747                splitAdd(y, LN_SPLIT_COEF[i], a);
    1748                y[0] = a[0];
    1749                y[1] = a[1];
    1750            }
    1751
    1752            splitMult(y, x, a);
    1753            y[0] = a[0];
    1754            y[1] = a[1];
    1755
    1756            return y;
    1757        }
    1758
    1759        /**
    1760         * For x between 0 and pi/4 compute sine.
    1761         * @param x number from which sine is requested
    1762         * @param result placeholder where to put the result in extended precision
    1763         * @return sin(x)
    1764         */
    1765        private static double slowSin(final double x, final double result[]) {
    1766            final double xs[] = new double[2];
    1767            final double ys[] = new double[2];
    1768            final double facts[] = new double[2];
    1769            final double as[] = new double[2];
    1770            split(x, xs);
    1771            ys[0] = ys[1] = 0.0;
    1772
    1773            for (int i = 19; i >= 0; i--) {
    1774                splitMult(xs, ys, as);
    1775                ys[0] = as[0]; ys[1] = as[1];
    1776
    1777                if ( (i & 1) == 0) {
    1778                    continue;
    1779                }
    1780
    1781                split(FACT[i], as);
    1782                splitReciprocal(as, facts);
    1783
    1784                if ( (i & 2) != 0 ) {
    1785                    facts[0] = -facts[0];
    1786                    facts[1] = -facts[1];
    1787                }
    1788
    1789                splitAdd(ys, facts, as);
    1790                ys[0] = as[0]; ys[1] = as[1];
    1791            }
    1792
    1793            if (result != null) {
    1794                result[0] = ys[0];
    1795                result[1] = ys[1];
    1796            }
    1797
    1798            return ys[0] + ys[1];
    1799        }
    1800
    1801        /**
    1802         *  For x between 0 and pi/4 compute cosine
    1803         * @param x number from which cosine is requested
    1804         * @param result placeholder where to put the result in extended precision
    1805         * @return cos(x)
    1806         */
    1807        private static double slowCos(final double x, final double result[]) {
    1808
    1809            final double xs[] = new double[2];
    1810            final double ys[] = new double[2];
    1811            final double facts[] = new double[2];
    1812            final double as[] = new double[2];
    1813            split(x, xs);
    1814            ys[0] = ys[1] = 0.0;
    1815
    1816            for (int i = 19; i >= 0; i--) {
    1817                splitMult(xs, ys, as);
    1818                ys[0] = as[0]; ys[1] = as[1];
    1819
    1820                if ( (i & 1) != 0) {
    1821                    continue;
    1822                }
    1823
    1824                split(FACT[i], as);
    1825                splitReciprocal(as, facts);
    1826
    1827                if ( (i & 2) != 0 ) {
    1828                    facts[0] = -facts[0];
    1829                    facts[1] = -facts[1];
    1830                }
    1831
    1832                splitAdd(ys, facts, as);
    1833                ys[0] = as[0]; ys[1] = as[1];
    1834            }
    1835
    1836            if (result != null) {
    1837                result[0] = ys[0];
    1838                result[1] = ys[1];
    1839            }
    1840
    1841            return ys[0] + ys[1];
    1842        }
    1843
    1844        /** Build the sine and cosine tables.
    1845         */
    1846        private static void buildSinCosTables() {
    1847            final double result[] = new double[2];
    1848
    1849            /* Use taylor series for 0 <= x <= 6/8 */
    1850            for (int i = 0; i < 7; i++) {
    1851                double x = i / 8.0;
    1852
    1853                slowSin(x, result);
    1854                SINE_TABLE_A[i] = result[0];
    1855                SINE_TABLE_B[i] = result[1];
    1856
    1857                slowCos(x, result);
    1858                COSINE_TABLE_A[i] = result[0];
    1859                COSINE_TABLE_B[i] = result[1];
    1860            }
    1861
    1862            /* Use angle addition formula to complete table to 13/8, just beyond pi/2 */
    1863            for (int i = 7; i < 14; i++) {
    1864                double xs[] = new double[2];
    1865                double ys[] = new double[2];
    1866                double as[] = new double[2];
    1867                double bs[] = new double[2];
    1868                double temps[] = new double[2];
    1869
    1870                if ( (i & 1) == 0) {
    1871                    // Even, use double angle
    1872                    xs[0] = SINE_TABLE_A[i/2];
    1873                    xs[1] = SINE_TABLE_B[i/2];
    1874                    ys[0] = COSINE_TABLE_A[i/2];
    1875                    ys[1] = COSINE_TABLE_B[i/2];
    1876
    1877                    /* compute sine */
    1878                    splitMult(xs, ys, result);
    1879                    SINE_TABLE_A[i] = result[0] * 2.0;
    1880                    SINE_TABLE_B[i] = result[1] * 2.0;
    1881
    1882                    /* Compute cosine */
    1883                    splitMult(ys, ys, as);
    1884                    splitMult(xs, xs, temps);
    1885                    temps[0] = -temps[0];
    1886                    temps[1] = -temps[1];
    1887                    splitAdd(as, temps, result);
    1888                    COSINE_TABLE_A[i] = result[0];
    1889                    COSINE_TABLE_B[i] = result[1];
    1890                } else {
    1891                    xs[0] = SINE_TABLE_A[i/2];
    1892                    xs[1] = SINE_TABLE_B[i/2];
    1893                    ys[0] = COSINE_TABLE_A[i/2];
    1894                    ys[1] = COSINE_TABLE_B[i/2];
    1895                    as[0] = SINE_TABLE_A[i/2+1];
    1896                    as[1] = SINE_TABLE_B[i/2+1];
    1897                    bs[0] = COSINE_TABLE_A[i/2+1];
    1898                    bs[1] = COSINE_TABLE_B[i/2+1];
    1899
    1900                    /* compute sine */
    1901                    splitMult(xs, bs, temps);
    1902                    splitMult(ys, as, result);
    1903                    splitAdd(result, temps, result);
    1904                    SINE_TABLE_A[i] = result[0];
    1905                    SINE_TABLE_B[i] = result[1];
    1906
    1907                    /* Compute cosine */
    1908                    splitMult(ys, bs, result);
    1909                    splitMult(xs, as, temps);
    1910                    temps[0] = -temps[0];
    1911                    temps[1] = -temps[1];
    1912                    splitAdd(result, temps, result);
    1913                    COSINE_TABLE_A[i] = result[0];
    1914                    COSINE_TABLE_B[i] = result[1];
    1915                }
    1916            }
    1917
    1918            /* Compute tangent = sine/cosine */
    1919            for (int i = 0; i < 14; i++) {
    1920                double xs[] = new double[2];
    1921                double ys[] = new double[2];
    1922                double as[] = new double[2];
    1923
    1924                as[0] = COSINE_TABLE_A[i];
    1925                as[1] = COSINE_TABLE_B[i];
    1926
    1927                splitReciprocal(as, ys);
    1928
    1929                xs[0] = SINE_TABLE_A[i];
    1930                xs[1] = SINE_TABLE_B[i];
    1931
    1932                splitMult(xs, ys, as);
    1933
    1934                TANGENT_TABLE_A[i] = as[0];
    1935                TANGENT_TABLE_B[i] = as[1];
    1936            }
    1937
    1938        }
    1939
    1940        /**
    1941         *  Computes sin(x) - x, where |x| < 1/16.
    1942         *  Use a Remez polynomial approximation.
    1943         *  @param x a number smaller than 1/16
    1944         *  @return sin(x) - x
    1945         */
    1946        private static double polySine(final double x)
    1947        {
    1948            double x2 = x*x;
    1949
    1950            double p = 2.7553817452272217E-6;
    1951            p = p * x2 + -1.9841269659586505E-4;
    1952            p = p * x2 + 0.008333333333329196;
    1953            p = p * x2 + -0.16666666666666666;
    1954            //p *= x2;
    1955            //p *= x;
    1956            p = p * x2 * x;
    1957
    1958            return p;
    1959        }
    1960
    1961        /**
    1962         *  Computes cos(x) - 1, where |x| < 1/16.
    1963         *  Use a Remez polynomial approximation.
    1964         *  @param x a number smaller than 1/16
    1965         *  @return cos(x) - 1
    1966         */
    1967        private static double polyCosine(double x) {
    1968            double x2 = x*x;
    1969
    1970            double p = 2.479773539153719E-5;
    1971            p = p * x2 + -0.0013888888689039883;
    1972            p = p * x2 + 0.041666666666621166;
    1973            p = p * x2 + -0.49999999999999994;
    1974            p *= x2;
    1975
    1976            return p;
    1977        }
    1978
    1979        /**
    1980         *  Compute sine over the first quadrant (0 < x < pi/2).
    1981         *  Use combination of table lookup and rational polynomial expansion.
    1982         *  @param xa number from which sine is requested
    1983         *  @param xb extra bits for x (may be 0.0)
    1984         *  @return sin(xa + xb)
    1985         */
    1986        private static double sinQ(double xa, double xb) {
    1987            int idx = (int) ((xa * 8.0) + 0.5);
    1988            final double epsilon = xa - EIGHTHS[idx]; //idx*0.125;
    1989
    1990            // Table lookups
    1991            final double sintA = SINE_TABLE_A[idx];
    1992            final double sintB = SINE_TABLE_B[idx];
    1993            final double costA = COSINE_TABLE_A[idx];
    1994            final double costB = COSINE_TABLE_B[idx];
    1995
    1996            // Polynomial eval of sin(epsilon), cos(epsilon)
    1997            double sinEpsA = epsilon;
    1998            double sinEpsB = polySine(epsilon);
    1999            final double cosEpsA = 1.0;
    2000            final double cosEpsB = polyCosine(epsilon);
    2001
    2002            // Split epsilon   xa + xb = x
    2003            final double temp = sinEpsA * HEX_40000000;
    2004            double temp2 = (sinEpsA + temp) - temp;
    2005            sinEpsB +=  sinEpsA - temp2;
    2006            sinEpsA = temp2;
    2007
    2008            /* Compute sin(x) by angle addition formula */
    2009            double result;
    2010
    2011            /* Compute the following sum:
    2012             *
    2013             * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
    2014             *          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
    2015             *
    2016             * Ranges of elements
    2017             *
    2018             * xxxtA   0            PI/2
    2019             * xxxtB   -1.5e-9      1.5e-9
    2020             * sinEpsA -0.0625      0.0625
    2021             * sinEpsB -6e-11       6e-11
    2022             * cosEpsA  1.0
    2023             * cosEpsB  0           -0.0625
    2024             *
    2025             */
    2026
    2027            //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
    2028            //          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
    2029
    2030            //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
    2031            //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
    2032            double a = 0;
    2033            double b = 0;
    2034
    2035            double t = sintA;
    2036            double c = a + t;
    2037            double d = -(c - a - t);
    2038            a = c;
    2039            b = b + d;
    2040
    2041            t = costA * sinEpsA;
    2042            c = a + t;
    2043            d = -(c - a - t);
    2044            a = c;
    2045            b = b + d;
    2046
    2047            b = b + sintA * cosEpsB + costA * sinEpsB;
    2048            /*
    2049        t = sintA*cosEpsB;
    2050        c = a + t;
    2051        d = -(c - a - t);
    2052        a = c;
    2053        b = b + d;
    2054
    2055        t = costA*sinEpsB;
    2056        c = a + t;
    2057        d = -(c - a - t);
    2058        a = c;
    2059        b = b + d;
    2060             */
    2061
    2062            b = b + sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB;
    2063            /*
    2064        t = sintB;
    2065        c = a + t;
    2066        d = -(c - a - t);
    2067        a = c;
    2068        b = b + d;
    2069
    2070        t = costB*sinEpsA;
    2071        c = a + t;
    2072        d = -(c - a - t);
    2073        a = c;
    2074        b = b + d;
    2075
    2076        t = sintB*cosEpsB;
    2077        c = a + t;
    2078        d = -(c - a - t);
    2079        a = c;
    2080        b = b + d;
    2081
    2082        t = costB*sinEpsB;
    2083        c = a + t;
    2084        d = -(c - a - t);
    2085        a = c;
    2086        b = b + d;
    2087             */
    2088
    2089            if (xb != 0.0) {
    2090                t = ((costA + costB) * (cosEpsA + cosEpsB) -
    2091                     (sintA + sintB) * (sinEpsA + sinEpsB)) * xb;  // approximate cosine*xb
    2092                c = a + t;
    2093                d = -(c - a - t);
    2094                a = c;
    2095                b = b + d;
    2096            }
    2097
    2098            result = a + b;
    2099
    2100            return result;
    2101        }
    2102
    2103        /**
    2104         * Compute cosine in the first quadrant by subtracting input from PI/2 and
    2105         * then calling sinQ.  This is more accurate as the input approaches PI/2.
    2106         *  @param xa number from which cosine is requested
    2107         *  @param xb extra bits for x (may be 0.0)
    2108         *  @return cos(xa + xb)
    2109         */
    2110        private static double cosQ(double xa, double xb) {
    2111            final double pi2a = 1.5707963267948966;
    2112            final double pi2b = 6.123233995736766E-17;
    2113
    2114            final double a = pi2a - xa;
    2115            double b = -(a - pi2a + xa);
    2116            b += pi2b - xb;
    2117
    2118            return sinQ(a, b);
    2119        }
    2120
    2121        /**
    2122         *  Compute tangent (or cotangent) over the first quadrant.   0 < x < pi/2
    2123         *  Use combination of table lookup and rational polynomial expansion.
    2124         *  @param xa number from which sine is requested
    2125         *  @param xb extra bits for x (may be 0.0)
    2126         *  @param cotanFlag if true, compute the cotangent instead of the tangent
    2127         *  @return tan(xa+xb) (or cotangent, depending on cotanFlag)
    2128         */
    2129        private static double tanQ(double xa, double xb, boolean cotanFlag) {
    2130
    2131            int idx = (int) ((xa * 8.0) + 0.5);
    2132            final double epsilon = xa - EIGHTHS[idx]; //idx*0.125;
    2133
    2134            // Table lookups
    2135            final double sintA = SINE_TABLE_A[idx];
    2136            final double sintB = SINE_TABLE_B[idx];
    2137            final double costA = COSINE_TABLE_A[idx];
    2138            final double costB = COSINE_TABLE_B[idx];
    2139
    2140            // Polynomial eval of sin(epsilon), cos(epsilon)
    2141            double sinEpsA = epsilon;
    2142            double sinEpsB = polySine(epsilon);
    2143            final double cosEpsA = 1.0;
    2144            final double cosEpsB = polyCosine(epsilon);
    2145
    2146            // Split epsilon   xa + xb = x
    2147            double temp = sinEpsA * HEX_40000000;
    2148            double temp2 = (sinEpsA + temp) - temp;
    2149            sinEpsB +=  sinEpsA - temp2;
    2150            sinEpsA = temp2;
    2151
    2152            /* Compute sin(x) by angle addition formula */
    2153
    2154            /* Compute the following sum:
    2155             *
    2156             * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
    2157             *          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
    2158             *
    2159             * Ranges of elements
    2160             *
    2161             * xxxtA   0            PI/2
    2162             * xxxtB   -1.5e-9      1.5e-9
    2163             * sinEpsA -0.0625      0.0625
    2164             * sinEpsB -6e-11       6e-11
    2165             * cosEpsA  1.0
    2166             * cosEpsB  0           -0.0625
    2167             *
    2168             */
    2169
    2170            //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
    2171            //          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
    2172
    2173            //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
    2174            //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
    2175            double a = 0;
    2176            double b = 0;
    2177
    2178            // Compute sine
    2179            double t = sintA;
    2180            double c = a + t;
    2181            double d = -(c - a - t);
    2182            a = c;
    2183            b = b + d;
    2184
    2185            t = costA*sinEpsA;
    2186            c = a + t;
    2187            d = -(c - a - t);
    2188            a = c;
    2189            b = b + d;
    2190
    2191            b = b + sintA*cosEpsB + costA*sinEpsB;
    2192            b = b + sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
    2193
    2194            double sina = a + b;
    2195            double sinb = -(sina - a - b);
    2196
    2197            // Compute cosine
    2198
    2199            a = b = c = d = 0.0;
    2200
    2201            t = costA*cosEpsA;
    2202            c = a + t;
    2203            d = -(c - a - t);
    2204            a = c;
    2205            b = b + d;
    2206
    2207            t = -sintA*sinEpsA;
    2208            c = a + t;
    2209            d = -(c - a - t);
    2210            a = c;
    2211            b = b + d;
    2212
    2213            b = b + costB*cosEpsA + costA*cosEpsB + costB*cosEpsB;
    2214            b = b - (sintB*sinEpsA + sintA*sinEpsB + sintB*sinEpsB);
    2215
    2216            double cosa = a + b;
    2217            double cosb = -(cosa - a - b);
    2218
    2219            if (cotanFlag) {
    2220                double tmp;
    2221                tmp = cosa; cosa = sina; sina = tmp;
    2222                tmp = cosb; cosb = sinb; sinb = tmp;
    2223            }
    2224
    2225
    2226            /* estimate and correct, compute 1.0/(cosa+cosb) */
    2227            /*
    2228        double est = (sina+sinb)/(cosa+cosb);
    2229        double err = (sina - cosa*est) + (sinb - cosb*est);
    2230        est += err/(cosa+cosb);
    2231        err = (sina - cosa*est) + (sinb - cosb*est);
    2232             */
    2233
    2234            // f(x) = 1/x,   f'(x) = -1/x^2
    2235
    2236            double est = sina/cosa;
    2237
    2238            /* Split the estimate to get more accurate read on division rounding */
    2239            temp = est * HEX_40000000;
    2240            double esta = (est + temp) - temp;
    2241            double estb =  est - esta;
    2242
    2243            temp = cosa * HEX_40000000;
    2244            double cosaa = (cosa + temp) - temp;
    2245            double cosab =  cosa - cosaa;
    2246
    2247            //double err = (sina - est*cosa)/cosa;  // Correction for division rounding
    2248            double err = (sina - esta*cosaa - esta*cosab - estb*cosaa - estb*cosab)/cosa;  // Correction for division rounding
    2249            err += sinb/cosa;                     // Change in est due to sinb
    2250            err += -sina * cosb / cosa / cosa;    // Change in est due to cosb
    2251
    2252            if (xb != 0.0) {
    2253                // tan' = 1 + tan^2      cot' = -(1 + cot^2)
    2254                // Approximate impact of xb
    2255                double xbadj = xb + est*est*xb;
    2256                if (cotanFlag) {
    2257                    xbadj = -xbadj;
    2258                }
    2259
    2260                err += xbadj;
    2261            }
    2262
    2263            return est+err;
    2264        }
    2265
    2266        /** Reduce the input argument using the Payne and Hanek method.
    2267         *  This is good for all inputs 0.0 < x < inf
    2268         *  Output is remainder after dividing by PI/2
    2269         *  The result array should contain 3 numbers.
    2270         *  result[0] is the integer portion, so mod 4 this gives the quadrant.
    2271         *  result[1] is the upper bits of the remainder
    2272         *  result[2] is the lower bits of the remainder
    2273         *
    2274         * @param x number to reduce
    2275         * @param result placeholder where to put the result
    2276         */
    2277        private static void reducePayneHanek(double x, double result[])
    2278        {
    2279            /* Convert input double to bits */
    2280            long inbits = Double.doubleToLongBits(x);
    2281            int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
    2282
    2283            /* Convert to fixed point representation */
    2284            inbits &= 0x000fffffffffffffL;
    2285            inbits |= 0x0010000000000000L;
    2286
    2287            /* Normalize input to be between 0.5 and 1.0 */
    2288            exponent++;
    2289            inbits <<= 11;
    2290
    2291            /* Based on the exponent, get a shifted copy of recip2pi */
    2292            long shpi0;
    2293            long shpiA;
    2294            long shpiB;
    2295            int idx = exponent >> 6;
    2296            int shift = exponent - (idx << 6);
    2297
    2298            if (shift != 0) {
    2299                shpi0 = (idx == 0) ? 0 : (RECIP_2PI[idx-1] << shift);
    2300                shpi0 |= RECIP_2PI[idx] >>> (64-shift);
    2301                shpiA = (RECIP_2PI[idx] << shift) | (RECIP_2PI[idx+1] >>> (64-shift));
    2302                shpiB = (RECIP_2PI[idx+1] << shift) | (RECIP_2PI[idx+2] >>> (64-shift));
    2303            } else {
    2304                shpi0 = (idx == 0) ? 0 : RECIP_2PI[idx-1];
    2305                shpiA = RECIP_2PI[idx];
    2306                shpiB = RECIP_2PI[idx+1];
    2307            }
    2308
    2309            /* Multiply input by shpiA */
    2310            long a = inbits >>> 32;
    2311            long b = inbits & 0xffffffffL;
    2312
    2313            long c = shpiA >>> 32;
    2314            long d = shpiA & 0xffffffffL;
    2315
    2316            long ac = a * c;
    2317            long bd = b * d;
    2318            long bc = b * c;
    2319            long ad = a * d;
    2320
    2321            long prodB = bd + (ad << 32);
    2322            long prodA = ac + (ad >>> 32);
    2323
    2324            boolean bita = (bd & 0x8000000000000000L) != 0;
    2325            boolean bitb = (ad & 0x80000000L ) != 0;
    2326            boolean bitsum = (prodB & 0x8000000000000000L) != 0;
    2327
    2328            /* Carry */
    2329            if ( (bita && bitb) ||
    2330                    ((bita || bitb) && !bitsum) ) {
    2331                prodA++;
    2332            }
    2333
    2334            bita = (prodB & 0x8000000000000000L) != 0;
    2335            bitb = (bc & 0x80000000L ) != 0;
    2336
    2337            prodB = prodB + (bc << 32);
    2338            prodA = prodA + (bc >>> 32);
    2339
    2340            bitsum = (prodB & 0x8000000000000000L) != 0;
    2341
    2342            /* Carry */
    2343            if ( (bita && bitb) ||
    2344                    ((bita || bitb) && !bitsum) ) {
    2345                prodA++;
    2346            }
    2347
    2348            /* Multiply input by shpiB */
    2349            c = shpiB >>> 32;
    2350            d = shpiB & 0xffffffffL;
    2351            ac = a * c;
    2352            bc = b * c;
    2353            ad = a * d;
    2354
    2355            /* Collect terms */
    2356            ac = ac + ((bc + ad) >>> 32);
    2357
    2358            bita = (prodB & 0x8000000000000000L) != 0;
    2359            bitb = (ac & 0x8000000000000000L ) != 0;
    2360            prodB += ac;
    2361            bitsum = (prodB & 0x8000000000000000L) != 0;
    2362            /* Carry */
    2363            if ( (bita && bitb) ||
    2364                    ((bita || bitb) && !bitsum) ) {
    2365                prodA++;
    2366            }
    2367
    2368            /* Multiply by shpi0 */
    2369            c = shpi0 >>> 32;
    2370            d = shpi0 & 0xffffffffL;
    2371
    2372            bd = b * d;
    2373            bc = b * c;
    2374            ad = a * d;
    2375
    2376            prodA += bd + ((bc + ad) << 32);
    2377
    2378            /*
    2379             * prodA, prodB now contain the remainder as a fraction of PI.  We want this as a fraction of
    2380             * PI/2, so use the following steps:
    2381             * 1.) multiply by 4.
    2382             * 2.) do a fixed point muliply by PI/4.
    2383             * 3.) Convert to floating point.
    2384             * 4.) Multiply by 2
    2385             */
    2386
    2387            /* This identifies the quadrant */
    2388            int intPart = (int)(prodA >>> 62);
    2389
    2390            /* Multiply by 4 */
    2391            prodA <<= 2;
    2392            prodA |= prodB >>> 62;
    2393            prodB <<= 2;
    2394
    2395            /* Multiply by PI/4 */
    2396            a = prodA >>> 32;
    2397            b = prodA & 0xffffffffL;
    2398
    2399            c = PI_O_4_BITS[0] >>> 32;
    2400            d = PI_O_4_BITS[0] & 0xffffffffL;
    2401
    2402            ac = a * c;
    2403            bd = b * d;
    2404            bc = b * c;
    2405            ad = a * d;
    2406
    2407            long prod2B = bd + (ad << 32);
    2408            long prod2A = ac + (ad >>> 32);
    2409
    2410            bita = (bd & 0x8000000000000000L) != 0;
    2411            bitb = (ad & 0x80000000L ) != 0;
    2412            bitsum = (prod2B & 0x8000000000000000L) != 0;
    2413
    2414            /* Carry */
    2415            if ( (bita && bitb) ||
    2416                    ((bita || bitb) && !bitsum) ) {
    2417                prod2A++;
    2418            }
    2419
    2420            bita = (prod2B & 0x8000000000000000L) != 0;
    2421            bitb = (bc & 0x80000000L ) != 0;
    2422
    2423            prod2B = prod2B + (bc << 32);
    2424            prod2A = prod2A + (bc >>> 32);
    2425
    2426            bitsum = (prod2B & 0x8000000000000000L) != 0;
    2427
    2428            /* Carry */
    2429            if ( (bita && bitb) ||
    2430                    ((bita || bitb) && !bitsum) ) {
    2431                prod2A++;
    2432            }
    2433
    2434            /* Multiply input by pio4bits[1] */
    2435            c = PI_O_4_BITS[1] >>> 32;
    2436            d = PI_O_4_BITS[1] & 0xffffffffL;
    2437            ac = a * c;
    2438            bc = b * c;
    2439            ad = a * d;
    2440
    2441            /* Collect terms */
    2442            ac = ac + ((bc + ad) >>> 32);
    2443
    2444            bita = (prod2B & 0x8000000000000000L) != 0;
    2445            bitb = (ac & 0x8000000000000000L ) != 0;
    2446            prod2B += ac;
    2447            bitsum = (prod2B & 0x8000000000000000L) != 0;
    2448            /* Carry */
    2449            if ( (bita && bitb) ||
    2450                    ((bita || bitb) && !bitsum) ) {
    2451                prod2A++;
    2452            }
    2453
    2454            /* Multiply inputB by pio4bits[0] */
    2455            a = prodB >>> 32;
    2456            b = prodB & 0xffffffffL;
    2457            c = PI_O_4_BITS[0] >>> 32;
    2458            d = PI_O_4_BITS[0] & 0xffffffffL;
    2459            ac = a * c;
    2460            bc = b * c;
    2461            ad = a * d;
    2462
    2463            /* Collect terms */
    2464            ac = ac + ((bc + ad) >>> 32);
    2465
    2466            bita = (prod2B & 0x8000000000000000L) != 0;
    2467            bitb = (ac & 0x8000000000000000L ) != 0;
    2468            prod2B += ac;
    2469            bitsum = (prod2B & 0x8000000000000000L) != 0;
    2470            /* Carry */
    2471            if ( (bita && bitb) ||
    2472                    ((bita || bitb) && !bitsum) ) {
    2473                prod2A++;
    2474            }
    2475
    2476            /* Convert to double */
    2477            double tmpA = (prod2A >>> 12) / TWO_POWER_52;  // High order 52 bits
    2478            double tmpB = (((prod2A & 0xfffL) << 40) + (prod2B >>> 24)) / TWO_POWER_52 / TWO_POWER_52; // Low bits
    2479
    2480            double sumA = tmpA + tmpB;
    2481            double sumB = -(sumA - tmpA - tmpB);
    2482
    2483            /* Multiply by PI/2 and return */
    2484            result[0] = intPart;
    2485            result[1] = sumA * 2.0;
    2486            result[2] = sumB * 2.0;
    2487        }
    2488
    2489        /**
    2490         *  Sine function.
    2491         *  @param x a number
    2492         *  @return sin(x)
    2493         */
    2494        public static double sin(double x) {
    2495            boolean negative = false;
    2496            int quadrant = 0;
    2497            double xa;
    2498            double xb = 0.0;
    2499
    2500            /* Take absolute value of the input */
    2501            xa = x;
    2502            if (x < 0) {
    2503                negative = true;
    2504                xa = -xa;
    2505            }
    2506
    2507            /* Check for zero and negative zero */
    2508            if (xa == 0.0) {
    2509                long bits = Double.doubleToLongBits(x);
    2510                if (bits < 0) {
    2511                    return -0.0;
    2512                }
    2513                return 0.0;
    2514            }
    2515
    2516            if (xa != xa || xa == Double.POSITIVE_INFINITY) {
    2517                return Double.NaN;
    2518            }
    2519
    2520            /* Perform any argument reduction */
    2521            if (xa > 3294198.0) {
    2522                // PI * (2**20)
    2523                // Argument too big for CodyWaite reduction.  Must use
    2524                // PayneHanek.
    2525                double reduceResults[] = new double[3];
    2526                reducePayneHanek(xa, reduceResults);
    2527                quadrant = ((int) reduceResults[0]) & 3;
    2528                xa = reduceResults[1];
    2529                xb = reduceResults[2];
    2530            } else if (xa > 1.5707963267948966) {
    2531                /* Inline the Cody/Waite reduction for performance */
    2532
    2533                // Estimate k
    2534                //k = (int)(xa / 1.5707963267948966);
    2535                int k = (int)(xa * 0.6366197723675814);
    2536
    2537                // Compute remainder
    2538                double remA;
    2539                double remB;
    2540                while (true) {
    2541                    double a = -k * 1.570796251296997;
    2542                    remA = xa + a;
    2543                    remB = -(remA - xa - a);
    2544
    2545                    a = -k * 7.549789948768648E-8;
    2546                    double b = remA;
    2547                    remA = a + b;
    2548                    remB += -(remA - b - a);
    2549
    2550                    a = -k * 6.123233995736766E-17;
    2551                    b = remA;
    2552                    remA = a + b;
    2553                    remB += -(remA - b - a);
    2554
    2555                    if (remA > 0.0)
    2556                        break;
    2557
    2558                    // Remainder is negative, so decrement k and try again.
    2559                    // This should only happen if the input is very close
    2560                    // to an even multiple of pi/2
    2561                    k--;
    2562                }
    2563                quadrant = k & 3;
    2564                xa = remA;
    2565                xb = remB;
    2566            }
    2567
    2568            if (negative) {
    2569                quadrant ^= 2;  // Flip bit 1
    2570            }
    2571
    2572            switch (quadrant) {
    2573                case 0:
    2574                    return sinQ(xa, xb);
    2575                case 1:
    2576                    return cosQ(xa, xb);
    2577                case 2:
    2578                    return -sinQ(xa, xb);
    2579                case 3:
    2580                    return -cosQ(xa, xb);
    2581                default:
    2582                    return Double.NaN;
    2583            }
    2584        }
    2585
    2586        /**
    2587         *  Cosine function
    2588         *  @param x a number
    2589         *  @return cos(x)
    2590         */
    2591        public static double cos(double x) {
    2592            int quadrant = 0;
    2593
    2594            /* Take absolute value of the input */
    2595            double xa = x;
    2596            if (x < 0) {
    2597                xa = -xa;
    2598            }
    2599
    2600            if (xa != xa || xa == Double.POSITIVE_INFINITY) {
    2601                return Double.NaN;
    2602            }
    2603
    2604            /* Perform any argument reduction */
    2605            double xb = 0;
    2606            if (xa > 3294198.0) {
    2607                // PI * (2**20)
    2608                // Argument too big for CodyWaite reduction.  Must use
    2609                // PayneHanek.
    2610                double reduceResults[] = new double[3];
    2611                reducePayneHanek(xa, reduceResults);
    2612                quadrant = ((int) reduceResults[0]) & 3;
    2613                xa = reduceResults[1];
    2614                xb = reduceResults[2];
    2615            } else if (xa > 1.5707963267948966) {
    2616                /* Inline the Cody/Waite reduction for performance */
    2617
    2618                // Estimate k
    2619                //k = (int)(xa / 1.5707963267948966);
    2620                int k = (int)(xa * 0.6366197723675814);
    2621
    2622                // Compute remainder
    2623                double remA;
    2624                double remB;
    2625                while (true) {
    2626                    double a = -k * 1.570796251296997;
    2627                    remA = xa + a;
    2628                    remB = -(remA - xa - a);
    2629
    2630                    a = -k * 7.549789948768648E-8;
    2631                    double b = remA;
    2632                    remA = a + b;
    2633                    remB += -(remA - b - a);
    2634
    2635                    a = -k * 6.123233995736766E-17;
    2636                    b = remA;
    2637                    remA = a + b;
    2638                    remB += -(remA - b - a);
    2639
    2640                    if (remA > 0.0)
    2641                        break;
    2642
    2643                    // Remainder is negative, so decrement k and try again.
    2644                    // This should only happen if the input is very close
    2645                    // to an even multiple of pi/2
    2646                    k--;
    2647                }
    2648                quadrant = k & 3;
    2649                xa = remA;
    2650                xb = remB;
    2651            }
    2652
    2653            //if (negative)
    2654            //  quadrant = (quadrant + 2) % 4;
    2655
    2656            switch (quadrant) {
    2657                case 0:
    2658                    return cosQ(xa, xb);
    2659                case 1:
    2660                    return -sinQ(xa, xb);
    2661                case 2:
    2662                    return -cosQ(xa, xb);
    2663                case 3:
    2664                    return sinQ(xa, xb);
    2665                default:
    2666                    return Double.NaN;
    2667            }
    2668        }
    2669
    2670        /**
    2671         *   Tangent function
    2672         *  @param x a number
    2673         *  @return tan(x)
    2674         */
    2675        public static double tan(double x) {
    2676            boolean negative = false;
    2677            int quadrant = 0;
    2678
    2679            /* Take absolute value of the input */
    2680            double xa = x;
    2681            if (x < 0) {
    2682                negative = true;
    2683                xa = -xa;
    2684            }
    2685
    2686            /* Check for zero and negative zero */
    2687            if (xa == 0.0) {
    2688                long bits = Double.doubleToLongBits(x);
    2689                if (bits < 0) {
    2690                    return -0.0;
    2691                }
    2692                return 0.0;
    2693            }
    2694
    2695            if (xa != xa || xa == Double.POSITIVE_INFINITY) {
    2696                return Double.NaN;
    2697            }
    2698
    2699            /* Perform any argument reduction */
    2700            double xb = 0;
    2701            if (xa > 3294198.0) {
    2702                // PI * (2**20)
    2703                // Argument too big for CodyWaite reduction.  Must use
    2704                // PayneHanek.
    2705                double reduceResults[] = new double[3];
    2706                reducePayneHanek(xa, reduceResults);
    2707                quadrant = ((int) reduceResults[0]) & 3;
    2708                xa = reduceResults[1];
    2709                xb = reduceResults[2];
    2710            } else if (xa > 1.5707963267948966) {
    2711                /* Inline the Cody/Waite reduction for performance */
    2712
    2713                // Estimate k
    2714                //k = (int)(xa / 1.5707963267948966);
    2715                int k = (int)(xa * 0.6366197723675814);
    2716
    2717                // Compute remainder
    2718                double remA;
    2719                double remB;
    2720                while (true) {
    2721                    double a = -k * 1.570796251296997;
    2722                    remA = xa + a;
    2723                    remB = -(remA - xa - a);
    2724
    2725                    a = -k * 7.549789948768648E-8;
    2726                    double b = remA;
    2727                    remA = a + b;
    2728                    remB += -(remA - b - a);
    2729
    2730                    a = -k * 6.123233995736766E-17;
    2731                    b = remA;
    2732                    remA = a + b;
    2733                    remB += -(remA - b - a);
    2734
    2735                    if (remA > 0.0)
    2736                        break;
    2737
    2738                    // Remainder is negative, so decrement k and try again.
    2739                    // This should only happen if the input is very close
    2740                    // to an even multiple of pi/2
    2741                    k--;
    2742                }
    2743                quadrant = k & 3;
    2744                xa = remA;
    2745                xb = remB;
    2746            }
    2747
    2748            if (xa > 1.5) {
    2749                // Accurracy suffers between 1.5 and PI/2
    2750                final double pi2a = 1.5707963267948966;
    2751                final double pi2b = 6.123233995736766E-17;
    2752
    2753                final double a = pi2a - xa;
    2754                double b = -(a - pi2a + xa);
    2755                b += pi2b - xb;
    2756
    2757                xa = a + b;
    2758                xb = -(xa - a - b);
    2759                quadrant ^= 1;
    2760                negative ^= true;
    2761            }
    2762
    2763            double result;
    2764            if ((quadrant & 1) == 0) {
    2765                result = tanQ(xa, xb, false);
    2766            } else {
    2767                result = -tanQ(xa, xb, true);
    2768            }
    2769
    2770            if (negative) {
    2771                result = -result;
    2772            }
    2773
    2774            return result;
    2775        }
    2776
    2777        /**
    2778         * Arctangent function
    2779         *  @param x a number
    2780         *  @return atan(x)
    2781         */
    2782        public static double atan(double x) {
    2783            return atan(x, 0.0, false);
    2784        }
    2785
    2786        /** Internal helper function to compute arctangent.
    2787         * @param xa number from which arctangent is requested
    2788         * @param xb extra bits for x (may be 0.0)
    2789         * @param leftPlane if true, result angle must be put in the left half plane
    2790         * @return atan(xa + xb) (or angle shifted by {@code PI} if leftPlane is true)
    2791         */
    2792        private static double atan(double xa, double xb, boolean leftPlane) {
    2793            boolean negate = false;
    2794            int idx;
    2795
    2796            if (xa == 0.0) { // Matches +/- 0.0; return correct sign
    2797                return leftPlane ? copySign(Math.PI, xa) : xa;
    2798            }
    2799
    2800            if (xa < 0) {
    2801                // negative
    2802                xa = -xa;
    2803                xb = -xb;
    2804                negate = true;
    2805            }
    2806
    2807            if (xa > 1.633123935319537E16) { // Very large input
    2808                return (negate ^ leftPlane) ? (-Math.PI/2.0) : (Math.PI/2.0);
    2809            }
    2810
    2811            /* Estimate the closest tabulated arctan value, compute eps = xa-tangentTable */
    2812            if (xa < 1.0) {
    2813                idx = (int) (((-1.7168146928204136 * xa * xa + 8.0) * xa) + 0.5);
    2814            } else {
    2815                double temp = 1.0/xa;
    2816                idx = (int) (-((-1.7168146928204136 * temp * temp + 8.0) * temp) + 13.07);
    2817            }
    2818            double epsA = xa - TANGENT_TABLE_A[idx];
    2819            double epsB = -(epsA - xa + TANGENT_TABLE_A[idx]);
    2820            epsB += xb - TANGENT_TABLE_B[idx];
    2821
    2822            double temp = epsA + epsB;
    2823            epsB = -(temp - epsA - epsB);
    2824            epsA = temp;
    2825
    2826            /* Compute eps = eps / (1.0 + xa*tangent) */
    2827            temp = xa * HEX_40000000;
    2828            double ya = xa + temp - temp;
    2829            double yb = xb + xa - ya;
    2830            xa = ya;
    2831            xb += yb;
    2832
    2833            //if (idx > 8 || idx == 0)
    2834            if (idx == 0) {
    2835                /* If the slope of the arctan is gentle enough (< 0.45), this approximation will suffice */
    2836                //double denom = 1.0 / (1.0 + xa*tangentTableA[idx] + xb*tangentTableA[idx] + xa*tangentTableB[idx] + xb*tangentTableB[idx]);
    2837                double denom = 1.0 / (1.0 + (xa + xb) * (TANGENT_TABLE_A[idx] + TANGENT_TABLE_B[idx]));
    2838                //double denom = 1.0 / (1.0 + xa*tangentTableA[idx]);
    2839                ya = epsA * denom;
    2840                yb = epsB * denom;
    2841            } else {
    2842                double temp2 = xa * TANGENT_TABLE_A[idx];
    2843                double za = 1.0 + temp2;
    2844                double zb = -(za - 1.0 - temp2);
    2845                temp2 = xb * TANGENT_TABLE_A[idx] + xa * TANGENT_TABLE_B[idx];
    2846                temp = za + temp2;
    2847                zb += -(temp - za - temp2);
    2848                za = temp;
    2849
    2850                zb += xb * TANGENT_TABLE_B[idx];
    2851                ya = epsA / za;
    2852
    2853                temp = ya * HEX_40000000;
    2854                final double yaa = (ya + temp) - temp;
    2855                final double yab = ya - yaa;
    2856
    2857                temp = za * HEX_40000000;
    2858                final double zaa = (za + temp) - temp;
    2859                final double zab = za - zaa;
    2860
    2861                /* Correct for rounding in division */
    2862                yb = (epsA - yaa * zaa - yaa * zab - yab * zaa - yab * zab) / za;
    2863
    2864                yb += -epsA * zb / za / za;
    2865                yb += epsB / za;
    2866            }
    2867
    2868
    2869            epsA = ya;
    2870            epsB = yb;
    2871
    2872            /* Evaluate polynomial */
    2873            double epsA2 = epsA*epsA;
    2874
    2875            /*
    2876        yb = -0.09001346640161823;
    2877        yb = yb * epsA2 + 0.11110718400605211;
    2878        yb = yb * epsA2 + -0.1428571349122913;
    2879        yb = yb * epsA2 + 0.19999999999273194;
    2880        yb = yb * epsA2 + -0.33333333333333093;
    2881        yb = yb * epsA2 * epsA;
    2882             */
    2883
    2884            yb = 0.07490822288864472;
    2885            yb = yb * epsA2 + -0.09088450866185192;
    2886            yb = yb * epsA2 + 0.11111095942313305;
    2887            yb = yb * epsA2 + -0.1428571423679182;
    2888            yb = yb * epsA2 + 0.19999999999923582;
    2889            yb = yb * epsA2 + -0.33333333333333287;
    2890            yb = yb * epsA2 * epsA;
    2891
    2892
    2893            ya = epsA;
    2894
    2895            temp = ya + yb;
    2896            yb = -(temp - ya - yb);
    2897            ya = temp;
    2898
    2899            /* Add in effect of epsB.   atan'(x) = 1/(1+x^2) */
    2900            yb += epsB / (1.0 + epsA * epsA);
    2901
    2902            double result;
    2903            double resultb;
    2904
    2905            //result = yb + eighths[idx] + ya;
    2906            double za = EIGHTHS[idx] + ya;
    2907            double zb = -(za - EIGHTHS[idx] - ya);
    2908            temp = za + yb;
    2909            zb += -(temp - za - yb);
    2910            za = temp;
    2911
    2912            result = za + zb;
    2913            resultb = -(result - za - zb);
    2914
    2915            if (leftPlane) {
    2916                // Result is in the left plane
    2917                final double pia = 1.5707963267948966*2.0;
    2918                final double pib = 6.123233995736766E-17*2.0;
    2919
    2920                za = pia - result;
    2921                zb = -(za - pia + result);
    2922                zb += pib - resultb;
    2923
    2924                result = za + zb;
    2925                resultb = -(result - za - zb);
    2926            }
    2927
    2928
    2929            if (negate ^ leftPlane) {
    2930                result = -result;
    2931            }
    2932
    2933            return result;
    2934        }
    2935
    2936        /**
    2937         * Two arguments arctangent function
    2938         * @param y ordinate
    2939         * @param x abscissa
    2940         * @return phase angle of point (x,y) between {@code -PI} and {@code PI}
    2941         */
    2942        public static double atan2(double y, double x) {
    2943            if (x !=x || y != y) {
    2944                return Double.NaN;
    2945            }
    2946
    2947            if (y == 0.0) {
    2948                double result = x*y;
    2949                double invx = 1.0/x;
    2950                double invy = 1.0/y;
    2951
    2952                if (invx == 0.0) { // X is infinite
    2953                    if (x > 0) {
    2954                        return y; // return +/- 0.0
    2955                    } else {
    2956                        return copySign(Math.PI, y);
    2957                    }
    2958                }
    2959
    2960                if (x < 0.0 || invx < 0.0) {
    2961                    if (y < 0.0 || invy < 0.0) {
    2962                        return -Math.PI;
    2963                    } else {
    2964                        return Math.PI;
    2965                    }
    2966                } else {
    2967                    return result;
    2968                }
    2969            }
    2970
    2971            // y cannot now be zero
    2972
    2973            if (y == Double.POSITIVE_INFINITY) {
    2974                if (x == Double.POSITIVE_INFINITY) {
    2975                    return Math.PI/4.0;
    2976                }
    2977
    2978                if (x == Double.NEGATIVE_INFINITY) {
    2979                    return Math.PI*3.0/4.0;
    2980                }
    2981
    2982                return Math.PI/2.0;
    2983            }
    2984
    2985            if (y == Double.NEGATIVE_INFINITY) {
    2986                if (x == Double.POSITIVE_INFINITY) {
    2987                    return -Math.PI/4.0;
    2988                }
    2989
    2990                if (x == Double.NEGATIVE_INFINITY) {
    2991                    return -Math.PI*3.0/4.0;
    2992                }
    2993
    2994                return -Math.PI/2.0;
    2995            }
    2996
    2997            if (x == Double.POSITIVE_INFINITY) {
    2998                if (y > 0.0 || 1/y > 0.0) {
    2999                    return 0.0;
    3000                }
    3001
    3002                if (y < 0.0 || 1/y < 0.0) {
    3003                    return -0.0;
    3004                }
    3005            }
    3006
    3007            if (x == Double.NEGATIVE_INFINITY)
    3008            {
    3009                if (y > 0.0 || 1/y > 0.0) {
    3010                    return Math.PI;
    3011                }
    3012
    3013                if (y < 0.0 || 1/y < 0.0) {
    3014                    return -Math.PI;
    3015                }
    3016            }
    3017
    3018            // Neither y nor x can be infinite or NAN here
    3019
    3020            if (x == 0) {
    3021                if (y > 0.0 || 1/y > 0.0) {
    3022                    return Math.PI/2.0;
    3023                }
    3024
    3025                if (y < 0.0 || 1/y < 0.0) {
    3026                    return -Math.PI/2.0;
    3027                }
    3028            }
    3029
    3030            // Compute ratio r = y/x
    3031            final double r = y/x;
    3032            if (Double.isInfinite(r)) { // bypass calculations that can create NaN
    3033                return atan(r, 0, x < 0);
    3034            }
    3035
    3036            double ra = doubleHighPart(r);
    3037            double rb = r - ra;
    3038
    3039            // Split x
    3040            final double xa = doubleHighPart(x);
    3041            final double xb = x - xa;
    3042
    3043            rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x;
    3044
    3045            double temp = ra + rb;
    3046            rb = -(temp - ra - rb);
    3047            ra = temp;
    3048
    3049            if (ra == 0) { // Fix up the sign so atan works correctly
    3050                ra = copySign(0.0, y);
    3051            }
    3052
    3053            // Call atan
    3054            double result = atan(ra, rb, x < 0);
    3055
    3056            return result;
    3057        }
    3058
    3059        /** Compute the arc sine of a number.
    3060         * @param x number on which evaluation is done
    3061         * @return arc sine of x
    3062         */
    3063        public static double asin(double x) {
    3064          if (x != x) {
    3065              return Double.NaN;
    3066          }
    3067
    3068          if (x > 1.0 || x < -1.0) {
    3069              return Double.NaN;
    3070          }
    3071
    3072          if (x == 1.0) {
    3073              return Math.PI/2.0;
    3074          }
    3075
    3076          if (x == -1.0) {
    3077              return -Math.PI/2.0;
    3078          }
    3079
    3080          if (x == 0.0) { // Matches +/- 0.0; return correct sign
    3081              return x;
    3082          }
    3083
    3084          /* Compute asin(x) = atan(x/sqrt(1-x*x)) */
    3085
    3086          /* Split x */
    3087          double temp = x * HEX_40000000;
    3088          final double xa = x + temp - temp;
    3089          final double xb = x - xa;
    3090
    3091          /* Square it */
    3092          double ya = xa*xa;
    3093          double yb = xa*xb*2.0 + xb*xb;
    3094
    3095          /* Subtract from 1 */
    3096          ya = -ya;
    3097          yb = -yb;
    3098
    3099          double za = 1.0 + ya;
    3100          double zb = -(za - 1.0 - ya);
    3101
    3102          temp = za + yb;
    3103          zb += -(temp - za - yb);
    3104          za = temp;
    3105
    3106          /* Square root */
    3107          double y;
    3108          y = sqrt(za);
    3109          temp = y * HEX_40000000;
    3110          ya = y + temp - temp;
    3111          yb = y - ya;
    3112
    3113          /* Extend precision of sqrt */
    3114          yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y);
    3115
    3116          /* Contribution of zb to sqrt */
    3117          double dx = zb / (2.0*y);
    3118
    3119          // Compute ratio r = x/y
    3120          double r = x/y;
    3121          temp = r * HEX_40000000;
    3122          double ra = r + temp - temp;
    3123          double rb = r - ra;
    3124
    3125          rb += (x - ra*ya - ra*yb - rb*ya - rb*yb) / y;  // Correct for rounding in division
    3126          rb += -x * dx / y / y;  // Add in effect additional bits of sqrt.
    3127
    3128          temp = ra + rb;
    3129          rb = -(temp - ra - rb);
    3130          ra = temp;
    3131
    3132          return atan(ra, rb, false);
    3133        }
    3134
    3135        /** Compute the arc cosine of a number.
    3136         * @param x number on which evaluation is done
    3137         * @return arc cosine of x
    3138         */
    3139        public static double acos(double x) {
    3140          if (x != x) {
    3141              return Double.NaN;
    3142          }
    3143
    3144          if (x > 1.0 || x < -1.0) {
    3145              return Double.NaN;
    3146          }
    3147
    3148          if (x == -1.0) {
    3149              return Math.PI;
    3150          }
    3151
    3152          if (x == 1.0) {
    3153              return 0.0;
    3154          }
    3155
    3156          if (x == 0) {
    3157              return Math.PI/2.0;
    3158          }
    3159
    3160          /* Compute acos(x) = atan(sqrt(1-x*x)/x) */
    3161
    3162          /* Split x */
    3163          double temp = x * HEX_40000000;
    3164          final double xa = x + temp - temp;
    3165          final double xb = x - xa;
    3166
    3167          /* Square it */
    3168          double ya = xa*xa;
    3169          double yb = xa*xb*2.0 + xb*xb;
    3170
    3171          /* Subtract from 1 */
    3172          ya = -ya;
    3173          yb = -yb;
    3174
    3175          double za = 1.0 + ya;
    3176          double zb = -(za - 1.0 - ya);
    3177
    3178          temp = za + yb;
    3179          zb += -(temp - za - yb);
    3180          za = temp;
    3181
    3182          /* Square root */
    3183          double y = sqrt(za);
    3184          temp = y * HEX_40000000;
    3185          ya = y + temp - temp;
    3186          yb = y - ya;
    3187
    3188          /* Extend precision of sqrt */
    3189          yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y);
    3190
    3191          /* Contribution of zb to sqrt */
    3192          yb += zb / (2.0*y);
    3193          y = ya+yb;
    3194          yb = -(y - ya - yb);
    3195
    3196          // Compute ratio r = y/x
    3197          double r = y/x;
    3198
    3199          // Did r overflow?
    3200          if (Double.isInfinite(r)) { // x is effectively zero
    3201              return Math.PI/2; // so return the appropriate value
    3202          }
    3203
    3204          double ra = doubleHighPart(r);
    3205          double rb = r - ra;
    3206
    3207          rb += (y - ra*xa - ra*xb - rb*xa - rb*xb) / x;  // Correct for rounding in division
    3208          rb += yb / x;  // Add in effect additional bits of sqrt.
    3209
    3210          temp = ra + rb;
    3211          rb = -(temp - ra - rb);
    3212          ra = temp;
    3213
    3214          return atan(ra, rb, x<0);
    3215        }
    3216
    3217        /** Compute the cubic root of a number.
    3218         * @param x number on which evaluation is done
    3219         * @return cubic root of x
    3220         */
    3221        public static double cbrt(double x) {
    3222          /* Convert input double to bits */
    3223          long inbits = Double.doubleToLongBits(x);
    3224          int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
    3225          boolean subnormal = false;
    3226
    3227          if (exponent == -1023) {
    3228              if (x == 0) {
    3229                  return x;
    3230              }
    3231
    3232              /* Subnormal, so normalize */
    3233              subnormal = true;
    3234              x *= 1.8014398509481984E16;  // 2^54
    3235              inbits = Double.doubleToLongBits(x);
    3236              exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
    3237          }
    3238
    3239          if (exponent == 1024) {
    3240              // Nan or infinity.  Don't care which.
    3241              return x;
    3242          }
    3243
    3244          /* Divide the exponent by 3 */
    3245          int exp3 = exponent / 3;
    3246
    3247          /* p2 will be the nearest power of 2 to x with its exponent divided by 3 */
    3248          double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) |
    3249                                              (long)(((exp3 + 1023) & 0x7ff)) << 52);
    3250
    3251          /* This will be a number between 1 and 2 */
    3252          final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L);
    3253
    3254          /* Estimate the cube root of mant by polynomial */
    3255          double est = -0.010714690733195933;
    3256          est = est * mant + 0.0875862700108075;
    3257          est = est * mant + -0.3058015757857271;
    3258          est = est * mant + 0.7249995199969751;
    3259          est = est * mant + 0.5039018405998233;
    3260
    3261          est *= CBRTTWO[exponent % 3 + 2];
    3262
    3263          // est should now be good to about 15 bits of precision.   Do 2 rounds of
    3264          // Newton's method to get closer,  this should get us full double precision
    3265          // Scale down x for the purpose of doing newtons method.  This avoids over/under flows.
    3266          final double xs = x / (p2*p2*p2);
    3267          est += (xs - est*est*est) / (3*est*est);
    3268          est += (xs - est*est*est) / (3*est*est);
    3269
    3270          // Do one round of Newton's method in extended precision to get the last bit right.
    3271          double temp = est * HEX_40000000;
    3272          double ya = est + temp - temp;
    3273          double yb = est - ya;
    3274
    3275          double za = ya * ya;
    3276          double zb = ya * yb * 2.0 + yb * yb;
    3277          temp = za * HEX_40000000;
    3278          double temp2 = za + temp - temp;
    3279          zb += za - temp2;
    3280          za = temp2;
    3281
    3282          zb = za * yb + ya * zb + zb * yb;
    3283          za = za * ya;
    3284
    3285          double na = xs - za;
    3286          double nb = -(na - xs + za);
    3287          nb -= zb;
    3288
    3289          est += (na+nb)/(3*est*est);
    3290
    3291          /* Scale by a power of two, so this is exact. */
    3292          est *= p2;
    3293
    3294          if (subnormal) {
    3295              est *= 3.814697265625E-6;  // 2^-18
    3296          }
    3297
    3298          return est;
    3299        }
    3300
    3301        /**
    3302         *  Convert degrees to radians, with error of less than 0.5 ULP
    3303         *  @param x angle in degrees
    3304         *  @return x converted into radians
    3305         */
    3306        public static double toRadians(double x)
    3307        {
    3308            if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
    3309                return x;
    3310            }
    3311
    3312            // These are PI/180 split into high and low order bits
    3313            final double facta = 0.01745329052209854;
    3314            final double factb = 1.997844754509471E-9;
    3315
    3316            double xa = doubleHighPart(x);
    3317            double xb = x - xa;
    3318
    3319            double result = xb * factb + xb * facta + xa * factb + xa * facta;
    3320            if (result == 0) {
    3321                result = result * x; // ensure correct sign if calculation underflows
    3322            }
    3323            return result;
    3324        }
    3325
    3326        /**
    3327         *  Convert radians to degrees, with error of less than 0.5 ULP
    3328         *  @param x angle in radians
    3329         *  @return x converted into degrees
    3330         */
    3331        public static double toDegrees(double x)
    3332        {
    3333            if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
    3334                return x;
    3335            }
    3336
    3337            // These are 180/PI split into high and low order bits
    3338            final double facta = 57.2957763671875;
    3339            final double factb = 3.145894820876798E-6;
    3340
    3341            double xa = doubleHighPart(x);
    3342            double xb = x - xa;
    3343
    3344            return xb * factb + xb * facta + xa * factb + xa * facta;
    3345        }
    3346
    3347        /**
    3348         * Absolute value.
    3349         * @param x number from which absolute value is requested
    3350         * @return abs(x)
    3351         */
    3352        public static int abs(final int x) {
    3353            return (x < 0) ? -x : x;
    3354        }
    3355
    3356        /**
    3357         * Absolute value.
    3358         * @param x number from which absolute value is requested
    3359         * @return abs(x)
    3360         */
    3361        public static long abs(final long x) {
    3362            return (x < 0l) ? -x : x;
    3363        }
    3364
    3365        /**
    3366         * Absolute value.
    3367         * @param x number from which absolute value is requested
    3368         * @return abs(x)
    3369         */
    3370        public static float abs(final float x) {
    3371            return (x < 0.0f) ? -x : (x == 0.0f) ? 0.0f : x; // -0.0 => +0.0
    3372        }
    3373
    3374        /**
    3375         * Absolute value.
    3376         * @param x number from which absolute value is requested
    3377         * @return abs(x)
    3378         */
    3379        public static double abs(double x) {
    3380            return (x < 0.0) ? -x : (x == 0.0) ? 0.0 : x; // -0.0 => +0.0
    3381        }
    3382
    3383        /**
    3384         * Compute least significant bit (Unit in Last Position) for a number.
    3385         * @param x number from which ulp is requested
    3386         * @return ulp(x)
    3387         */
    3388        public static double ulp(double x) {
    3389            if (Double.isInfinite(x)) {
    3390                return Double.POSITIVE_INFINITY;
    3391            }
    3392            return abs(x - Double.longBitsToDouble(Double.doubleToLongBits(x) ^ 1));
    3393        }
    3394
    3395        /**
    3396         * Compute least significant bit (Unit in Last Position) for a number.
    3397         * @param x number from which ulp is requested
    3398         * @return ulp(x)
    3399         */
    3400        public static float ulp(float x) {
    3401            if (Float.isInfinite(x)) {
    3402                return Float.POSITIVE_INFINITY;
    3403            }
    3404            return abs(x - Float.intBitsToFloat(Float.floatToIntBits(x) ^ 1));
    3405        }
    3406
    3407        /**
    3408         * Multiply a double number by a power of 2.
    3409         * @param d number to multiply
    3410         * @param n power of 2
    3411         * @return d &times; 2<sup>n</sup>
    3412         */
    3413        public static double scalb(final double d, final int n) {
    3414
    3415            // first simple and fast handling when 2^n can be represented using normal numbers
    3416            if ((n > -1023) && (n < 1024)) {
    3417                return d * Double.longBitsToDouble(((long) (n + 1023)) << 52);
    3418            }
    3419
    3420            // handle special cases
    3421            if (Double.isNaN(d) || Double.isInfinite(d) || (d == 0)) {
    3422                return d;
    3423            }
    3424            if (n < -2098) {
    3425                return (d > 0) ? 0.0 : -0.0;
    3426            }
    3427            if (n > 2097) {
    3428                return (d > 0) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
    3429            }
    3430
    3431            // decompose d
    3432            final long bits = Double.doubleToLongBits(d);
    3433            final long sign = bits & 0x8000000000000000L;
    3434            int  exponent   = ((int) (bits >>> 52)) & 0x7ff;
    3435            long mantissa   = bits & 0x000fffffffffffffL;
    3436
    3437            // compute scaled exponent
    3438            int scaledExponent = exponent + n;
    3439
    3440            if (n < 0) {
    3441                // we are really in the case n <= -1023
    3442                if (scaledExponent > 0) {
    3443                    // both the input and the result are normal numbers, we only adjust the exponent
    3444                    return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
    3445                } else if (scaledExponent > -53) {
    3446                    // the input is a normal number and the result is a subnormal number
    3447
    3448                    // recover the hidden mantissa bit
    3449                    mantissa = mantissa | (1L << 52);
    3450
    3451                    // scales down complete mantissa, hence losing least significant bits
    3452                    final long mostSignificantLostBit = mantissa & (1L << (-scaledExponent));
    3453                    mantissa = mantissa >>> (1 - scaledExponent);
    3454                    if (mostSignificantLostBit != 0) {
    3455                        // we need to add 1 bit to round up the result
    3456                        mantissa++;
    3457                    }
    3458                    return Double.longBitsToDouble(sign | mantissa);
    3459
    3460                } else {
    3461                    // no need to compute the mantissa, the number scales down to 0
    3462                    return (sign == 0L) ? 0.0 : -0.0;
    3463                }
    3464            } else {
    3465                // we are really in the case n >= 1024
    3466                if (exponent == 0) {
    3467
    3468                    // the input number is subnormal, normalize it
    3469                    while ((mantissa >>> 52) != 1) {
    3470                        mantissa = mantissa << 1;
    3471                        --scaledExponent;
    3472                    }
    3473                    ++scaledExponent;
    3474                    mantissa = mantissa & 0x000fffffffffffffL;
    3475
    3476                    if (scaledExponent < 2047) {
    3477                        return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
    3478                    } else {
    3479                        return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
    3480                    }
    3481
    3482                } else if (scaledExponent < 2047) {
    3483                    return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
    3484                } else {
    3485                    return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
    3486                }
    3487            }
    3488
    3489        }
    3490
    3491        /**
    3492         * Multiply a float number by a power of 2.
    3493         * @param f number to multiply
    3494         * @param n power of 2
    3495         * @return f &times; 2<sup>n</sup>
    3496         */
    3497        public static float scalb(final float f, final int n) {
    3498
    3499            // first simple and fast handling when 2^n can be represented using normal numbers
    3500            if ((n > -127) && (n < 128)) {
    3501                return f * Float.intBitsToFloat((n + 127) << 23);
    3502            }
    3503
    3504            // handle special cases
    3505            if (Float.isNaN(f) || Float.isInfinite(f) || (f == 0f)) {
    3506                return f;
    3507            }
    3508            if (n < -277) {
    3509                return (f > 0) ? 0.0f : -0.0f;
    3510            }
    3511            if (n > 276) {
    3512                return (f > 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
    3513            }
    3514
    3515            // decompose f
    3516            final int bits = Float.floatToIntBits(f);
    3517            final int sign = bits & 0x80000000;
    3518            int  exponent  = (bits >>> 23) & 0xff;
    3519            int mantissa   = bits & 0x007fffff;
    3520
    3521            // compute scaled exponent
    3522            int scaledExponent = exponent + n;
    3523
    3524            if (n < 0) {
    3525                // we are really in the case n <= -127
    3526                if (scaledExponent > 0) {
    3527                    // both the input and the result are normal numbers, we only adjust the exponent
    3528                    return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
    3529                } else if (scaledExponent > -24) {
    3530                    // the input is a normal number and the result is a subnormal number
    3531
    3532                    // recover the hidden mantissa bit
    3533                    mantissa = mantissa | (1 << 23);
    3534
    3535                    // scales down complete mantissa, hence losing least significant bits
    3536                    final int mostSignificantLostBit = mantissa & (1 << (-scaledExponent));
    3537                    mantissa = mantissa >>> (1 - scaledExponent);
    3538                    if (mostSignificantLostBit != 0) {
    3539                        // we need to add 1 bit to round up the result
    3540                        mantissa++;
    3541                    }
    3542                    return Float.intBitsToFloat(sign | mantissa);
    3543
    3544                } else {
    3545                    // no need to compute the mantissa, the number scales down to 0
    3546                    return (sign == 0) ? 0.0f : -0.0f;
    3547                }
    3548            } else {
    3549                // we are really in the case n >= 128
    3550                if (exponent == 0) {
    3551
    3552                    // the input number is subnormal, normalize it
    3553                    while ((mantissa >>> 23) != 1) {
    3554                        mantissa = mantissa << 1;
    3555                        --scaledExponent;
    3556                    }
    3557                    ++scaledExponent;
    3558                    mantissa = mantissa & 0x007fffff;
    3559
    3560                    if (scaledExponent < 255) {
    3561                        return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
    3562                    } else {
    3563                        return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
    3564                    }
    3565
    3566                } else if (scaledExponent < 255) {
    3567                    return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
    3568                } else {
    3569                    return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
    3570                }
    3571            }
    3572
    3573        }
    3574
    3575        /**
    3576         * Get the next machine representable number after a number, moving
    3577         * in the direction of another number.
    3578         * <p>
    3579         * The ordering is as follows (increasing):
    3580         * <ul>
    3581         * <li>-INFINITY</li>
    3582         * <li>-MAX_VALUE</li>
    3583         * <li>-MIN_VALUE</li>
    3584         * <li>-0.0</li>
    3585         * <li>+0.0</li>
    3586         * <li>+MIN_VALUE</li>
    3587         * <li>+MAX_VALUE</li>
    3588         * <li>+INFINITY</li>
    3589         * <li></li>
    3590         * <p>
    3591         * If arguments compare equal, then the second argument is returned.
    3592         * <p>
    3593         * If {@code direction} is greater than {@code d},
    3594         * the smallest machine representable number strictly greater than
    3595         * {@code d} is returned; if less, then the largest representable number
    3596         * strictly less than {@code d} is returned.</p>
    3597         * <p>
    3598         * If {@code d} is infinite and direction does not
    3599         * bring it back to finite numbers, it is returned unchanged.</p>
    3600         *
    3601         * @param d base number
    3602         * @param direction (the only important thing is whether
    3603         * {@code direction} is greater or smaller than {@code d})
    3604         * @return the next machine representable number in the specified direction
    3605         */
    3606        public static double nextAfter(double d, double direction) {
    3607
    3608            // handling of some important special cases
    3609            if (Double.isNaN(d) || Double.isNaN(direction)) {
    3610                return Double.NaN;
    3611            } else if (d == direction) {
    3612                return direction;
    3613            } else if (Double.isInfinite(d)) {
    3614                return (d < 0) ? -Double.MAX_VALUE : Double.MAX_VALUE;
    3615            } else if (d == 0) {
    3616                return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
    3617            }
    3618            // special cases MAX_VALUE to infinity and  MIN_VALUE to 0
    3619            // are handled just as normal numbers
    3620
    3621            final long bits = Double.doubleToLongBits(d);
    3622            final long sign = bits & 0x8000000000000000L;
    3623            if ((direction < d) ^ (sign == 0L)) {
    3624                return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) + 1));
    3625            } else {
    3626                return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) - 1));
    3627            }
    3628
    3629        }
    3630
    3631        /**
    3632         * Get the next machine representable number after a number, moving
    3633         * in the direction of another number.
    3634         * <p>
    3635         * The ordering is as follows (increasing):
    3636         * <ul>
    3637         * <li>-INFINITY</li>
    3638         * <li>-MAX_VALUE</li>
    3639         * <li>-MIN_VALUE</li>
    3640         * <li>-0.0</li>
    3641         * <li>+0.0</li>
    3642         * <li>+MIN_VALUE</li>
    3643         * <li>+MAX_VALUE</li>
    3644         * <li>+INFINITY</li>
    3645         * <li></li>
    3646         * <p>
    3647         * If arguments compare equal, then the second argument is returned.
    3648         * <p>
    3649         * If {@code direction} is greater than {@code f},
    3650         * the smallest machine representable number strictly greater than
    3651         * {@code f} is returned; if less, then the largest representable number
    3652         * strictly less than {@code f} is returned.</p>
    3653         * <p>
    3654         * If {@code f} is infinite and direction does not
    3655         * bring it back to finite numbers, it is returned unchanged.</p>
    3656         *
    3657         * @param f base number
    3658         * @param direction (the only important thing is whether
    3659         * {@code direction} is greater or smaller than {@code f})
    3660         * @return the next machine representable number in the specified direction
    3661         */
    3662        public static float nextAfter(final float f, final double direction) {
    3663
    3664            // handling of some important special cases
    3665            if (Double.isNaN(f) || Double.isNaN(direction)) {
    3666                return Float.NaN;
    3667            } else if (f == direction) {
    3668                return (float) direction;
    3669            } else if (Float.isInfinite(f)) {
    3670                return (f < 0f) ? -Float.MAX_VALUE : Float.MAX_VALUE;
    3671            } else if (f == 0f) {
    3672                return (direction < 0) ? -Float.MIN_VALUE : Float.MIN_VALUE;
    3673            }
    3674            // special cases MAX_VALUE to infinity and  MIN_VALUE to 0
    3675            // are handled just as normal numbers
    3676
    3677            final int bits = Float.floatToIntBits(f);
    3678            final int sign = bits & 0x80000000;
    3679            if ((direction < f) ^ (sign == 0)) {
    3680                return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) + 1));
    3681            } else {
    3682                return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) - 1));
    3683            }
    3684
    3685        }
    3686
    3687        /** Get the largest whole number smaller than x.
    3688         * @param x number from which floor is requested
    3689         * @return a double number f such that f is an integer f <= x < f + 1.0
    3690         */
    3691        public static double floor(double x) {
    3692            long y;
    3693
    3694            if (x != x) { // NaN
    3695                return x;
    3696            }
    3697
    3698            if (x >= TWO_POWER_52 || x <= -TWO_POWER_52) {
    3699                return x;
    3700            }
    3701
    3702            y = (long) x;
    3703            if (x < 0 && y != x) {
    3704                y--;
    3705            }
    3706
    3707            if (y == 0) {
    3708                return x*y;
    3709            }
    3710
    3711            return y;
    3712        }
    3713
    3714        /** Get the smallest whole number larger than x.
    3715         * @param x number from which ceil is requested
    3716         * @return a double number c such that c is an integer c - 1.0 < x <= c
    3717         */
    3718        public static double ceil(double x) {
    3719            double y;
    3720
    3721            if (x != x) { // NaN
    3722                return x;
    3723            }
    3724
    3725            y = floor(x);
    3726            if (y == x) {
    3727                return y;
    3728            }
    3729
    3730            y += 1.0;
    3731
    3732            if (y == 0) {
    3733                return x*y;
    3734            }
    3735
    3736            return y;
    3737        }
    3738
    3739        /** Get the whole number that is the nearest to x, or the even one if x is exactly half way between two integers.
    3740         * @param x number from which nearest whole number is requested
    3741         * @return a double number r such that r is an integer r - 0.5 <= x <= r + 0.5
    3742         */
    3743        public static double rint(double x) {
    3744            double y = floor(x);
    3745            double d = x - y;
    3746
    3747            if (d > 0.5) {
    3748                if (y == -1.0) {
    3749                    return -0.0; // Preserve sign of operand
    3750                }
    3751                return y+1.0;
    3752            }
    3753            if (d < 0.5) {
    3754                return y;
    3755            }
    3756
    3757            /* half way, round to even */
    3758            long z = (long) y;
    3759            return (z & 1) == 0 ? y : y + 1.0;
    3760        }
    3761
    3762        /** Get the closest long to x.
    3763         * @param x number from which closest long is requested
    3764         * @return closest long to x
    3765         */
    3766        public static long round(double x) {
    3767            return (long) floor(x + 0.5);
    3768        }
    3769
    3770        /** Get the closest int to x.
    3771         * @param x number from which closest int is requested
    3772         * @return closest int to x
    3773         */
    3774        public static int round(final float x) {
    3775            return (int) floor(x + 0.5f);
    3776        }
    3777
    3778        /** Compute the minimum of two values
    3779         * @param a first value
    3780         * @param b second value
    3781         * @return a if a is lesser or equal to b, b otherwise
    3782         */
    3783        public static int min(final int a, final int b) {
    3784            return (a <= b) ? a : b;
    3785        }
    3786
    3787        /** Compute the minimum of two values
    3788         * @param a first value
    3789         * @param b second value
    3790         * @return a if a is lesser or equal to b, b otherwise
    3791         */
    3792        public static long min(final long a, final long b) {
    3793            return (a <= b) ? a : b;
    3794        }
    3795
    3796        /** Compute the minimum of two values
    3797         * @param a first value
    3798         * @param b second value
    3799         * @return a if a is lesser or equal to b, b otherwise
    3800         */
    3801        public static float min(final float a, final float b) {
    3802            if (a > b) {
    3803                return b;
    3804            }
    3805            if (a < b) {
    3806                return a;
    3807            }
    3808            /* if either arg is NaN, return NaN */
    3809            if (a != b) {
    3810                return Float.NaN;
    3811            }
    3812            /* min(+0.0,-0.0) == -0.0 */
    3813            /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
    3814            int bits = Float.floatToRawIntBits(a);
    3815            if (bits == 0x80000000) {
    3816                return a;
    3817            }
    3818            return b;
    3819        }
    3820
    3821        /** Compute the minimum of two values
    3822         * @param a first value
    3823         * @param b second value
    3824         * @return a if a is lesser or equal to b, b otherwise
    3825         */
    3826        public static double min(final double a, final double b) {
    3827            if (a > b) {
    3828                return b;
    3829            }
    3830            if (a < b) {
    3831                return a;
    3832            }
    3833            /* if either arg is NaN, return NaN */
    3834            if (a != b) {
    3835                return Double.NaN;
    3836            }
    3837            /* min(+0.0,-0.0) == -0.0 */
    3838            /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
    3839            long bits = Double.doubleToRawLongBits(a);
    3840            if (bits == 0x8000000000000000L) {
    3841                return a;
    3842            }
    3843            return b;
    3844        }
    3845
    3846        /** Compute the maximum of two values
    3847         * @param a first value
    3848         * @param b second value
    3849         * @return b if a is lesser or equal to b, a otherwise
    3850         */
    3851        public static int max(final int a, final int b) {
    3852            return (a <= b) ? b : a;
    3853        }
    3854
    3855        /** Compute the maximum of two values
    3856         * @param a first value
    3857         * @param b second value
    3858         * @return b if a is lesser or equal to b, a otherwise
    3859         */
    3860        public static long max(final long a, final long b) {
    3861            return (a <= b) ? b : a;
    3862        }
    3863
    3864        /** Compute the maximum of two values
    3865         * @param a first value
    3866         * @param b second value
    3867         * @return b if a is lesser or equal to b, a otherwise
    3868         */
    3869        public static float max(final float a, final float b) {
    3870            if (a > b) {
    3871                return a;
    3872            }
    3873            if (a < b) {
    3874                return b;
    3875            }
    3876            /* if either arg is NaN, return NaN */
    3877            if (a != b) {
    3878                return Float.NaN;
    3879            }
    3880            /* min(+0.0,-0.0) == -0.0 */
    3881            /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
    3882            int bits = Float.floatToRawIntBits(a);
    3883            if (bits == 0x80000000) {
    3884                return b;
    3885            }
    3886            return a;
    3887        }
    3888
    3889        /** Compute the maximum of two values
    3890         * @param a first value
    3891         * @param b second value
    3892         * @return b if a is lesser or equal to b, a otherwise
    3893         */
    3894        public static double max(final double a, final double b) {
    3895            if (a > b) {
    3896                return a;
    3897            }
    3898            if (a < b) {
    3899                return b;
    3900            }
    3901            /* if either arg is NaN, return NaN */
    3902            if (a != b) {
    3903                return Double.NaN;
    3904            }
    3905            /* min(+0.0,-0.0) == -0.0 */
    3906            /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
    3907            long bits = Double.doubleToRawLongBits(a);
    3908            if (bits == 0x8000000000000000L) {
    3909                return b;
    3910            }
    3911            return a;
    3912        }
    3913
    3914        /**
    3915         * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
    3916         * - sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)<br/>
    3917         * avoiding intermediate overflow or underflow.
    3918         *
    3919         * <ul>
    3920         * <li> If either argument is infinite, then the result is positive infinity.</li>
    3921         * <li> else, if either argument is NaN then the result is NaN.</li>
    3922         * </ul>
    3923         *
    3924         * @param x a value
    3925         * @param y a value
    3926         * @return sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
    3927         */
    3928        public static double hypot(final double x, final double y) {
    3929            if (Double.isInfinite(x) || Double.isInfinite(y)) {
    3930                return Double.POSITIVE_INFINITY;
    3931            } else if (Double.isNaN(x) || Double.isNaN(y)) {
    3932                return Double.NaN;
    3933            } else {
    3934
    3935                final int expX = getExponent(x);
    3936                final int expY = getExponent(y);
    3937                if (expX > expY + 27) {
    3938                    // y is neglectible with respect to x
    3939                    return abs(x);
    3940                } else if (expY > expX + 27) {
    3941                    // x is neglectible with respect to y
    3942                    return abs(y);
    3943                } else {
    3944
    3945                    // find an intermediate scale to avoid both overflow and underflow
    3946                    final int middleExp = (expX + expY) / 2;
    3947
    3948                    // scale parameters without losing precision
    3949                    final double scaledX = scalb(x, -middleExp);
    3950                    final double scaledY = scalb(y, -middleExp);
    3951
    3952                    // compute scaled hypotenuse
    3953                    final double scaledH = sqrt(scaledX * scaledX + scaledY * scaledY);
    3954
    3955                    // remove scaling
    3956                    return scalb(scaledH, middleExp);
    3957
    3958                }
    3959
    3960            }
    3961        }
    3962
    3963        /**
    3964         * Computes the remainder as prescribed by the IEEE 754 standard.
    3965         * The remainder value is mathematically equal to {@code x - y*n}
    3966         * where {@code n} is the mathematical integer closest to the exact mathematical value
    3967         * of the quotient {@code x/y}.
    3968         * If two mathematical integers are equally close to {@code x/y} then
    3969         * {@code n} is the integer that is even.
    3970         * <p>
    3971         * <ul>
    3972         * <li>If either operand is NaN, the result is NaN.</li>
    3973         * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
    3974         * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
    3975         * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
    3976         * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
    3977         * </ul>
    3978         * <p><b>Note:</b> this implementation currently delegates to {@link StrictMath#IEEEremainder}
    3979         * @param dividend the number to be divided
    3980         * @param divisor the number by which to divide
    3981         * @return the remainder, rounded
    3982         */
    3983        public static double IEEEremainder(double dividend, double divisor) {
    3984            return StrictMath.IEEEremainder(dividend, divisor); // TODO provide our own implementation
    3985        }
    3986
    3987        /**
    3988         * Returns the first argument with the sign of the second argument.
    3989         * A NaN {@code sign} argument is treated as positive.
    3990         *
    3991         * @param magnitude the value to return
    3992         * @param sign the sign for the returned value
    3993         * @return the magnitude with the same sign as the {@code sign} argument
    3994         */
    3995        public static double copySign(double magnitude, double sign){
    3996            long m = Double.doubleToLongBits(magnitude);
    3997            long s = Double.doubleToLongBits(sign);
    3998            if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
    3999                return magnitude;
    4000            }
    4001            return -magnitude; // flip sign
    4002        }
    4003
    4004        /**
    4005         * Returns the first argument with the sign of the second argument.
    4006         * A NaN {@code sign} argument is treated as positive.
    4007         *
    4008         * @param magnitude the value to return
    4009         * @param sign the sign for the returned value
    4010         * @return the magnitude with the same sign as the {@code sign} argument
    4011         */
    4012        public static float copySign(float magnitude, float sign){
    4013            int m = Float.floatToIntBits(magnitude);
    4014            int s = Float.floatToIntBits(sign);
    4015            if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
    4016                return magnitude;
    4017            }
    4018            return -magnitude; // flip sign
    4019        }
    4020
    4021        /**
    4022         * Return the exponent of a double number, removing the bias.
    4023         * <p>
    4024         * For double numbers of the form 2<sup>x</sup>, the unbiased
    4025         * exponent is exactly x.
    4026         * </p>
    4027         * @param d number from which exponent is requested
    4028         * @return exponent for d in IEEE754 representation, without bias
    4029         */
    4030        public static int getExponent(final double d) {
    4031            return (int) ((Double.doubleToLongBits(d) >>> 52) & 0x7ff) - 1023;
    4032        }
    4033
    4034        /**
    4035         * Return the exponent of a float number, removing the bias.
    4036         * <p>
    4037         * For float numbers of the form 2<sup>x</sup>, the unbiased
    4038         * exponent is exactly x.
    4039         * </p>
    4040         * @param f number from which exponent is requested
    4041         * @return exponent for d in IEEE754 representation, without bias
    4042         */
    4043        public static int getExponent(final float f) {
    4044            return ((Float.floatToIntBits(f) >>> 23) & 0xff) - 127;
    4045        }
    4046
    4047    }