FastMath.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.  *      https://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.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.util;

  24. import org.hipparchus.CalculusFieldElement;
  25. import org.hipparchus.exception.LocalizedCoreFormats;
  26. import org.hipparchus.exception.MathRuntimeException;

  28. /**
  29.  * Faster, more accurate, portable alternative to {@link Math} and
  30.  * {@link StrictMath} for large scale computation.
  31.  * <p>
  32.  * FastMath is a drop-in replacement for both Math and StrictMath. This
  33.  * means that for any method in Math (say {@code Math.sin(x)} or
  34.  * {@code Math.cbrt(y)}), user can directly change the class and use the
  35.  * methods as is (using {@code FastMath.sin(x)} or {@code FastMath.cbrt(y)}
  36.  * in the previous example).
  37.  * <p>
  38.  * FastMath speed is achieved by relying heavily on optimizing compilers
  39.  * to native code present in many JVMs today and use of large tables.
  40.  * The larger tables are lazily initialized on first use, so that the setup
  41.  * time does not penalize methods that don't need them.
  42.  * <p>
  43.  * Note that FastMath is
  44.  * extensively used inside Hipparchus, so by calling some algorithms,
  45.  * the overhead when the the tables need to be initialized will occur
  46.  * regardless of the end-user calling FastMath methods directly or not.
  47.  * Performance figures for a specific JVM and hardware can be evaluated by
  48.  * running the FastMathTestPerformance tests in the test directory of the source
  49.  * distribution.
  50.  * <p>
  51.  * FastMath accuracy should be mostly independent of the JVM as it relies only
  52.  * on IEEE-754 basic operations and on embedded tables. Almost all operations
  53.  * are accurate to about 0.5 ulp throughout the domain range. This statement,
  54.  * of course is only a rough global observed behavior, it is <em>not</em> a
  55.  * guarantee for <em>every</em> double numbers input (see William Kahan's <a
  56.  * href="http://en.wikipedia.org/wiki/Rounding#The_table-maker.27s_dilemma">Table
  57.  * Maker's Dilemma</a>).
  58.  * <p>
  59.  * FastMath additionally implements the following methods not found in Math/StrictMath:
  60.  * <ul>
  61.  * <li>{@link #asinh(double)}</li>
  62.  * <li>{@link #acosh(double)}</li>
  63.  * <li>{@link #atanh(double)}</li>
  64.  * </ul>
  65.  * The following methods are found in Math/StrictMath since 1.6 only, they are provided
  66.  * by FastMath even in 1.5 Java virtual machines
  67.  * <ul>
  68.  * <li>{@link #copySign(double, double)}</li>
  69.  * <li>{@link #getExponent(double)}</li>
  70.  * <li>{@link #nextAfter(double,double)}</li>
  71.  * <li>{@link #nextUp(double)}</li>
  72.  * <li>{@link #scalb(double, int)}</li>
  73.  * <li>{@link #copySign(float, float)}</li>
  74.  * <li>{@link #getExponent(float)}</li>
  75.  * <li>{@link #nextAfter(float,double)}</li>
  76.  * <li>{@link #nextUp(float)}</li>
  77.  * <li>{@link #scalb(float, int)}</li>
  78.  * </ul>
  79.  */
  80. public class FastMath {
  81.     /** Archimede's constant PI, ratio of circle circumference to diameter. */
  82.     public static final double PI = 105414357.0 / 33554432.0 + 1.984187159361080883e-9;

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

  87.     /** Index of exp(0) in the array of integer exponentials. */
  88.     static final int EXP_INT_TABLE_MAX_INDEX = 750;
  89.     /** Length of the array of integer exponentials. */
  90.     static final int EXP_INT_TABLE_LEN = EXP_INT_TABLE_MAX_INDEX * 2;
  91.     /** Logarithm table length. */
  92.     static final int LN_MANT_LEN = 1024;
  93.     /** Exponential fractions table length. */
  94.     static final int EXP_FRAC_TABLE_LEN = 1025; // 0, 1/1024, ... 1024/1024

  96.     /** StrictMath.log(Double.MAX_VALUE): {@value} */
  97.     private static final double LOG_MAX_VALUE = StrictMath.log(Double.MAX_VALUE);

  99.     /** Indicator for tables initialization.
  100.      * <p>
  101.      * This compile-time constant should be set to true only if one explicitly
  102.      * wants to compute the tables at class loading time instead of using the
  103.      * already computed ones provided as literal arrays below.
  104.      * </p>
  105.      */
  106.     private static final boolean RECOMPUTE_TABLES_AT_RUNTIME = false;

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

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

  114.     /** Coefficients for log, when input 0.99 < x < 1.01. */
  115.     private static final double[][] LN_QUICK_COEF = {
  116.         {1.0, 5.669184079525E-24},
  117.         {-0.25, -0.25},
  118.         {0.3333333134651184, 1.986821492305628E-8},
  119.         {-0.25, -6.663542893624021E-14},
  120.         {0.19999998807907104, 1.1921056801463227E-8},
  121.         {-0.1666666567325592, -7.800414592973399E-9},
  122.         {0.1428571343421936, 5.650007086920087E-9},
  123.         {-0.12502530217170715, -7.44321345601866E-11},
  124.         {0.11113807559013367, 9.219544613762692E-9},
  125.     };

  127.     /** Coefficients for log in the range of 1.0 < x < 1.0 + 2^-10. */
  128.     private static final double[][] LN_HI_PREC_COEF = {
  129.         {1.0, -6.032174644509064E-23},
  130.         {-0.25, -0.25},
  131.         {0.3333333134651184, 1.9868161777724352E-8},
  132.         {-0.2499999701976776, -2.957007209750105E-8},
  133.         {0.19999954104423523, 1.5830993332061267E-10},
  134.         {-0.16624879837036133, -2.6033824355191673E-8}
  135.     };

  137.     /** Sine, Cosine, Tangent tables are for 0, 1/8, 2/8, ... 13/8 = PI/2 approx. */

  139.     /** Sine table (high bits). */
  140.     private static final double[] SINE_TABLE_A =
  141.         {
  142.         +0.0d,
  143.         +0.1246747374534607d,
  144.         +0.24740394949913025d,
  145.         +0.366272509098053d,
  146.         +0.4794255495071411d,
  147.         +0.5850973129272461d,
  148.         +0.6816387176513672d,
  149.         +0.7675435543060303d,
  150.         +0.8414709568023682d,
  151.         +0.902267575263977d,
  152.         +0.9489846229553223d,
  153.         +0.9808930158615112d,
  154.         +0.9974949359893799d,
  155.         +0.9985313415527344d,
  156.     };

  158.     /** Sine table (low bits). */
  159.     private static final double[] SINE_TABLE_B =
  160.         {
  161.         +0.0d,
  162.         -4.068233003401932E-9d,
  163.         +9.755392680573412E-9d,
  164.         +1.9987994582857286E-8d,
  165.         -1.0902938113007961E-8d,
  166.         -3.9986783938944604E-8d,
  167.         +4.23719669792332E-8d,
  168.         -5.207000323380292E-8d,
  169.         +2.800552834259E-8d,
  170.         +1.883511811213715E-8d,
  171.         -3.5997360512765566E-9d,
  172.         +4.116164446561962E-8d,
  173.         +5.0614674548127384E-8d,
  174.         -1.0129027912496858E-9d,
  175.     };

  177.     /** Cosine table (high bits). */
  178.     private static final double[] COSINE_TABLE_A =
  179.         {
  180.         +1.0d,
  181.         +0.9921976327896118d,
  182.         +0.9689123630523682d,
  183.         +0.9305076599121094d,
  184.         +0.8775825500488281d,
  185.         +0.8109631538391113d,
  186.         +0.7316888570785522d,
  187.         +0.6409968137741089d,
  188.         +0.5403022766113281d,
  189.         +0.4311765432357788d,
  190.         +0.3153223395347595d,
  191.         +0.19454771280288696d,
  192.         +0.07073719799518585d,
  193.         -0.05417713522911072d,
  194.     };

  196.     /** Cosine table (low bits). */
  197.     private static final double[] COSINE_TABLE_B =
  198.         {
  199.         +0.0d,
  200.         +3.4439717236742845E-8d,
  201.         +5.865827662008209E-8d,
  202.         -3.7999795083850525E-8d,
  203.         +1.184154459111628E-8d,
  204.         -3.43338934259355E-8d,
  205.         +1.1795268640216787E-8d,
  206.         +4.438921624363781E-8d,
  207.         +2.925681159240093E-8d,
  208.         -2.6437112632041807E-8d,
  209.         +2.2860509143963117E-8d,
  210.         -4.813899778443457E-9d,
  211.         +3.6725170580355583E-9d,
  212.         +2.0217439756338078E-10d,
  213.     };


  216.     /** Tangent table, used by atan() (high bits). */
  217.     private static final double[] TANGENT_TABLE_A =
  218.         {
  219.         +0.0d,
  220.         +0.1256551444530487d,
  221.         +0.25534194707870483d,
  222.         +0.3936265707015991d,
  223.         +0.5463024377822876d,
  224.         +0.7214844226837158d,
  225.         +0.9315965175628662d,
  226.         +1.1974215507507324d,
  227.         +1.5574076175689697d,
  228.         +2.092571258544922d,
  229.         +3.0095696449279785d,
  230.         +5.041914939880371d,
  231.         +14.101419448852539d,
  232.         -18.430862426757812d,
  233.     };

  235.     /** Tangent table, used by atan() (low bits). */
  236.     private static final double[] TANGENT_TABLE_B =
  237.         {
  238.         +0.0d,
  239.         -7.877917738262007E-9d,
  240.         -2.5857668567479893E-8d,
  241.         +5.2240336371356666E-9d,
  242.         +5.206150291559893E-8d,
  243.         +1.8307188599677033E-8d,
  244.         -5.7618793749770706E-8d,
  245.         +7.848361555046424E-8d,
  246.         +1.0708593250394448E-7d,
  247.         +1.7827257129423813E-8d,
  248.         +2.893485277253286E-8d,
  249.         +3.1660099222737955E-7d,
  250.         +4.983191803254889E-7d,
  251.         -3.356118100840571E-7d,
  252.     };

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

  275.     /** Bits of pi/4, need for reducePayneHanek(). */
  276.     private static final long[] PI_O_4_BITS = {
  277.         (0xc90fdaa2L << 32) | 0x2168c234L,
  278.         (0xc4c6628bL << 32) | 0x80dc1cd1L };

  280.     /** Eighths.
  281.      * This is used by sinQ, because its faster to do a table lookup than
  282.      * a multiply in this time-critical routine
  283.      */
  284.     private static final double[]
  285.                     EIGHTHS = {0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625};

  287.     /** Table of 2^((n+2)/3) */
  288.     private static final double[] CBRTTWO = { 0.6299605249474366,
  289.                                             0.7937005259840998,
  290.                                             1.0,
  291.                                             1.2599210498948732,
  292.                                             1.5874010519681994 };

  294.     /*
  295.      *  There are 52 bits in the mantissa of a double.
  296.      *  For additional precision, the code splits double numbers into two parts,
  297.      *  by clearing the low order 30 bits if possible, and then performs the arithmetic
  298.      *  on each half separately.
  299.      */

  301.     /**
  302.      * 0x40000000 - used to split a double into two parts, both with the low order bits cleared.
  303.      * Equivalent to 2^30.
  304.      */
  305.     private static final long HEX_40000000 = 0x40000000L; // 1073741824L

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

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

  313.     /** Mask used to clear the non-sign part of a long. */
  314.     private static final long MASK_NON_SIGN_LONG = 0x7fffffffffffffffl;

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

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

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

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

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

  363.     /**
  364.      * Private Constructor
  365.      */
  366.     private FastMath() {}

  368.     // Generic helper methods

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

  386.     /** Compute the square root of a number.
  387.      * <p><b>Note:</b> this implementation currently delegates to {@link Math#sqrt}
  388.      * @param a number on which evaluation is done
  389.      * @return square root of a
  390.      */
  391.     public static double sqrt(final double a) {
  392.         return Math.sqrt(a);
  393.     }

  395.     /** Compute the hyperbolic cosine of a number.
  396.      * @param x number on which evaluation is done
  397.      * @return hyperbolic cosine of x
  398.      */
  399.     public static double cosh(double x) {
  400.       if (Double.isNaN(x)) {
  401.           return x;
  402.       }

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

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

  409.       if (x > 20) {
  410.           if (x >= LOG_MAX_VALUE) {
  411.               // Avoid overflow (MATH-905).
  412.               final double t = exp(0.5 * x);
  413.               return (0.5 * t) * t;
  414.           } else {
  415.               return 0.5 * exp(x);
  416.           }
  417.       } else if (x < -20) {
  418.           if (x <= -LOG_MAX_VALUE) {
  419.               // Avoid overflow (MATH-905).
  420.               final double t = exp(-0.5 * x);
  421.               return (0.5 * t) * t;
  422.           } else {
  423.               return 0.5 * exp(-x);
  424.           }
  425.       }

  427.       final double[] hiPrec = new double[2];
  428.       if (x < 0.0) {
  429.           x = -x;
  430.       }
  431.       exp(x, 0.0, hiPrec);

  433.       double ya = hiPrec[0] + hiPrec[1];
  434.       double yb = -(ya - hiPrec[0] - hiPrec[1]);

  436.       double temp = ya * HEX_40000000;
  437.       double yaa = ya + temp - temp;
  438.       double yab = ya - yaa;

  440.       // recip = 1/y
  441.       double recip = 1.0/ya;
  442.       temp = recip * HEX_40000000;
  443.       double recipa = recip + temp - temp;
  444.       double recipb = recip - recipa;

  446.       // Correct for rounding in division
  447.       recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
  448.       // Account for yb
  449.       recipb += -yb * recip * recip;

  451.       // y = y + 1/y
  452.       temp = ya + recipa;
  453.       yb += -(temp - ya - recipa);
  454.       ya = temp;
  455.       temp = ya + recipb;
  456.       yb += -(temp - ya - recipb);
  457.       ya = temp;

  459.       double result = ya + yb;
  460.       result *= 0.5;
  461.       return result;
  462.     }

  464.     /** Compute the hyperbolic sine of a number.
  465.      * @param x number on which evaluation is done
  466.      * @return hyperbolic sine of x
  467.      */
  468.     public static double sinh(double x) {
  469.       boolean negate = false;
  470.       if (Double.isNaN(x)) {
  471.           return x;
  472.       }

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

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

  479.       if (x > 20) {
  480.           if (x >= LOG_MAX_VALUE) {
  481.               // Avoid overflow (MATH-905).
  482.               final double t = exp(0.5 * x);
  483.               return (0.5 * t) * t;
  484.           } else {
  485.               return 0.5 * exp(x);
  486.           }
  487.       } else if (x < -20) {
  488.           if (x <= -LOG_MAX_VALUE) {
  489.               // Avoid overflow (MATH-905).
  490.               final double t = exp(-0.5 * x);
  491.               return (-0.5 * t) * t;
  492.           } else {
  493.               return -0.5 * exp(-x);
  494.           }
  495.       }

  497.       if (x == 0) {
  498.           return x;
  499.       }

  501.       if (x < 0.0) {
  502.           x = -x;
  503.           negate = true;
  504.       }

  506.       double[] hiPrec = new double[2];
  507.       double result;

  509.       if (x > 0.25) {
  510.           exp(x, 0.0, hiPrec);

  512.           double ya = hiPrec[0] + hiPrec[1];
  513.           double yb = -(ya - hiPrec[0] - hiPrec[1]);

  515.           double temp = ya * HEX_40000000;
  516.           double yaa = ya + temp - temp;
  517.           double yab = ya - yaa;

  519.           // recip = 1/y
  520.           double recip = 1.0/ya;
  521.           temp = recip * HEX_40000000;
  522.           double recipa = recip + temp - temp;
  523.           double recipb = recip - recipa;

  525.           // Correct for rounding in division
  526.           recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
  527.           // Account for yb
  528.           recipb += -yb * recip * recip;

  530.           recipa = -recipa;
  531.           recipb = -recipb;

  533.           // y = y - 1/y
  534.           temp = ya + recipa;
  535.           yb += -(temp - ya - recipa);
  536.           ya = temp;
  537.           temp = ya + recipb;
  538.           yb += -(temp - ya - recipb);
  539.           ya = temp;

  541.           result = ya + yb;
  542.           result *= 0.5;
  543.       } else {
  544.           expm1(x, hiPrec);

  546.           double ya = hiPrec[0] + hiPrec[1];
  547.           double yb = -(ya - hiPrec[0] - hiPrec[1]);

  549.           /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
  550.           double denom = 1.0 + ya;
  551.           double denomr = 1.0 / denom;
  552.           double denomb = -(denom - 1.0 - ya) + yb;
  553.           double ratio = ya * denomr;
  554.           double temp = ratio * HEX_40000000;
  555.           double ra = ratio + temp - temp;
  556.           double rb = ratio - ra;

  558.           temp = denom * HEX_40000000;
  559.           double za = denom + temp - temp;
  560.           double zb = denom - za;

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

  564.           // Adjust for yb
  565.           rb += yb*denomr;                        // numerator
  566.           rb += -ya * denomb * denomr * denomr;   // denominator

  568.           // y = y - 1/y
  569.           temp = ya + ra;
  570.           yb += -(temp - ya - ra);
  571.           ya = temp;
  572.           temp = ya + rb;
  573.           yb += -(temp - ya - rb);
  574.           ya = temp;

  576.           result = ya + yb;
  577.           result *= 0.5;
  578.       }

  580.       if (negate) {
  581.           result = -result;
  582.       }

  584.       return result;
  585.     }

  587.     /**
  588.      * Combined hyperbolic sine and hyperbolic cosine function.
  589.      *
  590.      * @param x Argument.
  591.      * @return [sinh(x), cosh(x)]
  592.      */
  593.     public static SinhCosh sinhCosh(double x) {
  594.       boolean negate = false;
  595.       if (Double.isNaN(x)) {
  596.           return new SinhCosh(x, x);
  597.       }

  599.       // sinh[z] = (exp(z) - exp(-z) / 2
  600.       // cosh[z] = (exp(z) + exp(-z))/2

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

  605.       if (x > 20) {
  606.           final double e;
  607.           if (x >= LOG_MAX_VALUE) {
  608.               // Avoid overflow (MATH-905).
  609.               final double t = exp(0.5 * x);
  610.               e = (0.5 * t) * t;
  611.           } else {
  612.               e = 0.5 * exp(x);
  613.           }
  614.           return new SinhCosh(e, e);
  615.       } else if (x < -20) {
  616.           final double e;
  617.           if (x <= -LOG_MAX_VALUE) {
  618.               // Avoid overflow (MATH-905).
  619.               final double t = exp(-0.5 * x);
  620.               e = (-0.5 * t) * t;
  621.           } else {
  622.               e = -0.5 * exp(-x);
  623.           }
  624.           return new SinhCosh(e, -e);
  625.       }

  627.       if (x == 0) {
  628.           return new SinhCosh(x, 1.0);
  629.       }

  631.       if (x < 0.0) {
  632.           x = -x;
  633.           negate = true;
  634.       }

  636.       double[] hiPrec = new double[2];
  637.       double resultM;
  638.       double resultP;

  640.       if (x > 0.25) {
  641.           exp(x, 0.0, hiPrec);

  643.           final double ya = hiPrec[0] + hiPrec[1];
  644.           final double yb = -(ya - hiPrec[0] - hiPrec[1]);

  646.           double temp = ya * HEX_40000000;
  647.           double yaa = ya + temp - temp;
  648.           double yab = ya - yaa;

  650.           // recip = 1/y
  651.           double recip = 1.0/ya;
  652.           temp = recip * HEX_40000000;
  653.           double recipa = recip + temp - temp;
  654.           double recipb = recip - recipa;

  656.           // Correct for rounding in division
  657.           recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
  658.           // Account for yb
  659.           recipb += -yb * recip * recip;

  661.           // y = y - 1/y
  662.           temp = ya - recipa;
  663.           double ybM = yb - (temp - ya + recipa);
  664.           double yaM = temp;
  665.           temp = yaM - recipb;
  666.           ybM += -(temp - yaM + recipb);
  667.           yaM = temp;
  668.           resultM = yaM + ybM;
  669.           resultM *= 0.5;

  671.           // y = y + 1/y
  672.           temp = ya + recipa;
  673.           double ybP = yb - (temp - ya - recipa);
  674.           double yaP = temp;
  675.           temp = yaP + recipb;
  676.           ybP += -(temp - yaP - recipb);
  677.           yaP = temp;
  678.           resultP = yaP + ybP;
  679.           resultP *= 0.5;

  681.       } else {
  682.           expm1(x, hiPrec);

  684.           final double ya = hiPrec[0] + hiPrec[1];
  685.           final double yb = -(ya - hiPrec[0] - hiPrec[1]);

  687.           /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
  688.           double denom = 1.0 + ya;
  689.           double denomr = 1.0 / denom;
  690.           double denomb = -(denom - 1.0 - ya) + yb;
  691.           double ratio = ya * denomr;
  692.           double temp = ratio * HEX_40000000;
  693.           double ra = ratio + temp - temp;
  694.           double rb = ratio - ra;

  696.           temp = denom * HEX_40000000;
  697.           double za = denom + temp - temp;
  698.           double zb = denom - za;

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

  702.           // Adjust for yb
  703.           rb += yb*denomr;                        // numerator
  704.           rb += -ya * denomb * denomr * denomr;   // denominator

  706.           // y = y - 1/y
  707.           temp = ya + ra;
  708.           double ybM = yb - (temp - ya - ra);
  709.           double yaM = temp;
  710.           temp = yaM + rb;
  711.           ybM += -(temp - yaM - rb);
  712.           yaM = temp;
  713.           resultM = yaM + ybM;
  714.           resultM *= 0.5;

  716.           // y = y + 1/y + 2
  717.           temp = ya - ra;
  718.           double ybP = yb - (temp - ya + ra);
  719.           double yaP = temp;
  720.           temp = yaP - rb;
  721.           ybP += -(temp - yaP + rb);
  722.           yaP = temp;
  723.           resultP = yaP + ybP + 2;
  724.           resultP *= 0.5;
  725.       }

  727.       if (negate) {
  728.           resultM = -resultM;
  729.       }

  731.       return new SinhCosh(resultM, resultP);

  733.     }

  735.     /**
  736.      * Combined hyperbolic sine and hyperbolic cosine function.
  737.      *
  738.      * @param x Argument.
  739.      * @param <T> the type of the field element
  740.      * @return [sinh(x), cosh(x)]
  741.      */
  742.     public static <T extends CalculusFieldElement<T>> FieldSinhCosh<T> sinhCosh(T x) {
  743.         return x.sinhCosh();
  744.     }

  746.     /** Compute the hyperbolic tangent of a number.
  747.      * @param x number on which evaluation is done
  748.      * @return hyperbolic tangent of x
  749.      */
  750.     public static double tanh(double x) {
  751.       boolean negate = false;

  753.       if (Double.isNaN(x)) {
  754.           return x;
  755.       }

  757.       // tanh[z] = sinh[z] / cosh[z]
  758.       // = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
  759.       // = (exp(2x) - 1) / (exp(2x) + 1)

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

  763.       if (x > 20.0) {
  764.           return 1.0;
  765.       }

  767.       if (x < -20) {
  768.           return -1.0;
  769.       }

  771.       if (x == 0) {
  772.           return x;
  773.       }

  775.       if (x < 0.0) {
  776.           x = -x;
  777.           negate = true;
  778.       }

  780.       double result;
  781.       if (x >= 0.5) {
  782.           double[] hiPrec = new double[2];
  783.           // tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
  784.           exp(x*2.0, 0.0, hiPrec);

  786.           double ya = hiPrec[0] + hiPrec[1];
  787.           double yb = -(ya - hiPrec[0] - hiPrec[1]);

  789.           /* Numerator */
  790.           double na = -1.0 + ya;
  791.           double nb = -(na + 1.0 - ya);
  792.           double temp = na + yb;
  793.           nb += -(temp - na - yb);
  794.           na = temp;

  796.           /* Denominator */
  797.           double da = 1.0 + ya;
  798.           double db = -(da - 1.0 - ya);
  799.           temp = da + yb;
  800.           db += -(temp - da - yb);
  801.           da = temp;

  803.           temp = da * HEX_40000000;
  804.           double daa = da + temp - temp;
  805.           double dab = da - daa;

  807.           // ratio = na/da
  808.           double ratio = na/da;
  809.           temp = ratio * HEX_40000000;
  810.           double ratioa = ratio + temp - temp;
  811.           double ratiob = ratio - ratioa;

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

  816.           // Account for nb
  817.           ratiob += nb / da;
  818.           // Account for db
  819.           ratiob += -db * na / da / da;

  821.           result = ratioa + ratiob;
  822.       }
  823.       else {
  824.           double[] hiPrec = new double[2];
  825.           // tanh(x) = expm1(2x) / (expm1(2x) + 2)
  826.           expm1(x*2.0, hiPrec);

  828.           double ya = hiPrec[0] + hiPrec[1];
  829.           double yb = -(ya - hiPrec[0] - hiPrec[1]);

  831.           /* Numerator */
  832.           double na = ya;
  833.           double nb = yb;

  835.           /* Denominator */
  836.           double da = 2.0 + ya;
  837.           double db = -(da - 2.0 - ya);
  838.           double temp = da + yb;
  839.           db += -(temp - da - yb);
  840.           da = temp;

  842.           temp = da * HEX_40000000;
  843.           double daa = da + temp - temp;
  844.           double dab = da - daa;

  846.           // ratio = na/da
  847.           double ratio = na/da;
  848.           temp = ratio * HEX_40000000;
  849.           double ratioa = ratio + temp - temp;
  850.           double ratiob = ratio - ratioa;

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

  855.           // Account for nb
  856.           ratiob += nb / da;
  857.           // Account for db
  858.           ratiob += -db * na / da / da;

  860.           result = ratioa + ratiob;
  861.       }

  863.       if (negate) {
  864.           result = -result;
  865.       }

  867.       return result;
  868.     }

  870.     /** Compute the inverse hyperbolic cosine of a number.
  871.      * @param a number on which evaluation is done
  872.      * @return inverse hyperbolic cosine of a
  873.      */
  874.     public static double acosh(final double a) {
  875.         return log(a + sqrt(a * a - 1));
  876.     }

  878.     /** Compute the inverse hyperbolic sine of a number.
  879.      * @param a number on which evaluation is done
  880.      * @return inverse hyperbolic sine of a
  881.      */
  882.     public static double asinh(double a) {
  883.         boolean negative = false;
  884.         if (a < 0) {
  885.             negative = true;
  886.             a = -a;
  887.         }

  889.         double absAsinh;
  890.         if (a > 0.167) {
  891.             absAsinh = log(sqrt(a * a + 1) + a);
  892.         } else {
  893.             final double a2 = a * a;
  894.             if (a > 0.097) {
  895.                 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);
  896.             } else if (a > 0.036) {
  897.                 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);
  898.             } else if (a > 0.0036) {
  899.                 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);
  900.             } else {
  901.                 absAsinh = a * (1 - a2 * (F_1_3 - a2 * F_1_5 * F_3_4) * F_1_2);
  902.             }
  903.         }

  905.         return negative ? -absAsinh : absAsinh;
  906.     }

  908.     /** Compute the inverse hyperbolic tangent of a number.
  909.      * @param a number on which evaluation is done
  910.      * @return inverse hyperbolic tangent of a
  911.      */
  912.     public static double atanh(double a) {
  913.         boolean negative = false;
  914.         if (a < 0) {
  915.             negative = true;
  916.             a = -a;
  917.         }

  919.         double absAtanh;
  920.         if (a > 0.15) {
  921.             absAtanh = 0.5 * log((1 + a) / (1 - a));
  922.         } else {
  923.             final double a2 = a * a;
  924.             if (a > 0.087) {
  925.                 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))))))));
  926.             } else if (a > 0.031) {
  927.                 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))))));
  928.             } else if (a > 0.003) {
  929.                 absAtanh = a * (1 + a2 * (F_1_3 + a2 * (F_1_5 + a2 * (F_1_7 + a2 * F_1_9))));
  930.             } else {
  931.                 absAtanh = a * (1 + a2 * (F_1_3 + a2 * F_1_5));
  932.             }
  933.         }

  935.         return negative ? -absAtanh : absAtanh;
  936.     }

  938.     /** Compute the signum of a number.
  939.      * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise
  940.      * @param a number on which evaluation is done
  941.      * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
  942.      */
  943.     public static double signum(final double a) {
  944.         return (a < 0.0) ? -1.0 : ((a > 0.0) ? 1.0 : a); // return +0.0/-0.0/NaN depending on a
  945.     }

  947.     /** Compute the signum of a number.
  948.      * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise
  949.      * @param a number on which evaluation is done
  950.      * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
  951.      */
  952.     public static float signum(final float a) {
  953.         return (a < 0.0f) ? -1.0f : ((a > 0.0f) ? 1.0f : a); // return +0.0/-0.0/NaN depending on a
  954.     }

  956.     /** Compute next number towards positive infinity.
  957.      * @param a number to which neighbor should be computed
  958.      * @return neighbor of a towards positive infinity
  959.      */
  960.     public static double nextUp(final double a) {
  961.         return nextAfter(a, Double.POSITIVE_INFINITY);
  962.     }

  964.     /** Compute next number towards positive infinity.
  965.      * @param a number to which neighbor should be computed
  966.      * @return neighbor of a towards positive infinity
  967.      */
  968.     public static float nextUp(final float a) {
  969.         return nextAfter(a, Float.POSITIVE_INFINITY);
  970.     }

  972.     /** Compute next number towards negative infinity.
  973.      * @param a number to which neighbor should be computed
  974.      * @return neighbor of a towards negative infinity
  975.      */
  976.     public static double nextDown(final double a) {
  977.         return nextAfter(a, Double.NEGATIVE_INFINITY);
  978.     }

  980.     /** Compute next number towards negative infinity.
  981.      * @param a number to which neighbor should be computed
  982.      * @return neighbor of a towards negative infinity
  983.      */
  984.     public static float nextDown(final float a) {
  985.         return nextAfter(a, Float.NEGATIVE_INFINITY);
  986.     }

  988.     /** Clamp a value within an interval.
  989.      * @param value value to clamp
  990.      * @param inf lower bound of the clamping interval
  991.      * @param sup upper bound of the clamping interval
  992.      * @return value clamped within [inf; sup], or value if already within bounds.
  993.      * @since 3.0
  994.      */
  995.     public static int clamp(final int value, final int inf, final int sup) {
  996.         return max(inf, min(value, sup));
  997.     }

  999.     /** Clamp a value within an interval.
  1000.      * @param value value to clamp
  1001.      * @param inf lower bound of the clamping interval
  1002.      * @param sup upper bound of the clamping interval
  1003.      * @return value clamped within [inf; sup], or value if already within bounds.
  1004.      * @since 3.0
  1005.      */
  1006.     public static long clamp(final long value, final long inf, final long sup) {
  1007.         return max(inf, min(value, sup));
  1008.     }

  1010.     /** Clamp a value within an interval.
  1011.      * @param value value to clamp
  1012.      * @param inf lower bound of the clamping interval
  1013.      * @param sup upper bound of the clamping interval
  1014.      * @return value clamped within [inf; sup], or value if already within bounds.
  1015.      * @since 3.0
  1016.      */
  1017.     public static int clamp(final long value, final int inf, final int sup) {
  1018.         return (int) max(inf, min(value, sup));
  1019.     }

  1021.     /** Clamp a value within an interval.
  1022.      * <p>
  1023.      * This method assumes -0.0 is below +0.0
  1024.      * </p>
  1025.      * @param value value to clamp
  1026.      * @param inf lower bound of the clamping interval
  1027.      * @param sup upper bound of the clamping interval
  1028.      * @return value clamped within [inf; sup], or value if already within bounds.
  1029.      * @since 3.0
  1030.      */
  1031.     public static float clamp(final float value, final float inf, final float sup) {
  1032.         return max(inf, min(value, sup));
  1033.     }

  1035.     /** Clamp a value within an interval.
  1036.      * <p>
  1037.      * This method assumes -0.0 is below +0.0
  1038.      * </p>
  1039.      * @param value value to clamp
  1040.      * @param inf lower bound of the clamping interval
  1041.      * @param sup upper bound of the clamping interval
  1042.      * @return value clamped within [inf; sup], or value if already within bounds.
  1043.      * @since 3.0
  1044.      */
  1045.     public static double clamp(final double value, final double inf, final double sup) {
  1046.         return max(inf, min(value, sup));
  1047.     }

  1049.     /** Returns a pseudo-random number between 0.0 and 1.0.
  1050.      * <p><b>Note:</b> this implementation currently delegates to {@link Math#random}
  1051.      * @return a random number between 0.0 and 1.0
  1052.      */
  1053.     public static double random() {
  1054.         return Math.random();
  1055.     }

  1057.     /**
  1058.      * Exponential function.
  1059.      *
  1060.      * Computes exp(x), function result is nearly rounded.   It will be correctly
  1061.      * rounded to the theoretical value for 99.9% of input values, otherwise it will
  1062.      * have a 1 ULP error.
  1063.      *
  1064.      * Method:
  1065.      *    Lookup intVal = exp(int(x))
  1066.      *    Lookup fracVal = exp(int(x-int(x) / 1024.0) * 1024.0 );
  1067.      *    Compute z as the exponential of the remaining bits by a polynomial minus one
  1068.      *    exp(x) = intVal * fracVal * (1 + z)
  1069.      *
  1070.      * Accuracy:
  1071.      *    Calculation is done with 63 bits of precision, so result should be correctly
  1072.      *    rounded for 99.9% of input values, with less than 1 ULP error otherwise.
  1073.      *
  1074.      * @param x   a double
  1075.      * @return double e<sup>x</sup>
  1076.      */
  1077.     public static double exp(double x) {
  1078.         return exp(x, 0.0, null);
  1079.     }

  1081.     /**
  1082.      * Internal helper method for exponential function.
  1083.      * @param x original argument of the exponential function
  1084.      * @param extra extra bits of precision on input (To Be Confirmed)
  1085.      * @param hiPrec extra bits of precision on output (To Be Confirmed)
  1086.      * @return exp(x)
  1087.      */
  1088.     private static double exp(double x, double extra, double[] hiPrec) {
  1089.         double intPartA;
  1090.         double intPartB;
  1091.         int intVal = (int) x;

  1093.         /* Lookup exp(floor(x)).
  1094.          * intPartA will have the upper 22 bits, intPartB will have the lower
  1095.          * 52 bits.
  1096.          */
  1097.         if (x < 0.0) {

  1099.             // We don't check against intVal here as conversion of large negative double values
  1100.             // may be affected by a JIT bug. Subsequent comparisons can safely use intVal
  1101.             if (x < -746d) {
  1102.                 if (hiPrec != null) {
  1103.                     hiPrec[0] = 0.0;
  1104.                     hiPrec[1] = 0.0;
  1105.                 }
  1106.                 return 0.0;
  1107.             }

  1109.             if (intVal < -709) {
  1110.                 /* This will produce a subnormal output */
  1111.                 final double result = exp(x+40.19140625, extra, hiPrec) / 285040095144011776.0;
  1112.                 if (hiPrec != null) {
  1113.                     hiPrec[0] /= 285040095144011776.0;
  1114.                     hiPrec[1] /= 285040095144011776.0;
  1115.                 }
  1116.                 return result;
  1117.             }

  1119.             if (intVal == -709) {
  1120.                 /* exp(1.494140625) is nearly a machine number... */
  1121.                 final double result = exp(x+1.494140625, extra, hiPrec) / 4.455505956692756620;
  1122.                 if (hiPrec != null) {
  1123.                     hiPrec[0] /= 4.455505956692756620;
  1124.                     hiPrec[1] /= 4.455505956692756620;
  1125.                 }
  1126.                 return result;
  1127.             }

  1129.             intVal--;

  1131.         } else {
  1132.             if (intVal > 709) {
  1133.                 if (hiPrec != null) {
  1134.                     hiPrec[0] = Double.POSITIVE_INFINITY;
  1135.                     hiPrec[1] = 0.0;
  1136.                 }
  1137.                 return Double.POSITIVE_INFINITY;
  1138.             }

  1140.         }

  1142.         intPartA = ExpIntTable.EXP_INT_TABLE_A[EXP_INT_TABLE_MAX_INDEX+intVal];
  1143.         intPartB = ExpIntTable.EXP_INT_TABLE_B[EXP_INT_TABLE_MAX_INDEX+intVal];

  1145.         /* Get the fractional part of x, find the greatest multiple of 2^-10 less than
  1146.          * x and look up the exp function of it.
  1147.          * fracPartA will have the upper 22 bits, fracPartB the lower 52 bits.
  1148.          */
  1149.         final int intFrac = (int) ((x - intVal) * 1024.0);
  1150.         final double fracPartA = ExpFracTable.EXP_FRAC_TABLE_A[intFrac];
  1151.         final double fracPartB = ExpFracTable.EXP_FRAC_TABLE_B[intFrac];

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

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

  1164.         /* Remez generated polynomial.  Converges on the interval [0, 2^-10], error
  1165.        is less than 0.5 ULP */
  1166.         double z = 0.04168701738764507;
  1167.         z = z * epsilon + 0.1666666505023083;
  1168.         z = z * epsilon + 0.5000000000042687;
  1169.         z = z * epsilon + 1.0;
  1170.         z = z * epsilon + -3.940510424527919E-20;

  1172.         /* Compute (intPartA+intPartB) * (fracPartA+fracPartB) by binomial
  1173.        expansion.
  1174.        tempA is exact since intPartA and intPartB only have 22 bits each.
  1175.        tempB will have 52 bits of precision.
  1176.          */
  1177.         double tempA = intPartA * fracPartA;
  1178.         double tempB = intPartA * fracPartB + intPartB * fracPartA + intPartB * fracPartB;

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

  1186.         // If tempC is positive infinite, the evaluation below could result in NaN,
  1187.         // because z could be negative at the same time.
  1188.         if (tempC == Double.POSITIVE_INFINITY) {
  1189.             if (hiPrec != null) {
  1190.                 hiPrec[0] = Double.POSITIVE_INFINITY;
  1191.                 hiPrec[1] = 0.0;
  1192.             }
  1193.             return Double.POSITIVE_INFINITY;
  1194.         }

  1196.         final double result;
  1197.         if (extra != 0.0) {
  1198.             result = tempC*extra*z + tempC*extra + tempC*z + tempB + tempA;
  1199.         } else {
  1200.             result = tempC*z + tempB + tempA;
  1201.         }

  1203.         if (hiPrec != null) {
  1204.             // If requesting high precision
  1205.             hiPrec[0] = tempA;
  1206.             hiPrec[1] = tempC*extra*z + tempC*extra + tempC*z + tempB;
  1207.         }

  1209.         return result;
  1210.     }

  1212.     /** Compute exp(x) - 1
  1213.      * @param x number to compute shifted exponential
  1214.      * @return exp(x) - 1
  1215.      */
  1216.     public static double expm1(double x) {
  1217.       return expm1(x, null);
  1218.     }

  1220.     /** Internal helper method for expm1
  1221.      * @param x number to compute shifted exponential
  1222.      * @param hiPrecOut receive high precision result for -1.0 < x < 1.0
  1223.      * @return exp(x) - 1
  1224.      */
  1225.     private static double expm1(double x, double[] hiPrecOut) {
  1226.         if (Double.isNaN(x) || x == 0.0) { // NaN or zero
  1227.             return x;
  1228.         }

  1230.         if (x <= -1.0 || x >= 1.0) {
  1231.             // If not between +/- 1.0
  1232.             //return exp(x) - 1.0;
  1233.             double[] hiPrec = new double[2];
  1234.             exp(x, 0.0, hiPrec);
  1235.             if (x > 0.0) {
  1236.                 return -1.0 + hiPrec[0] + hiPrec[1];
  1237.             } else {
  1238.                 final double ra = -1.0 + hiPrec[0];
  1239.                 double rb = -(ra + 1.0 - hiPrec[0]);
  1240.                 rb += hiPrec[1];
  1241.                 return ra + rb;
  1242.             }
  1243.         }

  1245.         double baseA;
  1246.         double baseB;
  1247.         double epsilon;
  1248.         boolean negative = false;

  1250.         if (x < 0.0) {
  1251.             x = -x;
  1252.             negative = true;
  1253.         }

  1255.         {
  1256.             int intFrac = (int) (x * 1024.0);
  1257.             double tempA = ExpFracTable.EXP_FRAC_TABLE_A[intFrac] - 1.0;
  1258.             double tempB = ExpFracTable.EXP_FRAC_TABLE_B[intFrac];

  1260.             double temp = tempA + tempB;
  1261.             tempB = -(temp - tempA - tempB);
  1262.             tempA = temp;

  1264.             temp = tempA * HEX_40000000;
  1265.             baseA = tempA + temp - temp;
  1266.             baseB = tempB + (tempA - baseA);

  1268.             epsilon = x - intFrac/1024.0;
  1269.         }


  1272.         /* Compute expm1(epsilon) */
  1273.         double zb = 0.008336750013465571;
  1274.         zb = zb * epsilon + 0.041666663879186654;
  1275.         zb = zb * epsilon + 0.16666666666745392;
  1276.         zb = zb * epsilon + 0.49999999999999994;
  1277.         zb *= epsilon;
  1278.         zb *= epsilon;

  1280.         double za = epsilon;
  1281.         double temp = za + zb;
  1282.         zb = -(temp - za - zb);
  1283.         za = temp;

  1285.         temp = za * HEX_40000000;
  1286.         temp = za + temp - temp;
  1287.         zb += za - temp;
  1288.         za = temp;

  1290.         /* Combine the parts.   expm1(a+b) = expm1(a) + expm1(b) + expm1(a)*expm1(b) */
  1291.         double ya = za * baseA;
  1292.         //double yb = za*baseB + zb*baseA + zb*baseB;
  1293.         temp = ya + za * baseB;
  1294.         double yb = -(temp - ya - za * baseB);
  1295.         ya = temp;

  1297.         temp = ya + zb * baseA;
  1298.         yb += -(temp - ya - zb * baseA);
  1299.         ya = temp;

  1301.         temp = ya + zb * baseB;
  1302.         yb += -(temp - ya - zb*baseB);
  1303.         ya = temp;

  1305.         //ya = ya + za + baseA;
  1306.         //yb = yb + zb + baseB;
  1307.         temp = ya + baseA;
  1308.         yb += -(temp - baseA - ya);
  1309.         ya = temp;

  1311.         temp = ya + za;
  1312.         //yb += (ya > za) ? -(temp - ya - za) : -(temp - za - ya);
  1313.         yb += -(temp - ya - za);
  1314.         ya = temp;

  1316.         temp = ya + baseB;
  1317.         //yb += (ya > baseB) ? -(temp - ya - baseB) : -(temp - baseB - ya);
  1318.         yb += -(temp - ya - baseB);
  1319.         ya = temp;

  1321.         temp = ya + zb;
  1322.         //yb += (ya > zb) ? -(temp - ya - zb) : -(temp - zb - ya);
  1323.         yb += -(temp - ya - zb);
  1324.         ya = temp;

  1326.         if (negative) {
  1327.             /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
  1328.             double denom = 1.0 + ya;
  1329.             double denomr = 1.0 / denom;
  1330.             double denomb = -(denom - 1.0 - ya) + yb;
  1331.             double ratio = ya * denomr;
  1332.             temp = ratio * HEX_40000000;
  1333.             final double ra = ratio + temp - temp;
  1334.             double rb = ratio - ra;

  1336.             temp = denom * HEX_40000000;
  1337.             za = denom + temp - temp;
  1338.             zb = denom - za;

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

  1342.             // f(x) = x/1+x
  1343.             // Compute f'(x)
  1344.             // Product rule:  d(uv) = du*v + u*dv
  1345.             // Chain rule:  d(f(g(x)) = f'(g(x))*f(g'(x))
  1346.             // d(1/x) = -1/(x*x)
  1347.             // d(1/1+x) = -1/( (1+x)^2) *  1 =  -1/((1+x)*(1+x))
  1348.             // d(x/1+x) = -x/((1+x)(1+x)) + 1/1+x = 1 / ((1+x)(1+x))

  1350.             // Adjust for yb
  1351.             rb += yb * denomr;                      // numerator
  1352.             rb += -ya * denomb * denomr * denomr;   // denominator

  1354.             // negate
  1355.             ya = -ra;
  1356.             yb = -rb;
  1357.         }

  1359.         if (hiPrecOut != null) {
  1360.             hiPrecOut[0] = ya;
  1361.             hiPrecOut[1] = yb;
  1362.         }

  1364.         return ya + yb;
  1365.     }

  1367.     /**
  1368.      * Natural logarithm.
  1369.      *
  1370.      * @param x   a double
  1371.      * @return log(x)
  1372.      */
  1373.     public static double log(final double x) {
  1374.         return log(x, null);
  1375.     }

  1377.     /**
  1378.      * Internal helper method for natural logarithm function.
  1379.      * @param x original argument of the natural logarithm function
  1380.      * @param hiPrec extra bits of precision on output (To Be Confirmed)
  1381.      * @return log(x)
  1382.      */
  1383.     private static double log(final double x, final double[] hiPrec) {
  1384.         if (x==0) { // Handle special case of +0/-0
  1385.             return Double.NEGATIVE_INFINITY;
  1386.         }
  1387.         long bits = Double.doubleToRawLongBits(x);

  1389.         /* Handle special cases of negative input, and NaN */
  1390.         if (((bits & 0x8000000000000000L) != 0 || Double.isNaN(x)) && x != 0.0) {
  1391.             if (hiPrec != null) {
  1392.                 hiPrec[0] = Double.NaN;
  1393.             }

  1395.             return Double.NaN;
  1396.         }

  1398.         /* Handle special cases of Positive infinity. */
  1399.         if (x == Double.POSITIVE_INFINITY) {
  1400.             if (hiPrec != null) {
  1401.                 hiPrec[0] = Double.POSITIVE_INFINITY;
  1402.             }

  1404.             return Double.POSITIVE_INFINITY;
  1405.         }

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

  1410.         if ((bits & 0x7ff0000000000000L) == 0) {
  1411.             // Subnormal!
  1412.             if (x == 0) {
  1413.                 // Zero
  1414.                 if (hiPrec != null) {
  1415.                     hiPrec[0] = Double.NEGATIVE_INFINITY;
  1416.                 }

  1418.                 return Double.NEGATIVE_INFINITY;
  1419.             }

  1421.             /* Normalize the subnormal number. */
  1422.             bits <<= 1;
  1423.             while ( (bits & 0x0010000000000000L) == 0) {
  1424.                 --exp;
  1425.                 bits <<= 1;
  1426.             }
  1427.         }


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

  1434.             /* Compute x - 1.0 and split it */
  1435.             double xa = x - 1.0;
  1436.             double tmp = xa * HEX_40000000;
  1437.             double aa = xa + tmp - tmp;
  1438.             double ab = xa - aa;
  1439.             xa = aa;
  1440.             double xb = ab;

  1442.             final double[] lnCoef_last = LN_QUICK_COEF[LN_QUICK_COEF.length - 1];
  1443.             double ya = lnCoef_last[0];
  1444.             double yb = lnCoef_last[1];

  1446.             for (int i = LN_QUICK_COEF.length - 2; i >= 0; i--) {
  1447.                 /* Multiply a = y * x */
  1448.                 aa = ya * xa;
  1449.                 ab = ya * xb + yb * xa + yb * xb;
  1450.                 /* split, so now y = a */
  1451.                 tmp = aa * HEX_40000000;
  1452.                 ya = aa + tmp - tmp;
  1453.                 yb = aa - ya + ab;

  1455.                 /* Add  a = y + lnQuickCoef */
  1456.                 final double[] lnCoef_i = LN_QUICK_COEF[i];
  1457.                 aa = ya + lnCoef_i[0];
  1458.                 ab = yb + lnCoef_i[1];
  1459.                 /* Split y = a */
  1460.                 tmp = aa * HEX_40000000;
  1461.                 ya = aa + tmp - tmp;
  1462.                 yb = aa - ya + ab;
  1463.             }

  1465.             /* Multiply a = y * x */
  1466.             aa = ya * xa;
  1467.             ab = ya * xb + yb * xa + yb * xb;
  1468.             /* split, so now y = a */
  1469.             tmp = aa * HEX_40000000;
  1470.             ya = aa + tmp - tmp;
  1471.             yb = aa - ya + ab;

  1473.             return ya + yb;
  1474.         }

  1476.         // lnm is a log of a number in the range of 1.0 - 2.0, so 0 <= lnm < ln(2)
  1477.         final double[] lnm = lnMant.LN_MANT[(int)((bits & 0x000ffc0000000000L) >> 42)];

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

  1482.     epsilon -= 1.0;
  1483.          */

  1485.         // y is the most significant 10 bits of the mantissa
  1486.         //double y = Double.longBitsToDouble(bits & 0xfffffc0000000000L);
  1487.         //double epsilon = (x - y) / y;
  1488.         final double epsilon = (bits & 0x3ffffffffffL) / (TWO_POWER_52 + (bits & 0x000ffc0000000000L));

  1490.         double lnza;
  1491.         double lnzb = 0.0;

  1493.         if (hiPrec != null) {
  1494.             /* split epsilon -> x */
  1495.             double tmp = epsilon * HEX_40000000;
  1496.             double aa = epsilon + tmp - tmp;
  1497.             double ab = epsilon - aa;
  1498.             double xa = aa;
  1499.             double xb = ab;

  1501.             /* Need a more accurate epsilon, so adjust the division. */
  1502.             final double numer = bits & 0x3ffffffffffL;
  1503.             final double denom = TWO_POWER_52 + (bits & 0x000ffc0000000000L);
  1504.             aa = numer - xa*denom - xb * denom;
  1505.             xb += aa / denom;

  1507.             /* Remez polynomial evaluation */
  1508.             final double[] lnCoef_last = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1];
  1509.             double ya = lnCoef_last[0];
  1510.             double yb = lnCoef_last[1];

  1512.             for (int i = LN_HI_PREC_COEF.length - 2; i >= 0; i--) {
  1513.                 /* Multiply a = y * x */
  1514.                 aa = ya * xa;
  1515.                 ab = ya * xb + yb * xa + yb * xb;
  1516.                 /* split, so now y = a */
  1517.                 tmp = aa * HEX_40000000;
  1518.                 ya = aa + tmp - tmp;
  1519.                 yb = aa - ya + ab;

  1521.                 /* Add  a = y + lnHiPrecCoef */
  1522.                 final double[] lnCoef_i = LN_HI_PREC_COEF[i];
  1523.                 aa = ya + lnCoef_i[0];
  1524.                 ab = yb + lnCoef_i[1];
  1525.                 /* Split y = a */
  1526.                 tmp = aa * HEX_40000000;
  1527.                 ya = aa + tmp - tmp;
  1528.                 yb = aa - ya + ab;
  1529.             }

  1531.             /* Multiply a = y * x */
  1532.             aa = ya * xa;
  1533.             ab = ya * xb + yb * xa + yb * xb;

  1535.             /* split, so now lnz = a */
  1536.             /*
  1537.       tmp = aa * 1073741824.0;
  1538.       lnza = aa + tmp - tmp;
  1539.       lnzb = aa - lnza + ab;
  1540.              */
  1541.             lnza = aa + ab;
  1542.             lnzb = -(lnza - aa - ab);
  1543.         } else {
  1544.             /* High precision not required.  Eval Remez polynomial
  1545.          using standard double precision */
  1546.             lnza = -0.16624882440418567;
  1547.             lnza = lnza * epsilon + 0.19999954120254515;
  1548.             lnza = lnza * epsilon + -0.2499999997677497;
  1549.             lnza = lnza * epsilon + 0.3333333333332802;
  1550.             lnza = lnza * epsilon + -0.5;
  1551.             lnza = lnza * epsilon + 1.0;
  1552.             lnza *= epsilon;
  1553.         }

  1555.         /* Relative sizes:
  1556.          * lnzb     [0, 2.33E-10]
  1557.          * lnm[1]   [0, 1.17E-7]
  1558.          * ln2B*exp [0, 1.12E-4]
  1559.          * lnza      [0, 9.7E-4]
  1560.          * lnm[0]   [0, 0.692]
  1561.          * ln2A*exp [0, 709]
  1562.          */

  1564.         /* Compute the following sum:
  1565.          * lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
  1566.          */

  1568.         //return lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
  1569.         double a = LN_2_A*exp;
  1570.         double b = 0.0;
  1571.         double c = a+lnm[0];
  1572.         double d = -(c-a-lnm[0]);
  1573.         a = c;
  1574.         b += d;

  1576.         c = a + lnza;
  1577.         d = -(c - a - lnza);
  1578.         a = c;
  1579.         b += d;

  1581.         c = a + LN_2_B*exp;
  1582.         d = -(c - a - LN_2_B*exp);
  1583.         a = c;
  1584.         b += d;

  1586.         c = a + lnm[1];
  1587.         d = -(c - a - lnm[1]);
  1588.         a = c;
  1589.         b += d;

  1591.         c = a + lnzb;
  1592.         d = -(c - a - lnzb);
  1593.         a = c;
  1594.         b += d;

  1596.         if (hiPrec != null) {
  1597.             hiPrec[0] = a;
  1598.             hiPrec[1] = b;
  1599.         }

  1601.         return a + b;
  1602.     }

  1604.     /**
  1605.      * Computes log(1 + x).
  1606.      *
  1607.      * @param x Number.
  1608.      * @return {@code log(1 + x)}.
  1609.      */
  1610.     public static double log1p(final double x) {
  1611.         if (x == -1) {
  1612.             return Double.NEGATIVE_INFINITY;
  1613.         }

  1615.         if (x == Double.POSITIVE_INFINITY) {
  1616.             return Double.POSITIVE_INFINITY;
  1617.         }

  1619.         if (x > 1e-6 ||
  1620.             x < -1e-6) {
  1621.             final double xpa = 1 + x;
  1622.             final double xpb = -(xpa - 1 - x);

  1624.             final double[] hiPrec = new double[2];
  1625.             final double lores = log(xpa, hiPrec);
  1626.             if (Double.isInfinite(lores)) { // Don't allow this to be converted to NaN
  1627.                 return lores;
  1628.             }

  1630.             // Do a taylor series expansion around xpa:
  1631.             //   f(x+y) = f(x) + f'(x) y + f''(x)/2 y^2
  1632.             final double fx1 = xpb / xpa;
  1633.             final double epsilon = 0.5 * fx1 + 1;
  1634.             return epsilon * fx1 + hiPrec[1] + hiPrec[0];
  1635.         } else {
  1636.             // Value is small |x| < 1e6, do a Taylor series centered on 1.
  1637.             final double y = (x * F_1_3 - F_1_2) * x + 1;
  1638.             return y * x;
  1639.         }
  1640.     }

  1642.     /** Compute the base 10 logarithm.
  1643.      * @param x a number
  1644.      * @return log10(x)
  1645.      */
  1646.     public static double log10(final double x) {
  1647.         final double[] hiPrec = new double[2];

  1649.         final double lores = log(x, hiPrec);
  1650.         if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN
  1651.             return lores;
  1652.         }

  1654.         final double tmp = hiPrec[0] * HEX_40000000;
  1655.         final double lna = hiPrec[0] + tmp - tmp;
  1656.         final double lnb = hiPrec[0] - lna + hiPrec[1];

  1658.         final double rln10a = 0.4342944622039795;
  1659.         final double rln10b = 1.9699272335463627E-8;

  1661.         return rln10b * lnb + rln10b * lna + rln10a * lnb + rln10a * lna;
  1662.     }

  1664.     /**
  1665.      * Computes the <a href="http://mathworld.wolfram.com/Logarithm.html">
  1666.      * logarithm</a> in a given base.
  1667.      *
  1668.      * Returns {@code NaN} if either argument is negative.
  1669.      * If {@code base} is 0 and {@code x} is positive, 0 is returned.
  1670.      * If {@code base} is positive and {@code x} is 0,
  1671.      * {@code Double.NEGATIVE_INFINITY} is returned.
  1672.      * If both arguments are 0, the result is {@code NaN}.
  1673.      *
  1674.      * @param base Base of the logarithm, must be greater than 0.
  1675.      * @param x Argument, must be greater than 0.
  1676.      * @return the value of the logarithm, i.e. the number {@code y} such that
  1677.      * <code>base<sup>y</sup> = x</code>.
  1678.      */
  1679.     public static double log(double base, double x) {
  1680.         return log(x) / log(base);
  1681.     }

  1683.     /**
  1684.      * Power function.  Compute x^y.
  1685.      *
  1686.      * @param x   a double
  1687.      * @param y   a double
  1688.      * @return double
  1689.      */
  1690.     public static double pow(final double x, final double y) {

  1692.         if (y == 0) {
  1693.             // y = -0 or y = +0
  1694.             return 1.0;
  1695.         } else {

  1697.             final long yBits        = Double.doubleToRawLongBits(y);
  1698.             final int  yRawExp      = (int) ((yBits & MASK_DOUBLE_EXPONENT) >> 52);
  1699.             final long yRawMantissa = yBits & MASK_DOUBLE_MANTISSA;
  1700.             final long xBits        = Double.doubleToRawLongBits(x);
  1701.             final int  xRawExp      = (int) ((xBits & MASK_DOUBLE_EXPONENT) >> 52);
  1702.             final long xRawMantissa = xBits & MASK_DOUBLE_MANTISSA;

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

  1707.                 if ((yRawExp == 2047 && yRawMantissa != 0) ||
  1708.                     (xRawExp == 2047 && xRawMantissa != 0)) {
  1709.                     // NaN
  1710.                     return Double.NaN;
  1711.                 } else if (xRawExp == 1023 && xRawMantissa == 0) {
  1712.                     // x = -1.0 or x = +1.0
  1713.                     if (yRawExp == 2047) {
  1714.                         // y is infinite
  1715.                         return Double.NaN;
  1716.                     } else {
  1717.                         // y is a large even integer
  1718.                         return 1.0;
  1719.                     }
  1720.                 } else {
  1721.                     // the absolute value of x is either greater or smaller than 1.0

  1723.                     // if yRawExp == 2047 and mantissa is 0, y = -infinity or y = +infinity
  1724.                     // if 1085 < yRawExp < 2047, y is simply a large number, however, due to limited
  1725.                     // accuracy, at this magnitude it behaves just like infinity with regards to x
  1726.                     if ((y > 0) ^ (xRawExp < 1023)) {
  1727.                         // either y = +infinity (or large engouh) and abs(x) > 1.0
  1728.                         // or     y = -infinity (or large engouh) and abs(x) < 1.0
  1729.                         return Double.POSITIVE_INFINITY;
  1730.                     } else {
  1731.                         // either y = +infinity (or large engouh) and abs(x) < 1.0
  1732.                         // or     y = -infinity (or large engouh) and abs(x) > 1.0
  1733.                         return +0.0;
  1734.                     }
  1735.                 }

  1737.             } else {
  1738.                 // y is a regular non-zero number

  1740.                 if (yRawExp >= 1023) {
  1741.                     // y may be an integral value, which should be handled specifically
  1742.                     final long yFullMantissa = IMPLICIT_HIGH_BIT | yRawMantissa;
  1743.                     if (yRawExp < 1075) {
  1744.                         // normal number with negative shift that may have a fractional part
  1745.                         final long integralMask = (-1L) << (1075 - yRawExp);
  1746.                         if ((yFullMantissa & integralMask) == yFullMantissa) {
  1747.                             // all fractional bits are 0, the number is really integral
  1748.                             final long l = yFullMantissa >> (1075 - yRawExp);
  1749.                             return pow(x, (y < 0) ? -l : l);
  1750.                         }
  1751.                     } else {
  1752.                         // normal number with positive shift, always an integral value
  1753.                         // we know it fits in a primitive long because yRawExp > 1085 has been handled above
  1754.                         final long l =  yFullMantissa << (yRawExp - 1075);
  1755.                         return pow(x, (y < 0) ? -l : l);
  1756.                     }
  1757.                 }

  1759.                 // y is a non-integral value

  1761.                 if (x == 0) {
  1762.                     // x = -0 or x = +0
  1763.                     // the integer powers have already been handled above
  1764.                     return y < 0 ? Double.POSITIVE_INFINITY : +0.0;
  1765.                 } else if (xRawExp == 2047) {
  1766.                     if (xRawMantissa == 0) {
  1767.                         // x = -infinity or x = +infinity
  1768.                         return (y < 0) ? +0.0 : Double.POSITIVE_INFINITY;
  1769.                     } else {
  1770.                         // NaN
  1771.                         return Double.NaN;
  1772.                     }
  1773.                 } else if (x < 0) {
  1774.                     // the integer powers have already been handled above
  1775.                     return Double.NaN;
  1776.                 } else {

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

  1780.                     // Split y into ya and yb such that y = ya+yb
  1781.                     final double tmp = y * HEX_40000000;
  1782.                     final double ya = (y + tmp) - tmp;
  1783.                     final double yb = y - ya;

  1785.                     /* Compute ln(x) */
  1786.                     final double[] lns = new double[2];
  1787.                     final double lores = log(x, lns);
  1788.                     if (Double.isInfinite(lores)) { // don't allow this to be converted to NaN
  1789.                         return lores;
  1790.                     }

  1792.                     double lna = lns[0];
  1793.                     double lnb = lns[1];

  1795.                     /* resplit lns */
  1796.                     final double tmp1 = lna * HEX_40000000;
  1797.                     final double tmp2 = (lna + tmp1) - tmp1;
  1798.                     lnb += lna - tmp2;
  1799.                     lna = tmp2;

  1801.                     // y*ln(x) = (aa+ab)
  1802.                     final double aa = lna * ya;
  1803.                     final double ab = lna * yb + lnb * ya + lnb * yb;

  1805.                     lna = aa+ab;
  1806.                     lnb = -(lna - aa - ab);

  1808.                     double z = 1.0 / 120.0;
  1809.                     z = z * lnb + (1.0 / 24.0);
  1810.                     z = z * lnb + (1.0 / 6.0);
  1811.                     z = z * lnb + 0.5;
  1812.                     z = z * lnb + 1.0;
  1813.                     z *= lnb;

  1815.                     return exp(lna, z, null);

  1817.                 }
  1818.             }

  1820.         }

  1822.     }

  1824.     /**
  1825.      * Raise a double to an int power.
  1826.      *
  1827.      * @param d Number to raise.
  1828.      * @param e Exponent.
  1829.      * @return d<sup>e</sup>
  1830.      */
  1831.     public static double pow(double d, int e) {
  1832.         return pow(d, (long) e);
  1833.     }

  1835.     /**
  1836.      * Raise a double to a long power.
  1837.      *
  1838.      * @param d Number to raise.
  1839.      * @param e Exponent.
  1840.      * @return d<sup>e</sup>
  1841.      */
  1842.     public static double pow(double d, long e) {
  1843.         if (e == 0) {
  1844.             return 1.0;
  1845.         } else if (e > 0) {
  1846.             return new Split(d).pow(e).full;
  1847.         } else {
  1848.             return new Split(d).reciprocal().pow(-e).full;
  1849.         }
  1850.     }

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

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

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

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

  1864.         /** Full number. */
  1865.         private final double full;

  1867.         /** High order bits. */
  1868.         private final double high;

  1870.         /** Low order bits. */
  1871.         private final double low;

  1873.         /** Simple constructor.
  1874.          * @param x number to split
  1875.          */
  1876.         Split(final double x) {
  1877.             full = x;
  1878.             high = Double.longBitsToDouble(Double.doubleToRawLongBits(x) & ((-1L) << 27));
  1879.             low  = x - high;
  1880.         }

  1882.         /** Simple constructor.
  1883.          * @param high high order bits
  1884.          * @param low low order bits
  1885.          */
  1886.         Split(final double high, final double low) {
  1887.             this(high == 0.0 ? (low == 0.0 && Double.doubleToRawLongBits(high) == Long.MIN_VALUE /* negative zero */ ? -0.0 : low) : high + low, high, low);
  1888.         }

  1890.         /** Simple constructor.
  1891.          * @param full full number
  1892.          * @param high high order bits
  1893.          * @param low low order bits
  1894.          */
  1895.         Split(final double full, final double high, final double low) {
  1896.             this.full = full;
  1897.             this.high = high;
  1898.             this.low  = low;
  1899.         }

  1901.         /** Multiply the instance by another one.
  1902.          * @param b other instance to multiply by
  1903.          * @return product
  1904.          */
  1905.         public Split multiply(final Split b) {
  1906.             // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
  1907.             final Split  mulBasic  = new Split(full * b.full);
  1908.             final double mulError  = low * b.low - (((mulBasic.full - high * b.high) - low * b.high) - high * b.low);
  1909.             return new Split(mulBasic.high, mulBasic.low + mulError);
  1910.         }

  1912.         /** Compute the reciprocal of the instance.
  1913.          * @return reciprocal of the instance
  1914.          */
  1915.         public Split reciprocal() {

  1917.             final double approximateInv = 1.0 / full;
  1918.             final Split  splitInv       = new Split(approximateInv);

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

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

  1929.         }

  1931.         /** Computes this^e.
  1932.          * @param e exponent (beware, here it MUST be > 0; the only exclusion is Long.MIN_VALUE)
  1933.          * @return d^e, split in high and low bits
  1934.          */
  1935.         private Split pow(final long e) {

  1937.             // prepare result
  1938.             Split result = new Split(1);

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

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

  1945.                 if ((p & 0x1) != 0) {
  1946.                     // accurate multiplication result = result * d^(2p) using Veltkamp TwoProduct algorithm
  1947.                     result = result.multiply(d2p);
  1948.                 }

  1950.                 // accurate squaring d^(2(p+1)) = d^(2p) * d^(2p) using Veltkamp TwoProduct algorithm
  1951.                 d2p = d2p.multiply(d2p);

  1953.             }

  1955.             if (Double.isNaN(result.full)) {
  1956.                 if (Double.isNaN(full)) {
  1957.                     return NAN;
  1958.                 } else {
  1959.                     // some intermediate numbers exceeded capacity,
  1960.                     // and the low order bits became NaN (because infinity - infinity = NaN)
  1961.                     if (abs(full) < 1) {
  1962.                         return new Split(copySign(0.0, full), 0.0);
  1963.                     } else if (full < 0 && (e & 0x1) == 1) {
  1964.                         return NEGATIVE_INFINITY;
  1965.                     } else {
  1966.                         return POSITIVE_INFINITY;
  1967.                     }
  1968.                 }
  1969.             } else {
  1970.                 return result;
  1971.             }

  1973.         }

  1975.     }

  1977.     /**
  1978.      *  Computes sin(x) - x, where |x| < 1/16.
  1979.      *  Use a Remez polynomial approximation.
  1980.      *  @param x a number smaller than 1/16
  1981.      *  @return sin(x) - x
  1982.      */
  1983.     private static double polySine(final double x)
  1984.     {
  1985.         double x2 = x*x;

  1987.         double p = 2.7553817452272217E-6;
  1988.         p = p * x2 + -1.9841269659586505E-4;
  1989.         p = p * x2 + 0.008333333333329196;
  1990.         p = p * x2 + -0.16666666666666666;
  1991.         //p *= x2;
  1992.         //p *= x;
  1993.         p = p * x2 * x;

  1995.         return p;
  1996.     }

  1998.     /**
  1999.      *  Computes cos(x) - 1, where |x| < 1/16.
  2000.      *  Use a Remez polynomial approximation.
  2001.      *  @param x a number smaller than 1/16
  2002.      *  @return cos(x) - 1
  2003.      */
  2004.     private static double polyCosine(double x) {
  2005.         double x2 = x*x;

  2007.         double p = 2.479773539153719E-5;
  2008.         p = p * x2 + -0.0013888888689039883;
  2009.         p = p * x2 + 0.041666666666621166;
  2010.         p = p * x2 + -0.49999999999999994;
  2011.         p *= x2;

  2013.         return p;
  2014.     }

  2016.     /**
  2017.      *  Compute sine over the first quadrant (0 < x < pi/2).
  2018.      *  Use combination of table lookup and rational polynomial expansion.
  2019.      *  @param xa number from which sine is requested
  2020.      *  @param xb extra bits for x (may be 0.0)
  2021.      *  @return sin(xa + xb)
  2022.      */
  2023.     private static double sinQ(double xa, double xb) {
  2024.         int idx = (int) ((xa * 8.0) + 0.5);
  2025.         final double epsilon = xa - EIGHTHS[idx]; //idx*0.125;

  2027.         // Table lookups
  2028.         final double sintA = SINE_TABLE_A[idx];
  2029.         final double sintB = SINE_TABLE_B[idx];
  2030.         final double costA = COSINE_TABLE_A[idx];
  2031.         final double costB = COSINE_TABLE_B[idx];

  2033.         // Polynomial eval of sin(epsilon), cos(epsilon)
  2034.         double sinEpsA = epsilon;
  2035.         double sinEpsB = polySine(epsilon);
  2036.         final double cosEpsA = 1.0;
  2037.         final double cosEpsB = polyCosine(epsilon);

  2039.         // Split epsilon   xa + xb = x
  2040.         final double temp = sinEpsA * HEX_40000000;
  2041.         double temp2 = (sinEpsA + temp) - temp;
  2042.         sinEpsB +=  sinEpsA - temp2;
  2043.         sinEpsA = temp2;

  2045.         /* Compute sin(x) by angle addition formula */
  2046.         double result;

  2048.         /* Compute the following sum:
  2049.          *
  2050.          * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
  2051.          *          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
  2052.          *
  2053.          * Ranges of elements
  2054.          *
  2055.          * xxxtA   0            PI/2
  2056.          * xxxtB   -1.5e-9      1.5e-9
  2057.          * sinEpsA -0.0625      0.0625
  2058.          * sinEpsB -6e-11       6e-11
  2059.          * cosEpsA  1.0
  2060.          * cosEpsB  0           -0.0625
  2061.          *
  2062.          */

  2064.         //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
  2065.         //          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;

  2067.         //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
  2068.         //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
  2069.         double a = 0;
  2070.         double b = 0;

  2072.         double t = sintA;
  2073.         double c = a + t;
  2074.         double d = -(c - a - t);
  2075.         a = c;
  2076.         b += d;

  2078.         t = costA * sinEpsA;
  2079.         c = a + t;
  2080.         d = -(c - a - t);
  2081.         a = c;
  2082.         b += d;

  2084.         b = b + sintA * cosEpsB + costA * sinEpsB;
  2085.         /*
  2086.     t = sintA*cosEpsB;
  2087.     c = a + t;
  2088.     d = -(c - a - t);
  2089.     a = c;
  2090.     b = b + d;

  2092.     t = costA*sinEpsB;
  2093.     c = a + t;
  2094.     d = -(c - a - t);
  2095.     a = c;
  2096.     b = b + d;
  2097.          */

  2099.         b = b + sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB;
  2100.         /*
  2101.     t = sintB;
  2102.     c = a + t;
  2103.     d = -(c - a - t);
  2104.     a = c;
  2105.     b = b + d;

  2107.     t = costB*sinEpsA;
  2108.     c = a + t;
  2109.     d = -(c - a - t);
  2110.     a = c;
  2111.     b = b + d;

  2113.     t = sintB*cosEpsB;
  2114.     c = a + t;
  2115.     d = -(c - a - t);
  2116.     a = c;
  2117.     b = b + d;

  2119.     t = costB*sinEpsB;
  2120.     c = a + t;
  2121.     d = -(c - a - t);
  2122.     a = c;
  2123.     b = b + d;
  2124.          */

  2126.         if (xb != 0.0) {
  2127.             t = ((costA + costB) * (cosEpsA + cosEpsB) -
  2128.                  (sintA + sintB) * (sinEpsA + sinEpsB)) * xb;  // approximate cosine*xb
  2129.             c = a + t;
  2130.             d = -(c - a - t);
  2131.             a = c;
  2132.             b += d;
  2133.         }

  2135.         result = a + b;

  2137.         return result;
  2138.     }

  2140.     /**
  2141.      * Compute cosine in the first quadrant by subtracting input from PI/2 and
  2142.      * then calling sinQ.  This is more accurate as the input approaches PI/2.
  2143.      *  @param xa number from which cosine is requested
  2144.      *  @param xb extra bits for x (may be 0.0)
  2145.      *  @return cos(xa + xb)
  2146.      */
  2147.     private static double cosQ(double xa, double xb) {
  2148.         final double pi2a = 1.5707963267948966;
  2149.         final double pi2b = 6.123233995736766E-17;

  2151.         final double a = pi2a - xa;
  2152.         double b = -(a - pi2a + xa);
  2153.         b += pi2b - xb;

  2155.         return sinQ(a, b);
  2156.     }

  2158.     /**
  2159.      *  Compute tangent (or cotangent) over the first quadrant.   0 < x < pi/2
  2160.      *  Use combination of table lookup and rational polynomial expansion.
  2161.      *  @param xa number from which sine is requested
  2162.      *  @param xb extra bits for x (may be 0.0)
  2163.      *  @param cotanFlag if true, compute the cotangent instead of the tangent
  2164.      *  @return tan(xa+xb) (or cotangent, depending on cotanFlag)
  2165.      */
  2166.     private static double tanQ(double xa, double xb, boolean cotanFlag) {

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

  2171.         // Table lookups
  2172.         final double sintA = SINE_TABLE_A[idx];
  2173.         final double sintB = SINE_TABLE_B[idx];
  2174.         final double costA = COSINE_TABLE_A[idx];
  2175.         final double costB = COSINE_TABLE_B[idx];

  2177.         // Polynomial eval of sin(epsilon), cos(epsilon)
  2178.         double sinEpsA = epsilon;
  2179.         double sinEpsB = polySine(epsilon);
  2180.         final double cosEpsA = 1.0;
  2181.         final double cosEpsB = polyCosine(epsilon);

  2183.         // Split epsilon   xa + xb = x
  2184.         double temp = sinEpsA * HEX_40000000;
  2185.         double temp2 = (sinEpsA + temp) - temp;
  2186.         sinEpsB +=  sinEpsA - temp2;
  2187.         sinEpsA = temp2;

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

  2191.         /* Compute the following sum:
  2192.          *
  2193.          * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
  2194.          *          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
  2195.          *
  2196.          * Ranges of elements
  2197.          *
  2198.          * xxxtA   0            PI/2
  2199.          * xxxtB   -1.5e-9      1.5e-9
  2200.          * sinEpsA -0.0625      0.0625
  2201.          * sinEpsB -6e-11       6e-11
  2202.          * cosEpsA  1.0
  2203.          * cosEpsB  0           -0.0625
  2204.          *
  2205.          */

  2207.         //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
  2208.         //          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;

  2210.         //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
  2211.         //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
  2212.         double a = 0;
  2213.         double b = 0;

  2215.         // Compute sine
  2216.         double t = sintA;
  2217.         double c = a + t;
  2218.         double d = -(c - a - t);
  2219.         a = c;
  2220.         b += d;

  2222.         t = costA*sinEpsA;
  2223.         c = a + t;
  2224.         d = -(c - a - t);
  2225.         a = c;
  2226.         b += d;

  2228.         b += sintA*cosEpsB + costA*sinEpsB;
  2229.         b += sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;

  2231.         double sina = a + b;
  2232.         double sinb = -(sina - a - b);

  2234.         // Compute cosine

  2236.         a = 0.0;
  2237.         b = 0.0;

  2239.         t = costA*cosEpsA;
  2240.         c = a + t;
  2241.         d = -(c - a - t);
  2242.         a = c;
  2243.         b += d;

  2245.         t = -sintA*sinEpsA;
  2246.         c = a + t;
  2247.         d = -(c - a - t);
  2248.         a = c;
  2249.         b += d;

  2251.         b += costB*cosEpsA + costA*cosEpsB + costB*cosEpsB;
  2252.         b -= sintB*sinEpsA + sintA*sinEpsB + sintB*sinEpsB;

  2254.         double cosa = a + b;
  2255.         double cosb = -(cosa - a - b);

  2257.         if (cotanFlag) {
  2258.             double tmp;
  2259.             tmp = cosa; cosa = sina; sina = tmp;
  2260.             tmp = cosb; cosb = sinb; sinb = tmp;
  2261.         }


  2264.         /* estimate and correct, compute 1.0/(cosa+cosb) */
  2265.         /*
  2266.     double est = (sina+sinb)/(cosa+cosb);
  2267.     double err = (sina - cosa*est) + (sinb - cosb*est);
  2268.     est += err/(cosa+cosb);
  2269.     err = (sina - cosa*est) + (sinb - cosb*est);
  2270.          */

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

  2274.         double est = sina/cosa;

  2276.         /* Split the estimate to get more accurate read on division rounding */
  2277.         temp = est * HEX_40000000;
  2278.         double esta = (est + temp) - temp;
  2279.         double estb =  est - esta;

  2281.         temp = cosa * HEX_40000000;
  2282.         double cosaa = (cosa + temp) - temp;
  2283.         double cosab =  cosa - cosaa;

  2285.         //double err = (sina - est*cosa)/cosa;  // Correction for division rounding
  2286.         double err = (sina - esta*cosaa - esta*cosab - estb*cosaa - estb*cosab)/cosa;  // Correction for division rounding
  2287.         err += sinb/cosa;                     // Change in est due to sinb
  2288.         err += -sina * cosb / cosa / cosa;    // Change in est due to cosb

  2290.         if (xb != 0.0) {
  2291.             // tan' = 1 + tan^2      cot' = -(1 + cot^2)
  2292.             // Approximate impact of xb
  2293.             double xbadj = xb + est*est*xb;
  2294.             if (cotanFlag) {
  2295.                 xbadj = -xbadj;
  2296.             }

  2298.             err += xbadj;
  2299.         }

  2301.         return est+err;
  2302.     }

  2304.     /** Reduce the input argument using the Payne and Hanek method.
  2305.      *  This is good for all inputs 0.0 < x < inf
  2306.      *  Output is remainder after dividing by PI/2
  2307.      *  The result array should contain 3 numbers.
  2308.      *  result[0] is the integer portion, so mod 4 this gives the quadrant.
  2309.      *  result[1] is the upper bits of the remainder
  2310.      *  result[2] is the lower bits of the remainder
  2311.      *
  2312.      * @param x number to reduce
  2313.      * @param result placeholder where to put the result
  2314.      */
  2315.     private static void reducePayneHanek(double x, double[] result)
  2316.     {
  2317.         /* Convert input double to bits */
  2318.         long inbits = Double.doubleToRawLongBits(x);
  2319.         int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;

  2321.         /* Convert to fixed point representation */
  2322.         inbits &= 0x000fffffffffffffL;
  2323.         inbits |= 0x0010000000000000L;

  2325.         /* Normalize input to be between 0.5 and 1.0 */
  2326.         exponent++;
  2327.         inbits <<= 11;

  2329.         /* Based on the exponent, get a shifted copy of recip2pi */
  2330.         long shpi0;
  2331.         long shpiA;
  2332.         long shpiB;
  2333.         int idx = exponent >> 6;
  2334.         int shift = exponent - (idx << 6);

  2336.         if (shift != 0) {
  2337.             shpi0 = (idx == 0) ? 0 : (RECIP_2PI[idx-1] << shift);
  2338.             shpi0 |= RECIP_2PI[idx] >>> (64-shift);
  2339.             shpiA = (RECIP_2PI[idx] << shift) | (RECIP_2PI[idx+1] >>> (64-shift));
  2340.             shpiB = (RECIP_2PI[idx+1] << shift) | (RECIP_2PI[idx+2] >>> (64-shift));
  2341.         } else {
  2342.             shpi0 = (idx == 0) ? 0 : RECIP_2PI[idx-1];
  2343.             shpiA = RECIP_2PI[idx];
  2344.             shpiB = RECIP_2PI[idx+1];
  2345.         }

  2347.         /* Multiply input by shpiA */
  2348.         long a = inbits >>> 32;
  2349.         long b = inbits & 0xffffffffL;

  2351.         long c = shpiA >>> 32;
  2352.         long d = shpiA & 0xffffffffL;

  2354.         long ac = a * c;
  2355.         long bd = b * d;
  2356.         long bc = b * c;
  2357.         long ad = a * d;

  2359.         long prodB = bd + (ad << 32);
  2360.         long prodA = ac + (ad >>> 32);

  2362.         boolean bita = (bd & 0x8000000000000000L) != 0;
  2363.         boolean bitb = (ad & 0x80000000L ) != 0;
  2364.         boolean bitsum = (prodB & 0x8000000000000000L) != 0;

  2366.         /* Carry */
  2367.         if ( (bita && bitb) ||
  2368.                 ((bita || bitb) && !bitsum) ) {
  2369.             prodA++;
  2370.         }

  2372.         bita = (prodB & 0x8000000000000000L) != 0;
  2373.         bitb = (bc & 0x80000000L ) != 0;

  2375.         prodB += bc << 32;
  2376.         prodA += bc >>> 32;

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

  2380.         /* Carry */
  2381.         if ( (bita && bitb) ||
  2382.                 ((bita || bitb) && !bitsum) ) {
  2383.             prodA++;
  2384.         }

  2386.         /* Multiply input by shpiB */
  2387.         c = shpiB >>> 32;
  2388.         d = shpiB & 0xffffffffL;
  2389.         ac = a * c;
  2390.         bc = b * c;
  2391.         ad = a * d;

  2393.         /* Collect terms */
  2394.         ac += (bc + ad) >>> 32;

  2396.         bita = (prodB & 0x8000000000000000L) != 0;
  2397.         bitb = (ac & 0x8000000000000000L ) != 0;
  2398.         prodB += ac;
  2399.         bitsum = (prodB & 0x8000000000000000L) != 0;
  2400.         /* Carry */
  2401.         if ( (bita && bitb) ||
  2402.                 ((bita || bitb) && !bitsum) ) {
  2403.             prodA++;
  2404.         }

  2406.         /* Multiply by shpi0 */
  2407.         c = shpi0 >>> 32;
  2408.         d = shpi0 & 0xffffffffL;

  2410.         bd = b * d;
  2411.         bc = b * c;
  2412.         ad = a * d;

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

  2416.         /*
  2417.          * prodA, prodB now contain the remainder as a fraction of PI.  We want this as a fraction of
  2418.          * PI/2, so use the following steps:
  2419.          * 1.) multiply by 4.
  2420.          * 2.) do a fixed point muliply by PI/4.
  2421.          * 3.) Convert to floating point.
  2422.          * 4.) Multiply by 2
  2423.          */

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

  2428.         /* Multiply by 4 */
  2429.         prodA <<= 2;
  2430.         prodA |= prodB >>> 62;
  2431.         prodB <<= 2;

  2433.         /* Multiply by PI/4 */
  2434.         a = prodA >>> 32;
  2435.         b = prodA & 0xffffffffL;

  2437.         c = PI_O_4_BITS[0] >>> 32;
  2438.         d = PI_O_4_BITS[0] & 0xffffffffL;

  2440.         ac = a * c;
  2441.         bd = b * d;
  2442.         bc = b * c;
  2443.         ad = a * d;

  2445.         long prod2B = bd + (ad << 32);
  2446.         long prod2A = ac + (ad >>> 32);

  2448.         bita = (bd & 0x8000000000000000L) != 0;
  2449.         bitb = (ad & 0x80000000L ) != 0;
  2450.         bitsum = (prod2B & 0x8000000000000000L) != 0;

  2452.         /* Carry */
  2453.         if ( (bita && bitb) ||
  2454.                 ((bita || bitb) && !bitsum) ) {
  2455.             prod2A++;
  2456.         }

  2458.         bita = (prod2B & 0x8000000000000000L) != 0;
  2459.         bitb = (bc & 0x80000000L ) != 0;

  2461.         prod2B += bc << 32;
  2462.         prod2A += bc >>> 32;

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

  2466.         /* Carry */
  2467.         if ( (bita && bitb) ||
  2468.                 ((bita || bitb) && !bitsum) ) {
  2469.             prod2A++;
  2470.         }

  2472.         /* Multiply input by pio4bits[1] */
  2473.         c = PI_O_4_BITS[1] >>> 32;
  2474.         d = PI_O_4_BITS[1] & 0xffffffffL;
  2475.         ac = a * c;
  2476.         bc = b * c;
  2477.         ad = a * d;

  2479.         /* Collect terms */
  2480.         ac += (bc + ad) >>> 32;

  2482.         bita = (prod2B & 0x8000000000000000L) != 0;
  2483.         bitb = (ac & 0x8000000000000000L ) != 0;
  2484.         prod2B += ac;
  2485.         bitsum = (prod2B & 0x8000000000000000L) != 0;
  2486.         /* Carry */
  2487.         if ( (bita && bitb) ||
  2488.                 ((bita || bitb) && !bitsum) ) {
  2489.             prod2A++;
  2490.         }

  2492.         /* Multiply inputB by pio4bits[0] */
  2493.         a = prodB >>> 32;
  2494.         b = prodB & 0xffffffffL;
  2495.         c = PI_O_4_BITS[0] >>> 32;
  2496.         d = PI_O_4_BITS[0] & 0xffffffffL;
  2497.         ac = a * c;
  2498.         bc = b * c;
  2499.         ad = a * d;

  2501.         /* Collect terms */
  2502.         ac += (bc + ad) >>> 32;

  2504.         bita = (prod2B & 0x8000000000000000L) != 0;
  2505.         bitb = (ac & 0x8000000000000000L ) != 0;
  2506.         prod2B += ac;
  2507.         bitsum = (prod2B & 0x8000000000000000L) != 0;
  2508.         /* Carry */
  2509.         if ( (bita && bitb) ||
  2510.                 ((bita || bitb) && !bitsum) ) {
  2511.             prod2A++;
  2512.         }

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

  2518.         double sumA = tmpA + tmpB;
  2519.         double sumB = -(sumA - tmpA - tmpB);

  2521.         /* Multiply by PI/2 and return */
  2522.         result[0] = intPart;
  2523.         result[1] = sumA * 2.0;
  2524.         result[2] = sumB * 2.0;
  2525.     }

  2527.     /**
  2528.      * Sine function.
  2529.      *
  2530.      * @param x Argument.
  2531.      * @return sin(x)
  2532.      */
  2533.     public static double sin(double x) {
  2534.         boolean negative = false;
  2535.         int quadrant = 0;
  2536.         double xa;
  2537.         double xb = 0.0;

  2539.         /* Take absolute value of the input */
  2540.         xa = x;
  2541.         if (x < 0) {
  2542.             negative = true;
  2543.             xa = -xa;
  2544.         }

  2546.         /* Check for zero and negative zero */
  2547.         if (xa == 0.0) {
  2548.             long bits = Double.doubleToRawLongBits(x);
  2549.             if (bits < 0) {
  2550.                 return -0.0;
  2551.             }
  2552.             return 0.0;
  2553.         }

  2555.         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
  2556.             return Double.NaN;
  2557.         }

  2559.         /* Perform any argument reduction */
  2560.         if (xa > 3294198.0) {
  2561.             // PI * (2**20)
  2562.             // Argument too big for CodyWaite reduction.  Must use
  2563.             // PayneHanek.
  2564.             double[] reduceResults = new double[3];
  2565.             reducePayneHanek(xa, reduceResults);
  2566.             quadrant = ((int) reduceResults[0]) & 3;
  2567.             xa = reduceResults[1];
  2568.             xb = reduceResults[2];
  2569.         } else if (xa > 1.5707963267948966) {
  2570.             final CodyWaite cw = new CodyWaite(xa);
  2571.             quadrant = cw.getK() & 3;
  2572.             xa = cw.getRemA();
  2573.             xb = cw.getRemB();
  2574.         }

  2576.         if (negative) {
  2577.             quadrant ^= 2;  // Flip bit 1
  2578.         }

  2580.         switch (quadrant) {
  2581.             case 0:
  2582.                 return sinQ(xa, xb);
  2583.             case 1:
  2584.                 return cosQ(xa, xb);
  2585.             case 2:
  2586.                 return -sinQ(xa, xb);
  2587.             case 3:
  2588.                 return -cosQ(xa, xb);
  2589.             default:
  2590.                 return Double.NaN;
  2591.         }
  2592.     }

  2594.     /**
  2595.      * Cosine function.
  2596.      *
  2597.      * @param x Argument.
  2598.      * @return cos(x)
  2599.      */
  2600.     public static double cos(double x) {
  2601.         int quadrant = 0;

  2603.         /* Take absolute value of the input */
  2604.         double xa = x;
  2605.         if (x < 0) {
  2606.             xa = -xa;
  2607.         }

  2609.         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
  2610.             return Double.NaN;
  2611.         }

  2613.         /* Perform any argument reduction */
  2614.         double xb = 0;
  2615.         if (xa > 3294198.0) {
  2616.             // PI * (2**20)
  2617.             // Argument too big for CodyWaite reduction.  Must use
  2618.             // PayneHanek.
  2619.             double[] reduceResults = new double[3];
  2620.             reducePayneHanek(xa, reduceResults);
  2621.             quadrant = ((int) reduceResults[0]) & 3;
  2622.             xa = reduceResults[1];
  2623.             xb = reduceResults[2];
  2624.         } else if (xa > 1.5707963267948966) {
  2625.             final CodyWaite cw = new CodyWaite(xa);
  2626.             quadrant = cw.getK() & 3;
  2627.             xa = cw.getRemA();
  2628.             xb = cw.getRemB();
  2629.         }

  2631.         //if (negative)
  2632.         //  quadrant = (quadrant + 2) % 4;

  2634.         switch (quadrant) {
  2635.             case 0:
  2636.                 return cosQ(xa, xb);
  2637.             case 1:
  2638.                 return -sinQ(xa, xb);
  2639.             case 2:
  2640.                 return -cosQ(xa, xb);
  2641.             case 3:
  2642.                 return sinQ(xa, xb);
  2643.             default:
  2644.                 return Double.NaN;
  2645.         }
  2646.     }

  2648.     /**
  2649.      * Combined Sine and Cosine function.
  2650.      *
  2651.      * @param x Argument.
  2652.      * @return [sin(x), cos(x)]
  2653.      */
  2654.     public static SinCos sinCos(double x) {
  2655.         boolean negative = false;
  2656.         int quadrant = 0;
  2657.         double xa;
  2658.         double xb = 0.0;

  2660.         /* Take absolute value of the input */
  2661.         xa = x;
  2662.         if (x < 0) {
  2663.             negative = true;
  2664.             xa = -xa;
  2665.         }

  2667.         /* Check for zero and negative zero */
  2668.         if (xa == 0.0) {
  2669.             long bits = Double.doubleToRawLongBits(x);
  2670.             if (bits < 0) {
  2671.                 return new SinCos(-0.0, 1.0);
  2672.             }
  2673.             return new SinCos(0.0, 1.0);
  2674.         }

  2676.         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
  2677.             return new SinCos(Double.NaN, Double.NaN);
  2678.         }

  2680.         /* Perform any argument reduction */
  2681.         if (xa > 3294198.0) {
  2682.             // PI * (2**20)
  2683.             // Argument too big for CodyWaite reduction.  Must use
  2684.             // PayneHanek.
  2685.             double[] reduceResults = new double[3];
  2686.             reducePayneHanek(xa, reduceResults);
  2687.             quadrant = ((int) reduceResults[0]) & 3;
  2688.             xa = reduceResults[1];
  2689.             xb = reduceResults[2];
  2690.         } else if (xa > 1.5707963267948966) {
  2691.             final CodyWaite cw = new CodyWaite(xa);
  2692.             quadrant = cw.getK() & 3;
  2693.             xa = cw.getRemA();
  2694.             xb = cw.getRemB();
  2695.         }

  2697.         switch (quadrant) {
  2698.             case 0:
  2699.                 return new SinCos(negative ? -sinQ(xa, xb) :  sinQ(xa, xb),  cosQ(xa, xb));
  2700.             case 1:
  2701.                 return new SinCos(negative ? -cosQ(xa, xb) :  cosQ(xa, xb), -sinQ(xa, xb));
  2702.             case 2:
  2703.                 return new SinCos(negative ?  sinQ(xa, xb) : -sinQ(xa, xb), -cosQ(xa, xb));
  2704.             case 3:
  2705.                 return new SinCos(negative ?  cosQ(xa, xb) : -cosQ(xa, xb),  sinQ(xa, xb));
  2706.             default:
  2707.                 return new SinCos(Double.NaN, Double.NaN);
  2708.         }
  2709.     }

  2711.     /**
  2712.      * Combined Sine and Cosine function.
  2713.      *
  2714.      * @param x Argument.
  2715.      * @param <T> the type of the field element
  2716.      * @return [sin(x), cos(x)]
  2717.      * @since 1.4
  2718.      */
  2719.     public static <T extends CalculusFieldElement<T>> FieldSinCos<T> sinCos(T x) {
  2720.         return x.sinCos();
  2721.     }

  2723.     /**
  2724.      * Tangent function.
  2725.      *
  2726.      * @param x Argument.
  2727.      * @return tan(x)
  2728.      */
  2729.     public static double tan(double x) {
  2730.         boolean negative = false;
  2731.         int quadrant = 0;

  2733.         /* Take absolute value of the input */
  2734.         double xa = x;
  2735.         if (x < 0) {
  2736.             negative = true;
  2737.             xa = -xa;
  2738.         }

  2740.         /* Check for zero and negative zero */
  2741.         if (xa == 0.0) {
  2742.             long bits = Double.doubleToRawLongBits(x);
  2743.             if (bits < 0) {
  2744.                 return -0.0;
  2745.             }
  2746.             return 0.0;
  2747.         }

  2749.         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
  2750.             return Double.NaN;
  2751.         }

  2753.         /* Perform any argument reduction */
  2754.         double xb = 0;
  2755.         if (xa > 3294198.0) {
  2756.             // PI * (2**20)
  2757.             // Argument too big for CodyWaite reduction.  Must use
  2758.             // PayneHanek.
  2759.             double[] reduceResults = new double[3];
  2760.             reducePayneHanek(xa, reduceResults);
  2761.             quadrant = ((int) reduceResults[0]) & 3;
  2762.             xa = reduceResults[1];
  2763.             xb = reduceResults[2];
  2764.         } else if (xa > 1.5707963267948966) {
  2765.             final CodyWaite cw = new CodyWaite(xa);
  2766.             quadrant = cw.getK() & 3;
  2767.             xa = cw.getRemA();
  2768.             xb = cw.getRemB();
  2769.         }

  2771.         if (xa > 1.5) {
  2772.             // Accuracy suffers between 1.5 and PI/2
  2773.             final double pi2a = 1.5707963267948966;
  2774.             final double pi2b = 6.123233995736766E-17;

  2776.             final double a = pi2a - xa;
  2777.             double b = -(a - pi2a + xa);
  2778.             b += pi2b - xb;

  2780.             xa = a + b;
  2781.             xb = -(xa - a - b);
  2782.             quadrant ^= 1;
  2783.             negative ^= true;
  2784.         }

  2786.         double result;
  2787.         if ((quadrant & 1) == 0) {
  2788.             result = tanQ(xa, xb, false);
  2789.         } else {
  2790.             result = -tanQ(xa, xb, true);
  2791.         }

  2793.         if (negative) {
  2794.             result = -result;
  2795.         }

  2797.         return result;
  2798.     }

  2800.     /**
  2801.      * Arctangent function
  2802.      *  @param x a number
  2803.      *  @return atan(x)
  2804.      */
  2805.     public static double atan(double x) {
  2806.         return atan(x, 0.0, false);
  2807.     }

  2809.     /** Internal helper function to compute arctangent.
  2810.      * @param xa number from which arctangent is requested
  2811.      * @param xb extra bits for x (may be 0.0)
  2812.      * @param leftPlane if true, result angle must be put in the left half plane
  2813.      * @return atan(xa + xb) (or angle shifted by {@code PI} if leftPlane is true)
  2814.      */
  2815.     private static double atan(double xa, double xb, boolean leftPlane) {
  2816.         if (xa == 0.0) { // Matches +/- 0.0; return correct sign
  2817.             return leftPlane ? copySign(Math.PI, xa) : xa;
  2818.         }

  2820.         final boolean negate;
  2821.         if (xa < 0) {
  2822.             // negative
  2823.             xa = -xa;
  2824.             xb = -xb;
  2825.             negate = true;
  2826.         } else {
  2827.             negate = false;
  2828.         }

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

  2834.         /* Estimate the closest tabulated arctan value, compute eps = xa-tangentTable */
  2835.         final int idx;
  2836.         if (xa < 1) {
  2837.             idx = (int) (((-1.7168146928204136 * xa * xa + 8.0) * xa) + 0.5);
  2838.         } else {
  2839.             final double oneOverXa = 1 / xa;
  2840.             idx = (int) (-((-1.7168146928204136 * oneOverXa * oneOverXa + 8.0) * oneOverXa) + 13.07);
  2841.         }

  2843.         final double ttA = TANGENT_TABLE_A[idx];
  2844.         final double ttB = TANGENT_TABLE_B[idx];

  2846.         double epsA = xa - ttA;
  2847.         double epsB = -(epsA - xa + ttA);
  2848.         epsB += xb - ttB;

  2850.         double temp = epsA + epsB;
  2851.         epsB = -(temp - epsA - epsB);
  2852.         epsA = temp;

  2854.         /* Compute eps = eps / (1.0 + xa*tangent) */
  2855.         temp = xa * HEX_40000000;
  2856.         double ya = xa + temp - temp;
  2857.         double yb = xb + xa - ya;
  2858.         xa = ya;
  2859.         xb += yb;

  2861.         //if (idx > 8 || idx == 0)
  2862.         if (idx == 0) {
  2863.             /* If the slope of the arctan is gentle enough (< 0.45), this approximation will suffice */
  2864.             //double denom = 1.0 / (1.0 + xa*tangentTableA[idx] + xb*tangentTableA[idx] + xa*tangentTableB[idx] + xb*tangentTableB[idx]);
  2865.             final double denom = 1d / (1d + (xa + xb) * (ttA + ttB));
  2866.             //double denom = 1.0 / (1.0 + xa*tangentTableA[idx]);
  2867.             ya = epsA * denom;
  2868.             yb = epsB * denom;
  2869.         } else {
  2870.             double temp2 = xa * ttA;
  2871.             double za = 1d + temp2;
  2872.             double zb = -(za - 1d - temp2);
  2873.             temp2 = xb * ttA + xa * ttB;
  2874.             temp = za + temp2;
  2875.             zb += -(temp - za - temp2);
  2876.             za = temp;

  2878.             zb += xb * ttB;
  2879.             ya = epsA / za;

  2881.             temp = ya * HEX_40000000;
  2882.             final double yaa = (ya + temp) - temp;
  2883.             final double yab = ya - yaa;

  2885.             temp = za * HEX_40000000;
  2886.             final double zaa = (za + temp) - temp;
  2887.             final double zab = za - zaa;

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

  2892.             yb += -epsA * zb / za / za;
  2893.             yb += epsB / za;
  2894.         }


  2897.         epsA = ya;
  2898.         epsB = yb;

  2900.         /* Evaluate polynomial */
  2901.         final double epsA2 = epsA * epsA;

  2903.         /*
  2904.     yb = -0.09001346640161823;
  2905.     yb = yb * epsA2 + 0.11110718400605211;
  2906.     yb = yb * epsA2 + -0.1428571349122913;
  2907.     yb = yb * epsA2 + 0.19999999999273194;
  2908.     yb = yb * epsA2 + -0.33333333333333093;
  2909.     yb = yb * epsA2 * epsA;
  2910.          */

  2912.         yb = 0.07490822288864472;
  2913.         yb = yb * epsA2 - 0.09088450866185192;
  2914.         yb = yb * epsA2 + 0.11111095942313305;
  2915.         yb = yb * epsA2 - 0.1428571423679182;
  2916.         yb = yb * epsA2 + 0.19999999999923582;
  2917.         yb = yb * epsA2 - 0.33333333333333287;
  2918.         yb = yb * epsA2 * epsA;


  2921.         ya = epsA;

  2923.         temp = ya + yb;
  2924.         yb = -(temp - ya - yb);
  2925.         ya = temp;

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

  2930.         final double eighths = EIGHTHS[idx];

  2932.         //result = yb + eighths[idx] + ya;
  2933.         double za = eighths + ya;
  2934.         double zb = -(za - eighths - ya);
  2935.         temp = za + yb;
  2936.         zb += -(temp - za - yb);
  2937.         za = temp;

  2939.         double result = za + zb;

  2941.         if (leftPlane) {
  2942.             // Result is in the left plane
  2943.             final double resultb = -(result - za - zb);
  2944.             final double pia = 1.5707963267948966 * 2;
  2945.             final double pib = 6.123233995736766E-17 * 2;

  2947.             za = pia - result;
  2948.             zb = -(za - pia + result);
  2949.             zb += pib - resultb;

  2951.             result = za + zb;
  2952.         }


  2955.         if (negate ^ leftPlane) {
  2956.             result = -result;
  2957.         }

  2959.         return result;
  2960.     }

  2962.     /**
  2963.      * Two arguments arctangent function
  2964.      * @param y ordinate
  2965.      * @param x abscissa
  2966.      * @return phase angle of point (x,y) between {@code -PI} and {@code PI}
  2967.      */
  2968.     public static double atan2(double y, double x) {
  2969.         if (Double.isNaN(x) || Double.isNaN(y)) {
  2970.             return Double.NaN;
  2971.         }

  2973.         if (y == 0) {
  2974.             final double result = x * y;
  2975.             final double invx = 1d / x;

  2977.             if (invx == 0) { // X is infinite
  2978.                 if (x > 0) {
  2979.                     return y; // return +/- 0.0
  2980.                 } else {
  2981.                     return copySign(Math.PI, y);
  2982.                 }
  2983.             }

  2985.             if (x < 0 || invx < 0) {
  2986.                 return copySign(Math.PI, y);
  2987.             } else {
  2988.                 return result;
  2989.             }
  2990.         }

  2992.         // y cannot now be zero

  2994.         if (y == Double.POSITIVE_INFINITY) {
  2995.             if (x == Double.POSITIVE_INFINITY) {
  2996.                 return Math.PI * F_1_4;
  2997.             }

  2999.             if (x == Double.NEGATIVE_INFINITY) {
  3000.                 return Math.PI * F_3_4;
  3001.             }

  3003.             return Math.PI * F_1_2;
  3004.         }

  3006.         if (y == Double.NEGATIVE_INFINITY) {
  3007.             if (x == Double.POSITIVE_INFINITY) {
  3008.                 return -Math.PI * F_1_4;
  3009.             }

  3011.             if (x == Double.NEGATIVE_INFINITY) {
  3012.                 return -Math.PI * F_3_4;
  3013.             }

  3015.             return -Math.PI * F_1_2;
  3016.         }

  3018.         if (x == Double.POSITIVE_INFINITY) {
  3019.             return copySign(0d, y);
  3020.         }

  3022.         if (x == Double.NEGATIVE_INFINITY)
  3023.         {
  3024.             return copySign(Math.PI, y);
  3025.         }

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

  3029.         if (x == 0) {
  3030.             return copySign(Math.PI * F_1_2, y);
  3031.         }

  3033.         // Compute ratio r = y/x
  3034.         final double r = y / x;
  3035.         if (Double.isInfinite(r)) { // bypass calculations that can create NaN
  3036.             return atan(r, 0, x < 0);
  3037.         }

  3039.         double ra = doubleHighPart(r);
  3040.         double rb = r - ra;

  3042.         // Split x
  3043.         final double xa = doubleHighPart(x);
  3044.         final double xb = x - xa;

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

  3048.         final double temp = ra + rb;
  3049.         rb = -(temp - ra - rb);
  3050.         ra = temp;

  3052.         if (ra == 0) { // Fix up the sign so atan works correctly
  3053.             ra = copySign(0d, y);
  3054.         }

  3056.         // Call atan
  3057.         return atan(ra, rb, x < 0);

  3059.     }

  3061.     /** Compute the arc sine of a number.
  3062.      * @param x number on which evaluation is done
  3063.      * @return arc sine of x
  3064.      */
  3065.     public static double asin(double x) {
  3066.       if (Double.isNaN(x)) {
  3067.           return Double.NaN;
  3068.       }

  3070.       if (x > 1.0 || x < -1.0) {
  3071.           return Double.NaN;
  3072.       }

  3074.       if (x == 1.0) {
  3075.           return Math.PI/2.0;
  3076.       }

  3078.       if (x == -1.0) {
  3079.           return -Math.PI/2.0;
  3080.       }

  3082.       if (x == 0.0) { // Matches +/- 0.0; return correct sign
  3083.           return x;
  3084.       }

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

  3088.       /* Split x */
  3089.       double temp = x * HEX_40000000;
  3090.       final double xa = x + temp - temp;
  3091.       final double xb = x - xa;

  3093.       /* Square it */
  3094.       double ya = xa*xa;
  3095.       double yb = xa*xb*2.0 + xb*xb;

  3097.       /* Subtract from 1 */
  3098.       ya = -ya;
  3099.       yb = -yb;

  3101.       double za = 1.0 + ya;
  3102.       double zb = -(za - 1.0 - ya);

  3104.       temp = za + yb;
  3105.       zb += -(temp - za - yb);
  3106.       za = temp;

  3108.       /* Square root */
  3109.       double y;
  3110.       y = sqrt(za);
  3111.       temp = y * HEX_40000000;
  3112.       ya = y + temp - temp;
  3113.       yb = y - ya;

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

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

  3121.       // Compute ratio r = x/y
  3122.       double r = x/y;
  3123.       temp = r * HEX_40000000;
  3124.       double ra = r + temp - temp;
  3125.       double rb = r - ra;

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

  3130.       temp = ra + rb;
  3131.       rb = -(temp - ra - rb);
  3132.       ra = temp;

  3134.       return atan(ra, rb, false);
  3135.     }

  3137.     /** Compute the arc cosine of a number.
  3138.      * @param x number on which evaluation is done
  3139.      * @return arc cosine of x
  3140.      */
  3141.     public static double acos(double x) {
  3142.       if (Double.isNaN(x)) {
  3143.           return Double.NaN;
  3144.       }

  3146.       if (x > 1.0 || x < -1.0) {
  3147.           return Double.NaN;
  3148.       }

  3150.       if (x == -1.0) {
  3151.           return Math.PI;
  3152.       }

  3154.       if (x == 1.0) {
  3155.           return 0.0;
  3156.       }

  3158.       if (x == 0) {
  3159.           return Math.PI/2.0;
  3160.       }

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

  3164.       /* Split x */
  3165.       double temp = x * HEX_40000000;
  3166.       final double xa = x + temp - temp;
  3167.       final double xb = x - xa;

  3169.       /* Square it */
  3170.       double ya = xa*xa;
  3171.       double yb = xa*xb*2.0 + xb*xb;

  3173.       /* Subtract from 1 */
  3174.       ya = -ya;
  3175.       yb = -yb;

  3177.       double za = 1.0 + ya;
  3178.       double zb = -(za - 1.0 - ya);

  3180.       temp = za + yb;
  3181.       zb += -(temp - za - yb);
  3182.       za = temp;

  3184.       /* Square root */
  3185.       double y = sqrt(za);
  3186.       temp = y * HEX_40000000;
  3187.       ya = y + temp - temp;
  3188.       yb = y - ya;

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

  3193.       /* Contribution of zb to sqrt */
  3194.       yb += zb / (2.0*y);
  3195.       y = ya+yb;
  3196.       yb = -(y - ya - yb);

  3198.       // Compute ratio r = y/x
  3199.       double r = y/x;

  3201.       // Did r overflow?
  3202.       if (Double.isInfinite(r)) { // x is effectively zero
  3203.           return Math.PI/2; // so return the appropriate value
  3204.       }

  3206.       double ra = doubleHighPart(r);
  3207.       double rb = r - ra;

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

  3212.       temp = ra + rb;
  3213.       rb = -(temp - ra - rb);
  3214.       ra = temp;

  3216.       return atan(ra, rb, x<0);
  3217.     }

  3219.     /** Compute the cubic root of a number.
  3220.      * @param x number on which evaluation is done
  3221.      * @return cubic root of x
  3222.      */
  3223.     public static double cbrt(double x) {
  3224.       /* Convert input double to bits */
  3225.       long inbits = Double.doubleToRawLongBits(x);
  3226.       int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
  3227.       boolean subnormal = false;

  3229.       if (exponent == -1023) {
  3230.           if (x == 0) {
  3231.               return x;
  3232.           }

  3234.           /* Subnormal, so normalize */
  3235.           subnormal = true;
  3236.           x *= 1.8014398509481984E16;  // 2^54
  3237.           inbits = Double.doubleToRawLongBits(x);
  3238.           exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
  3239.       }

  3241.       if (exponent == 1024) {
  3242.           // Nan or infinity.  Don't care which.
  3243.           return x;
  3244.       }

  3246.       /* Divide the exponent by 3 */
  3247.       int exp3 = exponent / 3;

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

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

  3256.       /* Estimate the cube root of mant by polynomial */
  3257.       double est = -0.010714690733195933;
  3258.       est = est * mant + 0.0875862700108075;
  3259.       est = est * mant + -0.3058015757857271;
  3260.       est = est * mant + 0.7249995199969751;
  3261.       est = est * mant + 0.5039018405998233;

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

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

  3272.       // Do one round of Newton's method in extended precision to get the last bit right.
  3273.       double temp = est * HEX_40000000;
  3274.       double ya = est + temp - temp;
  3275.       double yb = est - ya;

  3277.       double za = ya * ya;
  3278.       double zb = ya * yb * 2.0 + yb * yb;
  3279.       temp = za * HEX_40000000;
  3280.       double temp2 = za + temp - temp;
  3281.       zb += za - temp2;
  3282.       za = temp2;

  3284.       zb = za * yb + ya * zb + zb * yb;
  3285.       za *= ya;

  3287.       double na = xs - za;
  3288.       double nb = -(na - xs + za);
  3289.       nb -= zb;

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

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

  3296.       if (subnormal) {
  3297.           est *= 3.814697265625E-6;  // 2^-18
  3298.       }

  3300.       return est;
  3301.     }

  3303.     /**
  3304.      *  Convert degrees to radians, with error of less than 0.5 ULP
  3305.      *  @param x angle in degrees
  3306.      *  @return x converted into radians
  3307.      */
  3308.     public static double toRadians(double x)
  3309.     {
  3310.         if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
  3311.             return x;
  3312.         }

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

  3318.         double xa = doubleHighPart(x);
  3319.         double xb = x - xa;

  3321.         double result = xb * factb + xb * facta + xa * factb + xa * facta;
  3322.         if (result == 0) {
  3323.             result *= x; // ensure correct sign if calculation underflows
  3324.         }
  3325.         return result;
  3326.     }

  3328.     /**
  3329.      *  Convert radians to degrees, with error of less than 0.5 ULP
  3330.      *  @param x angle in radians
  3331.      *  @return x converted into degrees
  3332.      */
  3333.     public static double toDegrees(double x)
  3334.     {
  3335.         if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
  3336.             return x;
  3337.         }

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

  3343.         double xa = doubleHighPart(x);
  3344.         double xb = x - xa;

  3346.         return xb * factb + xb * facta + xa * factb + xa * facta;
  3347.     }

  3349.     /**
  3350.      * Absolute value.
  3351.      * @param x number from which absolute value is requested
  3352.      * @return abs(x)
  3353.      */
  3354.     public static int abs(final int x) {
  3355.         final int i = x >>> 31;
  3356.         return (x ^ (~i + 1)) + i;
  3357.     }

  3359.     /**
  3360.      * Absolute value.
  3361.      * @param x number from which absolute value is requested
  3362.      * @return abs(x)
  3363.      */
  3364.     public static long abs(final long x) {
  3365.         final long l = x >>> 63;
  3366.         // l is one if x negative zero else
  3367.         // ~l+1 is zero if x is positive, -1 if x is negative
  3368.         // x^(~l+1) is x is x is positive, ~x if x is negative
  3369.         // add around
  3370.         return (x ^ (~l + 1)) + l;
  3371.     }

  3373.     /**
  3374.      * Absolute value.
  3375.      * @param x number from which absolute value is requested
  3376.      * @return abs(x), or throws an exception for {@code Integer.MIN_VALUE}
  3377.      * @since 2.0
  3378.      */
  3379.     public static int absExact(final int x) {
  3380.         if (x == Integer.MIN_VALUE) {
  3381.             throw new ArithmeticException();
  3382.         }
  3383.         return abs(x);
  3384.     }

  3386.     /**
  3387.      * Absolute value.
  3388.      * @param x number from which absolute value is requested
  3389.      * @return abs(x), or throws an exception for {@code Long.MIN_VALUE}
  3390.      * @since 2.0
  3391.      */
  3392.     public static long absExact(final long x) {
  3393.         if (x == Long.MIN_VALUE) {
  3394.             throw new ArithmeticException();
  3395.         }
  3396.         return abs(x);
  3397.     }

  3399.     /**
  3400.      * Absolute value.
  3401.      * @param x number from which absolute value is requested
  3402.      * @return abs(x)
  3403.      * @since 2.0
  3404.      */
  3405.     public static float abs(final float x) {
  3406.         return Float.intBitsToFloat(MASK_NON_SIGN_INT & Float.floatToRawIntBits(x));
  3407.     }

  3409.     /**
  3410.      * Absolute value.
  3411.      * @param x number from which absolute value is requested
  3412.      * @return abs(x)
  3413.      */
  3414.     public static double abs(double x) {
  3415.         return Double.longBitsToDouble(MASK_NON_SIGN_LONG & Double.doubleToRawLongBits(x));
  3416.     }

  3418.     /**
  3419.      * Negates the argument.
  3420.      * @param x number from which opposite value is requested
  3421.      * @return -x, or throws an exception for {@code Integer.MIN_VALUE}
  3422.      * @since 2.0
  3423.      */
  3424.     public static int negateExact(final int x) {
  3425.         if (x == Integer.MIN_VALUE) {
  3426.             throw new ArithmeticException();
  3427.         }
  3428.         return -x;
  3429.     }

  3431.     /**
  3432.      * Negates the argument.
  3433.      * @param x number from which opposite value is requested
  3434.      * @return -x, or throws an exception for {@code Long.MIN_VALUE}
  3435.      * @since 2.0
  3436.      */
  3437.     public static long negateExact(final long x) {
  3438.         if (x == Long.MIN_VALUE) {
  3439.             throw new ArithmeticException();
  3440.         }
  3441.         return -x;
  3442.     }

  3444.     /**
  3445.      * Compute least significant bit (Unit in Last Position) for a number.
  3446.      * @param x number from which ulp is requested
  3447.      * @return ulp(x)
  3448.      */
  3449.     public static double ulp(double x) {
  3450.         if (Double.isInfinite(x)) {
  3451.             return Double.POSITIVE_INFINITY;
  3452.         }
  3453.         return abs(x - Double.longBitsToDouble(Double.doubleToRawLongBits(x) ^ 1));
  3454.     }

  3456.     /**
  3457.      * Compute least significant bit (Unit in Last Position) for a number.
  3458.      * @param x number from which ulp is requested
  3459.      * @return ulp(x)
  3460.      */
  3461.     public static float ulp(float x) {
  3462.         if (Float.isInfinite(x)) {
  3463.             return Float.POSITIVE_INFINITY;
  3464.         }
  3465.         return abs(x - Float.intBitsToFloat(Float.floatToRawIntBits(x) ^ 1));
  3466.     }

  3468.     /**
  3469.      * Multiply a double number by a power of 2.
  3470.      * @param d number to multiply
  3471.      * @param n power of 2
  3472.      * @return d &times; 2<sup>n</sup>
  3473.      */
  3474.     public static double scalb(final double d, final int n) {

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

  3481.         // handle special cases
  3482.         if (Double.isNaN(d) || Double.isInfinite(d) || (d == 0)) {
  3483.             return d;
  3484.         }
  3485.         if (n < -2098) {
  3486.             return (d > 0) ? 0.0 : -0.0;
  3487.         }
  3488.         if (n > 2097) {
  3489.             return (d > 0) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
  3490.         }

  3492.         // decompose d
  3493.         final long bits = Double.doubleToRawLongBits(d);
  3494.         final long sign = bits & 0x8000000000000000L;
  3495.         int  exponent   = ((int) (bits >>> 52)) & 0x7ff;
  3496.         long mantissa   = bits & 0x000fffffffffffffL;

  3498.         // compute scaled exponent
  3499.         int scaledExponent = exponent + n;

  3501.         if (n < 0) {
  3502.             // we are really in the case n <= -1023
  3503.             if (scaledExponent > 0) {
  3504.                 // both the input and the result are normal numbers, we only adjust the exponent
  3505.                 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
  3506.             } else if (scaledExponent > -53) {
  3507.                 // the input is a normal number and the result is a subnormal number

  3509.                 // recover the hidden mantissa bit
  3510.                 mantissa |= 1L << 52;

  3512.                 // scales down complete mantissa, hence losing least significant bits
  3513.                 final long mostSignificantLostBit = mantissa & (1L << (-scaledExponent));
  3514.                 mantissa >>>= 1 - scaledExponent;
  3515.                 if (mostSignificantLostBit != 0) {
  3516.                     // we need to add 1 bit to round up the result
  3517.                     mantissa++;
  3518.                 }
  3519.                 return Double.longBitsToDouble(sign | mantissa);

  3521.             } else {
  3522.                 // no need to compute the mantissa, the number scales down to 0
  3523.                 return (sign == 0L) ? 0.0 : -0.0;
  3524.             }
  3525.         } else {
  3526.             // we are really in the case n >= 1024
  3527.             if (exponent == 0) {

  3529.                 // the input number is subnormal, normalize it
  3530.                 while ((mantissa >>> 52) != 1) {
  3531.                     mantissa <<= 1;
  3532.                     --scaledExponent;
  3533.                 }
  3534.                 ++scaledExponent;
  3535.                 mantissa &= 0x000fffffffffffffL;

  3537.                 if (scaledExponent < 2047) {
  3538.                     return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
  3539.                 } else {
  3540.                     return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
  3541.                 }

  3543.             } else if (scaledExponent < 2047) {
  3544.                 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
  3545.             } else {
  3546.                 return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
  3547.             }
  3548.         }

  3550.     }

  3552.     /**
  3553.      * Multiply a float number by a power of 2.
  3554.      * @param f number to multiply
  3555.      * @param n power of 2
  3556.      * @return f &times; 2<sup>n</sup>
  3557.      */
  3558.     public static float scalb(final float f, final int n) {

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

  3565.         // handle special cases
  3566.         if (Float.isNaN(f) || Float.isInfinite(f) || (f == 0f)) {
  3567.             return f;
  3568.         }
  3569.         if (n < -277) {
  3570.             return (f > 0) ? 0.0f : -0.0f;
  3571.         }
  3572.         if (n > 276) {
  3573.             return (f > 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
  3574.         }

  3576.         // decompose f
  3577.         final int bits = Float.floatToIntBits(f);
  3578.         final int sign = bits & 0x80000000;
  3579.         int  exponent  = (bits >>> 23) & 0xff;
  3580.         int mantissa   = bits & 0x007fffff;

  3582.         // compute scaled exponent
  3583.         int scaledExponent = exponent + n;

  3585.         if (n < 0) {
  3586.             // we are really in the case n <= -127
  3587.             if (scaledExponent > 0) {
  3588.                 // both the input and the result are normal numbers, we only adjust the exponent
  3589.                 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
  3590.             } else if (scaledExponent > -24) {
  3591.                 // the input is a normal number and the result is a subnormal number

  3593.                 // recover the hidden mantissa bit
  3594.                 mantissa |= 1 << 23;

  3596.                 // scales down complete mantissa, hence losing least significant bits
  3597.                 final int mostSignificantLostBit = mantissa & (1 << (-scaledExponent));
  3598.                 mantissa >>>= 1 - scaledExponent;
  3599.                 if (mostSignificantLostBit != 0) {
  3600.                     // we need to add 1 bit to round up the result
  3601.                     mantissa++;
  3602.                 }
  3603.                 return Float.intBitsToFloat(sign | mantissa);

  3605.             } else {
  3606.                 // no need to compute the mantissa, the number scales down to 0
  3607.                 return (sign == 0) ? 0.0f : -0.0f;
  3608.             }
  3609.         } else {
  3610.             // we are really in the case n >= 128
  3611.             if (exponent == 0) {

  3613.                 // the input number is subnormal, normalize it
  3614.                 while ((mantissa >>> 23) != 1) {
  3615.                     mantissa <<= 1;
  3616.                     --scaledExponent;
  3617.                 }
  3618.                 ++scaledExponent;
  3619.                 mantissa &= 0x007fffff;

  3621.                 if (scaledExponent < 255) {
  3622.                     return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
  3623.                 } else {
  3624.                     return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
  3625.                 }

  3627.             } else if (scaledExponent < 255) {
  3628.                 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
  3629.             } else {
  3630.                 return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
  3631.             }
  3632.         }

  3634.     }

  3636.     /**
  3637.      * Get the next machine representable number after a number, moving
  3638.      * in the direction of another number.
  3639.      * <p>
  3640.      * The ordering is as follows (increasing):
  3641.      * </p>
  3642.      * <ul>
  3643.      * <li>-INFINITY</li>
  3644.      * <li>-MAX_VALUE</li>
  3645.      * <li>-MIN_VALUE</li>
  3646.      * <li>-0.0</li>
  3647.      * <li>+0.0</li>
  3648.      * <li>+MIN_VALUE</li>
  3649.      * <li>+MAX_VALUE</li>
  3650.      * <li>+INFINITY</li>
  3651.      * </ul>
  3652.      * <p>
  3653.      * If arguments compare equal, then the second argument is returned.
  3654.      * </p>
  3655.      * <p>
  3656.      * If {@code direction} is greater than {@code d},
  3657.      * the smallest machine representable number strictly greater than
  3658.      * {@code d} is returned; if less, then the largest representable number
  3659.      * strictly less than {@code d} is returned.
  3660.      * </p>
  3661.      * <p>
  3662.      * If {@code d} is infinite and direction does not
  3663.      * bring it back to finite numbers, it is returned unchanged.
  3664.      * </p>
  3665.      *
  3666.      * @param d base number
  3667.      * @param direction (the only important thing is whether
  3668.      * {@code direction} is greater or smaller than {@code d})
  3669.      * @return the next machine representable number in the specified direction
  3670.      */
  3671.     public static double nextAfter(double d, double direction) {

  3673.         // handling of some important special cases
  3674.         if (Double.isNaN(d) || Double.isNaN(direction)) {
  3675.             return Double.NaN;
  3676.         } else if (d == direction) {
  3677.             return direction;
  3678.         } else if (Double.isInfinite(d)) {
  3679.             return (d < 0) ? -Double.MAX_VALUE : Double.MAX_VALUE;
  3680.         } else if (d == 0) {
  3681.             return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
  3682.         }
  3683.         // special cases MAX_VALUE to infinity and  MIN_VALUE to 0
  3684.         // are handled just as normal numbers
  3685.         // can use raw bits since already dealt with infinity and NaN
  3686.         final long bits = Double.doubleToRawLongBits(d);
  3687.         final long sign = bits & 0x8000000000000000L;
  3688.         if ((direction < d) ^ (sign == 0L)) {
  3689.             return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) + 1));
  3690.         } else {
  3691.             return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) - 1));
  3692.         }

  3694.     }

  3696.     /**
  3697.      * Get the next machine representable number after a number, moving
  3698.      * in the direction of another number.
  3699.      * <p>* The ordering is as follows (increasing):</p>
  3700.      * <ul>
  3701.      * <li>-INFINITY</li>
  3702.      * <li>-MAX_VALUE</li>
  3703.      * <li>-MIN_VALUE</li>
  3704.      * <li>-0.0</li>
  3705.      * <li>+0.0</li>
  3706.      * <li>+MIN_VALUE</li>
  3707.      * <li>+MAX_VALUE</li>
  3708.      * <li>+INFINITY</li>
  3709.      * </ul>
  3710.      * <p>
  3711.      * If arguments compare equal, then the second argument is returned.
  3712.      * </p>
  3713.      * <p>
  3714.      * If {@code direction} is greater than {@code f},
  3715.      * the smallest machine representable number strictly greater than
  3716.      * {@code f} is returned; if less, then the largest representable number
  3717.      * strictly less than {@code f} is returned.
  3718.      * </p>
  3719.      * <p>
  3720.      * If {@code f} is infinite and direction does not
  3721.      * bring it back to finite numbers, it is returned unchanged.
  3722.      * </p>
  3723.      *
  3724.      * @param f base number
  3725.      * @param direction (the only important thing is whether
  3726.      * {@code direction} is greater or smaller than {@code f})
  3727.      * @return the next machine representable number in the specified direction
  3728.      */
  3729.     public static float nextAfter(final float f, final double direction) {

  3731.         // handling of some important special cases
  3732.         if (Double.isNaN(f) || Double.isNaN(direction)) {
  3733.             return Float.NaN;
  3734.         } else if (f == direction) {
  3735.             return (float) direction;
  3736.         } else if (Float.isInfinite(f)) {
  3737.             return (f < 0f) ? -Float.MAX_VALUE : Float.MAX_VALUE;
  3738.         } else if (f == 0f) {
  3739.             return (direction < 0) ? -Float.MIN_VALUE : Float.MIN_VALUE;
  3740.         }
  3741.         // special cases MAX_VALUE to infinity and  MIN_VALUE to 0
  3742.         // are handled just as normal numbers

  3744.         final int bits = Float.floatToIntBits(f);
  3745.         final int sign = bits & 0x80000000;
  3746.         if ((direction < f) ^ (sign == 0)) {
  3747.             return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) + 1));
  3748.         } else {
  3749.             return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) - 1));
  3750.         }

  3752.     }

  3754.     /** Get the largest whole number smaller than x.
  3755.      * @param x number from which floor is requested
  3756.      * @return a double number f such that f is an integer f &lt;= x &lt; f + 1.0
  3757.      */
  3758.     public static double floor(double x) {
  3759.         long y;

  3761.         if (Double.isNaN(x)) {
  3762.             return x;
  3763.         }

  3765.         if (x >= TWO_POWER_52 || x <= -TWO_POWER_52) {
  3766.             return x;
  3767.         }

  3769.         y = (long) x;
  3770.         if (x < 0 && y != x) {
  3771.             y--;
  3772.         }

  3774.         if (y == 0) {
  3775.             return x*y;
  3776.         }

  3778.         return y;
  3779.     }

  3781.     /** Get the smallest whole number larger than x.
  3782.      * @param x number from which ceil is requested
  3783.      * @return a double number c such that c is an integer c - 1.0 &lt; x &lt;= c
  3784.      */
  3785.     public static double ceil(double x) {
  3786.         double y;

  3788.         if (Double.isNaN(x)) {
  3789.             return x;
  3790.         }

  3792.         y = floor(x);
  3793.         if (y == x) {
  3794.             return y;
  3795.         }

  3797.         y += 1.0;

  3799.         if (y == 0) {
  3800.             return x*y;
  3801.         }

  3803.         return y;
  3804.     }

  3806.     /** Get the whole number that is the nearest to x, or the even one if x is exactly half way between two integers.
  3807.      * @param x number from which nearest whole number is requested
  3808.      * @return a double number r such that r is an integer r - 0.5 &lt;= x &lt;= r + 0.5
  3809.      */
  3810.     public static double rint(double x) {
  3811.         double y = floor(x);
  3812.         double d = x - y;

  3814.         if (d > 0.5) {
  3815.             if (y == -1.0) {
  3816.                 return -0.0; // Preserve sign of operand
  3817.             }
  3818.             return y+1.0;
  3819.         }
  3820.         if (d < 0.5) {
  3821.             return y;
  3822.         }

  3824.         /* half way, round to even */
  3825.         long z = (long) y;
  3826.         return (z & 1) == 0 ? y : y + 1.0;
  3827.     }

  3829.     /** Get the closest long to x.
  3830.      * @param x number from which closest long is requested
  3831.      * @return closest long to x
  3832.      */
  3833.     public static long round(double x) {
  3834.         final long bits = Double.doubleToRawLongBits(x);
  3835.         final int biasedExp = ((int)(bits>>52)) & 0x7ff;
  3836.         // Shift to get rid of bits past comma except first one: will need to
  3837.         // 1-shift to the right to end up with correct magnitude.
  3838.         final int shift = (52 - 1 + Double.MAX_EXPONENT) - biasedExp;
  3839.         if ((shift & -64) == 0) {
  3840.             // shift in [0,63], so unbiased exp in [-12,51].
  3841.             long extendedMantissa = 0x0010000000000000L | (bits & 0x000fffffffffffffL);
  3842.             if (bits < 0) {
  3843.                 extendedMantissa = -extendedMantissa;
  3844.             }
  3845.             // If value is positive and first bit past comma is 0, rounding
  3846.             // to lower integer, else to upper one, which is what "+1" and
  3847.             // then ">>1" do.
  3848.             return ((extendedMantissa >> shift) + 1L) >> 1;
  3849.         } else {
  3850.             // +-Infinity, NaN, or a mathematical integer.
  3851.             return (long) x;
  3852.         }
  3853.     }

  3855.     /** Get the closest int to x.
  3856.      * @param x number from which closest int is requested
  3857.      * @return closest int to x
  3858.      */
  3859.     public static int round(final float x) {
  3860.         final int bits = Float.floatToRawIntBits(x);
  3861.         final int biasedExp = (bits>>23) & 0xff;
  3862.         // Shift to get rid of bits past comma except first one: will need to
  3863.         // 1-shift to the right to end up with correct magnitude.
  3864.         final int shift = (23 - 1 + Float.MAX_EXPONENT) - biasedExp;
  3865.         if ((shift & -32) == 0) {
  3866.             // shift in [0,31], so unbiased exp in [-9,22].
  3867.             int extendedMantissa = 0x00800000 | (bits & 0x007fffff);
  3868.             if (bits < 0) {
  3869.                 extendedMantissa = -extendedMantissa;
  3870.             }
  3871.             // If value is positive and first bit past comma is 0, rounding
  3872.             // to lower integer, else to upper one, which is what "+1" and
  3873.             // then ">>1" do.
  3874.             return ((extendedMantissa >> shift) + 1) >> 1;
  3875.         } else {
  3876.             // +-Infinity, NaN, or a mathematical integer.
  3877.             return (int) x;
  3878.         }
  3879.     }

  3881.     /** Compute the minimum of two values
  3882.      * @param a first value
  3883.      * @param b second value
  3884.      * @return a if a is lesser or equal to b, b otherwise
  3885.      */
  3886.     public static int min(final int a, final int b) {
  3887.         return (a <= b) ? a : b;
  3888.     }

  3890.     /** Compute the minimum of two values
  3891.      * @param a first value
  3892.      * @param b second value
  3893.      * @return a if a is lesser or equal to b, b otherwise
  3894.      */
  3895.     public static long min(final long a, final long b) {
  3896.         return (a <= b) ? a : b;
  3897.     }

  3899.     /** Compute the minimum of two values
  3900.      * @param a first value
  3901.      * @param b second value
  3902.      * @return a if a is lesser or equal to b, b otherwise
  3903.      */
  3904.     public static float min(final float a, final float b) {
  3905.         if (a > b) {
  3906.             return b;
  3907.         }
  3908.         if (a < b) {
  3909.             return a;
  3910.         }
  3911.         /* if either arg is NaN, return NaN */
  3912.         if (a != b) {
  3913.             return Float.NaN;
  3914.         }
  3915.         /* min(+0.0,-0.0) == -0.0 */
  3916.         /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
  3917.         int bits = Float.floatToRawIntBits(a);
  3918.         if (bits == 0x80000000) {
  3919.             return a;
  3920.         }
  3921.         return b;
  3922.     }

  3924.     /** Compute the minimum of two values
  3925.      * @param a first value
  3926.      * @param b second value
  3927.      * @return a if a is lesser or equal to b, b otherwise
  3928.      */
  3929.     public static double min(final double a, final double b) {
  3930.         if (a > b) {
  3931.             return b;
  3932.         }
  3933.         if (a < b) {
  3934.             return a;
  3935.         }
  3936.         /* if either arg is NaN, return NaN */
  3937.         if (a != b) {
  3938.             return Double.NaN;
  3939.         }
  3940.         /* min(+0.0,-0.0) == -0.0 */
  3941.         /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
  3942.         long bits = Double.doubleToRawLongBits(a);
  3943.         if (bits == 0x8000000000000000L) {
  3944.             return a;
  3945.         }
  3946.         return b;
  3947.     }

  3949.     /** Compute the maximum of two values
  3950.      * @param a first value
  3951.      * @param b second value
  3952.      * @return b if a is lesser or equal to b, a otherwise
  3953.      */
  3954.     public static int max(final int a, final int b) {
  3955.         return (a <= b) ? b : a;
  3956.     }

  3958.     /** Compute the maximum of two values
  3959.      * @param a first value
  3960.      * @param b second value
  3961.      * @return b if a is lesser or equal to b, a otherwise
  3962.      */
  3963.     public static long max(final long a, final long b) {
  3964.         return (a <= b) ? b : a;
  3965.     }

  3967.     /** Compute the maximum of two values
  3968.      * @param a first value
  3969.      * @param b second value
  3970.      * @return b if a is lesser or equal to b, a otherwise
  3971.      */
  3972.     public static float max(final float a, final float b) {
  3973.         if (a > b) {
  3974.             return a;
  3975.         }
  3976.         if (a < b) {
  3977.             return b;
  3978.         }
  3979.         /* if either arg is NaN, return NaN */
  3980.         if (a != b) {
  3981.             return Float.NaN;
  3982.         }
  3983.         /* min(+0.0,-0.0) == -0.0 */
  3984.         /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
  3985.         int bits = Float.floatToRawIntBits(a);
  3986.         if (bits == 0x80000000) {
  3987.             return b;
  3988.         }
  3989.         return a;
  3990.     }

  3992.     /** Compute the maximum of two values
  3993.      * @param a first value
  3994.      * @param b second value
  3995.      * @return b if a is lesser or equal to b, a otherwise
  3996.      */
  3997.     public static double max(final double a, final double b) {
  3998.         if (a > b) {
  3999.             return a;
  4000.         }
  4001.         if (a < b) {
  4002.             return b;
  4003.         }
  4004.         /* if either arg is NaN, return NaN */
  4005.         if (a != b) {
  4006.             return Double.NaN;
  4007.         }
  4008.         /* min(+0.0,-0.0) == -0.0 */
  4009.         /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
  4010.         long bits = Double.doubleToRawLongBits(a);
  4011.         if (bits == 0x8000000000000000L) {
  4012.             return b;
  4013.         }
  4014.         return a;
  4015.     }

  4017.     /**
  4018.      * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
  4019.      * - sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)<br>
  4020.      * avoiding intermediate overflow or underflow.
  4021.      *
  4022.      * <ul>
  4023.      * <li> If either argument is infinite, then the result is positive infinity.</li>
  4024.      * <li> else, if either argument is NaN then the result is NaN.</li>
  4025.      * </ul>
  4026.      *
  4027.      * @param x a value
  4028.      * @param y a value
  4029.      * @return sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
  4030.      */
  4031.     public static double hypot(final double x, final double y) {
  4032.         if (Double.isInfinite(x) || Double.isInfinite(y)) {
  4033.             return Double.POSITIVE_INFINITY;
  4034.         } else if (Double.isNaN(x) || Double.isNaN(y)) {
  4035.             return Double.NaN;
  4036.         } else {

  4038.             final int expX = getExponent(x);
  4039.             final int expY = getExponent(y);
  4040.             if (expX > expY + 27) {
  4041.                 // y is neglectible with respect to x
  4042.                 return abs(x);
  4043.             } else if (expY > expX + 27) {
  4044.                 // x is neglectible with respect to y
  4045.                 return abs(y);
  4046.             } else {

  4048.                 // find an intermediate scale to avoid both overflow and underflow
  4049.                 final int middleExp = (expX + expY) / 2;

  4051.                 // scale parameters without losing precision
  4052.                 final double scaledX = scalb(x, -middleExp);
  4053.                 final double scaledY = scalb(y, -middleExp);

  4055.                 // compute scaled hypotenuse
  4056.                 final double scaledH = sqrt(scaledX * scaledX + scaledY * scaledY);

  4058.                 // remove scaling
  4059.                 return scalb(scaledH, middleExp);

  4061.             }

  4063.         }
  4064.     }

  4066.     /**
  4067.      * Computes the remainder as prescribed by the IEEE 754 standard.
  4068.      * <p>
  4069.      * The remainder value is mathematically equal to {@code x - y*n}
  4070.      * where {@code n} is the mathematical integer closest to the exact mathematical value
  4071.      * of the quotient {@code x/y}.
  4072.      * If two mathematical integers are equally close to {@code x/y} then
  4073.      * {@code n} is the integer that is even.
  4074.      * </p>
  4075.      * <ul>
  4076.      * <li>If either operand is NaN, the result is NaN.</li>
  4077.      * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
  4078.      * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
  4079.      * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
  4080.      * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
  4081.      * </ul>
  4082.      * @param dividend the number to be divided
  4083.      * @param divisor the number by which to divide
  4084.      * @return the remainder, rounded
  4085.      */
  4086.     public static double IEEEremainder(final double dividend, final double divisor) {
  4087.         if (getExponent(dividend) == 1024 || getExponent(divisor) == 1024 || divisor == 0.0) {
  4088.             // we are in one of the special cases
  4089.             if (Double.isInfinite(divisor) && !Double.isInfinite(dividend)) {
  4090.                 return dividend;
  4091.             } else {
  4092.                 return Double.NaN;
  4093.             }
  4094.         } else {
  4095.             // we are in the general case
  4096.             final double n         = rint(dividend / divisor);
  4097.             final double remainder = Double.isInfinite(n) ? 0.0 : dividend - divisor * n;
  4098.             return (remainder == 0) ? copySign(remainder, dividend) : remainder;
  4099.         }
  4100.     }

  4102.     /** Convert a long to interger, detecting overflows
  4103.      * @param n number to convert to int
  4104.      * @return integer with same valie as n if no overflows occur
  4105.      * @exception MathRuntimeException if n cannot fit into an int
  4106.      */
  4107.     public static int toIntExact(final long n) throws MathRuntimeException {
  4108.         if (n < Integer.MIN_VALUE || n > Integer.MAX_VALUE) {
  4109.             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW);
  4110.         }
  4111.         return (int) n;
  4112.     }

  4114.     /** Increment a number, detecting overflows.
  4115.      * @param n number to increment
  4116.      * @return n+1 if no overflows occur
  4117.      * @exception MathRuntimeException if an overflow occurs
  4118.      */
  4119.     public static int incrementExact(final int n) throws MathRuntimeException {

  4121.         if (n == Integer.MAX_VALUE) {
  4122.             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_ADDITION, n, 1);
  4123.         }

  4125.         return n + 1;

  4127.     }

  4129.     /** Increment a number, detecting overflows.
  4130.      * @param n number to increment
  4131.      * @return n+1 if no overflows occur
  4132.      * @exception MathRuntimeException if an overflow occurs
  4133.      */
  4134.     public static long incrementExact(final long n) throws MathRuntimeException {

  4136.         if (n == Long.MAX_VALUE) {
  4137.             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_ADDITION, n, 1);
  4138.         }

  4140.         return n + 1;

  4142.     }

  4144.     /** Decrement a number, detecting overflows.
  4145.      * @param n number to decrement
  4146.      * @return n-1 if no overflows occur
  4147.      * @exception MathRuntimeException if an overflow occurs
  4148.      */
  4149.     public static int decrementExact(final int n) throws MathRuntimeException {

  4151.         if (n == Integer.MIN_VALUE) {
  4152.             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_SUBTRACTION, n, 1);
  4153.         }

  4155.         return n - 1;

  4157.     }

  4159.     /** Decrement a number, detecting overflows.
  4160.      * @param n number to decrement
  4161.      * @return n-1 if no overflows occur
  4162.      * @exception MathRuntimeException if an overflow occurs
  4163.      */
  4164.     public static long decrementExact(final long n) throws MathRuntimeException {

  4166.         if (n == Long.MIN_VALUE) {
  4167.             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_SUBTRACTION, n, 1);
  4168.         }

  4170.         return n - 1;

  4172.     }

  4174.     /** Add two numbers, detecting overflows.
  4175.      * @param a first number to add
  4176.      * @param b second number to add
  4177.      * @return a+b if no overflows occur
  4178.      * @exception MathRuntimeException if an overflow occurs
  4179.      */
  4180.     public static int addExact(final int a, final int b) throws MathRuntimeException {

  4182.         // compute sum
  4183.         final int sum = a + b;

  4185.         // check for overflow
  4186.         if ((a ^ b) >= 0 && (sum ^ b) < 0) {
  4187.             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_ADDITION, a, b);
  4188.         }

  4190.         return sum;

  4192.     }

  4194.     /** Add two numbers, detecting overflows.
  4195.      * @param a first number to add
  4196.      * @param b second number to add
  4197.      * @return a+b if no overflows occur
  4198.      * @exception MathRuntimeException if an overflow occurs
  4199.      */
  4200.     public static long addExact(final long a, final long b) throws MathRuntimeException {

  4202.         // compute sum
  4203.         final long sum = a + b;

  4205.         // check for overflow
  4206.         if ((a ^ b) >= 0 && (sum ^ b) < 0) {
  4207.             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_ADDITION, a, b);
  4208.         }

  4210.         return sum;

  4212.     }

  4214.     /** Subtract two numbers, detecting overflows.
  4215.      * @param a first number
  4216.      * @param b second number to subtract from a
  4217.      * @return a-b if no overflows occur
  4218.      * @exception MathRuntimeException if an overflow occurs
  4219.      */
  4220.     public static int subtractExact(final int a, final int b) {

  4222.         // compute subtraction
  4223.         final int sub = a - b;

  4225.         // check for overflow
  4226.         if ((a ^ b) < 0 && (sub ^ b) >= 0) {
  4227.             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_SUBTRACTION, a, b);
  4228.         }

  4230.         return sub;

  4232.     }

  4234.     /** Subtract two numbers, detecting overflows.
  4235.      * @param a first number
  4236.      * @param b second number to subtract from a
  4237.      * @return a-b if no overflows occur
  4238.      * @exception MathRuntimeException if an overflow occurs
  4239.      */
  4240.     public static long subtractExact(final long a, final long b) {

  4242.         // compute subtraction
  4243.         final long sub = a - b;

  4245.         // check for overflow
  4246.         if ((a ^ b) < 0 && (sub ^ b) >= 0) {
  4247.             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_SUBTRACTION, a, b);
  4248.         }

  4250.         return sub;

  4252.     }

  4254.     /** Multiply two numbers, detecting overflows.
  4255.      * @param a first number to multiply
  4256.      * @param b second number to multiply
  4257.      * @return a*b if no overflows occur
  4258.      * @exception MathRuntimeException if an overflow occurs
  4259.      */
  4260.     public static int multiplyExact(final int a, final int b) {
  4261.         if (((b  >  0)  && (a > Integer.MAX_VALUE / b || a < Integer.MIN_VALUE / b)) ||
  4262.             ((b  < -1)  && (a > Integer.MIN_VALUE / b || a < Integer.MAX_VALUE / b)) ||
  4263.             ((b == -1)  && (a == Integer.MIN_VALUE))) {
  4264.             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_MULTIPLICATION, a, b);
  4265.         }
  4266.         return a * b;
  4267.     }

  4269.     /** Multiply two numbers, detecting overflows.
  4270.      * @param a first number to multiply
  4271.      * @param b second number to multiply
  4272.      * @return a*b if no overflows occur
  4273.      * @exception MathRuntimeException if an overflow occurs
  4274.      * @since 1.3
  4275.      */
  4276.     public static long multiplyExact(final long a, final int b) {
  4277.         return multiplyExact(a, (long) b);
  4278.     }

  4280.     /** Multiply two numbers, detecting overflows.
  4281.      * @param a first number to multiply
  4282.      * @param b second number to multiply
  4283.      * @return a*b if no overflows occur
  4284.      * @exception MathRuntimeException if an overflow occurs
  4285.      */
  4286.     public static long multiplyExact(final long a, final long b) {
  4287.         if (((b  >  0l)  && (a > Long.MAX_VALUE / b || a < Long.MIN_VALUE / b)) ||
  4288.             ((b  < -1l)  && (a > Long.MIN_VALUE / b || a < Long.MAX_VALUE / b)) ||
  4289.             ((b == -1l)  && (a == Long.MIN_VALUE))) {
  4290.                 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_MULTIPLICATION, a, b);
  4291.             }
  4292.             return a * b;
  4293.     }

  4295.     /** Multiply two integers and give an exact result without overflow.
  4296.      * @param a first factor
  4297.      * @param b second factor
  4298.      * @return a * b exactly
  4299.      * @since 1.3
  4300.      */
  4301.     public static long multiplyFull(final int a, final int b) {
  4302.         return ((long) a) * ((long) b); // NOPMD - casts are intentional here
  4303.     }

  4305.     /** Multiply two long integers and give the 64 most significant bits of the result.
  4306.      * <p>
  4307.      * Beware that as Java primitive long are always considered to be signed, there are some
  4308.      * intermediate values {@code a} and {@code b} for which {@code a * b} exceeds {@code Long.MAX_VALUE}
  4309.      * but this method will still return 0l. This happens for example for {@code a = 2³¹} and
  4310.      * {@code b = 2³²} as {@code a * b = 2⁶³ = Long.MAX_VALUE + 1}, so it exceeds the max value
  4311.      * for a long, but still fits in 64 bits, so this method correctly returns 0l in this case,
  4312.      * but multiplication result would be considered negative (and in fact equal to {@code Long.MIN_VALUE}
  4313.      * </p>
  4314.      * @param a first factor
  4315.      * @param b second factor
  4316.      * @return a * b / 2<sup>64</sup>
  4317.      * @since 1.3
  4318.      */
  4319.     public static long multiplyHigh(final long a, final long b) {

  4321.         // all computations below are performed on unsigned numbers because we start
  4322.         // by using logical shifts (and not arithmetic shifts). We will therefore
  4323.         // need to take care of sign before returning
  4324.         // a negative long n between -2⁶³ and -1, interpreted as an unsigned long
  4325.         // corresponds to 2⁶⁴ + n (which is between 2⁶³ and 2⁶⁴-1)
  4326.         // so if this number is multiplied by p, what we really compute
  4327.         // is (2⁶⁴ + n) * p = 2⁶⁴ * p + n * p, therefore the part above 64 bits
  4328.         // will have an extra term p that we will need to remove
  4329.         final long tobeRemoved = ((a < 0) ? b : 0) + ((b < 0) ? a : 0);

  4331.         return unsignedMultiplyHigh(a, b) - tobeRemoved;

  4333.     }

  4335.     /** Multiply two long unsigned integers and give the 64 most significant bits of the unsigned result.
  4336.      * <p>
  4337.      * Beware that as Java primitive long are always considered to be signed, there are some
  4338.      * intermediate values {@code a} and {@code b} for which {@code a * b} exceeds {@code Long.MAX_VALUE}
  4339.      * but this method will still return 0l. This happens for example for {@code a = 2³¹} and
  4340.      * {@code b = 2³²} as {@code a * b = 2⁶³ = Long.MAX_VALUE + 1}, so it exceeds the max value
  4341.      * for a long, but still fits in 64 bits, so this method correctly returns 0l in this case,
  4342.      * but multiplication result would be considered negative (and in fact equal to {@code Long.MIN_VALUE}
  4343.      * </p>
  4344.      * @param a first factor
  4345.      * @param b second factor
  4346.      * @return a * b / 2<sup>64</sup>
  4347.      * @since 3.0
  4348.      */
  4349.     public static long unsignedMultiplyHigh(final long a, final long b) {

  4351.         // split numbers in two 32 bits parts
  4352.         final long aHigh  = a >>> 32;
  4353.         final long aLow   = a & 0xFFFFFFFFl;
  4354.         final long bHigh  = b >>> 32;
  4355.         final long bLow   = b & 0xFFFFFFFFl;

  4357.         // ab = aHigh * bHigh * 2⁶⁴ + (aHigh * bLow + aLow * bHigh) * 2³² + aLow * bLow
  4358.         final long hh     = aHigh * bHigh;
  4359.         final long hl1    = aHigh * bLow;
  4360.         final long hl2    = aLow  * bHigh;
  4361.         final long ll     = aLow  * bLow;

  4363.         // adds up everything in the above 64 bit part, taking care to avoid overflow
  4364.         final long hlHigh = (hl1 >>> 32) + (hl2 >>> 32);
  4365.         final long hlLow  = (hl1 & 0xFFFFFFFFl) + (hl2 & 0xFFFFFFFFl);
  4366.         final long carry  = (hlLow + (ll >>> 32)) >>> 32;

  4368.         return hh + hlHigh + carry;

  4370.     }

  4372.     /** Divide two integers, checking for overflow.
  4373.      * @param x dividend
  4374.      * @param y divisor
  4375.      * @return x / y
  4376.      * @exception MathRuntimeException if an overflow occurs
  4377.      * @since 3.0
  4378.      */
  4379.     public static int divideExact(final int x, final int y) {
  4380.         if (y == 0) {
  4381.             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
  4382.         }
  4383.         if (y == -1 && x == Integer.MIN_VALUE) {
  4384.             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, x, y);
  4385.         }
  4386.         return x / y;
  4387.     }

  4389.     /** Divide two long integers, checking for overflow.
  4390.      * @param x dividend
  4391.      * @param y divisor
  4392.      * @return x / y
  4393.      * @exception MathRuntimeException if an overflow occurs
  4394.      * @since 3.0
  4395.      */
  4396.     public static long divideExact(final long x, final long y) {
  4397.         if (y == 0l) {
  4398.             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
  4399.         }
  4400.         if (y == -1l && x == Long.MIN_VALUE) {
  4401.             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, x, y);
  4402.         }
  4403.         return x / y;
  4404.     }

  4406.     /** Finds q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}.
  4407.      * <p>
  4408.      * This methods returns the same value as integer division when
  4409.      * a and b are opposite signs, but returns a different value when
  4410.      * they are same (i.e. q is positive).
  4411.      *
  4412.      * @param a dividend
  4413.      * @param b divisor
  4414.      * @return q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}
  4415.      * @exception MathRuntimeException if b == 0
  4416.      * @since 3.0
  4417.      */
  4418.     public static int ceilDiv(final int a, final int b) throws MathRuntimeException {

  4420.         if (b == 0) {
  4421.             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
  4422.         }

  4424.         final int m = a % b;
  4425.         if ((a ^ b) < 0 || m == 0) {
  4426.             // a and b have opposite signs, or division is exact
  4427.             return a / b;
  4428.         } else {
  4429.             // a and b have same signs and division is not exact
  4430.             return (a / b) + 1;
  4431.         }

  4433.     }

  4435.     /** Finds q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}.
  4436.      * <p>
  4437.      * This methods returns the same value as integer division when
  4438.      * a and b are opposite signs, but returns a different value when
  4439.      * they are same (i.e. q is positive).
  4440.      *
  4441.      * @param a dividend
  4442.      * @param b divisor
  4443.      * @return q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}
  4444.      * @exception MathRuntimeException if b == 0 or if a == {@code Integer.MIN_VALUE} and b = -1
  4445.      * @since 3.0
  4446.      */
  4447.     public static int ceilDivExact(final int a, final int b) throws MathRuntimeException {

  4449.         if (a == Integer.MIN_VALUE && b == -1) {
  4450.             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, a, b);
  4451.         }

  4453.         return ceilDiv(a, b);

  4455.     }

  4457.     /** Finds q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}.
  4458.      * <p>
  4459.      * This methods returns the same value as integer division when
  4460.      * a and b are opposite signs, but returns a different value when
  4461.      * they are same (i.e. q is positive).
  4462.      *
  4463.      * @param a dividend
  4464.      * @param b divisor
  4465.      * @return q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}
  4466.      * @exception MathRuntimeException if b == 0
  4467.      * @since 3.0
  4468.      */
  4469.     public static long ceilDiv(final long a, final long b) throws MathRuntimeException {

  4471.         if (b == 0l) {
  4472.             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
  4473.         }

  4475.         final long m = a % b;
  4476.         if ((a ^ b) < 0 || m == 0l) {
  4477.             // a and b have opposite signs, or division is exact
  4478.             return a / b;
  4479.         } else {
  4480.             // a and b have same signs and division is not exact
  4481.             return (a / b) + 1l;
  4482.         }

  4484.     }

  4486.     /** Finds q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}.
  4487.      * <p>
  4488.      * This methods returns the same value as integer division when
  4489.      * a and b are opposite signs, but returns a different value when
  4490.      * they are same (i.e. q is positive).
  4491.      *
  4492.      * @param a dividend
  4493.      * @param b divisor
  4494.      * @return q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}
  4495.      * @exception MathRuntimeException if b == 0 or if a == {@code Long.MIN_VALUE} and b = -1
  4496.      * @since 3.0
  4497.      */
  4498.     public static long ceilDivExact(final long a, final long b) throws MathRuntimeException {

  4500.         if (a == Long.MIN_VALUE && b == -1l) {
  4501.             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, a, b);
  4502.         }

  4504.         return ceilDiv(a, b);

  4506.     }

  4508.     /** Finds q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}.
  4509.      * <p>
  4510.      * This methods returns the same value as integer division when
  4511.      * a and b are opposite signs, but returns a different value when
  4512.      * they are same (i.e. q is positive).
  4513.      *
  4514.      * @param a dividend
  4515.      * @param b divisor
  4516.      * @return q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}
  4517.      * @exception MathRuntimeException if b == 0
  4518.      * @since 3.0
  4519.      */
  4520.     public static long ceilDiv(final long a, final int b) throws MathRuntimeException {
  4521.         return ceilDiv(a, (long) b);
  4522.     }

  4524.     /** Finds r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}.
  4525.      * <p>
  4526.      * This methods returns the same value as integer modulo when
  4527.      * a and b are opposite signs, but returns a different value when
  4528.      * they are same (i.e. q is positive).
  4529.      *
  4530.      * @param a dividend
  4531.      * @param b divisor
  4532.      * @return r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}
  4533.      * @exception MathRuntimeException if b == 0
  4534.      * @since 3.0
  4535.      */
  4536.     public static int ceilMod(final int a, final int b) throws MathRuntimeException {

  4538.         if (b == 0) {
  4539.             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
  4540.         }

  4542.         final int m = a % b;
  4543.         if ((a ^ b) < 0 || m == 0) {
  4544.             // a and b have opposite signs, or division is exact
  4545.             return m;
  4546.         } else {
  4547.             // a and b have same signs and division is not exact
  4548.             return m - b;
  4549.         }

  4551.     }

  4553.     /** Finds r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}.
  4554.      * <p>
  4555.      * This methods returns the same value as integer modulo when
  4556.      * a and b are opposite signs, but returns a different value when
  4557.      * they are same (i.e. q is positive).
  4558.      *
  4559.      * @param a dividend
  4560.      * @param b divisor
  4561.      * @return r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}
  4562.      * @exception MathRuntimeException if b == 0
  4563.      * @since 3.0
  4564.      */
  4565.     public static int ceilMod(final long a, final int b) throws MathRuntimeException {
  4566.         return (int) ceilMod(a, (long) b);
  4567.     }

  4569.     /** Finds r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}.
  4570.      * <p>
  4571.      * This methods returns the same value as integer modulo when
  4572.      * a and b are opposite signs, but returns a different value when
  4573.      * they are same (i.e. q is positive).
  4574.      *
  4575.      * @param a dividend
  4576.      * @param b divisor
  4577.      * @return r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}
  4578.      * @exception MathRuntimeException if b == 0
  4579.      * @since 3.0
  4580.      */
  4581.     public static long ceilMod(final long a, final long b) throws MathRuntimeException {

  4583.         if (b == 0l) {
  4584.             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
  4585.         }

  4587.         final long m = a % b;
  4588.         if ((a ^ b) < 0l || m == 0l) {
  4589.             // a and b have opposite signs, or division is exact
  4590.             return m;
  4591.         } else {
  4592.             // a and b have same signs and division is not exact
  4593.             return m - b;
  4594.         }

  4596.     }

  4598.     /** Finds q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code  b < 0}.
  4599.      * <p>
  4600.      * This methods returns the same value as integer division when
  4601.      * a and b are same signs, but returns a different value when
  4602.      * they are opposite (i.e. q is negative).
  4603.      *
  4604.      * @param a dividend
  4605.      * @param b divisor
  4606.      * @return q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code  b < 0}
  4607.      * @exception MathRuntimeException if b == 0
  4608.      * @see #floorMod(int, int)
  4609.      */
  4610.     public static int floorDiv(final int a, final int b) throws MathRuntimeException {

  4612.         if (b == 0) {
  4613.             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
  4614.         }

  4616.         final int m = a % b;
  4617.         if ((a ^ b) >= 0 || m == 0) {
  4618.             // a and b have same sign, or division is exact
  4619.             return a / b;
  4620.         } else {
  4621.             // a and b have opposite signs and division is not exact
  4622.             return (a / b) - 1;
  4623.         }

  4625.     }

  4627.     /** Finds q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code  b < 0}.
  4628.      * <p>
  4629.      * This methods returns the same value as integer division when
  4630.      * a and b are same signs, but returns a different value when
  4631.      * they are opposite (i.e. q is negative).
  4632.      *
  4633.      * @param a dividend
  4634.      * @param b divisor
  4635.      * @return q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code  b < 0}
  4636.      * @exception MathRuntimeException if b == 0 or if a == {@code Integer.MIN_VALUE} and b = -1
  4637.      * @see #floorMod(int, int)
  4638.      * @since 3.0
  4639.      */
  4640.     public static int floorDivExact(final int a, final int b) throws MathRuntimeException {

  4642.         if (a == Integer.MIN_VALUE && b == -1) {
  4643.             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, a, b);
  4644.         }

  4646.         return floorDiv(a, b);

  4648.     }

  4650.     /** Finds q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
  4651.      * <p>
  4652.      * This methods returns the same value as integer division when
  4653.      * a and b are same signs, but returns a different value when
  4654.      * they are opposite (i.e. q is negative).
  4655.      *
  4656.      * @param a dividend
  4657.      * @param b divisor
  4658.      * @return q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
  4659.      * @exception MathRuntimeException if b == 0
  4660.      * @see #floorMod(long, int)
  4661.      * @since 1.3
  4662.      */
  4663.     public static long floorDiv(final long a, final int b) throws MathRuntimeException {
  4664.         return floorDiv(a, (long) b);
  4665.     }

  4667.     /** Finds q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
  4668.      * <p>
  4669.      * This methods returns the same value as integer division when
  4670.      * a and b are same signs, but returns a different value when
  4671.      * they are opposite (i.e. q is negative).
  4672.      *
  4673.      * @param a dividend
  4674.      * @param b divisor
  4675.      * @return q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
  4676.      * @exception MathRuntimeException if b == 0
  4677.      * @see #floorMod(long, long)
  4678.      */
  4679.     public static long floorDiv(final long a, final long b) throws MathRuntimeException {

  4681.         if (b == 0l) {
  4682.             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
  4683.         }

  4685.         final long m = a % b;
  4686.         if ((a ^ b) >= 0l || m == 0l) {
  4687.             // a and b have same sign, or division is exact
  4688.             return a / b;
  4689.         } else {
  4690.             // a and b have opposite signs and division is not exact
  4691.             return (a / b) - 1l;
  4692.         }

  4694.     }

  4696.     /** Finds q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
  4697.      * <p>
  4698.      * This methods returns the same value as integer division when
  4699.      * a and b are same signs, but returns a different value when
  4700.      * they are opposite (i.e. q is negative).
  4701.      *
  4702.      * @param a dividend
  4703.      * @param b divisor
  4704.      * @return q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
  4705.      * @exception MathRuntimeException if b == 0 or if a == {@code Long.MIN_VALUE} and b = -1
  4706.      * @see #floorMod(long, long)
  4707.      * @since 3.0
  4708.      */
  4709.     public static long floorDivExact(final long a, final long b) throws MathRuntimeException {

  4711.         if (a == Long.MIN_VALUE && b == -1l) {
  4712.             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, a, b);
  4713.         }

  4715.         return floorDiv(a, b);

  4717.     }

  4719.     /** Finds r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
  4720.      * <p>
  4721.      * This methods returns the same value as integer modulo when
  4722.      * a and b are same signs, but returns a different value when
  4723.      * they are opposite (i.e. q is negative).
  4724.      * </p>
  4725.      * @param a dividend
  4726.      * @param b divisor
  4727.      * @return r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
  4728.      * @exception MathRuntimeException if b == 0
  4729.      * @see #floorDiv(int, int)
  4730.      */
  4731.     public static int floorMod(final int a, final int b) throws MathRuntimeException {

  4733.         if (b == 0) {
  4734.             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
  4735.         }

  4737.         final int m = a % b;
  4738.         if ((a ^ b) >= 0 || m == 0) {
  4739.             // a and b have same sign, or division is exact
  4740.             return m;
  4741.         } else {
  4742.             // a and b have opposite signs and division is not exact
  4743.             return b + m;
  4744.         }

  4746.     }

  4748.     /** Finds r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
  4749.      * <p>
  4750.      * This methods returns the same value as integer modulo when
  4751.      * a and b are same signs, but returns a different value when
  4752.      * they are opposite (i.e. q is negative).
  4753.      * </p>
  4754.      * @param a dividend
  4755.      * @param b divisor
  4756.      * @return r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
  4757.      * @exception MathRuntimeException if b == 0
  4758.      * @see #floorDiv(long, int)
  4759.      * @since 1.3
  4760.      */
  4761.     public static int floorMod(final long a, final int b) {
  4762.         return (int) floorMod(a, (long) b);
  4763.     }

  4765.     /** Finds r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
  4766.      * <p>
  4767.      * This methods returns the same value as integer modulo when
  4768.      * a and b are same signs, but returns a different value when
  4769.      * they are opposite (i.e. q is negative).
  4770.      * </p>
  4771.      * @param a dividend
  4772.      * @param b divisor
  4773.      * @return r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
  4774.      * @exception MathRuntimeException if b == 0
  4775.      * @see #floorDiv(long, long)
  4776.      */
  4777.     public static long floorMod(final long a, final long b) {

  4779.         if (b == 0l) {
  4780.             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
  4781.         }

  4783.         final long m = a % b;
  4784.         if ((a ^ b) >= 0l || m == 0l) {
  4785.             // a and b have same sign, or division is exact
  4786.             return m;
  4787.         } else {
  4788.             // a and b have opposite signs and division is not exact
  4789.             return b + m;
  4790.         }

  4792.     }

  4794.     /**
  4795.      * Returns the first argument with the sign of the second argument.
  4796.      * A NaN {@code sign} argument is treated as positive.
  4797.      *
  4798.      * @param magnitude the value to return
  4799.      * @param sign the sign for the returned value
  4800.      * @return the magnitude with the same sign as the {@code sign} argument
  4801.      */
  4802.     public static double copySign(double magnitude, double sign){
  4803.         // The highest order bit is going to be zero if the
  4804.         // highest order bit of m and s is the same and one otherwise.
  4805.         // So (m^s) will be positive if both m and s have the same sign
  4806.         // and negative otherwise.
  4807.         final long m = Double.doubleToRawLongBits(magnitude); // don't care about NaN
  4808.         final long s = Double.doubleToRawLongBits(sign);
  4809.         if ((m ^ s) >= 0) {
  4810.             return magnitude;
  4811.         }
  4812.         return -magnitude; // flip sign
  4813.     }

  4815.     /**
  4816.      * Returns the first argument with the sign of the second argument.
  4817.      * A NaN {@code sign} argument is treated as positive.
  4818.      *
  4819.      * @param magnitude the value to return
  4820.      * @param sign the sign for the returned value
  4821.      * @return the magnitude with the same sign as the {@code sign} argument
  4822.      */
  4823.     public static float copySign(float magnitude, float sign){
  4824.         // The highest order bit is going to be zero if the
  4825.         // highest order bit of m and s is the same and one otherwise.
  4826.         // So (m^s) will be positive if both m and s have the same sign
  4827.         // and negative otherwise.
  4828.         final int m = Float.floatToRawIntBits(magnitude);
  4829.         final int s = Float.floatToRawIntBits(sign);
  4830.         if ((m ^ s) >= 0) {
  4831.             return magnitude;
  4832.         }
  4833.         return -magnitude; // flip sign
  4834.     }

  4836.     /**
  4837.      * Return the exponent of a double number, removing the bias.
  4838.      * <p>
  4839.      * For double numbers of the form 2<sup>x</sup>, the unbiased
  4840.      * exponent is exactly x.
  4841.      * </p>
  4842.      * @param d number from which exponent is requested
  4843.      * @return exponent for d in IEEE754 representation, without bias
  4844.      */
  4845.     public static int getExponent(final double d) {
  4846.         // NaN and Infinite will return 1024 anyhow so can use raw bits
  4847.         return (int) ((Double.doubleToRawLongBits(d) >>> 52) & 0x7ff) - 1023;
  4848.     }

  4850.     /**
  4851.      * Return the exponent of a float number, removing the bias.
  4852.      * <p>
  4853.      * For float numbers of the form 2<sup>x</sup>, the unbiased
  4854.      * exponent is exactly x.
  4855.      * </p>
  4856.      * @param f number from which exponent is requested
  4857.      * @return exponent for d in IEEE754 representation, without bias
  4858.      */
  4859.     public static int getExponent(final float f) {
  4860.         // NaN and Infinite will return the same exponent anyhow so can use raw bits
  4861.         return ((Float.floatToRawIntBits(f) >>> 23) & 0xff) - 127;
  4862.     }

  4864.     /** Compute Fused-multiply-add operation a * b + c.
  4865.      * <p>
  4866.      * This method was introduced in the regular {@code Math} and {@code StrictMath}
  4867.      * methods with Java 9, and then added to Hipparchus for consistency. However,
  4868.      * a more general method was available in Hipparchus that also allow to repeat
  4869.      * this computation across several terms: {@link MathArrays#linearCombination(double[], double[])}.
  4870.      * The linear combination method should probably be preferred in most cases.
  4871.      * </p>
  4872.      * @param a first factor
  4873.      * @param b second factor
  4874.      * @param c additive term
  4875.      * @return a * b + c, using extended precision in the multiplication
  4876.      * @see MathArrays#linearCombination(double[], double[])
  4877.      * @see MathArrays#linearCombination(double, double, double, double)
  4878.      * @see MathArrays#linearCombination(double, double, double, double, double, double)
  4879.      * @see MathArrays#linearCombination(double, double, double, double, double, double, double, double)
  4880.      * @since 1.3
  4881.      */
  4882.     public static double fma(final double a, final double b, final double c) {
  4883.         return MathArrays.linearCombination(a, b, 1.0, c);
  4884.     }

  4886.     /** Compute Fused-multiply-add operation a * b + c.
  4887.      * <p>
  4888.      * This method was introduced in the regular {@code Math} and {@code StrictMath}
  4889.      * methods with Java 9, and then added to Hipparchus for consistency. However,
  4890.      * a more general method was available in Hipparchus that also allow to repeat
  4891.      * this computation across several terms: {@link MathArrays#linearCombination(double[], double[])}.
  4892.      * The linear combination method should probably be preferred in most cases.
  4893.      * </p>
  4894.      * @param a first factor
  4895.      * @param b second factor
  4896.      * @param c additive term
  4897.      * @return a * b + c, using extended precision in the multiplication
  4898.      * @see MathArrays#linearCombination(double[], double[])
  4899.      * @see MathArrays#linearCombination(double, double, double, double)
  4900.      * @see MathArrays#linearCombination(double, double, double, double, double, double)
  4901.      * @see MathArrays#linearCombination(double, double, double, double, double, double, double, double)
  4902.      */
  4903.     public static float fma(final float a, final float b, final float c) {
  4904.         return (float) MathArrays.linearCombination(a, b, 1.0, c);
  4905.     }

  4907.     /** Compute the square root of a number.
  4908.      * @param a number on which evaluation is done
  4909.      * @param <T> the type of the field element
  4910.      * @return square root of a
  4911.      * @since 1.3
  4912.      */
  4913.     public static <T extends CalculusFieldElement<T>> T sqrt(final T a) {
  4914.         return a.sqrt();
  4915.     }

  4917.     /** Compute the hyperbolic cosine of a number.
  4918.      * @param x number on which evaluation is done
  4919.      * @param <T> the type of the field element
  4920.      * @return hyperbolic cosine of x
  4921.      * @since 1.3
  4922.      */
  4923.     public static <T extends CalculusFieldElement<T>> T cosh(final T x) {
  4924.         return x.cosh();
  4925.     }

  4927.     /** Compute the hyperbolic sine of a number.
  4928.      * @param x number on which evaluation is done
  4929.      * @param <T> the type of the field element
  4930.      * @return hyperbolic sine of x
  4931.      * @since 1.3
  4932.      */
  4933.     public static <T extends CalculusFieldElement<T>> T sinh(final T x) {
  4934.         return x.sinh();
  4935.     }

  4937.     /** Compute the hyperbolic tangent of a number.
  4938.      * @param x number on which evaluation is done
  4939.      * @param <T> the type of the field element
  4940.      * @return hyperbolic tangent of x
  4941.      * @since 1.3
  4942.      */
  4943.     public static <T extends CalculusFieldElement<T>> T tanh(final T x) {
  4944.         return x.tanh();
  4945.     }

  4947.     /** Compute the inverse hyperbolic cosine of a number.
  4948.      * @param a number on which evaluation is done
  4949.      * @param <T> the type of the field element
  4950.      * @return inverse hyperbolic cosine of a
  4951.      * @since 1.3
  4952.      */
  4953.     public static <T extends CalculusFieldElement<T>> T acosh(final T a) {
  4954.         return a.acosh();
  4955.     }

  4957.     /** Compute the inverse hyperbolic sine of a number.
  4958.      * @param a number on which evaluation is done
  4959.      * @param <T> the type of the field element
  4960.      * @return inverse hyperbolic sine of a
  4961.      * @since 1.3
  4962.      */
  4963.     public static <T extends CalculusFieldElement<T>> T asinh(final T a) {
  4964.         return a.asinh();
  4965.     }

  4967.     /** Compute the inverse hyperbolic tangent of a number.
  4968.      * @param a number on which evaluation is done
  4969.      * @param <T> the type of the field element
  4970.      * @return inverse hyperbolic tangent of a
  4971.      * @since 1.3
  4972.      */
  4973.     public static <T extends CalculusFieldElement<T>> T atanh(final T a) {
  4974.         return a.atanh();
  4975.     }

  4977.     /** Compute the sign of a number.
  4978.      * The sign is -1 for negative numbers, +1 for positive numbers and 0 otherwise,
  4979.      * for Complex number, it is extended on the unit circle (equivalent to z/|z|,
  4980.      * with special handling for 0 and NaN)
  4981.      * @param a number on which evaluation is done
  4982.      * @param <T> the type of the field element
  4983.      * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
  4984.      * @since 2.0
  4985.      */
  4986.     public static <T extends CalculusFieldElement<T>> T sign(final T a) {
  4987.         return a.sign();
  4988.     }

  4990.     /**
  4991.      * Exponential function.
  4992.      *
  4993.      * Computes exp(x), function result is nearly rounded.   It will be correctly
  4994.      * rounded to the theoretical value for 99.9% of input values, otherwise it will
  4995.      * have a 1 ULP error.
  4996.      *
  4997.      * Method:
  4998.      *    Lookup intVal = exp(int(x))
  4999.      *    Lookup fracVal = exp(int(x-int(x) / 1024.0) * 1024.0 );
  5000.      *    Compute z as the exponential of the remaining bits by a polynomial minus one
  5001.      *    exp(x) = intVal * fracVal * (1 + z)
  5002.      *
  5003.      * Accuracy:
  5004.      *    Calculation is done with 63 bits of precision, so result should be correctly
  5005.      *    rounded for 99.9% of input values, with less than 1 ULP error otherwise.
  5006.      *
  5007.      * @param x   a double
  5008.      * @param <T> the type of the field element
  5009.      * @return double e<sup>x</sup>
  5010.      * @since 1.3
  5011.      */
  5012.     public static <T extends CalculusFieldElement<T>> T exp(final T x) {
  5013.         return x.exp();
  5014.     }

  5016.     /** Compute exp(x) - 1
  5017.      * @param x number to compute shifted exponential
  5018.      * @param <T> the type of the field element
  5019.      * @return exp(x) - 1
  5020.      * @since 1.3
  5021.      */
  5022.     public static <T extends CalculusFieldElement<T>> T expm1(final T x) {
  5023.         return x.expm1();
  5024.     }

  5026.     /**
  5027.      * Natural logarithm.
  5028.      *
  5029.      * @param x   a double
  5030.      * @param <T> the type of the field element
  5031.      * @return log(x)
  5032.      * @since 1.3
  5033.      */
  5034.     public static <T extends CalculusFieldElement<T>> T log(final T x) {
  5035.         return x.log();
  5036.     }

  5038.     /**
  5039.      * Computes log(1 + x).
  5040.      *
  5041.      * @param x Number.
  5042.      * @param <T> the type of the field element
  5043.      * @return {@code log(1 + x)}.
  5044.      * @since 1.3
  5045.      */
  5046.     public static <T extends CalculusFieldElement<T>> T log1p(final T x) {
  5047.         return x.log1p();
  5048.     }

  5050.     /** Compute the base 10 logarithm.
  5051.      * @param x a number
  5052.      * @param <T> the type of the field element
  5053.      * @return log10(x)
  5054.      * @since 1.3
  5055.      */
  5056.     public static <T extends CalculusFieldElement<T>> T log10(final T x) {
  5057.         return x.log10();
  5058.     }

  5060.     /**
  5061.      * Power function.  Compute x<sup>y</sup>.
  5062.      *
  5063.      * @param x   a double
  5064.      * @param y   a double
  5065.      * @param <T> the type of the field element
  5066.      * @return x<sup>y</sup>
  5067.      * @since 1.3
  5068.      */
  5069.     public static <T extends CalculusFieldElement<T>> T pow(final T x, final T y) {
  5070.         return x.pow(y);
  5071.     }

  5073.     /**
  5074.      * Power function.  Compute x<sup>y</sup>.
  5075.      *
  5076.      * @param x   a double
  5077.      * @param y   a double
  5078.      * @param <T> the type of the field element
  5079.      * @return x<sup>y</sup>
  5080.      * @since 1.7
  5081.      */
  5082.     public static <T extends CalculusFieldElement<T>> T pow(final T x, final double y) {
  5083.         return x.pow(y);
  5084.     }

  5086.     /**
  5087.      * Raise a double to an int power.
  5088.      *
  5089.      * @param d Number to raise.
  5090.      * @param e Exponent.
  5091.      * @param <T> the type of the field element
  5092.      * @return d<sup>e</sup>
  5093.      * @since 1.3
  5094.      */
  5095.     public static <T extends CalculusFieldElement<T>> T pow(T d, int e) {
  5096.         return d.pow(e);
  5097.     }

  5099.     /**
  5100.      * Sine function.
  5101.      *
  5102.      * @param x Argument.
  5103.      * @param <T> the type of the field element
  5104.      * @return sin(x)
  5105.      * @since 1.3
  5106.      */
  5107.     public static <T extends CalculusFieldElement<T>> T sin(final T x) {
  5108.         return x.sin();
  5109.     }

  5111.     /**
  5112.      * Cosine function.
  5113.      *
  5114.      * @param x Argument.
  5115.      * @param <T> the type of the field element
  5116.      * @return cos(x)
  5117.      * @since 1.3
  5118.      */
  5119.     public static <T extends CalculusFieldElement<T>> T cos(final T x) {
  5120.         return x.cos();
  5121.     }

  5123.     /**
  5124.      * Tangent function.
  5125.      *
  5126.      * @param x Argument.
  5127.      * @param <T> the type of the field element
  5128.      * @return tan(x)
  5129.      * @since 1.3
  5130.      */
  5131.     public static <T extends CalculusFieldElement<T>> T tan(final T x) {
  5132.         return x.tan();
  5133.     }

  5135.     /**
  5136.      * Arctangent function
  5137.      *  @param x a number
  5138.      * @param <T> the type of the field element
  5139.      *  @return atan(x)
  5140.      * @since 1.3
  5141.      */
  5142.     public static <T extends CalculusFieldElement<T>> T atan(final T x) {
  5143.         return x.atan();
  5144.     }

  5146.     /**
  5147.      * Two arguments arctangent function
  5148.      * @param y ordinate
  5149.      * @param x abscissa
  5150.      * @param <T> the type of the field element
  5151.      * @return phase angle of point (x,y) between {@code -PI} and {@code PI}
  5152.      * @since 1.3
  5153.      */
  5154.     public static <T extends CalculusFieldElement<T>> T atan2(final T y, final T x) {
  5155.         return y.atan2(x);
  5156.     }

  5158.     /** Compute the arc sine of a number.
  5159.      * @param x number on which evaluation is done
  5160.      * @param <T> the type of the field element
  5161.      * @return arc sine of x
  5162.      * @since 1.3
  5163.      */
  5164.     public static <T extends CalculusFieldElement<T>> T asin(final T x) {
  5165.         return x.asin();
  5166.     }

  5168.     /** Compute the arc cosine of a number.
  5169.      * @param x number on which evaluation is done
  5170.      * @param <T> the type of the field element
  5171.      * @return arc cosine of x
  5172.      * @since 1.3
  5173.      */
  5174.     public static <T extends CalculusFieldElement<T>> T acos(final T x) {
  5175.         return x.acos();
  5176.     }

  5178.     /** Compute the cubic root of a number.
  5179.      * @param x number on which evaluation is done
  5180.      * @param <T> the type of the field element
  5181.      * @return cubic root of x
  5182.      * @since 1.3
  5183.      */
  5184.     public static <T extends CalculusFieldElement<T>> T cbrt(final T x) {
  5185.         return x.cbrt();
  5186.     }

  5188.     /**
  5189.      * Norm.
  5190.      * @param x number from which norm is requested
  5191.      * @param <T> the type of the field element
  5192.      * @return norm(x)
  5193.      * @since 2.0
  5194.      */
  5195.     public static <T extends CalculusFieldElement<T>> double norm(final T x) {
  5196.         return x.norm();
  5197.     }

  5199.     /**
  5200.      * Absolute value.
  5201.      * @param x number from which absolute value is requested
  5202.      * @param <T> the type of the field element
  5203.      * @return abs(x)
  5204.      * @since 2.0
  5205.      */
  5206.     public static <T extends CalculusFieldElement<T>> T abs(final T x) {
  5207.         return x.abs();
  5208.     }

  5210.     /**
  5211.      *  Convert degrees to radians, with error of less than 0.5 ULP
  5212.      *  @param x angle in degrees
  5213.      *  @param <T> the type of the field element
  5214.      *  @return x converted into radians
  5215.      */
  5216.     public static <T extends CalculusFieldElement<T>> T toRadians(T x) {
  5217.         return x.toRadians();
  5218.     }

  5220.     /**
  5221.      *  Convert radians to degrees, with error of less than 0.5 ULP
  5222.      *  @param x angle in radians
  5223.      *  @param <T> the type of the field element
  5224.      *  @return x converted into degrees
  5225.      */
  5226.     public static <T extends CalculusFieldElement<T>> T toDegrees(T x) {
  5227.         return x.toDegrees();
  5228.     }

  5230.     /**
  5231.      * Multiply a double number by a power of 2.
  5232.      * @param d number to multiply
  5233.      * @param n power of 2
  5234.      * @param <T> the type of the field element
  5235.      * @return d &times; 2<sup>n</sup>
  5236.      * @since 1.3
  5237.      */
  5238.     public static <T extends CalculusFieldElement<T>> T scalb(final T d, final int n) {
  5239.         return d.scalb(n);
  5240.     }

  5242.     /**
  5243.      * Compute least significant bit (Unit in Last Position) for a number.
  5244.      * @param x number from which ulp is requested
  5245.      * @param <T> the type of the field element
  5246.      * @return ulp(x)
  5247.      * @since 2.0
  5248.      */
  5249.     public static <T extends CalculusFieldElement<T>> T ulp(final T x) {
  5250.         if (Double.isInfinite(x.getReal())) {
  5251.             return x.newInstance(Double.POSITIVE_INFINITY);
  5252.         }
  5253.         return x.ulp();
  5254.     }

  5256.     /** Get the largest whole number smaller than x.
  5257.      * @param x number from which floor is requested
  5258.      * @param <T> the type of the field element
  5259.      * @return a double number f such that f is an integer f &lt;= x &lt; f + 1.0
  5260.      * @since 1.3
  5261.      */
  5262.     public static <T extends CalculusFieldElement<T>> T floor(final T x) {
  5263.         return x.floor();
  5264.     }

  5266.     /** Get the smallest whole number larger than x.
  5267.      * @param x number from which ceil is requested
  5268.      * @param <T> the type of the field element
  5269.      * @return a double number c such that c is an integer c - 1.0 &lt; x &lt;= c
  5270.      * @since 1.3
  5271.      */
  5272.     public static <T extends CalculusFieldElement<T>> T ceil(final T x) {
  5273.         return x.ceil();
  5274.     }

  5276.     /** Get the whole number that is the nearest to x, or the even one if x is exactly half way between two integers.
  5277.      * @param x number from which nearest whole number is requested
  5278.      * @param <T> the type of the field element
  5279.      * @return a double number r such that r is an integer r - 0.5 &lt;= x &lt;= r + 0.5
  5280.      * @since 1.3
  5281.      */
  5282.     public static <T extends CalculusFieldElement<T>> T rint(final T x) {
  5283.         return x.rint();
  5284.     }

  5286.     /** Get the closest long to x.
  5287.      * @param x number from which closest long is requested
  5288.      * @param <T> the type of the field element
  5289.      * @return closest long to x
  5290.      * @since 1.3
  5291.      */
  5292.     public static <T extends CalculusFieldElement<T>> long round(final T x) {
  5293.         return x.round();
  5294.     }

  5296.     /** Compute the minimum of two values
  5297.      * @param a first value
  5298.      * @param b second value
  5299.      * @param <T> the type of the field element
  5300.      * @return a if a is lesser or equal to b, b otherwise
  5301.      * @since 1.3
  5302.      */
  5303.     public static <T extends CalculusFieldElement<T>> T min(final T a, final T b) {
  5304.         final double aR = a.getReal();
  5305.         final double bR = b.getReal();
  5306.         if (aR < bR) {
  5307.             return a;
  5308.         } else if (bR < aR) {
  5309.             return b;
  5310.         } else {
  5311.             // either the numbers are equal, or one of them is a NaN
  5312.             return Double.isNaN(aR) ? a : b;
  5313.         }
  5314.     }

  5316.     /** Compute the minimum of two values
  5317.      * @param a first value
  5318.      * @param b second value
  5319.      * @param <T> the type of the field element
  5320.      * @return a if a is lesser or equal to b, b otherwise
  5321.      * @since 1.3
  5322.      */
  5323.     public static <T extends CalculusFieldElement<T>> T min(final T a, final double b) {
  5324.         final double aR = a.getReal();
  5325.         if (aR < b) {
  5326.             return a;
  5327.         } else if (b < aR) {
  5328.             return a.getField().getZero().add(b);
  5329.         } else {
  5330.             // either the numbers are equal, or one of them is a NaN
  5331.             return Double.isNaN(aR) ? a : a.getField().getZero().add(b);
  5332.         }
  5333.     }

  5335.     /** Compute the maximum of two values
  5336.      * @param a first value
  5337.      * @param b second value
  5338.      * @param <T> the type of the field element
  5339.      * @return b if a is lesser or equal to b, a otherwise
  5340.      * @since 1.3
  5341.      */
  5342.     public static <T extends CalculusFieldElement<T>> T max(final T a, final T b) {
  5343.         final double aR = a.getReal();
  5344.         final double bR = b.getReal();
  5345.         if (aR < bR) {
  5346.             return b;
  5347.         } else if (bR < aR) {
  5348.             return a;
  5349.         } else {
  5350.             // either the numbers are equal, or one of them is a NaN
  5351.             return Double.isNaN(aR) ? a : b;
  5352.         }
  5353.     }

  5355.     /** Compute the maximum of two values
  5356.      * @param a first value
  5357.      * @param b second value
  5358.      * @param <T> the type of the field element
  5359.      * @return b if a is lesser or equal to b, a otherwise
  5360.      * @since 1.3
  5361.      */
  5362.     public static <T extends CalculusFieldElement<T>> T max(final T a, final double b) {
  5363.         final double aR = a.getReal();
  5364.         if (aR < b) {
  5365.             return a.getField().getZero().add(b);
  5366.         } else if (b < aR) {
  5367.             return a;
  5368.         } else {
  5369.             // either the numbers are equal, or one of them is a NaN
  5370.             return Double.isNaN(aR) ? a : a.getField().getZero().add(b);
  5371.         }
  5372.     }

  5374.     /**
  5375.      * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
  5376.      * - sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)<br>
  5377.      * avoiding intermediate overflow or underflow.
  5378.      *
  5379.      * <ul>
  5380.      * <li> If either argument is infinite, then the result is positive infinity.</li>
  5381.      * <li> else, if either argument is NaN then the result is NaN.</li>
  5382.      * </ul>
  5383.      *
  5384.      * @param x a value
  5385.      * @param y a value
  5386.      * @param <T> the type of the field element
  5387.      * @return sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
  5388.      * @since 1.3
  5389.      */
  5390.     public static <T extends CalculusFieldElement<T>> T hypot(final T x, final T y) {
  5391.         return x.hypot(y);
  5392.     }

  5394.     /**
  5395.      * Computes the remainder as prescribed by the IEEE 754 standard.
  5396.      * <p>
  5397.      * The remainder value is mathematically equal to {@code x - y*n}
  5398.      * where {@code n} is the mathematical integer closest to the exact mathematical value
  5399.      * of the quotient {@code x/y}.
  5400.      * If two mathematical integers are equally close to {@code x/y} then
  5401.      * {@code n} is the integer that is even.
  5402.      * </p>
  5403.      * <ul>
  5404.      * <li>If either operand is NaN, the result is NaN.</li>
  5405.      * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
  5406.      * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
  5407.      * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
  5408.      * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
  5409.      * </ul>
  5410.      * @param dividend the number to be divided
  5411.      * @param divisor the number by which to divide
  5412.      * @param <T> the type of the field element
  5413.      * @return the remainder, rounded
  5414.      * @since 1.3
  5415.      */
  5416.     public static <T extends CalculusFieldElement<T>> T IEEEremainder(final T dividend, final double divisor) {
  5417.         return dividend.remainder(divisor);
  5418.     }

  5420.     /**
  5421.      * Computes the remainder as prescribed by the IEEE 754 standard.
  5422.      * <p>
  5423.      * The remainder value is mathematically equal to {@code x - y*n}
  5424.      * where {@code n} is the mathematical integer closest to the exact mathematical value
  5425.      * of the quotient {@code x/y}.
  5426.      * If two mathematical integers are equally close to {@code x/y} then
  5427.      * {@code n} is the integer that is even.
  5428.      * </p>
  5429.      * <ul>
  5430.      * <li>If either operand is NaN, the result is NaN.</li>
  5431.      * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
  5432.      * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
  5433.      * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
  5434.      * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
  5435.      * </ul>
  5436.      * @param dividend the number to be divided
  5437.      * @param divisor the number by which to divide
  5438.      * @param <T> the type of the field element
  5439.      * @return the remainder, rounded
  5440.      * @since 1.3
  5441.      */
  5442.     public static <T extends CalculusFieldElement<T>> T IEEEremainder(final T dividend, final T divisor) {
  5443.         return dividend.remainder(divisor);
  5444.     }

  5446.     /**
  5447.      * Returns the first argument with the sign of the second argument.
  5448.      * A NaN {@code sign} argument is treated as positive.
  5449.      *
  5450.      * @param magnitude the value to return
  5451.      * @param sign the sign for the returned value
  5452.      * @param <T> the type of the field element
  5453.      * @return the magnitude with the same sign as the {@code sign} argument
  5454.      * @since 1.3
  5455.      */
  5456.     public static <T extends CalculusFieldElement<T>> T copySign(T magnitude, T sign) {
  5457.         return magnitude.copySign(sign);
  5458.     }

  5460.     /**
  5461.      * Returns the first argument with the sign of the second argument.
  5462.      * A NaN {@code sign} argument is treated as positive.
  5463.      *
  5464.      * @param magnitude the value to return
  5465.      * @param sign the sign for the returned value
  5466.      * @param <T> the type of the field element
  5467.      * @return the magnitude with the same sign as the {@code sign} argument
  5468.      * @since 1.3
  5469.      */
  5470.     public static <T extends CalculusFieldElement<T>> T copySign(T magnitude, double sign) {
  5471.         return magnitude.copySign(sign);
  5472.     }

  5474. //    /**
  5475. //     * Print out contents of arrays, and check the length.
  5476. //     * <p>used to generate the preset arrays originally.</p>
  5477. //     * @param a unused
  5478. //     */
  5479. //    public static void main(String[] a) {
  5480. //        FastMathCalc.printarray(System.out, "EXP_INT_TABLE_A", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_A);
  5481. //        FastMathCalc.printarray(System.out, "EXP_INT_TABLE_B", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_B);
  5482. //        FastMathCalc.printarray(System.out, "EXP_FRAC_TABLE_A", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_A);
  5483. //        FastMathCalc.printarray(System.out, "EXP_FRAC_TABLE_B", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_B);
  5484. //        FastMathCalc.printarray(System.out, "LN_MANT",LN_MANT_LEN, lnMant.LN_MANT);
  5485. //        FastMathCalc.printarray(System.out, "SINE_TABLE_A", SINE_TABLE_LEN, SINE_TABLE_A);
  5486. //        FastMathCalc.printarray(System.out, "SINE_TABLE_B", SINE_TABLE_LEN, SINE_TABLE_B);
  5487. //        FastMathCalc.printarray(System.out, "COSINE_TABLE_A", SINE_TABLE_LEN, COSINE_TABLE_A);
  5488. //        FastMathCalc.printarray(System.out, "COSINE_TABLE_B", SINE_TABLE_LEN, COSINE_TABLE_B);
  5489. //        FastMathCalc.printarray(System.out, "TANGENT_TABLE_A", SINE_TABLE_LEN, TANGENT_TABLE_A);
  5490. //        FastMathCalc.printarray(System.out, "TANGENT_TABLE_B", SINE_TABLE_LEN, TANGENT_TABLE_B);
  5491. //    }

  5493.     /** Enclose large data table in nested static class so it's only loaded on first access. */
  5494.     private static class ExpIntTable {
  5495.         /** Exponential evaluated at integer values,
  5496.          * exp(x) =  expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX].
  5497.          */
  5498.         private static final double[] EXP_INT_TABLE_A;
  5499.         /** Exponential evaluated at integer values,
  5500.          * exp(x) =  expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX]
  5501.          */
  5502.         private static final double[] EXP_INT_TABLE_B;

  5504.         static {
  5505.             if (RECOMPUTE_TABLES_AT_RUNTIME) {
  5506.                 EXP_INT_TABLE_A = new double[EXP_INT_TABLE_LEN];
  5507.                 EXP_INT_TABLE_B = new double[EXP_INT_TABLE_LEN];

  5509.                 final double[] tmp = new double[2];
  5510.                 final double[] recip = new double[2];

  5512.                 // Populate expIntTable
  5513.                 for (int i = 0; i < EXP_INT_TABLE_MAX_INDEX; i++) {
  5514.                     FastMathCalc.expint(i, tmp);
  5515.                     EXP_INT_TABLE_A[i + EXP_INT_TABLE_MAX_INDEX] = tmp[0];
  5516.                     EXP_INT_TABLE_B[i + EXP_INT_TABLE_MAX_INDEX] = tmp[1];

  5518.                     if (i != 0) {
  5519.                         // Negative integer powers
  5520.                         FastMathCalc.splitReciprocal(tmp, recip);
  5521.                         EXP_INT_TABLE_A[EXP_INT_TABLE_MAX_INDEX - i] = recip[0];
  5522.                         EXP_INT_TABLE_B[EXP_INT_TABLE_MAX_INDEX - i] = recip[1];
  5523.                     }
  5524.                 }
  5525.             } else {
  5526.                 EXP_INT_TABLE_A = FastMathLiteralArrays.loadExpIntA();
  5527.                 EXP_INT_TABLE_B = FastMathLiteralArrays.loadExpIntB();
  5528.             }
  5529.         }
  5530.     }

  5532.     /** Enclose large data table in nested static class so it's only loaded on first access. */
  5533.     private static class ExpFracTable {
  5534.         /** Exponential over the range of 0 - 1 in increments of 2^-10
  5535.          * exp(x/1024) =  expFracTableA[x] + expFracTableB[x].
  5536.          * 1024 = 2^10
  5537.          */
  5538.         private static final double[] EXP_FRAC_TABLE_A;
  5539.         /** Exponential over the range of 0 - 1 in increments of 2^-10
  5540.          * exp(x/1024) =  expFracTableA[x] + expFracTableB[x].
  5541.          */
  5542.         private static final double[] EXP_FRAC_TABLE_B;

  5544.         static {
  5545.             if (RECOMPUTE_TABLES_AT_RUNTIME) {
  5546.                 EXP_FRAC_TABLE_A = new double[EXP_FRAC_TABLE_LEN];
  5547.                 EXP_FRAC_TABLE_B = new double[EXP_FRAC_TABLE_LEN];

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

  5551.                 // Populate expFracTable
  5552.                 final double factor = 1d / (EXP_FRAC_TABLE_LEN - 1);
  5553.                 for (int i = 0; i < EXP_FRAC_TABLE_A.length; i++) {
  5554.                     FastMathCalc.slowexp(i * factor, tmp);
  5555.                     EXP_FRAC_TABLE_A[i] = tmp[0];
  5556.                     EXP_FRAC_TABLE_B[i] = tmp[1];
  5557.                 }
  5558.             } else {
  5559.                 EXP_FRAC_TABLE_A = FastMathLiteralArrays.loadExpFracA();
  5560.                 EXP_FRAC_TABLE_B = FastMathLiteralArrays.loadExpFracB();
  5561.             }
  5562.         }
  5563.     }

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

  5570.         static {
  5571.             if (RECOMPUTE_TABLES_AT_RUNTIME) {
  5572.                 LN_MANT = new double[LN_MANT_LEN][];

  5574.                 // Populate lnMant table
  5575.                 for (int i = 0; i < LN_MANT.length; i++) {
  5576.                     final double d = Double.longBitsToDouble( (((long) i) << 42) | 0x3ff0000000000000L );
  5577.                     LN_MANT[i] = FastMathCalc.slowLog(d);
  5578.                 }
  5579.             } else {
  5580.                 LN_MANT = FastMathLiteralArrays.loadLnMant();
  5581.             }
  5582.         }
  5583.     }

  5585.     /** Enclose the Cody/Waite reduction (used in "sin", "cos" and "tan"). */
  5586.     private static class CodyWaite {
  5587.         /** k */
  5588.         private final int finalK;
  5589.         /** remA */
  5590.         private final double finalRemA;
  5591.         /** remB */
  5592.         private final double finalRemB;

  5594.         /**
  5595.          * @param xa Argument.
  5596.          */
  5597.         CodyWaite(double xa) {
  5598.             // Estimate k.
  5599.             //k = (int)(xa / 1.5707963267948966);
  5600.             int k = (int)(xa * 0.6366197723675814);

  5602.             // Compute remainder.
  5603.             double remA;
  5604.             double remB;
  5605.             while (true) {
  5606.                 double a = -k * 1.570796251296997;
  5607.                 remA = xa + a;
  5608.                 remB = -(remA - xa - a);

  5610.                 a = -k * 7.549789948768648E-8;
  5611.                 double b = remA;
  5612.                 remA = a + b;
  5613.                 remB += -(remA - b - a);

  5615.                 a = -k * 6.123233995736766E-17;
  5616.                 b = remA;
  5617.                 remA = a + b;
  5618.                 remB += -(remA - b - a);

  5620.                 if (remA > 0) {
  5621.                     break;
  5622.                 }

  5624.                 // Remainder is negative, so decrement k and try again.
  5625.                 // This should only happen if the input is very close
  5626.                 // to an even multiple of pi/2.
  5627.                 --k;
  5628.             }

  5630.             this.finalK = k;
  5631.             this.finalRemA = remA;
  5632.             this.finalRemB = remB;
  5633.         }

  5635.         /**
  5636.          * @return k
  5637.          */
  5638.         int getK() {
  5639.             return finalK;
  5640.         }
  5641.         /**
  5642.          * @return remA
  5643.          */
  5644.         double getRemA() {
  5645.             return finalRemA;
  5646.         }
  5647.         /**
  5648.          * @return remB
  5649.          */
  5650.         double getRemB() {
  5651.             return finalRemB;
  5652.         }
  5653.     }

  5655. }
Created with JaCoCo 0.8.12.202403310830