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