AccurateMath.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math4.core.jdkmath;

  19. import java.io.PrintStream;

  21. import org.apache.commons.numbers.core.Precision;

  23. /**
  24.  * Portable alternative to {@link Math} and {@link StrictMath}.
  25.  * <p>
  26.  * Caveat: At the time of implementation, the {@code AccurateMath} functions
  27.  * were often faster and/or more accurate than their JDK equivalent.
  28.  * Nowadays, it would not be surprising that they are always slower (due
  29.  * to the various JVM optimizations that have appeared since Java 5).
  30.  * However, any change to this class should ensure that the current
  31.  * accuracy is not lost.
  32.  * </p>
  33.  * <p>
  34.  * AccurateMath speed is achieved by relying heavily on optimizing compilers
  35.  * to native code present in many JVMs today and use of large tables.
  36.  * The larger tables are lazily initialized on first use, so that the setup
  37.  * time does not penalize methods that don't need them.
  38.  * </p>
  39.  * <p>
  40.  * Note that AccurateMath is extensively used inside Apache Commons Math,
  41.  * so by calling some algorithms, the overhead when the tables need to be
  42.  * initialized will occur regardless of the end-user calling AccurateMath
  43.  * methods directly or not.
  44.  * Performance figures for a specific JVM and hardware can be evaluated by
  45.  * running the AccurateMathTestPerformance tests in the test directory of
  46.  * the source distribution.
  47.  * </p>
  48.  * <p>
  49.  * AccurateMath accuracy should be mostly independent of the JVM as it relies only
  50.  * on IEEE-754 basic operations and on embedded tables. Almost all operations
  51.  * are accurate to about 0.5 ulp throughout the domain range. This statement,
  52.  * of course is only a rough global observed behavior, it is <em>not</em> a
  53.  * guarantee for <em>every</em> double numbers input (see William Kahan's
  54.  * <a href="http://en.wikipedia.org/wiki/Rounding#The_table-maker.27s_dilemma">
  55.  * Table Maker's Dilemma</a>).
  56.  * </p>
  57.  * <p>
  58.  * AccurateMath implements the following methods not found in Math/StrictMath:
  59.  * <ul>
  60.  * <li>{@link #asinh(double)}</li>
  61.  * <li>{@link #acosh(double)}</li>
  62.  * <li>{@link #atanh(double)}</li>
  63.  * </ul>
  64.  *
  65.  * @since 2.2
  66.  */
  67. public final class AccurateMath {
  68.     /** Archimede's constant PI, ratio of circle circumference to diameter. */
  69.     public static final double PI = 105414357.0 / 33554432.0 + 1.984187159361080883e-9;

  71.     /** Napier's constant e, base of the natural logarithm. */
  72.     public static final double E = 2850325.0 / 1048576.0 + 8.254840070411028747e-8;

  74.     /** Index of exp(0) in the array of integer exponentials. */
  75.     static final int EXP_INT_TABLE_MAX_INDEX = 750;
  76.     /** Length of the array of integer exponentials. */
  77.     static final int EXP_INT_TABLE_LEN = EXP_INT_TABLE_MAX_INDEX * 2;
  78.     /** Logarithm table length. */
  79.     static final int LN_MANT_LEN = 1024;
  80.     /** Exponential fractions table length. */
  81.     static final int EXP_FRAC_TABLE_LEN = 1025; // 0, 1/1024, ... 1024/1024

  83.     /** Error message. */
  84.     private static final String ZERO_DENOMINATOR_MSG = "Division by zero";
  85.     /** Error message. */
  86.     private static final String OVERFLOW_MSG = "Overflow";
  87.     /** StrictMath.log(Double.MAX_VALUE): {@value}. */
  88.     private static final double LOG_MAX_VALUE = StrictMath.log(Double.MAX_VALUE);

  90.     /** Indicator for tables initialization.
  91.      * <p>
  92.      * This compile-time constant should be set to true only if one explicitly
  93.      * wants to compute the tables at class loading time instead of using the
  94.      * already computed ones provided as literal arrays below.
  95.      * </p>
  96.      */
  97.     private static final boolean RECOMPUTE_TABLES_AT_RUNTIME = false;

  99.     /** log(2) (high bits). */
  100.     private static final double LN_2_A = 0.693147063255310059;

  102.     /** log(2) (low bits). */
  103.     private static final double LN_2_B = 1.17304635250823482e-7;

  105.     /** Coefficients for log, when input 0.99 < x < 1.01. */
  106.     private static final double[][] LN_QUICK_COEF = {
  107.         {1.0, 5.669184079525E-24},
  108.         {-0.25, -0.25},
  109.         {0.3333333134651184, 1.986821492305628E-8},
  110.         {-0.25, -6.663542893624021E-14},
  111.         {0.19999998807907104, 1.1921056801463227E-8},
  112.         {-0.1666666567325592, -7.800414592973399E-9},
  113.         {0.1428571343421936, 5.650007086920087E-9},
  114.         {-0.12502530217170715, -7.44321345601866E-11},
  115.         {0.11113807559013367, 9.219544613762692E-9},
  116.     };

  118.     /** Coefficients for log in the range of 1.0 < x < 1.0 + 2^-10. */
  119.     private static final double[][] LN_HI_PREC_COEF = {
  120.         {1.0, -6.032174644509064E-23},
  121.         {-0.25, -0.25},
  122.         {0.3333333134651184, 1.9868161777724352E-8},
  123.         {-0.2499999701976776, -2.957007209750105E-8},
  124.         {0.19999954104423523, 1.5830993332061267E-10},
  125.         {-0.16624879837036133, -2.6033824355191673E-8}
  126.     };

  128.     /** Sine, Cosine, Tangent tables are for 0, 1/8, 2/8, ... 13/8 = PI/2 approx. */
  129.     private static final int SINE_TABLE_LEN = 14;

  131.     /** Sine table (high bits). */
  132.     private static final double[] SINE_TABLE_A =
  133.         {
  134.         +0.0d,
  135.         +0.1246747374534607d,
  136.         +0.24740394949913025d,
  137.         +0.366272509098053d,
  138.         +0.4794255495071411d,
  139.         +0.5850973129272461d,
  140.         +0.6816387176513672d,
  141.         +0.7675435543060303d,
  142.         +0.8414709568023682d,
  143.         +0.902267575263977d,
  144.         +0.9489846229553223d,
  145.         +0.9808930158615112d,
  146.         +0.9974949359893799d,
  147.         +0.9985313415527344d,
  148.     };

  150.     /** Sine table (low bits). */
  151.     private static final double[] SINE_TABLE_B =
  152.         {
  153.         +0.0d,
  154.         -4.068233003401932E-9d,
  155.         +9.755392680573412E-9d,
  156.         +1.9987994582857286E-8d,
  157.         -1.0902938113007961E-8d,
  158.         -3.9986783938944604E-8d,
  159.         +4.23719669792332E-8d,
  160.         -5.207000323380292E-8d,
  161.         +2.800552834259E-8d,
  162.         +1.883511811213715E-8d,
  163.         -3.5997360512765566E-9d,
  164.         +4.116164446561962E-8d,
  165.         +5.0614674548127384E-8d,
  166.         -1.0129027912496858E-9d,
  167.     };

  169.     /** Cosine table (high bits). */
  170.     private static final double[] COSINE_TABLE_A =
  171.         {
  172.         +1.0d,
  173.         +0.9921976327896118d,
  174.         +0.9689123630523682d,
  175.         +0.9305076599121094d,
  176.         +0.8775825500488281d,
  177.         +0.8109631538391113d,
  178.         +0.7316888570785522d,
  179.         +0.6409968137741089d,
  180.         +0.5403022766113281d,
  181.         +0.4311765432357788d,
  182.         +0.3153223395347595d,
  183.         +0.19454771280288696d,
  184.         +0.07073719799518585d,
  185.         -0.05417713522911072d,
  186.     };

  188.     /** Cosine table (low bits). */
  189.     private static final double[] COSINE_TABLE_B =
  190.         {
  191.         +0.0d,
  192.         +3.4439717236742845E-8d,
  193.         +5.865827662008209E-8d,
  194.         -3.7999795083850525E-8d,
  195.         +1.184154459111628E-8d,
  196.         -3.43338934259355E-8d,
  197.         +1.1795268640216787E-8d,
  198.         +4.438921624363781E-8d,
  199.         +2.925681159240093E-8d,
  200.         -2.6437112632041807E-8d,
  201.         +2.2860509143963117E-8d,
  202.         -4.813899778443457E-9d,
  203.         +3.6725170580355583E-9d,
  204.         +2.0217439756338078E-10d,
  205.     };


  208.     /** Tangent table, used by atan() (high bits). */
  209.     private static final double[] TANGENT_TABLE_A =
  210.         {
  211.         +0.0d,
  212.         +0.1256551444530487d,
  213.         +0.25534194707870483d,
  214.         +0.3936265707015991d,
  215.         +0.5463024377822876d,
  216.         +0.7214844226837158d,
  217.         +0.9315965175628662d,
  218.         +1.1974215507507324d,
  219.         +1.5574076175689697d,
  220.         +2.092571258544922d,
  221.         +3.0095696449279785d,
  222.         +5.041914939880371d,
  223.         +14.101419448852539d,
  224.         -18.430862426757812d,
  225.     };

  227.     /** Tangent table, used by atan() (low bits). */
  228.     private static final double[] TANGENT_TABLE_B =
  229.         {
  230.         +0.0d,
  231.         -7.877917738262007E-9d,
  232.         -2.5857668567479893E-8d,
  233.         +5.2240336371356666E-9d,
  234.         +5.206150291559893E-8d,
  235.         +1.8307188599677033E-8d,
  236.         -5.7618793749770706E-8d,
  237.         +7.848361555046424E-8d,
  238.         +1.0708593250394448E-7d,
  239.         +1.7827257129423813E-8d,
  240.         +2.893485277253286E-8d,
  241.         +3.1660099222737955E-7d,
  242.         +4.983191803254889E-7d,
  243.         -3.356118100840571E-7d,
  244.     };

  246.     /** Bits of 1/(2*pi), need for reducePayneHanek(). */
  247.     private static final long[] RECIP_2PI = new long[] {
  248.         (0x28be60dbL << 32) | 0x9391054aL,
  249.         (0x7f09d5f4L << 32) | 0x7d4d3770L,
  250.         (0x36d8a566L << 32) | 0x4f10e410L,
  251.         (0x7f9458eaL << 32) | 0xf7aef158L,
  252.         (0x6dc91b8eL << 32) | 0x909374b8L,
  253.         (0x01924bbaL << 32) | 0x82746487L,
  254.         (0x3f877ac7L << 32) | 0x2c4a69cfL,
  255.         (0xba208d7dL << 32) | 0x4baed121L,
  256.         (0x3a671c09L << 32) | 0xad17df90L,
  257.         (0x4e64758eL << 32) | 0x60d4ce7dL,
  258.         (0x272117e2L << 32) | 0xef7e4a0eL,
  259.         (0xc7fe25ffL << 32) | 0xf7816603L,
  260.         (0xfbcbc462L << 32) | 0xd6829b47L,
  261.         (0xdb4d9fb3L << 32) | 0xc9f2c26dL,
  262.         (0xd3d18fd9L << 32) | 0xa797fa8bL,
  263.         (0x5d49eeb1L << 32) | 0xfaf97c5eL,
  264.         (0xcf41ce7dL << 32) | 0xe294a4baL,
  265.         (0x9afed7ecL << 32)};

  267.     /** Bits of pi/4, need for reducePayneHanek(). */
  268.     private static final long[] PI_O_4_BITS = new long[] {
  269.         (0xc90fdaa2L << 32) | 0x2168c234L,
  270.         (0xc4c6628bL << 32) | 0x80dc1cd1L };

  272.     /** Eighths.
  273.      * This is used by sinQ, because its faster to do a table lookup than
  274.      * a multiply in this time-critical routine
  275.      */
  276.     private static final double[] EIGHTHS = {
  277.         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};

  279.     /** Table of 2^((n+2)/3). */
  280.     private static final double[] CBRTTWO = {0.6299605249474366,
  281.                                              0.7937005259840998,
  282.                                              1.0,
  283.                                              1.2599210498948732,
  284.                                              1.5874010519681994};

  286.     /*
  287.      *  There are 52 bits in the mantissa of a double.
  288.      *  For additional precision, the code splits double numbers into two parts,
  289.      *  by clearing the low order 30 bits if possible, and then performs the arithmetic
  290.      *  on each half separately.
  291.      */

  293.     /**
  294.      * 0x40000000 - used to split a double into two parts, both with the low order bits cleared.
  295.      * Equivalent to 2^30.
  296.      */
  297.     private static final long HEX_40000000 = 0x40000000L; // 1073741824L

  299.     /** Mask used to clear low order 30 bits. */
  300.     private static final long MASK_30BITS = -1L - (HEX_40000000 - 1); // 0xFFFFFFFFC0000000L;

  302.     /** Mask used to clear the non-sign part of an int. */
  303.     private static final int MASK_NON_SIGN_INT = 0x7fffffff;

  305.     /** Mask used to clear the non-sign part of a long. */
  306.     private static final long MASK_NON_SIGN_LONG = 0x7fffffffffffffffL;

  308.     /** Mask used to extract exponent from double bits. */
  309.     private static final long MASK_DOUBLE_EXPONENT = 0x7ff0000000000000L;

  311.     /** Mask used to extract mantissa from double bits. */
  312.     private static final long MASK_DOUBLE_MANTISSA = 0x000fffffffffffffL;

  314.     /** Mask used to add implicit high order bit for normalized double. */
  315.     private static final long IMPLICIT_HIGH_BIT = 0x0010000000000000L;

  317.     /** 2^52 - double numbers this large must be integral (no fraction) or NaN or Infinite. */
  318.     private static final double TWO_POWER_52 = 4503599627370496.0;

  320.     /** Constant: {@value}. */
  321.     private static final double F_1_3 = 1d / 3d;
  322.     /** Constant: {@value}. */
  323.     private static final double F_1_5 = 1d / 5d;
  324.     /** Constant: {@value}. */
  325.     private static final double F_1_7 = 1d / 7d;
  326.     /** Constant: {@value}. */
  327.     private static final double F_1_9 = 1d / 9d;
  328.     /** Constant: {@value}. */
  329.     private static final double F_1_11 = 1d / 11d;
  330.     /** Constant: {@value}. */
  331.     private static final double F_1_13 = 1d / 13d;
  332.     /** Constant: {@value}. */
  333.     private static final double F_1_15 = 1d / 15d;
  334.     /** Constant: {@value}. */
  335.     private static final double F_1_17 = 1d / 17d;
  336.     /** Constant: {@value}. */
  337.     private static final double F_3_4 = 3d / 4d;
  338.     /** Constant: {@value}. */
  339.     private static final double F_15_16 = 15d / 16d;
  340.     /** Constant: {@value}. */
  341.     private static final double F_13_14 = 13d / 14d;
  342.     /** Constant: {@value}. */
  343.     private static final double F_11_12 = 11d / 12d;
  344.     /** Constant: {@value}. */
  345.     private static final double F_9_10 = 9d / 10d;
  346.     /** Constant: {@value}. */
  347.     private static final double F_7_8 = 7d / 8d;
  348.     /** Constant: {@value}. */
  349.     private static final double F_5_6 = 5d / 6d;
  350.     /** Constant: {@value}. */
  351.     private static final double F_1_2 = 1d / 2d;
  352.     /** Constant: {@value}. */
  353.     private static final double F_1_4 = 1d / 4d;

  355.     /**
  356.      * Private Constructor.
  357.      */
  358.     private AccurateMath() {}

  360.     // Generic helper methods

  362.     /**
  363.      * Get the high order bits from the mantissa.
  364.      * Equivalent to adding and subtracting HEX_40000 but also works for very large numbers
  365.      *
  366.      * @param d the value to split
  367.      * @return the high order part of the mantissa
  368.      */
  369.     private static double doubleHighPart(double d) {
  370.         if (d > -Precision.SAFE_MIN && d < Precision.SAFE_MIN) {
  371.             return d; // These are un-normalised - don't try to convert
  372.         }
  373.         long xl = Double.doubleToRawLongBits(d); // can take raw bits because just gonna convert it back
  374.         xl &= MASK_30BITS; // Drop low order bits
  375.         return Double.longBitsToDouble(xl);
  376.     }

  378.     /** Compute the hyperbolic cosine of a number.
  379.      * @param x number on which evaluation is done
  380.      * @return hyperbolic cosine of x
  381.      */
  382.     public static double cosh(double x) {
  383.         if (Double.isNaN(x)) {
  384.             return x;
  385.         }

  387.         // cosh[z] = (exp(z) + exp(-z))/2

  389.         // for numbers with magnitude 20 or so,
  390.         // exp(-z) can be ignored in comparison with exp(z)

  392.         if (x > 20) {
  393.             if (x >= LOG_MAX_VALUE) {
  394.                 // Avoid overflow (MATH-905).
  395.                 final double t = exp(0.5 * x);
  396.                 return (0.5 * t) * t;
  397.             } else {
  398.                 return 0.5 * exp(x);
  399.             }
  400.         } else if (x < -20) {
  401.             if (x <= -LOG_MAX_VALUE) {
  402.                 // Avoid overflow (MATH-905).
  403.                 final double t = exp(-0.5 * x);
  404.                 return (0.5 * t) * t;
  405.             } else {
  406.                 return 0.5 * exp(-x);
  407.             }
  408.         }

  410.         final double[] hiPrec = new double[2];
  411.         if (x < 0.0) {
  412.             x = -x;
  413.         }
  414.         exp(x, 0.0, hiPrec);

  416.         double ya = hiPrec[0] + hiPrec[1];
  417.         double yb = -(ya - hiPrec[0] - hiPrec[1]);

  419.         double temp = ya * HEX_40000000;
  420.         double yaa = ya + temp - temp;
  421.         double yab = ya - yaa;

  423.         // recip = 1/y
  424.         double recip = 1.0 / ya;
  425.         temp = recip * HEX_40000000;
  426.         double recipa = recip + temp - temp;
  427.         double recipb = recip - recipa;

  429.         // Correct for rounding in division
  430.         recipb += (1.0 - yaa * recipa - yaa * recipb - yab * recipa - yab * recipb) * recip;
  431.         // Account for yb
  432.         recipb += -yb * recip * recip;

  434.         // y = y + 1/y
  435.         temp = ya + recipa;
  436.         yb += -(temp - ya - recipa);
  437.         ya = temp;
  438.         temp = ya + recipb;
  439.         yb += -(temp - ya - recipb);
  440.         ya = temp;

  442.         double result = ya + yb;
  443.         result *= 0.5;
  444.         return result;
  445.     }

  447.     /** Compute the hyperbolic sine of a number.
  448.      * @param x number on which evaluation is done
  449.      * @return hyperbolic sine of x
  450.      */
  451.     public static double sinh(double x) {
  452.         boolean negate = false;
  453.         if (Double.isNaN(x)) {
  454.             return x;
  455.         }

  457.         // sinh[z] = (exp(z) - exp(-z) / 2

  459.         // for values of z larger than about 20,
  460.         // exp(-z) can be ignored in comparison with exp(z)

  462.         if (x > 20) {
  463.             if (x >= LOG_MAX_VALUE) {
  464.                 // Avoid overflow (MATH-905).
  465.                 final double t = exp(0.5 * x);
  466.                 return (0.5 * t) * t;
  467.             } else {
  468.                 return 0.5 * exp(x);
  469.             }
  470.         } else if (x < -20) {
  471.             if (x <= -LOG_MAX_VALUE) {
  472.                 // Avoid overflow (MATH-905).
  473.                 final double t = exp(-0.5 * x);
  474.                 return (-0.5 * t) * t;
  475.             } else {
  476.                 return -0.5 * exp(-x);
  477.             }
  478.         }

  480.         if (x == 0) {
  481.             return x;
  482.         }

  484.         if (x < 0.0) {
  485.             x = -x;
  486.             negate = true;
  487.         }

  489.         double result;

  491.         if (x > 0.25) {
  492.             double[] hiPrec = new double[2];
  493.             exp(x, 0.0, hiPrec);

  495.             double ya = hiPrec[0] + hiPrec[1];
  496.             double yb = -(ya - hiPrec[0] - hiPrec[1]);

  498.             double temp = ya * HEX_40000000;
  499.             double yaa = ya + temp - temp;
  500.             double yab = ya - yaa;

  502.             // recip = 1/y
  503.             double recip = 1.0 / ya;
  504.             temp = recip * HEX_40000000;
  505.             double recipa = recip + temp - temp;
  506.             double recipb = recip - recipa;

  508.             // Correct for rounding in division
  509.             recipb += (1.0 - yaa * recipa - yaa * recipb - yab * recipa - yab * recipb) * recip;
  510.             // Account for yb
  511.             recipb += -yb * recip * recip;

  513.             recipa = -recipa;
  514.             recipb = -recipb;

  516.             // y = y + 1/y
  517.             temp = ya + recipa;
  518.             yb += -(temp - ya - recipa);
  519.             ya = temp;
  520.             temp = ya + recipb;
  521.             yb += -(temp - ya - recipb);
  522.             ya = temp;

  524.             result = ya + yb;
  525.             result *= 0.5;
  526.         } else {
  527.             double[] hiPrec = new double[2];
  528.             expm1(x, hiPrec);

  530.             double ya = hiPrec[0] + hiPrec[1];
  531.             double yb = -(ya - hiPrec[0] - hiPrec[1]);

  533.             /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
  534.             double denom = 1.0 + ya;
  535.             double denomr = 1.0 / denom;
  536.             double denomb = -(denom - 1.0 - ya) + yb;
  537.             double ratio = ya * denomr;
  538.             double temp = ratio * HEX_40000000;
  539.             double ra = ratio + temp - temp;
  540.             double rb = ratio - ra;

  542.             temp = denom * HEX_40000000;
  543.             double za = denom + temp - temp;
  544.             double zb = denom - za;

  546.             rb += (ya - za * ra - za * rb - zb * ra - zb * rb) * denomr;

  548.             // Adjust for yb
  549.             rb += yb * denomr;                      // numerator
  550.             rb += -ya * denomb * denomr * denomr;   // denominator

  552.             // y = y - 1/y
  553.             temp = ya + ra;
  554.             yb += -(temp - ya - ra);
  555.             ya = temp;
  556.             temp = ya + rb;
  557.             yb += -(temp - ya - rb);
  558.             ya = temp;

  560.             result = ya + yb;
  561.             result *= 0.5;
  562.         }

  564.         if (negate) {
  565.             result = -result;
  566.         }

  568.         return result;
  569.     }

  571.     /** Compute the hyperbolic tangent of a number.
  572.      * @param x number on which evaluation is done
  573.      * @return hyperbolic tangent of x
  574.      */
  575.     public static double tanh(double x) {
  576.         boolean negate = false;

  578.         if (Double.isNaN(x)) {
  579.             return x;
  580.         }

  582.         // tanh[z] = sinh[z] / cosh[z]
  583.         // = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
  584.         // = (exp(2x) - 1) / (exp(2x) + 1)

  586.         // for magnitude > 20, sinh[z] == cosh[z] in double precision

  588.         if (x > 20.0) {
  589.             return 1.0;
  590.         }

  592.         if (x < -20) {
  593.             return -1.0;
  594.         }

  596.         if (x == 0) {
  597.             return x;
  598.         }

  600.         if (x < 0.0) {
  601.             x = -x;
  602.             negate = true;
  603.         }

  605.         double result;
  606.         if (x >= 0.5) {
  607.             double[] hiPrec = new double[2];
  608.             // tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
  609.             exp(x * 2.0, 0.0, hiPrec);

  611.             double ya = hiPrec[0] + hiPrec[1];
  612.             double yb = -(ya - hiPrec[0] - hiPrec[1]);

  614.             /* Numerator */
  615.             double na = -1.0 + ya;
  616.             double nb = -(na + 1.0 - ya);
  617.             double temp = na + yb;
  618.             nb += -(temp - na - yb);
  619.             na = temp;

  621.             /* Denominator */
  622.             double da = 1.0 + ya;
  623.             double db = -(da - 1.0 - ya);
  624.             temp = da + yb;
  625.             db += -(temp - da - yb);
  626.             da = temp;

  628.             temp = da * HEX_40000000;
  629.             double daa = da + temp - temp;
  630.             double dab = da - daa;

  632.             // ratio = na/da
  633.             double ratio = na / da;
  634.             temp = ratio * HEX_40000000;
  635.             double ratioa = ratio + temp - temp;
  636.             double ratiob = ratio - ratioa;

  638.             // Correct for rounding in division
  639.             ratiob += (na - daa * ratioa - daa * ratiob - dab * ratioa - dab * ratiob) / da;

  641.             // Account for nb
  642.             ratiob += nb / da;
  643.             // Account for db
  644.             ratiob += -db * na / da / da;

  646.             result = ratioa + ratiob;
  647.         } else {
  648.             double[] hiPrec = new double[2];
  649.             // tanh(x) = expm1(2x) / (expm1(2x) + 2)
  650.             expm1(x * 2.0, hiPrec);

  652.             double ya = hiPrec[0] + hiPrec[1];
  653.             double yb = -(ya - hiPrec[0] - hiPrec[1]);

  655.             /* Numerator */
  656.             double na = ya;
  657.             double nb = yb;

  659.             /* Denominator */
  660.             double da = 2.0 + ya;
  661.             double db = -(da - 2.0 - ya);
  662.             double temp = da + yb;
  663.             db += -(temp - da - yb);
  664.             da = temp;

  666.             temp = da * HEX_40000000;
  667.             double daa = da + temp - temp;
  668.             double dab = da - daa;

  670.             // ratio = na/da
  671.             double ratio = na / da;
  672.             temp = ratio * HEX_40000000;
  673.             double ratioa = ratio + temp - temp;
  674.             double ratiob = ratio - ratioa;

  676.             // Correct for rounding in division
  677.             ratiob += (na - daa * ratioa - daa * ratiob - dab * ratioa - dab * ratiob) / da;

  679.             // Account for nb
  680.             ratiob += nb / da;
  681.             // Account for db
  682.             ratiob += -db * na / da / da;

  684.             result = ratioa + ratiob;
  685.         }

  687.         if (negate) {
  688.             result = -result;
  689.         }

  691.         return result;
  692.     }

  694.     /** Compute the inverse hyperbolic cosine of a number.
  695.      * @param a number on which evaluation is done
  696.      * @return inverse hyperbolic cosine of a
  697.      */
  698.     public static double acosh(final double a) {
  699.         return AccurateMath.log(a + Math.sqrt(a * a - 1));
  700.     }

  702.     /** Compute the inverse hyperbolic sine of a number.
  703.      * @param a number on which evaluation is done
  704.      * @return inverse hyperbolic sine of a
  705.      */
  706.     public static double asinh(double a) {
  707.         boolean negative = false;
  708.         if (a < 0) {
  709.             negative = true;
  710.             a = -a;
  711.         }

  713.         double absAsinh;
  714.         if (a > 0.167) {
  715.             absAsinh = AccurateMath.log(Math.sqrt(a * a + 1) + a);
  716.         } else {
  717.             final double a2 = a * a;
  718.             if (a > 0.097) {
  719.                 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);
  720.             } else if (a > 0.036) {
  721.                 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);
  722.             } else if (a > 0.0036) {
  723.                 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);
  724.             } else {
  725.                 absAsinh = a * (1 - a2 * (F_1_3 - a2 * F_1_5 * F_3_4) * F_1_2);
  726.             }
  727.         }

  729.         return negative ? -absAsinh : absAsinh;
  730.     }

  732.     /** Compute the inverse hyperbolic tangent of a number.
  733.      * @param a number on which evaluation is done
  734.      * @return inverse hyperbolic tangent of a
  735.      */
  736.     public static double atanh(double a) {
  737.         boolean negative = false;
  738.         if (a < 0) {
  739.             negative = true;
  740.             a = -a;
  741.         }

  743.         double absAtanh;
  744.         if (a > 0.15) {
  745.             absAtanh = 0.5 * AccurateMath.log((1 + a) / (1 - a));
  746.         } else {
  747.             final double a2 = a * a;
  748.             if (a > 0.087) {
  749.                 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))))))));
  750.             } else if (a > 0.031) {
  751.                 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))))));
  752.             } else if (a > 0.003) {
  753.                 absAtanh = a * (1 + a2 * (F_1_3 + a2 * (F_1_5 + a2 * (F_1_7 + a2 * F_1_9))));
  754.             } else {
  755.                 absAtanh = a * (1 + a2 * (F_1_3 + a2 * F_1_5));
  756.             }
  757.         }

  759.         return negative ? -absAtanh : absAtanh;
  760.     }

  762.     /** Compute the signum of a number.
  763.      * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise
  764.      * @param a number on which evaluation is done
  765.      * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
  766.      */
  767.     public static double signum(final double a) {
  768.         return (a < 0.0) ? -1.0 : ((a > 0.0) ? 1.0 : a); // return +0.0/-0.0/NaN depending on a
  769.     }

  771.     /** Compute the signum of a number.
  772.      * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise
  773.      * @param a number on which evaluation is done
  774.      * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
  775.      */
  776.     public static float signum(final float a) {
  777.         return (a < 0.0f) ? -1.0f : ((a > 0.0f) ? 1.0f : a); // return +0.0/-0.0/NaN depending on a
  778.     }

  780.     /** Compute next number towards positive infinity.
  781.      * @param a number to which neighbor should be computed
  782.      * @return neighbor of a towards positive infinity
  783.      */
  784.     public static double nextUp(final double a) {
  785.         return nextAfter(a, Double.POSITIVE_INFINITY);
  786.     }

  788.     /** Compute next number towards positive infinity.
  789.      * @param a number to which neighbor should be computed
  790.      * @return neighbor of a towards positive infinity
  791.      */
  792.     public static float nextUp(final float a) {
  793.         return nextAfter(a, Float.POSITIVE_INFINITY);
  794.     }

  796.     /** Compute next number towards negative infinity.
  797.      * @param a number to which neighbor should be computed
  798.      * @return neighbor of a towards negative infinity
  799.      * @since 3.4
  800.      */
  801.     public static double nextDown(final double a) {
  802.         return nextAfter(a, Double.NEGATIVE_INFINITY);
  803.     }

  805.     /** Compute next number towards negative infinity.
  806.      * @param a number to which neighbor should be computed
  807.      * @return neighbor of a towards negative infinity
  808.      * @since 3.4
  809.      */
  810.     public static float nextDown(final float a) {
  811.         return nextAfter(a, Float.NEGATIVE_INFINITY);
  812.     }

  814.     /**
  815.      * Exponential function.
  816.      *
  817.      * Computes exp(x), function result is nearly rounded.   It will be correctly
  818.      * rounded to the theoretical value for 99.9% of input values, otherwise it will
  819.      * have a 1 ULP error.
  820.      *
  821.      * Method:
  822.      *    Lookup intVal = exp(int(x))
  823.      *    Lookup fracVal = exp(int(x-int(x) / 1024.0) * 1024.0 );
  824.      *    Compute z as the exponential of the remaining bits by a polynomial minus one
  825.      *    exp(x) = intVal * fracVal * (1 + z)
  826.      *
  827.      * Accuracy:
  828.      *    Calculation is done with 63 bits of precision, so result should be correctly
  829.      *    rounded for 99.9% of input values, with less than 1 ULP error otherwise.
  830.      *
  831.      * @param x   a double
  832.      * @return double e<sup>x</sup>
  833.      */
  834.     public static double exp(double x) {
  835.         return exp(x, 0.0, null);
  836.     }

  838.     /**
  839.      * Internal helper method for exponential function.
  840.      * @param x original argument of the exponential function
  841.      * @param extra extra bits of precision on input (To Be Confirmed)
  842.      * @param hiPrec extra bits of precision on output (To Be Confirmed)
  843.      * @return exp(x)
  844.      */
  845.     private static double exp(double x, double extra, double[] hiPrec) {
  846.         double intPartA;
  847.         double intPartB;
  848.         int intVal = (int) x;

  850.         /* Lookup exp(floor(x)).
  851.          * intPartA will have the upper 22 bits, intPartB will have the lower
  852.          * 52 bits.
  853.          */
  854.         if (x < 0.0) {

  856.             // We don't check against intVal here as conversion of large negative double values
  857.             // may be affected by a JIT bug. Subsequent comparisons can safely use intVal
  858.             if (x < -746d) {
  859.                 if (hiPrec != null) {
  860.                     hiPrec[0] = 0.0;
  861.                     hiPrec[1] = 0.0;
  862.                 }
  863.                 return 0.0;
  864.             }

  866.             if (intVal < -709) {
  867.                 /* This will produce a subnormal output */
  868.                 final double result = exp(x + 40.19140625, extra, hiPrec) / 285040095144011776.0;
  869.                 if (hiPrec != null) {
  870.                     hiPrec[0] /= 285040095144011776.0;
  871.                     hiPrec[1] /= 285040095144011776.0;
  872.                 }
  873.                 return result;
  874.             }

  876.             if (intVal == -709) {
  877.                 /* exp(1.494140625) is nearly a machine number... */
  878.                 final double result = exp(x + 1.494140625, extra, hiPrec) / 4.455505956692756620;
  879.                 if (hiPrec != null) {
  880.                     hiPrec[0] /= 4.455505956692756620;
  881.                     hiPrec[1] /= 4.455505956692756620;
  882.                 }
  883.                 return result;
  884.             }

  886.             intVal--;
  887.         } else {
  888.             if (intVal > 709) {
  889.                 if (hiPrec != null) {
  890.                     hiPrec[0] = Double.POSITIVE_INFINITY;
  891.                     hiPrec[1] = 0.0;
  892.                 }
  893.                 return Double.POSITIVE_INFINITY;
  894.             }
  895.         }

  897.         intPartA = ExpIntTable.EXP_INT_TABLE_A[EXP_INT_TABLE_MAX_INDEX + intVal];
  898.         intPartB = ExpIntTable.EXP_INT_TABLE_B[EXP_INT_TABLE_MAX_INDEX + intVal];

  900.         /* Get the fractional part of x, find the greatest multiple of 2^-10 less than
  901.          * x and look up the exp function of it.
  902.          * fracPartA will have the upper 22 bits, fracPartB the lower 52 bits.
  903.          */
  904.         final int intFrac = (int) ((x - intVal) * 1024.0);
  905.         final double fracPartA = ExpFracTable.EXP_FRAC_TABLE_A[intFrac];
  906.         final double fracPartB = ExpFracTable.EXP_FRAC_TABLE_B[intFrac];

  908.         /* epsilon is the difference in x from the nearest multiple of 2^-10.  It
  909.          * has a value in the range 0 <= epsilon < 2^-10.
  910.          * Do the subtraction from x as the last step to avoid possible loss of precision.
  911.          */
  912.         final double epsilon = x - (intVal + intFrac / 1024.0);

  914.         /* Compute z = exp(epsilon) - 1.0 via a minimax polynomial.  z has
  915.        full double precision (52 bits).  Since z < 2^-10, we will have
  916.        62 bits of precision when combined with the constant 1.  This will be
  917.        used in the last addition below to get proper rounding. */

  919.         /* Remez generated polynomial.  Converges on the interval [0, 2^-10], error
  920.        is less than 0.5 ULP */
  921.         double z = 0.04168701738764507;
  922.         z = z * epsilon + 0.1666666505023083;
  923.         z = z * epsilon + 0.5000000000042687;
  924.         z = z * epsilon + 1.0;
  925.         z = z * epsilon + -3.940510424527919E-20;

  927.         /* Compute (intPartA+intPartB) * (fracPartA+fracPartB) by binomial
  928.        expansion.
  929.        tempA is exact since intPartA and intPartB only have 22 bits each.
  930.        tempB will have 52 bits of precision.
  931.          */
  932.         double tempA = intPartA * fracPartA;
  933.         double tempB = intPartA * fracPartB + intPartB * fracPartA + intPartB * fracPartB;

  935.         /* Compute the result.  (1+z)(tempA+tempB).  Order of operations is
  936.        important.  For accuracy add by increasing size.  tempA is exact and
  937.        much larger than the others.  If there are extra bits specified from the
  938.        pow() function, use them. */
  939.         final double tempC = tempB + tempA;

  941.         // If tempC is positive infinite, the evaluation below could result in NaN,
  942.         // because z could be negative at the same time.
  943.         if (tempC == Double.POSITIVE_INFINITY) {
  944.             return Double.POSITIVE_INFINITY;
  945.         }

  947.         final double result;
  948.         if (extra != 0.0) {
  949.             result = tempC * extra * z + tempC * extra + tempC * z + tempB + tempA;
  950.         } else {
  951.             result = tempC * z + tempB + tempA;
  952.         }

  954.         if (hiPrec != null) {
  955.             // If requesting high precision
  956.             hiPrec[0] = tempA;
  957.             hiPrec[1] = tempC * extra * z + tempC * extra + tempC * z + tempB;
  958.         }

  960.         return result;
  961.     }

  963.     /**Compute exp(x) - 1.
  964.      *
  965.      * @param x number to compute shifted exponential
  966.      * @return exp(x) - 1
  967.      */
  968.     public static double expm1(double x) {
  969.         return expm1(x, null);
  970.     }

  972.     /** Internal helper method for expm1.
  973.      * @param x number to compute shifted exponential
  974.      * @param hiPrecOut receive high precision result for -1.0 < x < 1.0
  975.      * @return exp(x) - 1
  976.      */
  977.     private static double expm1(double x, double[] hiPrecOut) {
  978.         if (Double.isNaN(x) || x == 0.0) { // NaN or zero
  979.             return x;
  980.         }

  982.         if (x <= -1.0 || x >= 1.0) {
  983.             // If not between +/- 1.0
  984.             //return exp(x) - 1.0;
  985.             double[] hiPrec = new double[2];
  986.             exp(x, 0.0, hiPrec);
  987.             if (x > 0.0) {
  988.                 return -1.0 + hiPrec[0] + hiPrec[1];
  989.             } else {
  990.                 final double ra = -1.0 + hiPrec[0];
  991.                 double rb = -(ra + 1.0 - hiPrec[0]);
  992.                 rb += hiPrec[1];
  993.                 return ra + rb;
  994.             }
  995.         }

  997.         double baseA;
  998.         double baseB;
  999.         double epsilon;
  1000.         boolean negative = false;

  1002.         if (x < 0.0) {
  1003.             x = -x;
  1004.             negative = true;
  1005.         }

  1007.         int intFrac = (int) (x * 1024.0);
  1008.         double tempA = ExpFracTable.EXP_FRAC_TABLE_A[intFrac] - 1.0;
  1009.         double tempB = ExpFracTable.EXP_FRAC_TABLE_B[intFrac];

  1011.         double temp = tempA + tempB;
  1012.         tempB = -(temp - tempA - tempB);
  1013.         tempA = temp;

  1015.         temp = tempA * HEX_40000000;
  1016.         baseA = tempA + temp - temp;
  1017.         baseB = tempB + (tempA - baseA);

  1019.         epsilon = x - intFrac / 1024.0;


  1022.         /* Compute expm1(epsilon) */
  1023.         double zb = 0.008336750013465571;
  1024.         zb = zb * epsilon + 0.041666663879186654;
  1025.         zb = zb * epsilon + 0.16666666666745392;
  1026.         zb = zb * epsilon + 0.49999999999999994;
  1027.         zb *= epsilon;
  1028.         zb *= epsilon;

  1030.         double za = epsilon;
  1031.         temp = za + zb;
  1032.         zb = -(temp - za - zb);
  1033.         za = temp;

  1035.         temp = za * HEX_40000000;
  1036.         temp = za + temp - temp;
  1037.         zb += za - temp;
  1038.         za = temp;

  1040.         /* Combine the parts.   expm1(a+b) = expm1(a) + expm1(b) + expm1(a)*expm1(b) */
  1041.         double ya = za * baseA;
  1042.         //double yb = za*baseB + zb*baseA + zb*baseB;
  1043.         temp = ya + za * baseB;
  1044.         double yb = -(temp - ya - za * baseB);
  1045.         ya = temp;

  1047.         temp = ya + zb * baseA;
  1048.         yb += -(temp - ya - zb * baseA);
  1049.         ya = temp;

  1051.         temp = ya + zb * baseB;
  1052.         yb += -(temp - ya - zb * baseB);
  1053.         ya = temp;

  1055.         //ya = ya + za + baseA;
  1056.         //yb = yb + zb + baseB;
  1057.         temp = ya + baseA;
  1058.         yb += -(temp - baseA - ya);
  1059.         ya = temp;

  1061.         temp = ya + za;
  1062.         //yb += (ya > za) ? -(temp - ya - za) : -(temp - za - ya);
  1063.         yb += -(temp - ya - za);
  1064.         ya = temp;

  1066.         temp = ya + baseB;
  1067.         //yb += (ya > baseB) ? -(temp - ya - baseB) : -(temp - baseB - ya);
  1068.         yb += -(temp - ya - baseB);
  1069.         ya = temp;

  1071.         temp = ya + zb;
  1072.         //yb += (ya > zb) ? -(temp - ya - zb) : -(temp - zb - ya);
  1073.         yb += -(temp - ya - zb);
  1074.         ya = temp;

  1076.         if (negative) {
  1077.             /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
  1078.             double denom = 1.0 + ya;
  1079.             double denomr = 1.0 / denom;
  1080.             double denomb = -(denom - 1.0 - ya) + yb;
  1081.             double ratio = ya * denomr;
  1082.             temp = ratio * HEX_40000000;
  1083.             final double ra = ratio + temp - temp;
  1084.             double rb = ratio - ra;

  1086.             temp = denom * HEX_40000000;
  1087.             za = denom + temp - temp;
  1088.             zb = denom - za;

  1090.             rb += (ya - za * ra - za * rb - zb * ra - zb * rb) * denomr;

  1092.             // f(x) = x/1+x
  1093.             // Compute f'(x)
  1094.             // Product rule:  d(uv) = du*v + u*dv
  1095.             // Chain rule:  d(f(g(x)) = f'(g(x))*f(g'(x))
  1096.             // d(1/x) = -1/(x*x)
  1097.             // d(1/1+x) = -1/( (1+x)^2) *  1 =  -1/((1+x)*(1+x))
  1098.             // d(x/1+x) = -x/((1+x)(1+x)) + 1/1+x = 1 / ((1+x)(1+x))

  1100.             // Adjust for yb
  1101.             rb += yb * denomr;                      // numerator
  1102.             rb += -ya * denomb * denomr * denomr;   // denominator

  1104.             // negate
  1105.             ya = -ra;
  1106.             yb = -rb;
  1107.         }

  1109.         if (hiPrecOut != null) {
  1110.             hiPrecOut[0] = ya;
  1111.             hiPrecOut[1] = yb;
  1112.         }

  1114.         return ya + yb;
  1115.     }

  1117.     /**
  1118.      * Natural logarithm.
  1119.      *
  1120.      * @param x   a double
  1121.      * @return log(x)
  1122.      */
  1123.     public static double log(final double x) {
  1124.         return log(x, null);
  1125.     }

  1127.     /**
  1128.      * Internal helper method for natural logarithm function.
  1129.      * @param x original argument of the natural logarithm function
  1130.      * @param hiPrec extra bits of precision on output (To Be Confirmed)
  1131.      * @return log(x)
  1132.      */
  1133.     private static double log(final double x, final double[] hiPrec) {
  1134.         if (x == 0) { // Handle special case of +0/-0
  1135.             return Double.NEGATIVE_INFINITY;
  1136.         }
  1137.         long bits = Double.doubleToRawLongBits(x);

  1139.         /* Handle special cases of negative input, and NaN */
  1140.         if (((bits & 0x8000000000000000L) != 0 || Double.isNaN(x)) && x != 0.0) {
  1141.             if (hiPrec != null) {
  1142.                 hiPrec[0] = Double.NaN;
  1143.             }

  1145.             return Double.NaN;
  1146.         }

  1148.         /* Handle special cases of Positive infinity. */
  1149.         if (x == Double.POSITIVE_INFINITY) {
  1150.             if (hiPrec != null) {
  1151.                 hiPrec[0] = Double.POSITIVE_INFINITY;
  1152.             }

  1154.             return Double.POSITIVE_INFINITY;
  1155.         }

  1157.         /* Extract the exponent */
  1158.         int exp = (int) (bits >> 52) - 1023;

  1160.         if ((bits & 0x7ff0000000000000L) == 0) {
  1161.             // Subnormal!
  1162.             if (x == 0) {
  1163.                 // Zero
  1164.                 if (hiPrec != null) {
  1165.                     hiPrec[0] = Double.NEGATIVE_INFINITY;
  1166.                 }

  1168.                 return Double.NEGATIVE_INFINITY;
  1169.             }

  1171.             /* Normalize the subnormal number. */
  1172.             bits <<= 1;
  1173.             while ((bits & 0x0010000000000000L) == 0) {
  1174.                 --exp;
  1175.                 bits <<= 1;
  1176.             }
  1177.         }


  1180.         if ((exp == -1 || exp == 0) && x < 1.01 && x > 0.99 && hiPrec == null) {
  1181.             /* The normal method doesn't work well in the range [0.99, 1.01], so call do a straight
  1182.            polynomial expansion in higer precision. */

  1184.             /* Compute x - 1.0 and split it */
  1185.             double xa = x - 1.0;
  1186.             double xb = xa - x + 1.0;
  1187.             double tmp = xa * HEX_40000000;
  1188.             double aa = xa + tmp - tmp;
  1189.             double ab = xa - aa;
  1190.             xa = aa;
  1191.             xb = ab;

  1193.             final double[] lnCoef_last = LN_QUICK_COEF[LN_QUICK_COEF.length - 1];
  1194.             double ya = lnCoef_last[0];
  1195.             double yb = lnCoef_last[1];

  1197.             for (int i = LN_QUICK_COEF.length - 2; i >= 0; i--) {
  1198.                 /* Multiply a = y * x */
  1199.                 aa = ya * xa;
  1200.                 ab = ya * xb + yb * xa + yb * xb;
  1201.                 /* split, so now y = a */
  1202.                 tmp = aa * HEX_40000000;
  1203.                 ya = aa + tmp - tmp;
  1204.                 yb = aa - ya + ab;

  1206.                 /* Add  a = y + lnQuickCoef */
  1207.                 final double[] lnCoef_i = LN_QUICK_COEF[i];
  1208.                 aa = ya + lnCoef_i[0];
  1209.                 ab = yb + lnCoef_i[1];
  1210.                 /* Split y = a */
  1211.                 tmp = aa * HEX_40000000;
  1212.                 ya = aa + tmp - tmp;
  1213.                 yb = aa - ya + ab;
  1214.             }

  1216.             /* Multiply a = y * x */
  1217.             aa = ya * xa;
  1218.             ab = ya * xb + yb * xa + yb * xb;
  1219.             /* split, so now y = a */
  1220.             tmp = aa * HEX_40000000;
  1221.             ya = aa + tmp - tmp;
  1222.             yb = aa - ya + ab;

  1224.             return ya + yb;
  1225.         }

  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)];

  1230.         /*
  1231.     double epsilon = x / Double.longBitsToDouble(bits & 0xfffffc0000000000L);

  1233.     epsilon -= 1.0;
  1234.          */

  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));

  1241.         double lnza = 0.0;
  1242.         double lnzb = 0.0;

  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;

  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;

  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];

  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;

  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.             }

  1282.             /* Multiply a = y * x */
  1283.             aa = ya * xa;
  1284.             ab = ya * xb + yb * xa + yb * xb;

  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 *= epsilon;
  1304.         }

  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.          */

  1315.         /* Compute the following sum:
  1316.          * lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
  1317.          */

  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 += d;

  1327.         c = a + lnza;
  1328.         d = -(c - a - lnza);
  1329.         a = c;
  1330.         b += d;

  1332.         c = a + LN_2_B * exp;
  1333.         d = -(c - a - LN_2_B * exp);
  1334.         a = c;
  1335.         b += d;

  1337.         c = a + lnm[1];
  1338.         d = -(c - a - lnm[1]);
  1339.         a = c;
  1340.         b += d;

  1342.         c = a + lnzb;
  1343.         d = -(c - a - lnzb);
  1344.         a = c;
  1345.         b += d;

  1347.         if (hiPrec != null) {
  1348.             hiPrec[0] = a;
  1349.             hiPrec[1] = b;
  1350.         }

  1352.         return a + b;
  1353.     }

  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.         }

  1366.         if (x == Double.POSITIVE_INFINITY) {
  1367.             return Double.POSITIVE_INFINITY;
  1368.         }

  1370.         if (x > 1e-6 ||
  1371.             x < -1e-6) {
  1372.             final double xpa = 1 + x;
  1373.             final double xpb = -(xpa - 1 - x);

  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.             }

  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.     }

  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];

  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.         }

  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];

  1409.         final double rln10a = 0.4342944622039795;
  1410.         final double rln10b = 1.9699272335463627E-8;

  1412.         return rln10b * lnb + rln10b * lna + rln10a * lnb + rln10a * lna;
  1413.     }

  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 3.0
  1430.      */
  1431.     public static double log(double base, double x) {
  1432.         return log(x) / log(base);
  1433.     }

  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(final double x, final double y) {

  1444.         if (y == 0) {
  1445.             // y = -0 or y = +0
  1446.             return 1.0;
  1447.         } else {

  1449.             final long yBits        = Double.doubleToRawLongBits(y);
  1450.             final int  yRawExp      = (int) ((yBits & MASK_DOUBLE_EXPONENT) >> 52);
  1451.             final long yRawMantissa = yBits & MASK_DOUBLE_MANTISSA;
  1452.             final long xBits        = Double.doubleToRawLongBits(x);
  1453.             final int  xRawExp      = (int) ((xBits & MASK_DOUBLE_EXPONENT) >> 52);
  1454.             final long xRawMantissa = xBits & MASK_DOUBLE_MANTISSA;

  1456.             if (yRawExp > 1085) {
  1457.                 // y is either a very large integral value that does not fit in a long or it is a special number

  1459.                 if ((yRawExp == 2047 && yRawMantissa != 0) ||
  1460.                     (xRawExp == 2047 && xRawMantissa != 0)) {
  1461.                     // NaN
  1462.                     return Double.NaN;
  1463.                 } else if (xRawExp == 1023 && xRawMantissa == 0) {
  1464.                     // x = -1.0 or x = +1.0
  1465.                     if (yRawExp == 2047) {
  1466.                         // y is infinite
  1467.                         return Double.NaN;
  1468.                     } else {
  1469.                         // y is a large even integer
  1470.                         return 1.0;
  1471.                     }
  1472.                 } else {
  1473.                     // the absolute value of x is either greater or smaller than 1.0

  1475.                     // if yRawExp == 2047 and mantissa is 0, y = -infinity or y = +infinity
  1476.                     // if 1085 < yRawExp < 2047, y is simply a large number, however, due to limited
  1477.                     // accuracy, at this magnitude it behaves just like infinity with regards to x
  1478.                     if ((y > 0) ^ (xRawExp < 1023)) {
  1479.                         // either y = +infinity (or large engouh) and abs(x) > 1.0
  1480.                         // or     y = -infinity (or large engouh) and abs(x) < 1.0
  1481.                         return Double.POSITIVE_INFINITY;
  1482.                     } else {
  1483.                         // either y = +infinity (or large engouh) and abs(x) < 1.0
  1484.                         // or     y = -infinity (or large engouh) and abs(x) > 1.0
  1485.                         return +0.0;
  1486.                     }
  1487.                 }
  1488.             } else {
  1489.                 // y is a regular non-zero number

  1491.                 if (yRawExp >= 1023) {
  1492.                     // y may be an integral value, which should be handled specifically
  1493.                     final long yFullMantissa = IMPLICIT_HIGH_BIT | yRawMantissa;
  1494.                     if (yRawExp < 1075) {
  1495.                         // normal number with negative shift that may have a fractional part
  1496.                         final long integralMask = (-1L) << (1075 - yRawExp);
  1497.                         if ((yFullMantissa & integralMask) == yFullMantissa) {
  1498.                             // all fractional bits are 0, the number is really integral
  1499.                             final long l = yFullMantissa >> (1075 - yRawExp);
  1500.                             return AccurateMath.pow(x, (y < 0) ? -l : l);
  1501.                         }
  1502.                     } else {
  1503.                         // normal number with positive shift, always an integral value
  1504.                         // we know it fits in a primitive long because yRawExp > 1085 has been handled above
  1505.                         final long l =  yFullMantissa << (yRawExp - 1075);
  1506.                         return AccurateMath.pow(x, (y < 0) ? -l : l);
  1507.                     }
  1508.                 }

  1510.                 // y is a non-integral value

  1512.                 if (x == 0) {
  1513.                     // x = -0 or x = +0
  1514.                     // the integer powers have already been handled above
  1515.                     return y < 0 ? Double.POSITIVE_INFINITY : +0.0;
  1516.                 } else if (xRawExp == 2047) {
  1517.                     if (xRawMantissa == 0) {
  1518.                         // x = -infinity or x = +infinity
  1519.                         return (y < 0) ? +0.0 : Double.POSITIVE_INFINITY;
  1520.                     } else {
  1521.                         // NaN
  1522.                         return Double.NaN;
  1523.                     }
  1524.                 } else if (x < 0) {
  1525.                     // the integer powers have already been handled above
  1526.                     return Double.NaN;
  1527.                 } else {

  1529.                     // this is the general case, for regular fractional numbers x and y

  1531.                     // Split y into ya and yb such that y = ya+yb
  1532.                     final double tmp = y * HEX_40000000;
  1533.                     final double ya = (y + tmp) - tmp;
  1534.                     final double yb = y - ya;

  1536.                     /* Compute ln(x) */
  1537.                     final double[] lns = new double[2];
  1538.                     final double lores = log(x, lns);
  1539.                     if (Double.isInfinite(lores)) { // don't allow this to be converted to NaN
  1540.                         return lores;
  1541.                     }

  1543.                     double lna = lns[0];
  1544.                     double lnb = lns[1];

  1546.                     /* resplit lns */
  1547.                     final double tmp1 = lna * HEX_40000000;
  1548.                     final double tmp2 = (lna + tmp1) - tmp1;
  1549.                     lnb += lna - tmp2;
  1550.                     lna = tmp2;

  1552.                     // y*ln(x) = (aa+ab)
  1553.                     final double aa = lna * ya;
  1554.                     final double ab = lna * yb + lnb * ya + lnb * yb;

  1556.                     lna = aa + ab;
  1557.                     lnb = -(lna - aa - ab);

  1559.                     double z = 1.0 / 120.0;
  1560.                     z = z * lnb + (1.0 / 24.0);
  1561.                     z = z * lnb + (1.0 / 6.0);
  1562.                     z = z * lnb + 0.5;
  1563.                     z = z * lnb + 1.0;
  1564.                     z *= lnb;

  1566.                     final double result = exp(lna, z, null);
  1567.                     //result = result + result * z;
  1568.                     return result;
  1569.                 }
  1570.             }
  1571.         }
  1572.     }

  1574.     /**
  1575.      * Raise a double to an int power.
  1576.      *
  1577.      * @param d Number to raise.
  1578.      * @param e Exponent.
  1579.      * @return d<sup>e</sup>
  1580.      * @since 3.1
  1581.      */
  1582.     public static double pow(double d, int e) {
  1583.         return pow(d, (long) e);
  1584.     }


  1587.     /**
  1588.      * Raise a double to a long power.
  1589.      *
  1590.      * @param d Number to raise.
  1591.      * @param e Exponent.
  1592.      * @return d<sup>e</sup>
  1593.      * @since 3.6
  1594.      */
  1595.     public static double pow(double d, long e) {
  1596.         if (e == 0) {
  1597.             return 1.0;
  1598.         } else if (e > 0) {
  1599.             return new Split(d).pow(e).full;
  1600.         } else {
  1601.             return new Split(d).reciprocal().pow(-e).full;
  1602.         }
  1603.     }

  1605.     /** Class operator on double numbers split into one 26 bits number and one 27 bits number. */
  1606.     private static class Split {

  1608.         /** Split version of NaN. */
  1609.         public static final Split NAN = new Split(Double.NaN, 0);

  1611.         /** Split version of positive infinity. */
  1612.         public static final Split POSITIVE_INFINITY = new Split(Double.POSITIVE_INFINITY, 0);

  1614.         /** Split version of negative infinity. */
  1615.         public static final Split NEGATIVE_INFINITY = new Split(Double.NEGATIVE_INFINITY, 0);

  1617.         /** Full number. */
  1618.         private final double full;

  1620.         /** High order bits. */
  1621.         private final double high;

  1623.         /** Low order bits. */
  1624.         private final double low;

  1626.         /** Simple constructor.
  1627.          * @param x number to split
  1628.          */
  1629.         Split(final double x) {
  1630.             full = x;
  1631.             high = Double.longBitsToDouble(Double.doubleToRawLongBits(x) & ((-1L) << 27));
  1632.             low  = x - high;
  1633.         }

  1635.         /** Simple constructor.
  1636.          * @param high high order bits
  1637.          * @param low low order bits
  1638.          */
  1639.         Split(final double high, final double low) {
  1640.             this(high == 0.0 ? (low == 0.0 && Double.doubleToRawLongBits(high) == Long.MIN_VALUE /* negative zero */ ? -0.0 : low) : high + low, high, low);
  1641.         }

  1643.         /** Simple constructor.
  1644.          * @param full full number
  1645.          * @param high high order bits
  1646.          * @param low low order bits
  1647.          */
  1648.         Split(final double full, final double high, final double low) {
  1649.             this.full = full;
  1650.             this.high = high;
  1651.             this.low  = low;
  1652.         }

  1654.         /** Multiply the instance by another one.
  1655.          * @param b other instance to multiply by
  1656.          * @return product
  1657.          */
  1658.         public Split multiply(final Split b) {
  1659.             // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
  1660.             final Split  mulBasic  = new Split(full * b.full);
  1661.             final double mulError  = low * b.low - (((mulBasic.full - high * b.high) - low * b.high) - high * b.low);
  1662.             return new Split(mulBasic.high, mulBasic.low + mulError);
  1663.         }

  1665.         /** Compute the reciprocal of the instance.
  1666.          * @return reciprocal of the instance
  1667.          */
  1668.         public Split reciprocal() {

  1670.             final double approximateInv = 1.0 / full;
  1671.             final Split  splitInv       = new Split(approximateInv);

  1673.             // if 1.0/d were computed perfectly, remultiplying it by d should give 1.0
  1674.             // we want to estimate the error so we can fix the low order bits of approximateInvLow
  1675.             // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
  1676.             final Split product = multiply(splitInv);
  1677.             final double error  = (product.high - 1) + product.low;

  1679.             // better accuracy estimate of reciprocal
  1680.             return Double.isNaN(error) ? splitInv : new Split(splitInv.high, splitInv.low - error / full);
  1681.         }

  1683.         /** Computes this^e.
  1684.          * @param e exponent (beware, here it MUST be > 0; the only exclusion is Long.MIN_VALUE)
  1685.          * @return d^e, split in high and low bits
  1686.          * @since 3.6
  1687.          */
  1688.         private Split pow(final long e) {

  1690.             // prepare result
  1691.             Split result = new Split(1);

  1693.             // d^(2p)
  1694.             Split d2p = new Split(full, high, low);

  1696.             for (long p = e; p != 0; p >>>= 1) {

  1698.                 if ((p & 0x1) != 0) {
  1699.                     // accurate multiplication result = result * d^(2p) using Veltkamp TwoProduct algorithm
  1700.                     result = result.multiply(d2p);
  1701.                 }

  1703.                 // accurate squaring d^(2(p+1)) = d^(2p) * d^(2p) using Veltkamp TwoProduct algorithm
  1704.                 d2p = d2p.multiply(d2p);
  1705.             }

  1707.             if (Double.isNaN(result.full)) {
  1708.                 if (Double.isNaN(full)) {
  1709.                     return Split.NAN;
  1710.                 } else {
  1711.                     // some intermediate numbers exceeded capacity,
  1712.                     // and the low order bits became NaN (because infinity - infinity = NaN)
  1713.                     if (AccurateMath.abs(full) < 1) {
  1714.                         return new Split(AccurateMath.copySign(0.0, full), 0.0);
  1715.                     } else if (full < 0 && (e & 0x1) == 1) {
  1716.                         return Split.NEGATIVE_INFINITY;
  1717.                     } else {
  1718.                         return Split.POSITIVE_INFINITY;
  1719.                     }
  1720.                 }
  1721.             } else {
  1722.                 return result;
  1723.             }
  1724.         }
  1725.     }

  1727.     /**
  1728.      *  Computes sin(x) - x, where |x| < 1/16.
  1729.      *  Use a Remez polynomial approximation.
  1730.      *  @param x a number smaller than 1/16
  1731.      *  @return sin(x) - x
  1732.      */
  1733.     private static double polySine(final double x) {
  1734.         double x2 = x * x;

  1736.         double p = 2.7553817452272217E-6;
  1737.         p = p * x2 + -1.9841269659586505E-4;
  1738.         p = p * x2 + 0.008333333333329196;
  1739.         p = p * x2 + -0.16666666666666666;
  1740.         //p *= x2;
  1741.         //p *= x;
  1742.         p = p * x2 * x;

  1744.         return p;
  1745.     }

  1747.     /**
  1748.      *  Computes cos(x) - 1, where |x| < 1/16.
  1749.      *  Use a Remez polynomial approximation.
  1750.      *  @param x a number smaller than 1/16
  1751.      *  @return cos(x) - 1
  1752.      */
  1753.     private static double polyCosine(double x) {
  1754.         double x2 = x * x;

  1756.         double p = 2.479773539153719E-5;
  1757.         p = p * x2 + -0.0013888888689039883;
  1758.         p = p * x2 + 0.041666666666621166;
  1759.         p = p * x2 + -0.49999999999999994;
  1760.         p *= x2;

  1762.         return p;
  1763.     }

  1765.     /**
  1766.      *  Compute sine over the first quadrant (0 < x < pi/2).
  1767.      *  Use combination of table lookup and rational polynomial expansion.
  1768.      *  @param xa number from which sine is requested
  1769.      *  @param xb extra bits for x (may be 0.0)
  1770.      *  @return sin(xa + xb)
  1771.      */
  1772.     private static double sinQ(double xa, double xb) {
  1773.         int idx = (int) ((xa * 8.0) + 0.5);
  1774.         final double epsilon = xa - EIGHTHS[idx]; //idx*0.125;

  1776.         // Table lookups
  1777.         final double sintA = SINE_TABLE_A[idx];
  1778.         final double sintB = SINE_TABLE_B[idx];
  1779.         final double costA = COSINE_TABLE_A[idx];
  1780.         final double costB = COSINE_TABLE_B[idx];

  1782.         // Polynomial eval of sin(epsilon), cos(epsilon)
  1783.         double sinEpsA = epsilon;
  1784.         double sinEpsB = polySine(epsilon);
  1785.         final double cosEpsA = 1.0;
  1786.         final double cosEpsB = polyCosine(epsilon);

  1788.         // Split epsilon   xa + xb = x
  1789.         final double temp = sinEpsA * HEX_40000000;
  1790.         double temp2 = (sinEpsA + temp) - temp;
  1791.         sinEpsB +=  sinEpsA - temp2;
  1792.         sinEpsA = temp2;

  1794.         /* Compute sin(x) by angle addition formula */
  1795.         double result;

  1797.         /* Compute the following sum:
  1798.          *
  1799.          * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
  1800.          *          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
  1801.          *
  1802.          * Ranges of elements
  1803.          *
  1804.          * xxxtA   0            PI/2
  1805.          * xxxtB   -1.5e-9      1.5e-9
  1806.          * sinEpsA -0.0625      0.0625
  1807.          * sinEpsB -6e-11       6e-11
  1808.          * cosEpsA  1.0
  1809.          * cosEpsB  0           -0.0625
  1810.          *
  1811.          */

  1813.         //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
  1814.         //          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;

  1816.         //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
  1817.         //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
  1818.         double a = 0;
  1819.         double b = 0;

  1821.         double t = sintA;
  1822.         double c = a + t;
  1823.         double d = -(c - a - t);
  1824.         a = c;
  1825.         b += d;

  1827.         t = costA * sinEpsA;
  1828.         c = a + t;
  1829.         d = -(c - a - t);
  1830.         a = c;
  1831.         b += d;

  1833.         b = b + sintA * cosEpsB + costA * sinEpsB;
  1834.         /*
  1835.         t = sintA*cosEpsB;
  1836.         c = a + t;
  1837.         d = -(c - a - t);
  1838.         a = c;
  1839.         b = b + d;

  1841.         t = costA*sinEpsB;
  1842.         c = a + t;
  1843.         d = -(c - a - t);
  1844.         a = c;
  1845.         b = b + d;
  1846.          */

  1848.         b = b + sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB;
  1849.         /*
  1850.         t = sintB;
  1851.         c = a + t;
  1852.         d = -(c - a - t);
  1853.         a = c;
  1854.         b = b + d;

  1856.         t = costB*sinEpsA;
  1857.         c = a + t;
  1858.         d = -(c - a - t);
  1859.         a = c;
  1860.         b = b + d;

  1862.         t = sintB*cosEpsB;
  1863.         c = a + t;
  1864.         d = -(c - a - t);
  1865.         a = c;
  1866.         b = b + d;

  1868.         t = costB*sinEpsB;
  1869.         c = a + t;
  1870.         d = -(c - a - t);
  1871.         a = c;
  1872.         b = b + d;
  1873.          */

  1875.         if (xb != 0.0) {
  1876.             t = ((costA + costB) * (cosEpsA + cosEpsB) -
  1877.                  (sintA + sintB) * (sinEpsA + sinEpsB)) * xb;  // approximate cosine*xb
  1878.             c = a + t;
  1879.             d = -(c - a - t);
  1880.             a = c;
  1881.             b += d;
  1882.         }

  1884.         result = a + b;

  1886.         return result;
  1887.     }

  1889.     /**
  1890.      * Compute cosine in the first quadrant by subtracting input from PI/2 and
  1891.      * then calling sinQ.  This is more accurate as the input approaches PI/2.
  1892.      *  @param xa number from which cosine is requested
  1893.      *  @param xb extra bits for x (may be 0.0)
  1894.      *  @return cos(xa + xb)
  1895.      */
  1896.     private static double cosQ(double xa, double xb) {
  1897.         final double pi2a = 1.5707963267948966;
  1898.         final double pi2b = 6.123233995736766E-17;

  1900.         final double a = pi2a - xa;
  1901.         double b = -(a - pi2a + xa);
  1902.         b += pi2b - xb;

  1904.         return sinQ(a, b);
  1905.     }

  1907.     /**
  1908.      *  Compute tangent (or cotangent) over the first quadrant.   0 < x < pi/2
  1909.      *  Use combination of table lookup and rational polynomial expansion.
  1910.      *  @param xa number from which sine is requested
  1911.      *  @param xb extra bits for x (may be 0.0)
  1912.      *  @param cotanFlag if true, compute the cotangent instead of the tangent
  1913.      *  @return tan(xa+xb) (or cotangent, depending on cotanFlag)
  1914.      */
  1915.     private static double tanQ(double xa, double xb, boolean cotanFlag) {

  1917.         int idx = (int) ((xa * 8.0) + 0.5);
  1918.         final double epsilon = xa - EIGHTHS[idx]; //idx*0.125;

  1920.         // Table lookups
  1921.         final double sintA = SINE_TABLE_A[idx];
  1922.         final double sintB = SINE_TABLE_B[idx];
  1923.         final double costA = COSINE_TABLE_A[idx];
  1924.         final double costB = COSINE_TABLE_B[idx];

  1926.         // Polynomial eval of sin(epsilon), cos(epsilon)
  1927.         double sinEpsA = epsilon;
  1928.         double sinEpsB = polySine(epsilon);
  1929.         final double cosEpsA = 1.0;
  1930.         final double cosEpsB = polyCosine(epsilon);

  1932.         // Split epsilon   xa + xb = x
  1933.         double temp = sinEpsA * HEX_40000000;
  1934.         double temp2 = (sinEpsA + temp) - temp;
  1935.         sinEpsB +=  sinEpsA - temp2;
  1936.         sinEpsA = temp2;

  1938.         /* Compute sin(x) by angle addition formula */

  1940.         /* Compute the following sum:
  1941.          *
  1942.          * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
  1943.          *          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
  1944.          *
  1945.          * Ranges of elements
  1946.          *
  1947.          * xxxtA   0            PI/2
  1948.          * xxxtB   -1.5e-9      1.5e-9
  1949.          * sinEpsA -0.0625      0.0625
  1950.          * sinEpsB -6e-11       6e-11
  1951.          * cosEpsA  1.0
  1952.          * cosEpsB  0           -0.0625
  1953.          *
  1954.          */

  1956.         //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
  1957.         //          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;

  1959.         //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
  1960.         //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
  1961.         double a = 0;
  1962.         double b = 0;

  1964.         // Compute sine
  1965.         double t = sintA;
  1966.         double c = a + t;
  1967.         double d = -(c - a - t);
  1968.         a = c;
  1969.         b += d;

  1971.         t = costA * sinEpsA;
  1972.         c = a + t;
  1973.         d = -(c - a - t);
  1974.         a = c;
  1975.         b += d;

  1977.         b += sintA * cosEpsB + costA * sinEpsB;
  1978.         b += sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB;

  1980.         double sina = a + b;
  1981.         double sinb = -(sina - a - b);

  1983.         // Compute cosine

  1985.         a = b = c = d = 0.0;

  1987.         t = costA * cosEpsA;
  1988.         c = a + t;
  1989.         d = -(c - a - t);
  1990.         a = c;
  1991.         b += d;

  1993.         t = -sintA * sinEpsA;
  1994.         c = a + t;
  1995.         d = -(c - a - t);
  1996.         a = c;
  1997.         b += d;

  1999.         b += costB * cosEpsA + costA * cosEpsB + costB * cosEpsB;
  2000.         b -= sintB * sinEpsA + sintA * sinEpsB + sintB * sinEpsB;

  2002.         double cosa = a + b;
  2003.         double cosb = -(cosa - a - b);

  2005.         if (cotanFlag) {
  2006.             double tmp;
  2007.             tmp = cosa; cosa = sina; sina = tmp;
  2008.             tmp = cosb; cosb = sinb; sinb = tmp;
  2009.         }


  2012.         /* estimate and correct, compute 1.0/(cosa+cosb) */
  2013.         /*
  2014.     double est = (sina+sinb)/(cosa+cosb);
  2015.     double err = (sina - cosa*est) + (sinb - cosb*est);
  2016.     est += err/(cosa+cosb);
  2017.     err = (sina - cosa*est) + (sinb - cosb*est);
  2018.          */

  2020.         // f(x) = 1/x,   f'(x) = -1/x^2

  2022.         double est = sina / cosa;

  2024.         /* Split the estimate to get more accurate read on division rounding */
  2025.         temp = est * HEX_40000000;
  2026.         double esta = (est + temp) - temp;
  2027.         double estb =  est - esta;

  2029.         temp = cosa * HEX_40000000;
  2030.         double cosaa = (cosa + temp) - temp;
  2031.         double cosab =  cosa - cosaa;

  2033.         //double err = (sina - est*cosa)/cosa;  // Correction for division rounding
  2034.         double err = (sina - esta * cosaa - esta * cosab - estb * cosaa - estb * cosab) / cosa;  // Correction for division rounding
  2035.         err += sinb / cosa;                   // Change in est due to sinb
  2036.         err += -sina * cosb / cosa / cosa;    // Change in est due to cosb

  2038.         if (xb != 0.0) {
  2039.             // tan' = 1 + tan^2      cot' = -(1 + cot^2)
  2040.             // Approximate impact of xb
  2041.             double xbadj = xb + est * est * xb;
  2042.             if (cotanFlag) {
  2043.                 xbadj = -xbadj;
  2044.             }

  2046.             err += xbadj;
  2047.         }

  2049.         return est + err;
  2050.     }

  2052.     /** Reduce the input argument using the Payne and Hanek method.
  2053.      *  This is good for all inputs 0.0 < x < inf
  2054.      *  Output is remainder after dividing by PI/2
  2055.      *  The result array should contain 3 numbers.
  2056.      *  result[0] is the integer portion, so mod 4 this gives the quadrant.
  2057.      *  result[1] is the upper bits of the remainder
  2058.      *  result[2] is the lower bits of the remainder
  2059.      *
  2060.      * @param x number to reduce
  2061.      * @param result placeholder where to put the result
  2062.      */
  2063.     private static void reducePayneHanek(double x, double[] result) {
  2064.         /* Convert input double to bits */
  2065.         long inbits = Double.doubleToRawLongBits(x);
  2066.         int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;

  2068.         /* Convert to fixed point representation */
  2069.         inbits &= 0x000fffffffffffffL;
  2070.         inbits |= 0x0010000000000000L;

  2072.         /* Normalize input to be between 0.5 and 1.0 */
  2073.         exponent++;
  2074.         inbits <<= 11;

  2076.         /* Based on the exponent, get a shifted copy of recip2pi */
  2077.         long shpi0;
  2078.         long shpiA;
  2079.         long shpiB;
  2080.         int idx = exponent >> 6;
  2081.         int shift = exponent - (idx << 6);

  2083.         if (shift != 0) {
  2084.             shpi0 = (idx == 0) ? 0 : (RECIP_2PI[idx - 1] << shift);
  2085.             shpi0 |= RECIP_2PI[idx] >>> (64 - shift);
  2086.             shpiA = (RECIP_2PI[idx] << shift) | (RECIP_2PI[idx + 1] >>> (64 - shift));
  2087.             shpiB = (RECIP_2PI[idx + 1] << shift) | (RECIP_2PI[idx + 2] >>> (64 - shift));
  2088.         } else {
  2089.             shpi0 = (idx == 0) ? 0 : RECIP_2PI[idx - 1];
  2090.             shpiA = RECIP_2PI[idx];
  2091.             shpiB = RECIP_2PI[idx + 1];
  2092.         }

  2094.         /* Multiply input by shpiA */
  2095.         long a = inbits >>> 32;
  2096.         long b = inbits & 0xffffffffL;

  2098.         long c = shpiA >>> 32;
  2099.         long d = shpiA & 0xffffffffL;

  2101.         long ac = a * c;
  2102.         long bd = b * d;
  2103.         long bc = b * c;
  2104.         long ad = a * d;

  2106.         long prodB = bd + (ad << 32);
  2107.         long prodA = ac + (ad >>> 32);

  2109.         boolean bita = (bd & 0x8000000000000000L) != 0;
  2110.         boolean bitb = (ad & 0x80000000L) != 0;
  2111.         boolean bitsum = (prodB & 0x8000000000000000L) != 0;

  2113.         /* Carry */
  2114.         if ((bita && bitb) ||
  2115.                 ((bita || bitb) && !bitsum)) {
  2116.             prodA++;
  2117.         }

  2119.         bita = (prodB & 0x8000000000000000L) != 0;
  2120.         bitb = (bc & 0x80000000L) != 0;

  2122.         prodB += bc << 32;
  2123.         prodA += bc >>> 32;

  2125.         bitsum = (prodB & 0x8000000000000000L) != 0;

  2127.         /* Carry */
  2128.         if ((bita && bitb) ||
  2129.                 ((bita || bitb) && !bitsum)) {
  2130.             prodA++;
  2131.         }

  2133.         /* Multiply input by shpiB */
  2134.         c = shpiB >>> 32;
  2135.         d = shpiB & 0xffffffffL;
  2136.         ac = a * c;
  2137.         bc = b * c;
  2138.         ad = a * d;

  2140.         /* Collect terms */
  2141.         ac += (bc + ad) >>> 32;

  2143.         bita = (prodB & 0x8000000000000000L) != 0;
  2144.         bitb = (ac & 0x8000000000000000L) != 0;
  2145.         prodB += ac;
  2146.         bitsum = (prodB & 0x8000000000000000L) != 0;
  2147.         /* Carry */
  2148.         if ((bita && bitb) ||
  2149.                 ((bita || bitb) && !bitsum)) {
  2150.             prodA++;
  2151.         }

  2153.         /* Multiply by shpi0 */
  2154.         c = shpi0 >>> 32;
  2155.         d = shpi0 & 0xffffffffL;

  2157.         bd = b * d;
  2158.         bc = b * c;
  2159.         ad = a * d;

  2161.         prodA += bd + ((bc + ad) << 32);

  2163.         /*
  2164.          * prodA, prodB now contain the remainder as a fraction of PI.  We want this as a fraction of
  2165.          * PI/2, so use the following steps:
  2166.          * 1.) multiply by 4.
  2167.          * 2.) do a fixed point muliply by PI/4.
  2168.          * 3.) Convert to floating point.
  2169.          * 4.) Multiply by 2
  2170.          */

  2172.         /* This identifies the quadrant */
  2173.         int intPart = (int)(prodA >>> 62);

  2175.         /* Multiply by 4 */
  2176.         prodA <<= 2;
  2177.         prodA |= prodB >>> 62;
  2178.         prodB <<= 2;

  2180.         /* Multiply by PI/4 */
  2181.         a = prodA >>> 32;
  2182.         b = prodA & 0xffffffffL;

  2184.         c = PI_O_4_BITS[0] >>> 32;
  2185.         d = PI_O_4_BITS[0] & 0xffffffffL;

  2187.         ac = a * c;
  2188.         bd = b * d;
  2189.         bc = b * c;
  2190.         ad = a * d;

  2192.         long prod2B = bd + (ad << 32);
  2193.         long prod2A = ac + (ad >>> 32);

  2195.         bita = (bd & 0x8000000000000000L) != 0;
  2196.         bitb = (ad & 0x80000000L) != 0;
  2197.         bitsum = (prod2B & 0x8000000000000000L) != 0;

  2199.         /* Carry */
  2200.         if ((bita && bitb) ||
  2201.                 ((bita || bitb) && !bitsum)) {
  2202.             prod2A++;
  2203.         }

  2205.         bita = (prod2B & 0x8000000000000000L) != 0;
  2206.         bitb = (bc & 0x80000000L) != 0;

  2208.         prod2B += bc << 32;
  2209.         prod2A += bc >>> 32;

  2211.         bitsum = (prod2B & 0x8000000000000000L) != 0;

  2213.         /* Carry */
  2214.         if ((bita && bitb) ||
  2215.                 ((bita || bitb) && !bitsum)) {
  2216.             prod2A++;
  2217.         }

  2219.         /* Multiply input by pio4bits[1] */
  2220.         c = PI_O_4_BITS[1] >>> 32;
  2221.         d = PI_O_4_BITS[1] & 0xffffffffL;
  2222.         ac = a * c;
  2223.         bc = b * c;
  2224.         ad = a * d;

  2226.         /* Collect terms */
  2227.         ac += (bc + ad) >>> 32;

  2229.         bita = (prod2B & 0x8000000000000000L) != 0;
  2230.         bitb = (ac & 0x8000000000000000L) != 0;
  2231.         prod2B += ac;
  2232.         bitsum = (prod2B & 0x8000000000000000L) != 0;
  2233.         /* Carry */
  2234.         if ((bita && bitb) ||
  2235.                 ((bita || bitb) && !bitsum)) {
  2236.             prod2A++;
  2237.         }

  2239.         /* Multiply inputB by pio4bits[0] */
  2240.         a = prodB >>> 32;
  2241.         b = prodB & 0xffffffffL;
  2242.         c = PI_O_4_BITS[0] >>> 32;
  2243.         d = PI_O_4_BITS[0] & 0xffffffffL;
  2244.         ac = a * c;
  2245.         bc = b * c;
  2246.         ad = a * d;

  2248.         /* Collect terms */
  2249.         ac += (bc + ad) >>> 32;

  2251.         bita = (prod2B & 0x8000000000000000L) != 0;
  2252.         bitb = (ac & 0x8000000000000000L) != 0;
  2253.         prod2B += ac;
  2254.         bitsum = (prod2B & 0x8000000000000000L) != 0;
  2255.         /* Carry */
  2256.         if ((bita && bitb) ||
  2257.                 ((bita || bitb) && !bitsum)) {
  2258.             prod2A++;
  2259.         }

  2261.         /* Convert to double */
  2262.         double tmpA = (prod2A >>> 12) / TWO_POWER_52;  // High order 52 bits
  2263.         double tmpB = (((prod2A & 0xfffL) << 40) + (prod2B >>> 24)) / TWO_POWER_52 / TWO_POWER_52; // Low bits

  2265.         double sumA = tmpA + tmpB;
  2266.         double sumB = -(sumA - tmpA - tmpB);

  2268.         /* Multiply by PI/2 and return */
  2269.         result[0] = intPart;
  2270.         result[1] = sumA * 2.0;
  2271.         result[2] = sumB * 2.0;
  2272.     }

  2274.     /**
  2275.      * Sine function.
  2276.      *
  2277.      * @param x Argument.
  2278.      * @return sin(x)
  2279.      */
  2280.     public static double sin(double x) {
  2281.         boolean negative = false;
  2282.         int quadrant = 0;
  2283.         double xa;
  2284.         double xb = 0.0;

  2286.         /* Take absolute value of the input */
  2287.         xa = x;
  2288.         if (x < 0) {
  2289.             negative = true;
  2290.             xa = -xa;
  2291.         }

  2293.         /* Check for zero and negative zero */
  2294.         if (xa == 0.0) {
  2295.             long bits = Double.doubleToRawLongBits(x);
  2296.             if (bits < 0) {
  2297.                 return -0.0;
  2298.             }
  2299.             return 0.0;
  2300.         }

  2302.         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
  2303.             return Double.NaN;
  2304.         }

  2306.         /* Perform any argument reduction */
  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.         }

  2323.         if (negative) {
  2324.             quadrant ^= 2;  // Flip bit 1
  2325.         }

  2327.         switch (quadrant) {
  2328.         case 0:
  2329.             return sinQ(xa, xb);
  2330.         case 1:
  2331.             return cosQ(xa, xb);
  2332.         case 2:
  2333.             return -sinQ(xa, xb);
  2334.         case 3:
  2335.             return -cosQ(xa, xb);
  2336.         default:
  2337.             return Double.NaN;
  2338.         }
  2339.     }

  2341.     /**
  2342.      * Cosine function.
  2343.      *
  2344.      * @param x Argument.
  2345.      * @return cos(x)
  2346.      */
  2347.     public static double cos(double x) {
  2348.         int quadrant = 0;

  2350.         /* Take absolute value of the input */
  2351.         double xa = x;
  2352.         if (x < 0) {
  2353.             xa = -xa;
  2354.         }

  2356.         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
  2357.             return Double.NaN;
  2358.         }

  2360.         /* Perform any argument reduction */
  2361.         double xb = 0;
  2362.         if (xa > 3294198.0) {
  2363.             // PI * (2**20)
  2364.             // Argument too big for CodyWaite reduction.  Must use
  2365.             // PayneHanek.
  2366.             double[] reduceResults = new double[3];
  2367.             reducePayneHanek(xa, reduceResults);
  2368.             quadrant = ((int) reduceResults[0]) & 3;
  2369.             xa = reduceResults[1];
  2370.             xb = reduceResults[2];
  2371.         } else if (xa > 1.5707963267948966) {
  2372.             final CodyWaite cw = new CodyWaite(xa);
  2373.             quadrant = cw.getK() & 3;
  2374.             xa = cw.getRemA();
  2375.             xb = cw.getRemB();
  2376.         }

  2378.         //if (negative)
  2379.         //  quadrant = (quadrant + 2) % 4;

  2381.         switch (quadrant) {
  2382.         case 0:
  2383.             return cosQ(xa, xb);
  2384.         case 1:
  2385.             return -sinQ(xa, xb);
  2386.         case 2:
  2387.             return -cosQ(xa, xb);
  2388.         case 3:
  2389.             return sinQ(xa, xb);
  2390.         default:
  2391.             return Double.NaN;
  2392.         }
  2393.     }

  2395.     /**
  2396.      * Tangent function.
  2397.      *
  2398.      * @param x Argument.
  2399.      * @return tan(x)
  2400.      */
  2401.     public static double tan(double x) {
  2402.         boolean negative = false;
  2403.         int quadrant = 0;

  2405.         /* Take absolute value of the input */
  2406.         double xa = x;
  2407.         if (x < 0) {
  2408.             negative = true;
  2409.             xa = -xa;
  2410.         }

  2412.         /* Check for zero and negative zero */
  2413.         if (xa == 0.0) {
  2414.             long bits = Double.doubleToRawLongBits(x);
  2415.             if (bits < 0) {
  2416.                 return -0.0;
  2417.             }
  2418.             return 0.0;
  2419.         }

  2421.         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
  2422.             return Double.NaN;
  2423.         }

  2425.         /* Perform any argument reduction */
  2426.         double xb = 0;
  2427.         if (xa > 3294198.0) {
  2428.             // PI * (2**20)
  2429.             // Argument too big for CodyWaite reduction.  Must use
  2430.             // PayneHanek.
  2431.             double[] reduceResults = new double[3];
  2432.             reducePayneHanek(xa, reduceResults);
  2433.             quadrant = ((int) reduceResults[0]) & 3;
  2434.             xa = reduceResults[1];
  2435.             xb = reduceResults[2];
  2436.         } else if (xa > 1.5707963267948966) {
  2437.             final CodyWaite cw = new CodyWaite(xa);
  2438.             quadrant = cw.getK() & 3;
  2439.             xa = cw.getRemA();
  2440.             xb = cw.getRemB();
  2441.         }

  2443.         if (xa > 1.5) {
  2444.             // Accuracy suffers between 1.5 and PI/2
  2445.             final double pi2a = 1.5707963267948966;
  2446.             final double pi2b = 6.123233995736766E-17;

  2448.             final double a = pi2a - xa;
  2449.             double b = -(a - pi2a + xa);
  2450.             b += pi2b - xb;

  2452.             xa = a + b;
  2453.             xb = -(xa - a - b);
  2454.             quadrant ^= 1;
  2455.             negative ^= true;
  2456.         }

  2458.         double result;
  2459.         if ((quadrant & 1) == 0) {
  2460.             result = tanQ(xa, xb, false);
  2461.         } else {
  2462.             result = -tanQ(xa, xb, true);
  2463.         }

  2465.         if (negative) {
  2466.             result = -result;
  2467.         }

  2469.         return result;
  2470.     }

  2472.     /**
  2473.      * Arctangent function.
  2474.      *  @param x a number
  2475.      *  @return atan(x)
  2476.      */
  2477.     public static double atan(double x) {
  2478.         return atan(x, 0.0, false);
  2479.     }

  2481.     /** Internal helper function to compute arctangent.
  2482.      * @param xa number from which arctangent is requested
  2483.      * @param xb extra bits for x (may be 0.0)
  2484.      * @param leftPlane if true, result angle must be put in the left half plane
  2485.      * @return atan(xa + xb) (or angle shifted by {@code PI} if leftPlane is true)
  2486.      */
  2487.     private static double atan(double xa, double xb, boolean leftPlane) {
  2488.         if (xa == 0.0) { // Matches +/- 0.0; return correct sign
  2489.             return leftPlane ? copySign(Math.PI, xa) : xa;
  2490.         }

  2492.         final boolean negate;
  2493.         if (xa < 0) {
  2494.             // negative
  2495.             xa = -xa;
  2496.             xb = -xb;
  2497.             negate = true;
  2498.         } else {
  2499.             negate = false;
  2500.         }

  2502.         if (xa > 1.633123935319537E16) { // Very large input
  2503.             return (negate ^ leftPlane) ? (-Math.PI * F_1_2) : (Math.PI * F_1_2);
  2504.         }

  2506.         /* Estimate the closest tabulated arctan value, compute eps = xa-tangentTable */
  2507.         final int idx;
  2508.         if (xa < 1) {
  2509.             idx = (int) (((-1.7168146928204136 * xa * xa + 8.0) * xa) + 0.5);
  2510.         } else {
  2511.             final double oneOverXa = 1 / xa;
  2512.             idx = (int) (-((-1.7168146928204136 * oneOverXa * oneOverXa + 8.0) * oneOverXa) + 13.07);
  2513.         }

  2515.         final double ttA = TANGENT_TABLE_A[idx];
  2516.         final double ttB = TANGENT_TABLE_B[idx];

  2518.         double epsA = xa - ttA;
  2519.         double epsB = -(epsA - xa + ttA);
  2520.         epsB += xb - ttB;

  2522.         double temp = epsA + epsB;
  2523.         epsB = -(temp - epsA - epsB);
  2524.         epsA = temp;

  2526.         /* Compute eps = eps / (1.0 + xa*tangent) */
  2527.         temp = xa * HEX_40000000;
  2528.         double ya = xa + temp - temp;
  2529.         double yb = xb + xa - ya;
  2530.         xa = ya;
  2531.         xb += yb;

  2533.         //if (idx > 8 || idx == 0)
  2534.         if (idx == 0) {
  2535.             /* If the slope of the arctan is gentle enough (< 0.45), this approximation will suffice */
  2536.             //double denom = 1.0 / (1.0 + xa*tangentTableA[idx] + xb*tangentTableA[idx] + xa*tangentTableB[idx] + xb*tangentTableB[idx]);
  2537.             final double denom = 1d / (1d + (xa + xb) * (ttA + ttB));
  2538.             //double denom = 1.0 / (1.0 + xa*tangentTableA[idx]);
  2539.             ya = epsA * denom;
  2540.             yb = epsB * denom;
  2541.         } else {
  2542.             double temp2 = xa * ttA;
  2543.             double za = 1d + temp2;
  2544.             double zb = -(za - 1d - temp2);
  2545.             temp2 = xb * ttA + xa * ttB;
  2546.             temp = za + temp2;
  2547.             zb += -(temp - za - temp2);
  2548.             za = temp;

  2550.             zb += xb * ttB;
  2551.             ya = epsA / za;

  2553.             temp = ya * HEX_40000000;
  2554.             final double yaa = (ya + temp) - temp;
  2555.             final double yab = ya - yaa;

  2557.             temp = za * HEX_40000000;
  2558.             final double zaa = (za + temp) - temp;
  2559.             final double zab = za - zaa;

  2561.             /* Correct for rounding in division */
  2562.             yb = (epsA - yaa * zaa - yaa * zab - yab * zaa - yab * zab) / za;

  2564.             yb += -epsA * zb / za / za;
  2565.             yb += epsB / za;
  2566.         }


  2569.         epsA = ya;
  2570.         epsB = yb;

  2572.         /* Evaluate polynomial */
  2573.         final double epsA2 = epsA * epsA;

  2575.         /*
  2576.     yb = -0.09001346640161823;
  2577.     yb = yb * epsA2 + 0.11110718400605211;
  2578.     yb = yb * epsA2 + -0.1428571349122913;
  2579.     yb = yb * epsA2 + 0.19999999999273194;
  2580.     yb = yb * epsA2 + -0.33333333333333093;
  2581.     yb = yb * epsA2 * epsA;
  2582.          */

  2584.         yb = 0.07490822288864472;
  2585.         yb = yb * epsA2 - 0.09088450866185192;
  2586.         yb = yb * epsA2 + 0.11111095942313305;
  2587.         yb = yb * epsA2 - 0.1428571423679182;
  2588.         yb = yb * epsA2 + 0.19999999999923582;
  2589.         yb = yb * epsA2 - 0.33333333333333287;
  2590.         yb = yb * epsA2 * epsA;


  2593.         ya = epsA;

  2595.         temp = ya + yb;
  2596.         yb = -(temp - ya - yb);
  2597.         ya = temp;

  2599.         /* Add in effect of epsB.   atan'(x) = 1/(1+x^2) */
  2600.         yb += epsB / (1d + epsA * epsA);

  2602.         final double eighths = EIGHTHS[idx];

  2604.         //result = yb + eighths[idx] + ya;
  2605.         double za = eighths + ya;
  2606.         double zb = -(za - eighths - ya);
  2607.         temp = za + yb;
  2608.         zb += -(temp - za - yb);
  2609.         za = temp;

  2611.         double result = za + zb;

  2613.         if (leftPlane) {
  2614.             // Result is in the left plane
  2615.             final double resultb = -(result - za - zb);
  2616.             final double pia = 1.5707963267948966 * 2;
  2617.             final double pib = 6.123233995736766E-17 * 2;

  2619.             za = pia - result;
  2620.             zb = -(za - pia + result);
  2621.             zb += pib - resultb;

  2623.             result = za + zb;
  2624.         }


  2627.         if (negate ^ leftPlane) {
  2628.             result = -result;
  2629.         }

  2631.         return result;
  2632.     }

  2634.     /**
  2635.      * Two arguments arctangent function.
  2636.      * @param y ordinate
  2637.      * @param x abscissa
  2638.      * @return phase angle of point (x,y) between {@code -PI} and {@code PI}
  2639.      */
  2640.     public static double atan2(double y, double x) {
  2641.         if (Double.isNaN(x) || Double.isNaN(y)) {
  2642.             return Double.NaN;
  2643.         }

  2645.         if (y == 0) {
  2646.             final double invx = 1d / x;

  2648.             if (invx == 0) { // X is infinite
  2649.                 if (x > 0) {
  2650.                     return y; // return +/- 0.0
  2651.                 }
  2652.                 return copySign(Math.PI, y);
  2653.             }

  2655.             if (x < 0 || invx < 0) {
  2656.                 return copySign(Math.PI, y);
  2657.             }
  2658.             return x * y;
  2659.         }

  2661.         // y cannot now be zero

  2663.         if (y == Double.POSITIVE_INFINITY) {
  2664.             if (x == Double.POSITIVE_INFINITY) {
  2665.                 return Math.PI * F_1_4;
  2666.             }

  2668.             if (x == Double.NEGATIVE_INFINITY) {
  2669.                 return Math.PI * F_3_4;
  2670.             }

  2672.             return Math.PI * F_1_2;
  2673.         }

  2675.         if (y == Double.NEGATIVE_INFINITY) {
  2676.             if (x == Double.POSITIVE_INFINITY) {
  2677.                 return -Math.PI * F_1_4;
  2678.             }

  2680.             if (x == Double.NEGATIVE_INFINITY) {
  2681.                 return -Math.PI * F_3_4;
  2682.             }

  2684.             return -Math.PI * F_1_2;
  2685.         }

  2687.         if (x == Double.POSITIVE_INFINITY) {
  2688.             return y > 0 ? 0d : -0d;
  2689.         }

  2691.         if (x == Double.NEGATIVE_INFINITY) {
  2692.             return y > 0 ? Math.PI : -Math.PI;
  2693.         }

  2695.         // Neither y nor x can be infinite or NAN here

  2697.         if (x == 0) {
  2698.             return y > 0 ? Math.PI * F_1_2 : -Math.PI * F_1_2;
  2699.         }

  2701.         // Compute ratio r = y/x
  2702.         final double r = y / x;
  2703.         if (Double.isInfinite(r)) { // bypass calculations that can create NaN
  2704.             return atan(r, 0, x < 0);
  2705.         }

  2707.         double ra = doubleHighPart(r);
  2708.         double rb = r - ra;

  2710.         // Split x
  2711.         final double xa = doubleHighPart(x);
  2712.         final double xb = x - xa;

  2714.         rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x;

  2716.         final double temp = ra + rb;
  2717.         rb = -(temp - ra - rb);
  2718.         ra = temp;

  2720.         if (ra == 0) { // Fix up the sign so atan works correctly
  2721.             ra = copySign(0d, y);
  2722.         }

  2724.         // Call atan
  2725.         return atan(ra, rb, x < 0);
  2726.     }

  2728.     /** Compute the arc sine of a number.
  2729.      * @param x number on which evaluation is done
  2730.      * @return arc sine of x
  2731.      */
  2732.     public static double asin(double x) {
  2733.         if (Double.isNaN(x)) {
  2734.             return Double.NaN;
  2735.         }

  2737.         if (x > 1.0 || x < -1.0) {
  2738.             return Double.NaN;
  2739.         }

  2741.         if (x == 1.0) {
  2742.             return Math.PI / 2.0;
  2743.         }

  2745.         if (x == -1.0) {
  2746.             return -Math.PI / 2.0;
  2747.         }

  2749.         if (x == 0.0) { // Matches +/- 0.0; return correct sign
  2750.             return x;
  2751.         }

  2753.         /* Compute asin(x) = atan(x/sqrt(1-x*x)) */

  2755.         /* Split x */
  2756.         double temp = x * HEX_40000000;
  2757.         final double xa = x + temp - temp;
  2758.         final double xb = x - xa;

  2760.         /* Square it */
  2761.         double ya = xa * xa;
  2762.         double yb = xa * xb * 2.0 + xb * xb;

  2764.         /* Subtract from 1 */
  2765.         ya = -ya;
  2766.         yb = -yb;

  2768.         double za = 1.0 + ya;
  2769.         double zb = -(za - 1.0 - ya);

  2771.         temp = za + yb;
  2772.         zb += -(temp - za - yb);
  2773.         za = temp;

  2775.         /* Square root */
  2776.         double y;
  2777.         y = Math.sqrt(za);
  2778.         temp = y * HEX_40000000;
  2779.         ya = y + temp - temp;
  2780.         yb = y - ya;

  2782.         /* Extend precision of sqrt */
  2783.         yb += (za - ya * ya - 2 * ya * yb - yb * yb) / (2.0 * y);

  2785.         /* Contribution of zb to sqrt */
  2786.         double dx = zb / (2.0 * y);

  2788.         // Compute ratio r = x/y
  2789.         double r = x / y;
  2790.         temp = r * HEX_40000000;
  2791.         double ra = r + temp - temp;
  2792.         double rb = r - ra;

  2794.         rb += (x - ra * ya - ra * yb - rb * ya - rb * yb) / y; // Correct for rounding in division
  2795.         rb += -x * dx / y / y; // Add in effect additional bits of sqrt.

  2797.         temp = ra + rb;
  2798.         rb = -(temp - ra - rb);
  2799.         ra = temp;

  2801.         return atan(ra, rb, false);
  2802.     }

  2804.     /** Compute the arc cosine of a number.
  2805.      * @param x number on which evaluation is done
  2806.      * @return arc cosine of x
  2807.      */
  2808.     public static double acos(double x) {
  2809.         if (Double.isNaN(x)) {
  2810.             return Double.NaN;
  2811.         }

  2813.         if (x > 1.0 || x < -1.0) {
  2814.             return Double.NaN;
  2815.         }

  2817.         if (x == -1.0) {
  2818.             return Math.PI;
  2819.         }

  2821.         if (x == 1.0) {
  2822.             return 0.0;
  2823.         }

  2825.         if (x == 0) {
  2826.             return Math.PI / 2.0;
  2827.         }

  2829.         /* Compute acos(x) = atan(sqrt(1-x*x)/x) */

  2831.         /* Split x */
  2832.         double temp = x * HEX_40000000;
  2833.         final double xa = x + temp - temp;
  2834.         final double xb = x - xa;

  2836.         /* Square it */
  2837.         double ya = xa * xa;
  2838.         double yb = xa * xb * 2.0 + xb * xb;

  2840.         /* Subtract from 1 */
  2841.         ya = -ya;
  2842.         yb = -yb;

  2844.         double za = 1.0 + ya;
  2845.         double zb = -(za - 1.0 - ya);

  2847.         temp = za + yb;
  2848.         zb += -(temp - za - yb);
  2849.         za = temp;

  2851.         /* Square root */
  2852.         double y = Math.sqrt(za);
  2853.         temp = y * HEX_40000000;
  2854.         ya = y + temp - temp;
  2855.         yb = y - ya;

  2857.         /* Extend precision of sqrt */
  2858.         yb += (za - ya * ya - 2 * ya * yb - yb * yb) / (2.0 * y);

  2860.         /* Contribution of zb to sqrt */
  2861.         yb += zb / (2.0 * y);
  2862.         y = ya + yb;
  2863.         yb = -(y - ya - yb);

  2865.         // Compute ratio r = y/x
  2866.         double r = y / x;

  2868.         // Did r overflow?
  2869.         if (Double.isInfinite(r)) { // x is effectively zero
  2870.             return Math.PI / 2; // so return the appropriate value
  2871.         }

  2873.         double ra = doubleHighPart(r);
  2874.         double rb = r - ra;

  2876.         rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x; // Correct for rounding in division
  2877.         rb += yb / x; // Add in effect additional bits of sqrt.

  2879.         temp = ra + rb;
  2880.         rb = -(temp - ra - rb);
  2881.         ra = temp;

  2883.         return atan(ra, rb, x < 0);
  2884.     }

  2886.     /** Compute the cubic root of a number.
  2887.      * @param x number on which evaluation is done
  2888.      * @return cubic root of x
  2889.      */
  2890.     public static double cbrt(double x) {
  2891.         /* Convert input double to bits */
  2892.         long inbits = Double.doubleToRawLongBits(x);
  2893.         int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
  2894.         boolean subnormal = false;

  2896.         if (exponent == -1023) {
  2897.             if (x == 0) {
  2898.                 return x;
  2899.             }

  2901.             /* Subnormal, so normalize */
  2902.             subnormal = true;
  2903.             x *= 1.8014398509481984E16; // 2^54
  2904.             inbits = Double.doubleToRawLongBits(x);
  2905.             exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
  2906.         }

  2908.         if (exponent == 1024) {
  2909.             // Nan or infinity. Don't care which.
  2910.             return x;
  2911.         }

  2913.         /* Divide the exponent by 3 */
  2914.         int exp3 = exponent / 3;

  2916.         /* p2 will be the nearest power of 2 to x with its exponent divided by 3 */
  2917.         double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) | (long) (((exp3 + 1023) & 0x7ff)) << 52);

  2919.         /* This will be a number between 1 and 2 */
  2920.         final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L);

  2922.         /* Estimate the cube root of mant by polynomial */
  2923.         double est = -0.010714690733195933;
  2924.         est = est * mant + 0.0875862700108075;
  2925.         est = est * mant + -0.3058015757857271;
  2926.         est = est * mant + 0.7249995199969751;
  2927.         est = est * mant + 0.5039018405998233;

  2929.         est *= CBRTTWO[exponent % 3 + 2];

  2931.         // est should now be good to about 15 bits of precision.   Do 2 rounds of
  2932.         // Newton's method to get closer,  this should get us full double precision
  2933.         // Scale down x for the purpose of doing newtons method.  This avoids over/under flows.
  2934.         final double xs = x / (p2 * p2 * p2);
  2935.         est += (xs - est * est * est) / (3 * est * est);
  2936.         est += (xs - est * est * est) / (3 * est * est);

  2938.         // Do one round of Newton's method in extended precision to get the last bitright.
  2939.         double temp = est * HEX_40000000;
  2940.         double ya = est + temp - temp;
  2941.         double yb = est - ya;

  2943.         double za = ya * ya;
  2944.         double zb = ya * yb * 2.0 + yb * yb;
  2945.         temp = za * HEX_40000000;
  2946.         double temp2 = za + temp - temp;
  2947.         zb += za - temp2;
  2948.         za = temp2;

  2950.         zb = za * yb + ya * zb + zb * yb;
  2951.         za *= ya;

  2953.         double na = xs - za;
  2954.         double nb = -(na - xs + za);
  2955.         nb -= zb;

  2957.         est += (na + nb) / (3 * est * est);

  2959.         /* Scale by a power of two, so this is exact. */
  2960.         est *= p2;

  2962.         if (subnormal) {
  2963.             est *= 3.814697265625E-6; // 2^-18
  2964.         }

  2966.         return est;
  2967.     }

  2969.     /**
  2970.      *  Convert degrees to radians, with error of less than 0.5 ULP.
  2971.      *  @param x angle in degrees
  2972.      *  @return x converted into radians
  2973.      */
  2974.     public static double toRadians(double x) {
  2975.         if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
  2976.             return x;
  2977.         }

  2979.         // These are PI/180 split into high and low order bits
  2980.         final double facta = 0.01745329052209854;
  2981.         final double factb = 1.997844754509471E-9;

  2983.         double xa = doubleHighPart(x);
  2984.         double xb = x - xa;

  2986.         double result = xb * factb + xb * facta + xa * factb + xa * facta;
  2987.         if (result == 0) {
  2988.             result *= x; // ensure correct sign if calculation underflows
  2989.         }
  2990.         return result;
  2991.     }

  2993.     /**
  2994.      *  Convert radians to degrees, with error of less than 0.5 ULP.
  2995.      *  @param x angle in radians
  2996.      *  @return x converted into degrees
  2997.      */
  2998.     public static double toDegrees(double x) {
  2999.         if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
  3000.             return x;
  3001.         }

  3003.         // These are 180/PI split into high and low order bits
  3004.         final double facta = 57.2957763671875;
  3005.         final double factb = 3.145894820876798E-6;

  3007.         double xa = doubleHighPart(x);
  3008.         double xb = x - xa;

  3010.         return xb * factb + xb * facta + xa * factb + xa * facta;
  3011.     }

  3013.     /**
  3014.      * Absolute value.
  3015.      * @param x number from which absolute value is requested
  3016.      * @return abs(x)
  3017.      */
  3018.     public static int abs(final int x) {
  3019.         final int i = x >>> 31;
  3020.         return (x ^ (~i + 1)) + i;
  3021.     }

  3023.     /**
  3024.      * Absolute value.
  3025.      * @param x number from which absolute value is requested
  3026.      * @return abs(x)
  3027.      */
  3028.     public static long abs(final long x) {
  3029.         final long l = x >>> 63;
  3030.         // l is one if x negative zero else
  3031.         // ~l+1 is zero if x is positive, -1 if x is negative
  3032.         // x^(~l+1) is x is x is positive, ~x if x is negative
  3033.         // add around
  3034.         return (x ^ (~l + 1)) + l;
  3035.     }

  3037.     /**
  3038.      * Absolute value.
  3039.      * @param x number from which absolute value is requested
  3040.      * @return abs(x)
  3041.      */
  3042.     public static float abs(final float x) {
  3043.         return Float.intBitsToFloat(MASK_NON_SIGN_INT & Float.floatToRawIntBits(x));
  3044.     }

  3046.     /**
  3047.      * Absolute value.
  3048.      * @param x number from which absolute value is requested
  3049.      * @return abs(x)
  3050.      */
  3051.     public static double abs(double x) {
  3052.         return Double.longBitsToDouble(MASK_NON_SIGN_LONG & Double.doubleToRawLongBits(x));
  3053.     }

  3055.     /**
  3056.      * Compute least significant bit (Unit in Last Position) for a number.
  3057.      * @param x number from which ulp is requested
  3058.      * @return ulp(x)
  3059.      */
  3060.     public static double ulp(double x) {
  3061.         if (Double.isInfinite(x)) {
  3062.             return Double.POSITIVE_INFINITY;
  3063.         }
  3064.         return abs(x - Double.longBitsToDouble(Double.doubleToRawLongBits(x) ^ 1));
  3065.     }

  3067.     /**
  3068.      * Compute least significant bit (Unit in Last Position) for a number.
  3069.      * @param x number from which ulp is requested
  3070.      * @return ulp(x)
  3071.      */
  3072.     public static float ulp(float x) {
  3073.         if (Float.isInfinite(x)) {
  3074.             return Float.POSITIVE_INFINITY;
  3075.         }
  3076.         return abs(x - Float.intBitsToFloat(Float.floatToIntBits(x) ^ 1));
  3077.     }

  3079.     /**
  3080.      * Multiply a double number by a power of 2.
  3081.      * @param d number to multiply
  3082.      * @param n power of 2
  3083.      * @return d &times; 2<sup>n</sup>
  3084.      */
  3085.     public static double scalb(final double d, final int n) {

  3087.         // first simple and fast handling when 2^n can be represented using normal numbers
  3088.         if ((n > -1023) && (n < 1024)) {
  3089.             return d * Double.longBitsToDouble(((long) (n + 1023)) << 52);
  3090.         }

  3092.         // handle special cases
  3093.         if (Double.isNaN(d) || Double.isInfinite(d) || (d == 0)) {
  3094.             return d;
  3095.         }
  3096.         if (n < -2098) {
  3097.             return (d > 0) ? 0.0 : -0.0;
  3098.         }
  3099.         if (n > 2097) {
  3100.             return (d > 0) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
  3101.         }

  3103.         // decompose d
  3104.         final long bits = Double.doubleToRawLongBits(d);
  3105.         final long sign = bits & 0x8000000000000000L;
  3106.         int  exponent   = ((int) (bits >>> 52)) & 0x7ff;
  3107.         long mantissa   = bits & 0x000fffffffffffffL;

  3109.         // compute scaled exponent
  3110.         int scaledExponent = exponent + n;

  3112.         if (n < 0) {
  3113.             // we are really in the case n <= -1023
  3114.             if (scaledExponent > 0) {
  3115.                 // both the input and the result are normal numbers, we only adjust the exponent
  3116.                 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
  3117.             } else if (scaledExponent > -53) {
  3118.                 // the input is a normal number and the result is a subnormal number

  3120.                 // recover the hidden mantissa bit
  3121.                 mantissa |= 1L << 52;

  3123.                 // scales down complete mantissa, hence losing least significant bits
  3124.                 final long mostSignificantLostBit = mantissa & (1L << (-scaledExponent));
  3125.                 mantissa >>>= 1 - scaledExponent;
  3126.                 if (mostSignificantLostBit != 0) {
  3127.                     // we need to add 1 bit to round up the result
  3128.                     mantissa++;
  3129.                 }
  3130.                 return Double.longBitsToDouble(sign | mantissa);
  3131.             } else {
  3132.                 // no need to compute the mantissa, the number scales down to 0
  3133.                 return (sign == 0L) ? 0.0 : -0.0;
  3134.             }
  3135.         } else {
  3136.             // we are really in the case n >= 1024
  3137.             if (exponent == 0) {

  3139.                 // the input number is subnormal, normalize it
  3140.                 while ((mantissa >>> 52) != 1) {
  3141.                     mantissa <<= 1;
  3142.                     --scaledExponent;
  3143.                 }
  3144.                 ++scaledExponent;
  3145.                 mantissa &= 0x000fffffffffffffL;

  3147.                 if (scaledExponent < 2047) {
  3148.                     return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
  3149.                 } else {
  3150.                     return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
  3151.                 }
  3152.             } else if (scaledExponent < 2047) {
  3153.                 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
  3154.             } else {
  3155.                 return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
  3156.             }
  3157.         }
  3158.     }

  3160.     /**
  3161.      * Multiply a float number by a power of 2.
  3162.      * @param f number to multiply
  3163.      * @param n power of 2
  3164.      * @return f &times; 2<sup>n</sup>
  3165.      */
  3166.     public static float scalb(final float f, final int n) {

  3168.         // first simple and fast handling when 2^n can be represented using normal numbers
  3169.         if ((n > -127) && (n < 128)) {
  3170.             return f * Float.intBitsToFloat((n + 127) << 23);
  3171.         }

  3173.         // handle special cases
  3174.         if (Float.isNaN(f) || Float.isInfinite(f) || (f == 0f)) {
  3175.             return f;
  3176.         }
  3177.         if (n < -277) {
  3178.             return (f > 0) ? 0.0f : -0.0f;
  3179.         }
  3180.         if (n > 276) {
  3181.             return (f > 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
  3182.         }

  3184.         // decompose f
  3185.         final int bits = Float.floatToIntBits(f);
  3186.         final int sign = bits & 0x80000000;
  3187.         int  exponent  = (bits >>> 23) & 0xff;
  3188.         int mantissa   = bits & 0x007fffff;

  3190.         // compute scaled exponent
  3191.         int scaledExponent = exponent + n;

  3193.         if (n < 0) {
  3194.             // we are really in the case n <= -127
  3195.             if (scaledExponent > 0) {
  3196.                 // both the input and the result are normal numbers, we only adjust the exponent
  3197.                 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
  3198.             } else if (scaledExponent > -24) {
  3199.                 // the input is a normal number and the result is a subnormal number

  3201.                 // recover the hidden mantissa bit
  3202.                 mantissa |= 1 << 23;

  3204.                 // scales down complete mantissa, hence losing least significant bits
  3205.                 final int mostSignificantLostBit = mantissa & (1 << (-scaledExponent));
  3206.                 mantissa >>>= 1 - scaledExponent;
  3207.                 if (mostSignificantLostBit != 0) {
  3208.                     // we need to add 1 bit to round up the result
  3209.                     mantissa++;
  3210.                 }
  3211.                 return Float.intBitsToFloat(sign | mantissa);
  3212.             } else {
  3213.                 // no need to compute the mantissa, the number scales down to 0
  3214.                 return (sign == 0) ? 0.0f : -0.0f;
  3215.             }
  3216.         } else {
  3217.             // we are really in the case n >= 128
  3218.             if (exponent == 0) {

  3220.                 // the input number is subnormal, normalize it
  3221.                 while ((mantissa >>> 23) != 1) {
  3222.                     mantissa <<= 1;
  3223.                     --scaledExponent;
  3224.                 }
  3225.                 ++scaledExponent;
  3226.                 mantissa &= 0x007fffff;

  3228.                 if (scaledExponent < 255) {
  3229.                     return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
  3230.                 } else {
  3231.                     return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
  3232.                 }
  3233.             } else if (scaledExponent < 255) {
  3234.                 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
  3235.             } else {
  3236.                 return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
  3237.             }
  3238.         }
  3239.     }

  3241.     /**
  3242.      * Get the next machine representable number after a number, moving
  3243.      * in the direction of another number.
  3244.      * <p>
  3245.      * The ordering is as follows (increasing):
  3246.      * <ul>
  3247.      * <li>-INFINITY</li>
  3248.      * <li>-MAX_VALUE</li>
  3249.      * <li>-MIN_VALUE</li>
  3250.      * <li>-0.0</li>
  3251.      * <li>+0.0</li>
  3252.      * <li>+MIN_VALUE</li>
  3253.      * <li>+MAX_VALUE</li>
  3254.      * <li>+INFINITY</li>
  3255.      * <li></li>
  3256.      * </ul>
  3257.      * <p>
  3258.      * If arguments compare equal, then the second argument is returned.
  3259.      * <p>
  3260.      * If {@code direction} is greater than {@code d},
  3261.      * the smallest machine representable number strictly greater than
  3262.      * {@code d} is returned; if less, then the largest representable number
  3263.      * strictly less than {@code d} is returned.</p>
  3264.      * <p>
  3265.      * If {@code d} is infinite and direction does not
  3266.      * bring it back to finite numbers, it is returned unchanged.</p>
  3267.      *
  3268.      * @param d base number
  3269.      * @param direction (the only important thing is whether
  3270.      * {@code direction} is greater or smaller than {@code d})
  3271.      * @return the next machine representable number in the specified direction
  3272.      */
  3273.     public static double nextAfter(double d, double direction) {
  3274.         // handling of some important special cases
  3275.         if (Double.isNaN(d) || Double.isNaN(direction)) {
  3276.             return Double.NaN;
  3277.         } else if (d == direction) {
  3278.             return direction;
  3279.         } else if (Double.isInfinite(d)) {
  3280.             return (d < 0) ? -Double.MAX_VALUE : Double.MAX_VALUE;
  3281.         } else if (d == 0) {
  3282.             return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
  3283.         }
  3284.         // special cases MAX_VALUE to infinity and  MIN_VALUE to 0
  3285.         // are handled just as normal numbers
  3286.         // can use raw bits since already dealt with infinity and NaN
  3287.         final long bits = Double.doubleToRawLongBits(d);
  3288.         final long sign = bits & 0x8000000000000000L;
  3289.         if ((direction < d) ^ (sign == 0L)) {
  3290.             return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) + 1));
  3291.         } else {
  3292.             return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) - 1));
  3293.         }
  3294.     }

  3296.     /**
  3297.      * Get the next machine representable number after a number, moving
  3298.      * in the direction of another number.
  3299.      * <p>
  3300.      * The ordering is as follows (increasing):
  3301.      * <ul>
  3302.      * <li>-INFINITY</li>
  3303.      * <li>-MAX_VALUE</li>
  3304.      * <li>-MIN_VALUE</li>
  3305.      * <li>-0.0</li>
  3306.      * <li>+0.0</li>
  3307.      * <li>+MIN_VALUE</li>
  3308.      * <li>+MAX_VALUE</li>
  3309.      * <li>+INFINITY</li>
  3310.      * <li></li>
  3311.      * </ul>
  3312.      * <p>
  3313.      * If arguments compare equal, then the second argument is returned.
  3314.      * <p>
  3315.      * If {@code direction} is greater than {@code f},
  3316.      * the smallest machine representable number strictly greater than
  3317.      * {@code f} is returned; if less, then the largest representable number
  3318.      * strictly less than {@code f} is returned.</p>
  3319.      * <p>
  3320.      * If {@code f} is infinite and direction does not
  3321.      * bring it back to finite numbers, it is returned unchanged.</p>
  3322.      *
  3323.      * @param f base number
  3324.      * @param direction (the only important thing is whether
  3325.      * {@code direction} is greater or smaller than {@code f})
  3326.      * @return the next machine representable number in the specified direction
  3327.      */
  3328.     public static float nextAfter(final float f, final double direction) {

  3330.         // handling of some important special cases
  3331.         if (Double.isNaN(f) || Double.isNaN(direction)) {
  3332.             return Float.NaN;
  3333.         } else if (f == direction) {
  3334.             return (float) direction;
  3335.         } else if (Float.isInfinite(f)) {
  3336.             return (f < 0f) ? -Float.MAX_VALUE : Float.MAX_VALUE;
  3337.         } else if (f == 0f) {
  3338.             return (direction < 0) ? -Float.MIN_VALUE : Float.MIN_VALUE;
  3339.         }
  3340.         // special cases MAX_VALUE to infinity and  MIN_VALUE to 0
  3341.         // are handled just as normal numbers

  3343.         final int bits = Float.floatToIntBits(f);
  3344.         final int sign = bits & 0x80000000;
  3345.         if ((direction < f) ^ (sign == 0)) {
  3346.             return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) + 1));
  3347.         } else {
  3348.             return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) - 1));
  3349.         }
  3350.     }

  3352.     /** Get the largest whole number smaller than x.
  3353.      * @param x number from which floor is requested
  3354.      * @return a double number f such that f is an integer f &lt;= x &lt; f + 1.0
  3355.      */
  3356.     public static double floor(double x) {
  3357.         long y;

  3359.         if (Double.isNaN(x)) {
  3360.             return x;
  3361.         }

  3363.         if (x >= TWO_POWER_52 || x <= -TWO_POWER_52) {
  3364.             return x;
  3365.         }

  3367.         y = (long) x;
  3368.         if (x < 0 && y != x) {
  3369.             y--;
  3370.         }

  3372.         if (y == 0) {
  3373.             return x * y;
  3374.         }

  3376.         return y;
  3377.     }

  3379.     /** Get the smallest whole number larger than x.
  3380.      * @param x number from which ceil is requested
  3381.      * @return a double number c such that c is an integer c - 1.0 &lt; x &lt;= c
  3382.      */
  3383.     public static double ceil(double x) {
  3384.         double y;

  3386.         if (Double.isNaN(x)) {
  3387.             return x;
  3388.         }

  3390.         y = floor(x);
  3391.         if (y == x) {
  3392.             return y;
  3393.         }

  3395.         y += 1.0;

  3397.         if (y == 0) {
  3398.             return x * y;
  3399.         }

  3401.         return y;
  3402.     }

  3404.     /** Get the whole number that is the nearest to x, or the even one if x is exactly half way between two integers.
  3405.      * @param x number from which nearest whole number is requested
  3406.      * @return a double number r such that r is an integer r - 0.5 &lt;= x &lt;= r + 0.5
  3407.      */
  3408.     public static double rint(double x) {
  3409.         double y = floor(x);
  3410.         double d = x - y;

  3412.         if (d > 0.5) {
  3413.             if (y == -1.0) {
  3414.                 return -0.0; // Preserve sign of operand
  3415.             }
  3416.             return y + 1.0;
  3417.         }
  3418.         if (d < 0.5) {
  3419.             return y;
  3420.         }

  3422.         /* half way, round to even */
  3423.         long z = (long) y;
  3424.         return (z & 1) == 0 ? y : y + 1.0;
  3425.     }

  3427.     /** Get the closest long to x.
  3428.      * @param x number from which closest long is requested
  3429.      * @return closest long to x
  3430.      */
  3431.     public static long round(double x) {
  3432.         final long bits = Double.doubleToRawLongBits(x);
  3433.         final int biasedExp = ((int) (bits >> 52)) & 0x7ff;
  3434.         // Shift to get rid of bits past comma except first one: will need to
  3435.         // 1-shift to the right to end up with correct magnitude.
  3436.         final int shift = (52 - 1 + Double.MAX_EXPONENT) - biasedExp;
  3437.         if ((shift & -64) == 0) {
  3438.             // shift in [0,63], so unbiased exp in [-12,51].
  3439.             long extendedMantissa = 0x0010000000000000L | (bits & 0x000fffffffffffffL);
  3440.             if (bits < 0) {
  3441.                 extendedMantissa = -extendedMantissa;
  3442.             }
  3443.             // If value is positive and first bit past comma is 0, rounding
  3444.             // to lower integer, else to upper one, which is what "+1" and
  3445.             // then ">>1" do.
  3446.             return ((extendedMantissa >> shift) + 1L) >> 1;
  3447.         } else {
  3448.             // +-Infinity, NaN, or a mathematical integer.
  3449.             return (long) x;
  3450.         }
  3451.     }

  3453.     /** Get the closest int to x.
  3454.      * @param x number from which closest int is requested
  3455.      * @return closest int to x
  3456.      */
  3457.     public static int round(final float x) {
  3458.         final int bits = Float.floatToRawIntBits(x);
  3459.         final int biasedExp = (bits >> 23) & 0xff;
  3460.         // Shift to get rid of bits past comma except first one: will need to
  3461.         // 1-shift to the right to end up with correct magnitude.
  3462.         final int shift = (23 - 1 + Float.MAX_EXPONENT) - biasedExp;
  3463.         if ((shift & -32) == 0) {
  3464.             // shift in [0,31], so unbiased exp in [-9,22].
  3465.             int extendedMantissa = 0x00800000 | (bits & 0x007fffff);
  3466.             if (bits < 0) {
  3467.                 extendedMantissa = -extendedMantissa;
  3468.             }
  3469.             // If value is positive and first bit past comma is 0, rounding
  3470.             // to lower integer, else to upper one, which is what "+1" and
  3471.             // then ">>1" do.
  3472.             return ((extendedMantissa >> shift) + 1) >> 1;
  3473.         } else {
  3474.             // +-Infinity, NaN, or a mathematical integer.
  3475.             return (int) x;
  3476.         }
  3477.     }

  3479.     /** Compute the minimum of two values.
  3480.      * @param a first value
  3481.      * @param b second value
  3482.      * @return a if a is lesser or equal to b, b otherwise
  3483.      */
  3484.     public static int min(final int a, final int b) {
  3485.         return (a <= b) ? a : b;
  3486.     }

  3488.     /** Compute the minimum of two values.
  3489.      * @param a first value
  3490.      * @param b second value
  3491.      * @return a if a is lesser or equal to b, b otherwise
  3492.      */
  3493.     public static long min(final long a, final long b) {
  3494.         return (a <= b) ? a : b;
  3495.     }

  3497.     /** Compute the minimum of two values.
  3498.      * @param a first value
  3499.      * @param b second value
  3500.      * @return a if a is lesser or equal to b, b otherwise
  3501.      */
  3502.     public static float min(final float a, final float b) {
  3503.         if (a > b) {
  3504.             return b;
  3505.         }
  3506.         if (a < b) {
  3507.             return a;
  3508.         }
  3509.         /* if either arg is NaN, return NaN */
  3510.         if (a != b) {
  3511.             return Float.NaN;
  3512.         }
  3513.         /* min(+0.0,-0.0) == -0.0 */
  3514.         /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
  3515.         int bits = Float.floatToRawIntBits(a);
  3516.         if (bits == 0x80000000) {
  3517.             return a;
  3518.         }
  3519.         return b;
  3520.     }

  3522.     /** Compute the minimum of two values.
  3523.      * @param a first value
  3524.      * @param b second value
  3525.      * @return a if a is lesser or equal to b, b otherwise
  3526.      */
  3527.     public static double min(final double a, final double b) {
  3528.         if (a > b) {
  3529.             return b;
  3530.         }
  3531.         if (a < b) {
  3532.             return a;
  3533.         }
  3534.         /* if either arg is NaN, return NaN */
  3535.         if (a != b) {
  3536.             return Double.NaN;
  3537.         }
  3538.         /* min(+0.0,-0.0) == -0.0 */
  3539.         /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
  3540.         long bits = Double.doubleToRawLongBits(a);
  3541.         if (bits == 0x8000000000000000L) {
  3542.             return a;
  3543.         }
  3544.         return b;
  3545.     }

  3547.     /** Compute the maximum of two values.
  3548.      * @param a first value
  3549.      * @param b second value
  3550.      * @return b if a is lesser or equal to b, a otherwise
  3551.      */
  3552.     public static int max(final int a, final int b) {
  3553.         return (a <= b) ? b : a;
  3554.     }

  3556.     /** Compute the maximum of two values.
  3557.      * @param a first value
  3558.      * @param b second value
  3559.      * @return b if a is lesser or equal to b, a otherwise
  3560.      */
  3561.     public static long max(final long a, final long b) {
  3562.         return (a <= b) ? b : a;
  3563.     }

  3565.     /** Compute the maximum of two values.
  3566.      * @param a first value
  3567.      * @param b second value
  3568.      * @return b if a is lesser or equal to b, a otherwise
  3569.      */
  3570.     public static float max(final float a, final float b) {
  3571.         if (a > b) {
  3572.             return a;
  3573.         }
  3574.         if (a < b) {
  3575.             return b;
  3576.         }
  3577.         /* if either arg is NaN, return NaN */
  3578.         if (a != b) {
  3579.             return Float.NaN;
  3580.         }
  3581.         /* min(+0.0,-0.0) == -0.0 */
  3582.         /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
  3583.         int bits = Float.floatToRawIntBits(a);
  3584.         if (bits == 0x80000000) {
  3585.             return b;
  3586.         }
  3587.         return a;
  3588.     }

  3590.     /** Compute the maximum of two values.
  3591.      * @param a first value
  3592.      * @param b second value
  3593.      * @return b if a is lesser or equal to b, a otherwise
  3594.      */
  3595.     public static double max(final double a, final double b) {
  3596.         if (a > b) {
  3597.             return a;
  3598.         }
  3599.         if (a < b) {
  3600.             return b;
  3601.         }
  3602.         /* if either arg is NaN, return NaN */
  3603.         if (a != b) {
  3604.             return Double.NaN;
  3605.         }
  3606.         /* min(+0.0,-0.0) == -0.0 */
  3607.         /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
  3608.         long bits = Double.doubleToRawLongBits(a);
  3609.         if (bits == 0x8000000000000000L) {
  3610.             return b;
  3611.         }
  3612.         return a;
  3613.     }

  3615.     /**
  3616.      * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
  3617.      * - sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)<br>
  3618.      * avoiding intermediate overflow or underflow.
  3619.      *
  3620.      * <ul>
  3621.      * <li> If either argument is infinite, then the result is positive infinity.</li>
  3622.      * <li> else, if either argument is NaN then the result is NaN.</li>
  3623.      * </ul>
  3624.      *
  3625.      * @param x a value
  3626.      * @param y a value
  3627.      * @return sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
  3628.      */
  3629.     public static double hypot(final double x, final double y) {
  3630.         if (Double.isInfinite(x) || Double.isInfinite(y)) {
  3631.             return Double.POSITIVE_INFINITY;
  3632.         } else if (Double.isNaN(x) || Double.isNaN(y)) {
  3633.             return Double.NaN;
  3634.         } else {
  3635.             final int expX = getExponent(x);
  3636.             final int expY = getExponent(y);
  3637.             if (expX > expY + 27) {
  3638.                 // y is neglectible with respect to x
  3639.                 return abs(x);
  3640.             } else if (expY > expX + 27) {
  3641.                 // x is neglectible with respect to y
  3642.                 return abs(y);
  3643.             } else {
  3644.                 // find an intermediate scale to avoid both overflow and underflow
  3645.                 final int middleExp = (expX + expY) / 2;

  3647.                 // scale parameters without losing precision
  3648.                 final double scaledX = scalb(x, -middleExp);
  3649.                 final double scaledY = scalb(y, -middleExp);

  3651.                 // compute scaled hypotenuse
  3652.                 final double scaledH = Math.sqrt(scaledX * scaledX + scaledY * scaledY);

  3654.                 // remove scaling
  3655.                 return scalb(scaledH, middleExp);
  3656.             }
  3657.         }
  3658.     }

  3660.     /**
  3661.      * Computes the remainder as prescribed by the IEEE 754 standard.
  3662.      * The remainder value is mathematically equal to {@code x - y*n}
  3663.      * where {@code n} is the mathematical integer closest to the exact mathematical value
  3664.      * of the quotient {@code x/y}.
  3665.      * If two mathematical integers are equally close to {@code x/y} then
  3666.      * {@code n} is the integer that is even.
  3667.      * <ul>
  3668.      * <li>If either operand is NaN, the result is NaN.</li>
  3669.      * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
  3670.      * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
  3671.      * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
  3672.      * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
  3673.      * </ul>
  3674.      * @param dividend the number to be divided
  3675.      * @param divisor the number by which to divide
  3676.      * @return the remainder, rounded
  3677.      */
  3678.     public static double IEEEremainder(final double dividend, final double divisor) {
  3679.         if (getExponent(dividend) == 1024 || getExponent(divisor) == 1024 || divisor == 0.0) {
  3680.             // we are in one of the special cases
  3681.             if (Double.isInfinite(divisor) && !Double.isInfinite(dividend)) {
  3682.                 return dividend;
  3683.             } else {
  3684.                 return Double.NaN;
  3685.             }
  3686.         } else {
  3687.             // we are in the general case
  3688.             final double n         = AccurateMath.rint(dividend / divisor);
  3689.             final double remainder = Double.isInfinite(n) ? 0.0 : dividend - divisor * n;
  3690.             return (remainder == 0) ? AccurateMath.copySign(remainder, dividend) : remainder;
  3691.         }
  3692.     }

  3694.     /** Convert a long to integer, detecting overflows.
  3695.      * @param n number to convert to int
  3696.      * @return integer with same value as n if no overflows occur
  3697.      * @exception ArithmeticException if n cannot fit into an int
  3698.      * @since 3.4
  3699.      */
  3700.     public static int toIntExact(final long n) {
  3701.         if (n < Integer.MIN_VALUE ||
  3702.             n > Integer.MAX_VALUE) {
  3703.             throw new ArithmeticException(OVERFLOW_MSG);
  3704.         }
  3705.         return (int) n;
  3706.     }

  3708.     /** Increment a number, detecting overflows.
  3709.      * @param n number to increment
  3710.      * @return n+1 if no overflows occur
  3711.      * @exception ArithmeticException if an overflow occurs
  3712.      * @since 3.4
  3713.      */
  3714.     public static int incrementExact(final int n) {
  3715.         if (n == Integer.MAX_VALUE) {
  3716.             throw new ArithmeticException(OVERFLOW_MSG);
  3717.         }
  3718.         return n + 1;
  3719.     }

  3721.     /** Increment a number, detecting overflows.
  3722.      * @param n number to increment
  3723.      * @return n+1 if no overflows occur
  3724.      * @exception ArithmeticException if an overflow occurs
  3725.      * @since 3.4
  3726.      */
  3727.     public static long incrementExact(final long n) {
  3728.         if (n == Long.MAX_VALUE) {
  3729.             throw new ArithmeticException(OVERFLOW_MSG);
  3730.         }
  3731.         return n + 1;
  3732.     }

  3734.     /** Decrement a number, detecting overflows.
  3735.      * @param n number to decrement
  3736.      * @return n-1 if no overflows occur
  3737.      * @exception ArithmeticException if an overflow occurs
  3738.      * @since 3.4
  3739.      */
  3740.     public static int decrementExact(final int n) {
  3741.         if (n == Integer.MIN_VALUE) {
  3742.             throw new ArithmeticException(OVERFLOW_MSG);
  3743.         }
  3744.         return n - 1;
  3745.     }

  3747.     /** Decrement a number, detecting overflows.
  3748.      * @param n number to decrement
  3749.      * @return n-1 if no overflows occur
  3750.      * @exception ArithmeticException if an overflow occurs
  3751.      * @since 3.4
  3752.      */
  3753.     public static long decrementExact(final long n) {
  3754.         if (n == Long.MIN_VALUE) {
  3755.             throw new ArithmeticException(OVERFLOW_MSG);
  3756.         }
  3757.         return n - 1;
  3758.     }

  3760.     /** Add two numbers, detecting overflows.
  3761.      * @param a first number to add
  3762.      * @param b second number to add
  3763.      * @return a+b if no overflows occur
  3764.      * @exception ArithmeticException if an overflow occurs
  3765.      * @since 3.4
  3766.      */
  3767.     public static int addExact(final int a, final int b) {
  3768.         // compute sum
  3769.         final int sum = a + b;

  3771.         // check for overflow
  3772.         if ((a ^ b) >= 0 &&
  3773.             (sum ^ b) < 0) {
  3774.             throw new ArithmeticException(OVERFLOW_MSG);
  3775.         }

  3777.         return sum;
  3778.     }

  3780.     /** Add two numbers, detecting overflows.
  3781.      * @param a first number to add
  3782.      * @param b second number to add
  3783.      * @return a+b if no overflows occur
  3784.      * @exception ArithmeticException if an overflow occurs
  3785.      * @since 3.4
  3786.      */
  3787.     public static long addExact(final long a, final long b) {
  3788.         // compute sum
  3789.         final long sum = a + b;

  3791.         // check for overflow
  3792.         if ((a ^ b) >= 0 &&
  3793.             (sum ^ b) < 0) {
  3794.             throw new ArithmeticException(OVERFLOW_MSG);
  3795.         }

  3797.         return sum;
  3798.     }

  3800.     /** Subtract two numbers, detecting overflows.
  3801.      * @param a first number
  3802.      * @param b second number to subtract from a
  3803.      * @return a-b if no overflows occur
  3804.      * @exception ArithmeticException if an overflow occurs
  3805.      * @since 3.4
  3806.      */
  3807.     public static int subtractExact(final int a, final int b) {
  3808.         // compute subtraction
  3809.         final int sub = a - b;

  3811.         // check for overflow
  3812.         if ((a ^ b) < 0 &&
  3813.             (sub ^ b) >= 0) {
  3814.             throw new ArithmeticException(OVERFLOW_MSG);
  3815.         }

  3817.         return sub;
  3818.     }

  3820.     /** Subtract two numbers, detecting overflows.
  3821.      * @param a first number
  3822.      * @param b second number to subtract from a
  3823.      * @return a-b if no overflows occur
  3824.      * @exception ArithmeticException if an overflow occurs
  3825.      * @since 3.4
  3826.      */
  3827.     public static long subtractExact(final long a, final long b) {
  3828.         // compute subtraction
  3829.         final long sub = a - b;

  3831.         // check for overflow
  3832.         if ((a ^ b) < 0 &&
  3833.             (sub ^ b) >= 0) {
  3834.             throw new ArithmeticException(OVERFLOW_MSG);
  3835.         }

  3837.         return sub;
  3838.     }

  3840.     /** Multiply two numbers, detecting overflows.
  3841.      * @param a first number to multiply
  3842.      * @param b second number to multiply
  3843.      * @return a*b if no overflows occur
  3844.      * @exception ArithmeticException if an overflow occurs
  3845.      * @since 3.4
  3846.      */
  3847.     public static int multiplyExact(final int a, final int b) {
  3848.         if (((b > 0) &&
  3849.              (a > Integer.MAX_VALUE / b ||
  3850.               a < Integer.MIN_VALUE / b)) ||
  3851.             ((b < -1) &&
  3852.              (a > Integer.MIN_VALUE / b ||
  3853.               a < Integer.MAX_VALUE / b)) ||
  3854.             ((b == -1) &&
  3855.              (a == Integer.MIN_VALUE))) {
  3856.             throw new ArithmeticException(OVERFLOW_MSG);
  3857.         }
  3858.         return a * b;
  3859.     }

  3861.     /** Multiply two numbers, detecting overflows.
  3862.      * @param a first number to multiply
  3863.      * @param b second number to multiply
  3864.      * @return a*b if no overflows occur
  3865.      * @exception ArithmeticException if an overflow occurs
  3866.      * @since 3.4
  3867.      */
  3868.     public static long multiplyExact(final long a, final long b) {
  3869.         if (((b  >  0L) &&
  3870.              (a > Long.MAX_VALUE / b ||
  3871.               a < Long.MIN_VALUE / b)) ||
  3872.             ((b  < -1L) &&
  3873.              (a > Long.MIN_VALUE / b ||
  3874.               a < Long.MAX_VALUE / b)) ||
  3875.             ((b == -1L) &&
  3876.              (a == Long.MIN_VALUE))) {
  3877.             throw new ArithmeticException(OVERFLOW_MSG);
  3878.         }
  3879.         return a * b;
  3880.     }

  3882.     /** Finds q such that a = q b + r with 0 &lt;= r &lt; b if b &gt; 0 and b &lt; r &lt;= 0 if b &lt; 0.
  3883.      * <p>
  3884.      * This methods returns the same value as integer division when
  3885.      * a and b are same signs, but returns a different value when
  3886.      * they are opposite (i.e. q is negative).
  3887.      * </p>
  3888.      * @param a dividend
  3889.      * @param b divisor
  3890.      * @return q such that a = q b + r with 0 &lt;= r &lt; b if b &gt; 0 and b &lt; r &lt;= 0 if b &lt; 0
  3891.      * @exception ArithmeticException if b == 0
  3892.      * @see #floorMod(int, int)
  3893.      * @since 3.4
  3894.      */
  3895.     public static int floorDiv(final int a, final int b) {
  3896.         if (b == 0) {
  3897.             throw new ArithmeticException(ZERO_DENOMINATOR_MSG);
  3898.         }

  3900.         final int m = a % b;
  3901.         if ((a ^ b) >= 0 || m == 0) {
  3902.             // a an b have same sign, or division is exact
  3903.             return a / b;
  3904.         } else {
  3905.             // a and b have opposite signs and division is not exact
  3906.             return (a / b) - 1;
  3907.         }
  3908.     }

  3910.     /** Finds q such that a = q b + r with 0 &lt;= r &lt; b if b &gt; 0 and b &lt; r &lt;= 0 if b &lt; 0.
  3911.      * <p>
  3912.      * This methods returns the same value as integer division when
  3913.      * a and b are same signs, but returns a different value when
  3914.      * they are opposite (i.e. q is negative).
  3915.      * </p>
  3916.      * @param a dividend
  3917.      * @param b divisor
  3918.      * @return q such that a = q b + r with 0 &lt;= r &lt; b if b &gt; 0 and b &lt; r &lt;= 0 if b &lt; 0
  3919.      * @exception ArithmeticException if b == 0
  3920.      * @see #floorMod(long, long)
  3921.      * @since 3.4
  3922.      */
  3923.     public static long floorDiv(final long a, final long b) {
  3924.         if (b == 0L) {
  3925.             throw new ArithmeticException(ZERO_DENOMINATOR_MSG);
  3926.         }

  3928.         final long m = a % b;
  3929.         if ((a ^ b) >= 0L || m == 0L) {
  3930.             // a an b have same sign, or division is exact
  3931.             return a / b;
  3932.         } else {
  3933.             // a and b have opposite signs and division is not exact
  3934.             return (a / b) - 1L;
  3935.         }
  3936.     }

  3938.     /** Finds r such that a = q b + r with 0 &lt;= r &lt; b if b &gt; 0 and b &lt; r &lt;= 0 if b &lt; 0.
  3939.      * <p>
  3940.      * This methods returns the same value as integer modulo when
  3941.      * a and b are same signs, but returns a different value when
  3942.      * they are opposite (i.e. q is negative).
  3943.      * </p>
  3944.      * @param a dividend
  3945.      * @param b divisor
  3946.      * @return r such that a = q b + r with 0 &lt;= r &lt; b if b &gt; 0 and b &lt; r &lt;= 0 if b &lt; 0
  3947.      * @exception ArithmeticException if b == 0
  3948.      * @see #floorDiv(int, int)
  3949.      * @since 3.4
  3950.      */
  3951.     public static int floorMod(final int a, final int b) {
  3952.         if (b == 0) {
  3953.             throw new ArithmeticException(ZERO_DENOMINATOR_MSG);
  3954.         }

  3956.         final int m = a % b;
  3957.         if ((a ^ b) >= 0 || m == 0) {
  3958.             // a an b have same sign, or division is exact
  3959.             return m;
  3960.         } else {
  3961.             // a and b have opposite signs and division is not exact
  3962.             return b + m;
  3963.         }
  3964.     }

  3966.     /** Finds r such that a = q b + r with 0 &lt;= r &lt; b if b &gt; 0 and b &lt; r &lt;= 0 if b &lt; 0.
  3967.      * <p>
  3968.      * This methods returns the same value as integer modulo when
  3969.      * a and b are same signs, but returns a different value when
  3970.      * they are opposite (i.e. q is negative).
  3971.      * </p>
  3972.      * @param a dividend
  3973.      * @param b divisor
  3974.      * @return r such that a = q b + r with 0 &lt;= r &lt; b if b &gt; 0 and b &lt; r &lt;= 0 if b &lt; 0
  3975.      * @exception ArithmeticException if b == 0
  3976.      * @see #floorDiv(long, long)
  3977.      * @since 3.4
  3978.      */
  3979.     public static long floorMod(final long a, final long b) {
  3980.         if (b == 0L) {
  3981.             throw new ArithmeticException(ZERO_DENOMINATOR_MSG);
  3982.         }

  3984.         final long m = a % b;
  3985.         if ((a ^ b) >= 0L || m == 0L) {
  3986.             // a an b have same sign, or division is exact
  3987.             return m;
  3988.         } else {
  3989.             // a and b have opposite signs and division is not exact
  3990.             return b + m;
  3991.         }
  3992.     }

  3994.     /**
  3995.      * Returns the first argument with the sign of the second argument.
  3996.      * A NaN {@code sign} argument is treated as positive.
  3997.      *
  3998.      * @param magnitude the value to return
  3999.      * @param sign the sign for the returned value
  4000.      * @return the magnitude with the same sign as the {@code sign} argument
  4001.      */
  4002.     public static double copySign(double magnitude, double sign) {
  4003.         // The highest order bit is going to be zero if the
  4004.         // highest order bit of m and s is the same and one otherwise.
  4005.         // So (m^s) will be positive if both m and s have the same sign
  4006.         // and negative otherwise.
  4007.         final long m = Double.doubleToRawLongBits(magnitude); // don't care about NaN
  4008.         final long s = Double.doubleToRawLongBits(sign);
  4009.         if ((m ^ s) >= 0) {
  4010.             return magnitude;
  4011.         }
  4012.         return -magnitude; // flip sign
  4013.     }

  4015.     /**
  4016.      * Returns the first argument with the sign of the second argument.
  4017.      * A NaN {@code sign} argument is treated as positive.
  4018.      *
  4019.      * @param magnitude the value to return
  4020.      * @param sign the sign for the returned value
  4021.      * @return the magnitude with the same sign as the {@code sign} argument
  4022.      */
  4023.     public static float copySign(float magnitude, float sign) {
  4024.         // The highest order bit is going to be zero if the
  4025.         // highest order bit of m and s is the same and one otherwise.
  4026.         // So (m^s) will be positive if both m and s have the same sign
  4027.         // and negative otherwise.
  4028.         final int m = Float.floatToRawIntBits(magnitude);
  4029.         final int s = Float.floatToRawIntBits(sign);
  4030.         if ((m ^ s) >= 0) {
  4031.             return magnitude;
  4032.         }
  4033.         return -magnitude; // flip sign
  4034.     }

  4036.     /**
  4037.      * Return the exponent of a double number, removing the bias.
  4038.      * <p>
  4039.      * For double numbers of the form 2<sup>x</sup>, the unbiased
  4040.      * exponent is exactly x.
  4041.      * </p>
  4042.      * @param d number from which exponent is requested
  4043.      * @return exponent for d in IEEE754 representation, without bias
  4044.      */
  4045.     public static int getExponent(final double d) {
  4046.         // NaN and Infinite will return 1024 anywho so can use raw bits
  4047.         return (int) ((Double.doubleToRawLongBits(d) >>> 52) & 0x7ff) - 1023;
  4048.     }

  4050.     /**
  4051.      * Return the exponent of a float number, removing the bias.
  4052.      * <p>
  4053.      * For float numbers of the form 2<sup>x</sup>, the unbiased
  4054.      * exponent is exactly x.
  4055.      * </p>
  4056.      * @param f number from which exponent is requested
  4057.      * @return exponent for d in IEEE754 representation, without bias
  4058.      */
  4059.     public static int getExponent(final float f) {
  4060.         // NaN and Infinite will return the same exponent anywho so can use raw bits
  4061.         return ((Float.floatToRawIntBits(f) >>> 23) & 0xff) - 127;
  4062.     }

  4064.     /**
  4065.      * Print out contents of arrays, and check the length.
  4066.      * <p>used to generate the preset arrays originally.</p>
  4067.      * @param a unused
  4068.      */
  4069.     public static void main(String[] a) {
  4070.         PrintStream out = System.out;
  4071.         AccurateMathCalc.printarray(out, "EXP_INT_TABLE_A", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_A);
  4072.         AccurateMathCalc.printarray(out, "EXP_INT_TABLE_B", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_B);
  4073.         AccurateMathCalc.printarray(out, "EXP_FRAC_TABLE_A", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_A);
  4074.         AccurateMathCalc.printarray(out, "EXP_FRAC_TABLE_B", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_B);
  4075.         AccurateMathCalc.printarray(out, "LN_MANT", LN_MANT_LEN, lnMant.LN_MANT);
  4076.         AccurateMathCalc.printarray(out, "SINE_TABLE_A", SINE_TABLE_LEN, SINE_TABLE_A);
  4077.         AccurateMathCalc.printarray(out, "SINE_TABLE_B", SINE_TABLE_LEN, SINE_TABLE_B);
  4078.         AccurateMathCalc.printarray(out, "COSINE_TABLE_A", SINE_TABLE_LEN, COSINE_TABLE_A);
  4079.         AccurateMathCalc.printarray(out, "COSINE_TABLE_B", SINE_TABLE_LEN, COSINE_TABLE_B);
  4080.         AccurateMathCalc.printarray(out, "TANGENT_TABLE_A", SINE_TABLE_LEN, TANGENT_TABLE_A);
  4081.         AccurateMathCalc.printarray(out, "TANGENT_TABLE_B", SINE_TABLE_LEN, TANGENT_TABLE_B);
  4082.     }

  4084.     /** Enclose large data table in nested static class so it's only loaded on first access. */
  4085.     private static final class ExpIntTable {
  4086.         /** Exponential evaluated at integer values,
  4087.          * exp(x) =  expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX].
  4088.          */
  4089.         private static final double[] EXP_INT_TABLE_A;
  4090.         /** Exponential evaluated at integer values.
  4091.          * exp(x) =  expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX]
  4092.          */
  4093.         private static final double[] EXP_INT_TABLE_B;

  4095.         static {
  4096.             if (RECOMPUTE_TABLES_AT_RUNTIME) {
  4097.                 EXP_INT_TABLE_A = new double[AccurateMath.EXP_INT_TABLE_LEN];
  4098.                 EXP_INT_TABLE_B = new double[AccurateMath.EXP_INT_TABLE_LEN];

  4100.                 final double[] tmp = new double[2];
  4101.                 final double[] recip = new double[2];

  4103.                 // Populate expIntTable
  4104.                 for (int i = 0; i < AccurateMath.EXP_INT_TABLE_MAX_INDEX; i++) {
  4105.                     AccurateMathCalc.expint(i, tmp);
  4106.                     EXP_INT_TABLE_A[i + AccurateMath.EXP_INT_TABLE_MAX_INDEX] = tmp[0];
  4107.                     EXP_INT_TABLE_B[i + AccurateMath.EXP_INT_TABLE_MAX_INDEX] = tmp[1];

  4109.                     if (i != 0) {
  4110.                         // Negative integer powers
  4111.                         AccurateMathCalc.splitReciprocal(tmp, recip);
  4112.                         EXP_INT_TABLE_A[AccurateMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[0];
  4113.                         EXP_INT_TABLE_B[AccurateMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[1];
  4114.                     }
  4115.                 }
  4116.             } else {
  4117.                 EXP_INT_TABLE_A = AccurateMathLiteralArrays.loadExpIntA();
  4118.                 EXP_INT_TABLE_B = AccurateMathLiteralArrays.loadExpIntB();
  4119.             }
  4120.         }
  4121.     }

  4123.     /** Enclose large data table in nested static class so it's only loaded on first access. */
  4124.     private static final class ExpFracTable {
  4125.         /** Exponential over the range of 0 - 1 in increments of 2^-10
  4126.          * exp(x/1024) =  expFracTableA[x] + expFracTableB[x].
  4127.          * 1024 = 2^10
  4128.          */
  4129.         private static final double[] EXP_FRAC_TABLE_A;
  4130.         /** Exponential over the range of 0 - 1 in increments of 2^-10
  4131.          * exp(x/1024) =  expFracTableA[x] + expFracTableB[x].
  4132.          */
  4133.         private static final double[] EXP_FRAC_TABLE_B;

  4135.         static {
  4136.             if (RECOMPUTE_TABLES_AT_RUNTIME) {
  4137.                 EXP_FRAC_TABLE_A = new double[AccurateMath.EXP_FRAC_TABLE_LEN];
  4138.                 EXP_FRAC_TABLE_B = new double[AccurateMath.EXP_FRAC_TABLE_LEN];

  4140.                 final double[] tmp = new double[2];

  4142.                 // Populate expFracTable
  4143.                 final double factor = 1d / (EXP_FRAC_TABLE_LEN - 1);
  4144.                 for (int i = 0; i < EXP_FRAC_TABLE_A.length; i++) {
  4145.                     AccurateMathCalc.slowexp(i * factor, tmp);
  4146.                     EXP_FRAC_TABLE_A[i] = tmp[0];
  4147.                     EXP_FRAC_TABLE_B[i] = tmp[1];
  4148.                 }
  4149.             } else {
  4150.                 EXP_FRAC_TABLE_A = AccurateMathLiteralArrays.loadExpFracA();
  4151.                 EXP_FRAC_TABLE_B = AccurateMathLiteralArrays.loadExpFracB();
  4152.             }
  4153.         }
  4154.     }

  4156.     /** Enclose large data table in nested static class so it's only loaded on first access. */
  4157.     private static final class lnMant {
  4158.         /** Extended precision logarithm table over the range 1 - 2 in increments of 2^-10. */
  4159.         private static final double[][] LN_MANT;

  4161.         static {
  4162.             if (RECOMPUTE_TABLES_AT_RUNTIME) {
  4163.                 LN_MANT = new double[AccurateMath.LN_MANT_LEN][];

  4165.                 // Populate lnMant table
  4166.                 for (int i = 0; i < LN_MANT.length; i++) {
  4167.                     final double d = Double.longBitsToDouble((((long) i) << 42) | 0x3ff0000000000000L);
  4168.                     LN_MANT[i] = AccurateMathCalc.slowLog(d);
  4169.                 }
  4170.             } else {
  4171.                 LN_MANT = AccurateMathLiteralArrays.loadLnMant();
  4172.             }
  4173.         }
  4174.     }

  4176.     /** Enclose the Cody/Waite reduction (used in "sin", "cos" and "tan"). */
  4177.     private static class CodyWaite {
  4178.         /** k. */
  4179.         private final int finalK;
  4180.         /** remA. */
  4181.         private final double finalRemA;
  4182.         /** remB. */
  4183.         private final double finalRemB;

  4185.         /**
  4186.          * @param xa Argument.
  4187.          */
  4188.         CodyWaite(double xa) {
  4189.             // Estimate k.
  4190.             //k = (int)(xa / 1.5707963267948966);
  4191.             int k = (int)(xa * 0.6366197723675814);

  4193.             // Compute remainder.
  4194.             double remA;
  4195.             double remB;
  4196.             while (true) {
  4197.                 double a = -k * 1.570796251296997;
  4198.                 remA = xa + a;
  4199.                 remB = -(remA - xa - a);

  4201.                 a = -k * 7.549789948768648E-8;
  4202.                 double b = remA;
  4203.                 remA = a + b;
  4204.                 remB += -(remA - b - a);

  4206.                 a = -k * 6.123233995736766E-17;
  4207.                 b = remA;
  4208.                 remA = a + b;
  4209.                 remB += -(remA - b - a);

  4211.                 if (remA > 0) {
  4212.                     break;
  4213.                 }

  4215.                 // Remainder is negative, so decrement k and try again.
  4216.                 // This should only happen if the input is very close
  4217.                 // to an even multiple of pi/2.
  4218.                 --k;
  4219.             }

  4221.             this.finalK = k;
  4222.             this.finalRemA = remA;
  4223.             this.finalRemB = remB;
  4224.         }

  4226.         /**
  4227.          * @return k
  4228.          */
  4229.         int getK() {
  4230.             return finalK;
  4231.         }
  4232.         /**
  4233.          * @return remA
  4234.          */
  4235.         double getRemA() {
  4236.             return finalRemA;
  4237.         }
  4238.         /**
  4239.          * @return remB
  4240.          */
  4241.         double getRemB() {
  4242.             return finalRemB;
  4243.         }
  4244.     }
  4245. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT