diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math/util/FastMath.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math/util/FastMath.java | 4047 |
1 files changed, 4047 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math/util/FastMath.java b/src/main/java/org/apache/commons/math/util/FastMath.java new file mode 100644 index 0000000..1907b32 --- /dev/null +++ b/src/main/java/org/apache/commons/math/util/FastMath.java @@ -0,0 +1,4047 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math.util; + +/** + * Faster, more accurate, portable alternative to {@link StrictMath}. + * <p> + * Additionally implements the following methods not found in StrictMath: + * <ul> + * <li>{@link #asinh(double)}</li> + * <li>{@link #acosh(double)}</li> + * <li>{@link #atanh(double)}</li> + * </ul> + * The following methods are found in StrictMath since 1.6 only + * <ul> + * <li>{@link #copySign(double, double)}</li> + * <li>{@link #getExponent(double)}</li> + * <li>{@link #nextAfter(double,double)}</li> + * <li>{@link #nextUp(double)}</li> + * <li>{@link #scalb(double, int)}</li> + * <li>{@link #copySign(float, float)}</li> + * <li>{@link #getExponent(float)}</li> + * <li>{@link #nextAfter(float,double)}</li> + * <li>{@link #nextUp(float)}</li> + * <li>{@link #scalb(float, int)}</li> + * </ul> + * @version $Revision: 1074294 $ $Date: 2011-02-24 22:18:59 +0100 (jeu. 24 févr. 2011) $ + * @since 2.2 + */ +public class FastMath { + + /** Archimede's constant PI, ratio of circle circumference to diameter. */ + public static final double PI = 105414357.0 / 33554432.0 + 1.984187159361080883e-9; + + /** Napier's constant e, base of the natural logarithm. */ + public static final double E = 2850325.0 / 1048576.0 + 8.254840070411028747e-8; + + /** Exponential evaluated at integer values, + * exp(x) = expIntTableA[x + 750] + expIntTableB[x+750]. + */ + private static final double EXP_INT_TABLE_A[] = new double[1500]; + + /** Exponential evaluated at integer values, + * exp(x) = expIntTableA[x + 750] + expIntTableB[x+750] + */ + private static final double EXP_INT_TABLE_B[] = new double[1500]; + + /** Exponential over the range of 0 - 1 in increments of 2^-10 + * exp(x/1024) = expFracTableA[x] + expFracTableB[x]. + */ + private static final double EXP_FRAC_TABLE_A[] = new double[1025]; + + /** Exponential over the range of 0 - 1 in increments of 2^-10 + * exp(x/1024) = expFracTableA[x] + expFracTableB[x]. + */ + private static final double EXP_FRAC_TABLE_B[] = new double[1025]; + + /** Factorial table, for Taylor series expansions. */ + private static final double FACT[] = new double[20]; + + /** Extended precision logarithm table over the range 1 - 2 in increments of 2^-10. */ + private static final double LN_MANT[][] = new double[1024][]; + + /** log(2) (high bits). */ + private static final double LN_2_A = 0.693147063255310059; + + /** log(2) (low bits). */ + private static final double LN_2_B = 1.17304635250823482e-7; + + /** Coefficients for slowLog. */ + private static final double LN_SPLIT_COEF[][] = { + {2.0, 0.0}, + {0.6666666269302368, 3.9736429850260626E-8}, + {0.3999999761581421, 2.3841857910019882E-8}, + {0.2857142686843872, 1.7029898543501842E-8}, + {0.2222222089767456, 1.3245471311735498E-8}, + {0.1818181574344635, 2.4384203044354907E-8}, + {0.1538461446762085, 9.140260083262505E-9}, + {0.13333332538604736, 9.220590270857665E-9}, + {0.11764700710773468, 1.2393345855018391E-8}, + {0.10526403784751892, 8.251545029714408E-9}, + {0.0952233225107193, 1.2675934823758863E-8}, + {0.08713622391223907, 1.1430250008909141E-8}, + {0.07842259109020233, 2.404307984052299E-9}, + {0.08371849358081818, 1.176342548272881E-8}, + {0.030589580535888672, 1.2958646899018938E-9}, + {0.14982303977012634, 1.225743062930824E-8}, + }; + + /** Coefficients for log, when input 0.99 < x < 1.01. */ + private static final double LN_QUICK_COEF[][] = { + {1.0, 5.669184079525E-24}, + {-0.25, -0.25}, + {0.3333333134651184, 1.986821492305628E-8}, + {-0.25, -6.663542893624021E-14}, + {0.19999998807907104, 1.1921056801463227E-8}, + {-0.1666666567325592, -7.800414592973399E-9}, + {0.1428571343421936, 5.650007086920087E-9}, + {-0.12502530217170715, -7.44321345601866E-11}, + {0.11113807559013367, 9.219544613762692E-9}, + }; + + /** Coefficients for log in the range of 1.0 < x < 1.0 + 2^-10. */ + private static final double LN_HI_PREC_COEF[][] = { + {1.0, -6.032174644509064E-23}, + {-0.25, -0.25}, + {0.3333333134651184, 1.9868161777724352E-8}, + {-0.2499999701976776, -2.957007209750105E-8}, + {0.19999954104423523, 1.5830993332061267E-10}, + {-0.16624879837036133, -2.6033824355191673E-8} + }; + + /** Sine table (high bits). */ + private static final double SINE_TABLE_A[] = new double[14]; + + /** Sine table (low bits). */ + private static final double SINE_TABLE_B[] = new double[14]; + + /** Cosine table (high bits). */ + private static final double COSINE_TABLE_A[] = new double[14]; + + /** Cosine table (low bits). */ + private static final double COSINE_TABLE_B[] = new double[14]; + + /** Tangent table, used by atan() (high bits). */ + private static final double TANGENT_TABLE_A[] = new double[14]; + + /** Tangent table, used by atan() (low bits). */ + private static final double TANGENT_TABLE_B[] = new double[14]; + + /** Bits of 1/(2*pi), need for reducePayneHanek(). */ + private static final long RECIP_2PI[] = new long[] { + (0x28be60dbL << 32) | 0x9391054aL, + (0x7f09d5f4L << 32) | 0x7d4d3770L, + (0x36d8a566L << 32) | 0x4f10e410L, + (0x7f9458eaL << 32) | 0xf7aef158L, + (0x6dc91b8eL << 32) | 0x909374b8L, + (0x01924bbaL << 32) | 0x82746487L, + (0x3f877ac7L << 32) | 0x2c4a69cfL, + (0xba208d7dL << 32) | 0x4baed121L, + (0x3a671c09L << 32) | 0xad17df90L, + (0x4e64758eL << 32) | 0x60d4ce7dL, + (0x272117e2L << 32) | 0xef7e4a0eL, + (0xc7fe25ffL << 32) | 0xf7816603L, + (0xfbcbc462L << 32) | 0xd6829b47L, + (0xdb4d9fb3L << 32) | 0xc9f2c26dL, + (0xd3d18fd9L << 32) | 0xa797fa8bL, + (0x5d49eeb1L << 32) | 0xfaf97c5eL, + (0xcf41ce7dL << 32) | 0xe294a4baL, + 0x9afed7ecL << 32 }; + + /** Bits of pi/4, need for reducePayneHanek(). */ + private static final long PI_O_4_BITS[] = new long[] { + (0xc90fdaa2L << 32) | 0x2168c234L, + (0xc4c6628bL << 32) | 0x80dc1cd1L }; + + /** Eighths. + * This is used by sinQ, because its faster to do a table lookup than + * a multiply in this time-critical routine + */ + private static final double EIGHTHS[] = {0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625}; + + /** Table of 2^((n+2)/3) */ + private static final double CBRTTWO[] = { 0.6299605249474366, + 0.7937005259840998, + 1.0, + 1.2599210498948732, + 1.5874010519681994 }; + + /* + * There are 52 bits in the mantissa of a double. + * For additional precision, the code splits double numbers into two parts, + * by clearing the low order 30 bits if possible, and then performs the arithmetic + * on each half separately. + */ + + /** + * 0x40000000 - used to split a double into two parts, both with the low order bits cleared. + * Equivalent to 2^30. + */ + private static final long HEX_40000000 = 0x40000000L; // 1073741824L + + /** Mask used to clear low order 30 bits */ + private static final long MASK_30BITS = -1L - (HEX_40000000 -1); // 0xFFFFFFFFC0000000L; + + /** 2^52 - double numbers this large must be integral (no fraction) or NaN or Infinite */ + private static final double TWO_POWER_52 = 4503599627370496.0; + + // Initialize tables + static { + int i; + + // Generate an array of factorials + FACT[0] = 1.0; + for (i = 1; i < FACT.length; i++) { + FACT[i] = FACT[i-1] * i; + } + + double tmp[] = new double[2]; + double recip[] = new double[2]; + + // Populate expIntTable + for (i = 0; i < 750; i++) { + expint(i, tmp); + EXP_INT_TABLE_A[i+750] = tmp[0]; + EXP_INT_TABLE_B[i+750] = tmp[1]; + + if (i != 0) { + // Negative integer powers + splitReciprocal(tmp, recip); + EXP_INT_TABLE_A[750-i] = recip[0]; + EXP_INT_TABLE_B[750-i] = recip[1]; + } + } + + // Populate expFracTable + for (i = 0; i < EXP_FRAC_TABLE_A.length; i++) { + slowexp(i/1024.0, tmp); + EXP_FRAC_TABLE_A[i] = tmp[0]; + EXP_FRAC_TABLE_B[i] = tmp[1]; + } + + // Populate lnMant table + for (i = 0; i < LN_MANT.length; i++) { + double d = Double.longBitsToDouble( (((long) i) << 42) | 0x3ff0000000000000L ); + LN_MANT[i] = slowLog(d); + } + + // Build the sine and cosine tables + buildSinCosTables(); + } + + /** + * Private Constructor + */ + private FastMath() { + } + + // Generic helper methods + + /** + * Get the high order bits from the mantissa. + * Equivalent to adding and subtracting HEX_40000 but also works for very large numbers + * + * @param d the value to split + * @return the high order part of the mantissa + */ + private static double doubleHighPart(double d) { + if (d > -MathUtils.SAFE_MIN && d < MathUtils.SAFE_MIN){ + return d; // These are un-normalised - don't try to convert + } + long xl = Double.doubleToLongBits(d); + xl = xl & MASK_30BITS; // Drop low order bits + return Double.longBitsToDouble(xl); + } + + /** Compute the square root of a number. + * <p><b>Note:</b> this implementation currently delegates to {@link Math#sqrt} + * @param a number on which evaluation is done + * @return square root of a + */ + public static double sqrt(final double a) { + return Math.sqrt(a); + } + + /** Compute the hyperbolic cosine of a number. + * @param x number on which evaluation is done + * @return hyperbolic cosine of x + */ + public static double cosh(double x) { + if (x != x) { + return x; + } + + if (x > 20.0) { + return exp(x)/2.0; + } + + if (x < -20) { + return exp(-x)/2.0; + } + + double hiPrec[] = new double[2]; + if (x < 0.0) { + x = -x; + } + exp(x, 0.0, hiPrec); + + double ya = hiPrec[0] + hiPrec[1]; + double yb = -(ya - hiPrec[0] - hiPrec[1]); + + double temp = ya * HEX_40000000; + double yaa = ya + temp - temp; + double yab = ya - yaa; + + // recip = 1/y + double recip = 1.0/ya; + temp = recip * HEX_40000000; + double recipa = recip + temp - temp; + double recipb = recip - recipa; + + // Correct for rounding in division + recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip; + // Account for yb + recipb += -yb * recip * recip; + + // y = y + 1/y + temp = ya + recipa; + yb += -(temp - ya - recipa); + ya = temp; + temp = ya + recipb; + yb += -(temp - ya - recipb); + ya = temp; + + double result = ya + yb; + result *= 0.5; + return result; + } + + /** Compute the hyperbolic sine of a number. + * @param x number on which evaluation is done + * @return hyperbolic sine of x + */ + public static double sinh(double x) { + boolean negate = false; + if (x != x) { + return x; + } + + if (x > 20.0) { + return exp(x)/2.0; + } + + if (x < -20) { + return -exp(-x)/2.0; + } + + if (x == 0) { + return x; + } + + if (x < 0.0) { + x = -x; + negate = true; + } + + double result; + + if (x > 0.25) { + double hiPrec[] = new double[2]; + exp(x, 0.0, hiPrec); + + double ya = hiPrec[0] + hiPrec[1]; + double yb = -(ya - hiPrec[0] - hiPrec[1]); + + double temp = ya * HEX_40000000; + double yaa = ya + temp - temp; + double yab = ya - yaa; + + // recip = 1/y + double recip = 1.0/ya; + temp = recip * HEX_40000000; + double recipa = recip + temp - temp; + double recipb = recip - recipa; + + // Correct for rounding in division + recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip; + // Account for yb + recipb += -yb * recip * recip; + + recipa = -recipa; + recipb = -recipb; + + // y = y + 1/y + temp = ya + recipa; + yb += -(temp - ya - recipa); + ya = temp; + temp = ya + recipb; + yb += -(temp - ya - recipb); + ya = temp; + + result = ya + yb; + result *= 0.5; + } + else { + double hiPrec[] = new double[2]; + expm1(x, hiPrec); + + double ya = hiPrec[0] + hiPrec[1]; + double yb = -(ya - hiPrec[0] - hiPrec[1]); + + /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */ + double denom = 1.0 + ya; + double denomr = 1.0 / denom; + double denomb = -(denom - 1.0 - ya) + yb; + double ratio = ya * denomr; + double temp = ratio * HEX_40000000; + double ra = ratio + temp - temp; + double rb = ratio - ra; + + temp = denom * HEX_40000000; + double za = denom + temp - temp; + double zb = denom - za; + + rb += (ya - za*ra - za*rb - zb*ra - zb*rb) * denomr; + + // Adjust for yb + rb += yb*denomr; // numerator + rb += -ya * denomb * denomr * denomr; // denominator + + // y = y - 1/y + temp = ya + ra; + yb += -(temp - ya - ra); + ya = temp; + temp = ya + rb; + yb += -(temp - ya - rb); + ya = temp; + + result = ya + yb; + result *= 0.5; + } + + if (negate) { + result = -result; + } + + return result; + } + + /** Compute the hyperbolic tangent of a number. + * @param x number on which evaluation is done + * @return hyperbolic tangent of x + */ + public static double tanh(double x) { + boolean negate = false; + + if (x != x) { + return x; + } + + if (x > 20.0) { + return 1.0; + } + + if (x < -20) { + return -1.0; + } + + if (x == 0) { + return x; + } + + if (x < 0.0) { + x = -x; + negate = true; + } + + double result; + if (x >= 0.5) { + double hiPrec[] = new double[2]; + // tanh(x) = (exp(2x) - 1) / (exp(2x) + 1) + exp(x*2.0, 0.0, hiPrec); + + double ya = hiPrec[0] + hiPrec[1]; + double yb = -(ya - hiPrec[0] - hiPrec[1]); + + /* Numerator */ + double na = -1.0 + ya; + double nb = -(na + 1.0 - ya); + double temp = na + yb; + nb += -(temp - na - yb); + na = temp; + + /* Denominator */ + double da = 1.0 + ya; + double db = -(da - 1.0 - ya); + temp = da + yb; + db += -(temp - da - yb); + da = temp; + + temp = da * HEX_40000000; + double daa = da + temp - temp; + double dab = da - daa; + + // ratio = na/da + double ratio = na/da; + temp = ratio * HEX_40000000; + double ratioa = ratio + temp - temp; + double ratiob = ratio - ratioa; + + // Correct for rounding in division + ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da; + + // Account for nb + ratiob += nb / da; + // Account for db + ratiob += -db * na / da / da; + + result = ratioa + ratiob; + } + else { + double hiPrec[] = new double[2]; + // tanh(x) = expm1(2x) / (expm1(2x) + 2) + expm1(x*2.0, hiPrec); + + double ya = hiPrec[0] + hiPrec[1]; + double yb = -(ya - hiPrec[0] - hiPrec[1]); + + /* Numerator */ + double na = ya; + double nb = yb; + + /* Denominator */ + double da = 2.0 + ya; + double db = -(da - 2.0 - ya); + double temp = da + yb; + db += -(temp - da - yb); + da = temp; + + temp = da * HEX_40000000; + double daa = da + temp - temp; + double dab = da - daa; + + // ratio = na/da + double ratio = na/da; + temp = ratio * HEX_40000000; + double ratioa = ratio + temp - temp; + double ratiob = ratio - ratioa; + + // Correct for rounding in division + ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da; + + // Account for nb + ratiob += nb / da; + // Account for db + ratiob += -db * na / da / da; + + result = ratioa + ratiob; + } + + if (negate) { + result = -result; + } + + return result; + } + + /** Compute the inverse hyperbolic cosine of a number. + * @param a number on which evaluation is done + * @return inverse hyperbolic cosine of a + */ + public static double acosh(final double a) { + return FastMath.log(a + FastMath.sqrt(a * a - 1)); + } + + /** Compute the inverse hyperbolic sine of a number. + * @param a number on which evaluation is done + * @return inverse hyperbolic sine of a + */ + public static double asinh(double a) { + + boolean negative = false; + if (a < 0) { + negative = true; + a = -a; + } + + double absAsinh; + if (a > 0.167) { + absAsinh = FastMath.log(FastMath.sqrt(a * a + 1) + a); + } else { + final double a2 = a * a; + if (a > 0.097) { + absAsinh = a * (1 - a2 * (1 / 3.0 - a2 * (1 / 5.0 - a2 * (1 / 7.0 - a2 * (1 / 9.0 - a2 * (1.0 / 11.0 - a2 * (1.0 / 13.0 - a2 * (1.0 / 15.0 - a2 * (1.0 / 17.0) * 15.0 / 16.0) * 13.0 / 14.0) * 11.0 / 12.0) * 9.0 / 10.0) * 7.0 / 8.0) * 5.0 / 6.0) * 3.0 / 4.0) / 2.0); + } else if (a > 0.036) { + absAsinh = a * (1 - a2 * (1 / 3.0 - a2 * (1 / 5.0 - a2 * (1 / 7.0 - a2 * (1 / 9.0 - a2 * (1.0 / 11.0 - a2 * (1.0 / 13.0) * 11.0 / 12.0) * 9.0 / 10.0) * 7.0 / 8.0) * 5.0 / 6.0) * 3.0 / 4.0) / 2.0); + } else if (a > 0.0036) { + absAsinh = a * (1 - a2 * (1 / 3.0 - a2 * (1 / 5.0 - a2 * (1 / 7.0 - a2 * (1 / 9.0) * 7.0 / 8.0) * 5.0 / 6.0) * 3.0 / 4.0) / 2.0); + } else { + absAsinh = a * (1 - a2 * (1 / 3.0 - a2 * (1 / 5.0) * 3.0 / 4.0) / 2.0); + } + } + + return negative ? -absAsinh : absAsinh; + + } + + /** Compute the inverse hyperbolic tangent of a number. + * @param a number on which evaluation is done + * @return inverse hyperbolic tangent of a + */ + public static double atanh(double a) { + + boolean negative = false; + if (a < 0) { + negative = true; + a = -a; + } + + double absAtanh; + if (a > 0.15) { + absAtanh = 0.5 * FastMath.log((1 + a) / (1 - a)); + } else { + final double a2 = a * a; + if (a > 0.087) { + absAtanh = a * (1 + a2 * (1.0 / 3.0 + a2 * (1.0 / 5.0 + a2 * (1.0 / 7.0 + a2 * (1.0 / 9.0 + a2 * (1.0 / 11.0 + a2 * (1.0 / 13.0 + a2 * (1.0 / 15.0 + a2 * (1.0 / 17.0))))))))); + } else if (a > 0.031) { + absAtanh = a * (1 + a2 * (1.0 / 3.0 + a2 * (1.0 / 5.0 + a2 * (1.0 / 7.0 + a2 * (1.0 / 9.0 + a2 * (1.0 / 11.0 + a2 * (1.0 / 13.0))))))); + } else if (a > 0.003) { + absAtanh = a * (1 + a2 * (1.0 / 3.0 + a2 * (1.0 / 5.0 + a2 * (1.0 / 7.0 + a2 * (1.0 / 9.0))))); + } else { + absAtanh = a * (1 + a2 * (1.0 / 3.0 + a2 * (1.0 / 5.0))); + } + } + + return negative ? -absAtanh : absAtanh; + + } + + /** Compute the signum of a number. + * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise + * @param a number on which evaluation is done + * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a + */ + public static double signum(final double a) { + return (a < 0.0) ? -1.0 : ((a > 0.0) ? 1.0 : a); // return +0.0/-0.0/NaN depending on a + } + + /** Compute the signum of a number. + * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise + * @param a number on which evaluation is done + * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a + */ + public static float signum(final float a) { + return (a < 0.0f) ? -1.0f : ((a > 0.0f) ? 1.0f : a); // return +0.0/-0.0/NaN depending on a + } + + /** Compute next number towards positive infinity. + * @param a number to which neighbor should be computed + * @return neighbor of a towards positive infinity + */ + public static double nextUp(final double a) { + return nextAfter(a, Double.POSITIVE_INFINITY); + } + + /** Compute next number towards positive infinity. + * @param a number to which neighbor should be computed + * @return neighbor of a towards positive infinity + */ + public static float nextUp(final float a) { + return nextAfter(a, Float.POSITIVE_INFINITY); + } + + /** Returns a pseudo-random number between 0.0 and 1.0. + * <p><b>Note:</b> this implementation currently delegates to {@link Math#random} + * @return a random number between 0.0 and 1.0 + */ + public static double random() { + return Math.random(); + } + + /** + * Exponential function. + * + * Computes exp(x), function result is nearly rounded. It will be correctly + * rounded to the theoretical value for 99.9% of input values, otherwise it will + * have a 1 UPL error. + * + * Method: + * Lookup intVal = exp(int(x)) + * Lookup fracVal = exp(int(x-int(x) / 1024.0) * 1024.0 ); + * Compute z as the exponential of the remaining bits by a polynomial minus one + * exp(x) = intVal * fracVal * (1 + z) + * + * Accuracy: + * Calculation is done with 63 bits of precision, so result should be correctly + * rounded for 99.9% of input values, with less than 1 ULP error otherwise. + * + * @param x a double + * @return double e<sup>x</sup> + */ + public static double exp(double x) { + return exp(x, 0.0, null); + } + + /** + * Internal helper method for exponential function. + * @param x original argument of the exponential function + * @param extra extra bits of precision on input (To Be Confirmed) + * @param hiPrec extra bits of precision on output (To Be Confirmed) + * @return exp(x) + */ + private static double exp(double x, double extra, double[] hiPrec) { + double intPartA; + double intPartB; + int intVal; + + /* Lookup exp(floor(x)). + * intPartA will have the upper 22 bits, intPartB will have the lower + * 52 bits. + */ + if (x < 0.0) { + intVal = (int) -x; + + if (intVal > 746) { + if (hiPrec != null) { + hiPrec[0] = 0.0; + hiPrec[1] = 0.0; + } + return 0.0; + } + + if (intVal > 709) { + /* This will produce a subnormal output */ + final double result = exp(x+40.19140625, extra, hiPrec) / 285040095144011776.0; + if (hiPrec != null) { + hiPrec[0] /= 285040095144011776.0; + hiPrec[1] /= 285040095144011776.0; + } + return result; + } + + if (intVal == 709) { + /* exp(1.494140625) is nearly a machine number... */ + final double result = exp(x+1.494140625, extra, hiPrec) / 4.455505956692756620; + if (hiPrec != null) { + hiPrec[0] /= 4.455505956692756620; + hiPrec[1] /= 4.455505956692756620; + } + return result; + } + + intVal++; + + intPartA = EXP_INT_TABLE_A[750-intVal]; + intPartB = EXP_INT_TABLE_B[750-intVal]; + + intVal = -intVal; + } else { + intVal = (int) x; + + if (intVal > 709) { + if (hiPrec != null) { + hiPrec[0] = Double.POSITIVE_INFINITY; + hiPrec[1] = 0.0; + } + return Double.POSITIVE_INFINITY; + } + + intPartA = EXP_INT_TABLE_A[750+intVal]; + intPartB = EXP_INT_TABLE_B[750+intVal]; + } + + /* Get the fractional part of x, find the greatest multiple of 2^-10 less than + * x and look up the exp function of it. + * fracPartA will have the upper 22 bits, fracPartB the lower 52 bits. + */ + final int intFrac = (int) ((x - intVal) * 1024.0); + final double fracPartA = EXP_FRAC_TABLE_A[intFrac]; + final double fracPartB = EXP_FRAC_TABLE_B[intFrac]; + + /* epsilon is the difference in x from the nearest multiple of 2^-10. It + * has a value in the range 0 <= epsilon < 2^-10. + * Do the subtraction from x as the last step to avoid possible loss of percison. + */ + final double epsilon = x - (intVal + intFrac / 1024.0); + + /* Compute z = exp(epsilon) - 1.0 via a minimax polynomial. z has + full double precision (52 bits). Since z < 2^-10, we will have + 62 bits of precision when combined with the contant 1. This will be + used in the last addition below to get proper rounding. */ + + /* Remez generated polynomial. Converges on the interval [0, 2^-10], error + is less than 0.5 ULP */ + double z = 0.04168701738764507; + z = z * epsilon + 0.1666666505023083; + z = z * epsilon + 0.5000000000042687; + z = z * epsilon + 1.0; + z = z * epsilon + -3.940510424527919E-20; + + /* Compute (intPartA+intPartB) * (fracPartA+fracPartB) by binomial + expansion. + tempA is exact since intPartA and intPartB only have 22 bits each. + tempB will have 52 bits of precision. + */ + double tempA = intPartA * fracPartA; + double tempB = intPartA * fracPartB + intPartB * fracPartA + intPartB * fracPartB; + + /* Compute the result. (1+z)(tempA+tempB). Order of operations is + important. For accuracy add by increasing size. tempA is exact and + much larger than the others. If there are extra bits specified from the + pow() function, use them. */ + final double tempC = tempB + tempA; + final double result; + if (extra != 0.0) { + result = tempC*extra*z + tempC*extra + tempC*z + tempB + tempA; + } else { + result = tempC*z + tempB + tempA; + } + + if (hiPrec != null) { + // If requesting high precision + hiPrec[0] = tempA; + hiPrec[1] = tempC*extra*z + tempC*extra + tempC*z + tempB; + } + + return result; + } + + /** Compute exp(x) - 1 + * @param x number to compute shifted exponential + * @return exp(x) - 1 + */ + public static double expm1(double x) { + return expm1(x, null); + } + + /** Internal helper method for expm1 + * @param x number to compute shifted exponential + * @param hiPrecOut receive high precision result for -1.0 < x < 1.0 + * @return exp(x) - 1 + */ + private static double expm1(double x, double hiPrecOut[]) { + if (x != x || x == 0.0) { // NaN or zero + return x; + } + + if (x <= -1.0 || x >= 1.0) { + // If not between +/- 1.0 + //return exp(x) - 1.0; + double hiPrec[] = new double[2]; + exp(x, 0.0, hiPrec); + if (x > 0.0) { + return -1.0 + hiPrec[0] + hiPrec[1]; + } else { + final double ra = -1.0 + hiPrec[0]; + double rb = -(ra + 1.0 - hiPrec[0]); + rb += hiPrec[1]; + return ra + rb; + } + } + + double baseA; + double baseB; + double epsilon; + boolean negative = false; + + if (x < 0.0) { + x = -x; + negative = true; + } + + { + int intFrac = (int) (x * 1024.0); + double tempA = EXP_FRAC_TABLE_A[intFrac] - 1.0; + double tempB = EXP_FRAC_TABLE_B[intFrac]; + + double temp = tempA + tempB; + tempB = -(temp - tempA - tempB); + tempA = temp; + + temp = tempA * HEX_40000000; + baseA = tempA + temp - temp; + baseB = tempB + (tempA - baseA); + + epsilon = x - intFrac/1024.0; + } + + + /* Compute expm1(epsilon) */ + double zb = 0.008336750013465571; + zb = zb * epsilon + 0.041666663879186654; + zb = zb * epsilon + 0.16666666666745392; + zb = zb * epsilon + 0.49999999999999994; + zb = zb * epsilon; + zb = zb * epsilon; + + double za = epsilon; + double temp = za + zb; + zb = -(temp - za - zb); + za = temp; + + temp = za * HEX_40000000; + temp = za + temp - temp; + zb += za - temp; + za = temp; + + /* Combine the parts. expm1(a+b) = expm1(a) + expm1(b) + expm1(a)*expm1(b) */ + double ya = za * baseA; + //double yb = za*baseB + zb*baseA + zb*baseB; + temp = ya + za * baseB; + double yb = -(temp - ya - za * baseB); + ya = temp; + + temp = ya + zb * baseA; + yb += -(temp - ya - zb * baseA); + ya = temp; + + temp = ya + zb * baseB; + yb += -(temp - ya - zb*baseB); + ya = temp; + + //ya = ya + za + baseA; + //yb = yb + zb + baseB; + temp = ya + baseA; + yb += -(temp - baseA - ya); + ya = temp; + + temp = ya + za; + //yb += (ya > za) ? -(temp - ya - za) : -(temp - za - ya); + yb += -(temp - ya - za); + ya = temp; + + temp = ya + baseB; + //yb += (ya > baseB) ? -(temp - ya - baseB) : -(temp - baseB - ya); + yb += -(temp - ya - baseB); + ya = temp; + + temp = ya + zb; + //yb += (ya > zb) ? -(temp - ya - zb) : -(temp - zb - ya); + yb += -(temp - ya - zb); + ya = temp; + + if (negative) { + /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */ + double denom = 1.0 + ya; + double denomr = 1.0 / denom; + double denomb = -(denom - 1.0 - ya) + yb; + double ratio = ya * denomr; + temp = ratio * HEX_40000000; + final double ra = ratio + temp - temp; + double rb = ratio - ra; + + temp = denom * HEX_40000000; + za = denom + temp - temp; + zb = denom - za; + + rb += (ya - za * ra - za * rb - zb * ra - zb * rb) * denomr; + + // f(x) = x/1+x + // Compute f'(x) + // Product rule: d(uv) = du*v + u*dv + // Chain rule: d(f(g(x)) = f'(g(x))*f(g'(x)) + // d(1/x) = -1/(x*x) + // d(1/1+x) = -1/( (1+x)^2) * 1 = -1/((1+x)*(1+x)) + // d(x/1+x) = -x/((1+x)(1+x)) + 1/1+x = 1 / ((1+x)(1+x)) + + // Adjust for yb + rb += yb * denomr; // numerator + rb += -ya * denomb * denomr * denomr; // denominator + + // negate + ya = -ra; + yb = -rb; + } + + if (hiPrecOut != null) { + hiPrecOut[0] = ya; + hiPrecOut[1] = yb; + } + + return ya + yb; + } + + /** + * For x between 0 and 1, returns exp(x), uses extended precision + * @param x argument of exponential + * @param result placeholder where to place exp(x) split in two terms + * for extra precision (i.e. exp(x) = result[0] ° result[1] + * @return exp(x) + */ + private static double slowexp(final double x, final double result[]) { + final double xs[] = new double[2]; + final double ys[] = new double[2]; + final double facts[] = new double[2]; + final double as[] = new double[2]; + split(x, xs); + ys[0] = ys[1] = 0.0; + + for (int i = 19; i >= 0; i--) { + splitMult(xs, ys, as); + ys[0] = as[0]; + ys[1] = as[1]; + + split(FACT[i], as); + splitReciprocal(as, facts); + + splitAdd(ys, facts, as); + ys[0] = as[0]; + ys[1] = as[1]; + } + + if (result != null) { + result[0] = ys[0]; + result[1] = ys[1]; + } + + return ys[0] + ys[1]; + } + + /** Compute split[0], split[1] such that their sum is equal to d, + * and split[0] has its 30 least significant bits as zero. + * @param d number to split + * @param split placeholder where to place the result + */ + private static void split(final double d, final double split[]) { + if (d < 8e298 && d > -8e298) { + final double a = d * HEX_40000000; + split[0] = (d + a) - a; + split[1] = d - split[0]; + } else { + final double a = d * 9.31322574615478515625E-10; + split[0] = (d + a - d) * HEX_40000000; + split[1] = d - split[0]; + } + } + + /** Recompute a split. + * @param a input/out array containing the split, changed + * on output + */ + private static void resplit(final double a[]) { + final double c = a[0] + a[1]; + final double d = -(c - a[0] - a[1]); + + if (c < 8e298 && c > -8e298) { + double z = c * HEX_40000000; + a[0] = (c + z) - z; + a[1] = c - a[0] + d; + } else { + double z = c * 9.31322574615478515625E-10; + a[0] = (c + z - c) * HEX_40000000; + a[1] = c - a[0] + d; + } + } + + /** Multiply two numbers in split form. + * @param a first term of multiplication + * @param b second term of multiplication + * @param ans placeholder where to put the result + */ + private static void splitMult(double a[], double b[], double ans[]) { + ans[0] = a[0] * b[0]; + ans[1] = a[0] * b[1] + a[1] * b[0] + a[1] * b[1]; + + /* Resplit */ + resplit(ans); + } + + /** Add two numbers in split form. + * @param a first term of addition + * @param b second term of addition + * @param ans placeholder where to put the result + */ + private static void splitAdd(final double a[], final double b[], final double ans[]) { + ans[0] = a[0] + b[0]; + ans[1] = a[1] + b[1]; + + resplit(ans); + } + + /** Compute the reciprocal of in. Use the following algorithm. + * in = c + d. + * want to find x + y such that x+y = 1/(c+d) and x is much + * larger than y and x has several zero bits on the right. + * + * Set b = 1/(2^22), a = 1 - b. Thus (a+b) = 1. + * Use following identity to compute (a+b)/(c+d) + * + * (a+b)/(c+d) = a/c + (bc - ad) / (c^2 + cd) + * set x = a/c and y = (bc - ad) / (c^2 + cd) + * This will be close to the right answer, but there will be + * some rounding in the calculation of X. So by carefully + * computing 1 - (c+d)(x+y) we can compute an error and + * add that back in. This is done carefully so that terms + * of similar size are subtracted first. + * @param in initial number, in split form + * @param result placeholder where to put the result + */ + private static void splitReciprocal(final double in[], final double result[]) { + final double b = 1.0/4194304.0; + final double a = 1.0 - b; + + if (in[0] == 0.0) { + in[0] = in[1]; + in[1] = 0.0; + } + + result[0] = a / in[0]; + result[1] = (b*in[0]-a*in[1]) / (in[0]*in[0] + in[0]*in[1]); + + if (result[1] != result[1]) { // can happen if result[1] is NAN + result[1] = 0.0; + } + + /* Resplit */ + resplit(result); + + for (int i = 0; i < 2; i++) { + /* this may be overkill, probably once is enough */ + double err = 1.0 - result[0] * in[0] - result[0] * in[1] - + result[1] * in[0] - result[1] * in[1]; + /*err = 1.0 - err; */ + err = err * (result[0] + result[1]); + /*printf("err = %16e\n", err); */ + result[1] += err; + } + } + + /** Compute (a[0] + a[1]) * (b[0] + b[1]) in extended precision. + * @param a first term of the multiplication + * @param b second term of the multiplication + * @param result placeholder where to put the result + */ + private static void quadMult(final double a[], final double b[], final double result[]) { + final double xs[] = new double[2]; + final double ys[] = new double[2]; + final double zs[] = new double[2]; + + /* a[0] * b[0] */ + split(a[0], xs); + split(b[0], ys); + splitMult(xs, ys, zs); + + result[0] = zs[0]; + result[1] = zs[1]; + + /* a[0] * b[1] */ + split(b[1], ys); + splitMult(xs, ys, zs); + + double tmp = result[0] + zs[0]; + result[1] = result[1] - (tmp - result[0] - zs[0]); + result[0] = tmp; + tmp = result[0] + zs[1]; + result[1] = result[1] - (tmp - result[0] - zs[1]); + result[0] = tmp; + + /* a[1] * b[0] */ + split(a[1], xs); + split(b[0], ys); + splitMult(xs, ys, zs); + + tmp = result[0] + zs[0]; + result[1] = result[1] - (tmp - result[0] - zs[0]); + result[0] = tmp; + tmp = result[0] + zs[1]; + result[1] = result[1] - (tmp - result[0] - zs[1]); + result[0] = tmp; + + /* a[1] * b[0] */ + split(a[1], xs); + split(b[1], ys); + splitMult(xs, ys, zs); + + tmp = result[0] + zs[0]; + result[1] = result[1] - (tmp - result[0] - zs[0]); + result[0] = tmp; + tmp = result[0] + zs[1]; + result[1] = result[1] - (tmp - result[0] - zs[1]); + result[0] = tmp; + } + + /** Compute exp(p) for a integer p in extended precision. + * @param p integer whose exponential is requested + * @param result placeholder where to put the result in extended precision + * @return exp(p) in standard precision (equal to result[0] + result[1]) + */ + private static double expint(int p, final double result[]) { + //double x = M_E; + final double xs[] = new double[2]; + final double as[] = new double[2]; + final double ys[] = new double[2]; + //split(x, xs); + //xs[1] = (double)(2.7182818284590452353602874713526625L - xs[0]); + //xs[0] = 2.71827697753906250000; + //xs[1] = 4.85091998273542816811e-06; + //xs[0] = Double.longBitsToDouble(0x4005bf0800000000L); + //xs[1] = Double.longBitsToDouble(0x3ed458a2bb4a9b00L); + + /* E */ + xs[0] = 2.718281828459045; + xs[1] = 1.4456468917292502E-16; + + split(1.0, ys); + + while (p > 0) { + if ((p & 1) != 0) { + quadMult(ys, xs, as); + ys[0] = as[0]; ys[1] = as[1]; + } + + quadMult(xs, xs, as); + xs[0] = as[0]; xs[1] = as[1]; + + p >>= 1; + } + + if (result != null) { + result[0] = ys[0]; + result[1] = ys[1]; + + resplit(result); + } + + return ys[0] + ys[1]; + } + + + /** + * Natural logarithm. + * + * @param x a double + * @return log(x) + */ + public static double log(final double x) { + return log(x, null); + } + + /** + * Internal helper method for natural logarithm function. + * @param x original argument of the natural logarithm function + * @param hiPrec extra bits of precision on output (To Be Confirmed) + * @return log(x) + */ + private static double log(final double x, final double[] hiPrec) { + if (x==0) { // Handle special case of +0/-0 + return Double.NEGATIVE_INFINITY; + } + long bits = Double.doubleToLongBits(x); + + /* Handle special cases of negative input, and NaN */ + if ((bits & 0x8000000000000000L) != 0 || x != x) { + if (x != 0.0) { + if (hiPrec != null) { + hiPrec[0] = Double.NaN; + } + + return Double.NaN; + } + } + + /* Handle special cases of Positive infinity. */ + if (x == Double.POSITIVE_INFINITY) { + if (hiPrec != null) { + hiPrec[0] = Double.POSITIVE_INFINITY; + } + + return Double.POSITIVE_INFINITY; + } + + /* Extract the exponent */ + int exp = (int)(bits >> 52)-1023; + + if ((bits & 0x7ff0000000000000L) == 0) { + // Subnormal! + if (x == 0) { + // Zero + if (hiPrec != null) { + hiPrec[0] = Double.NEGATIVE_INFINITY; + } + + return Double.NEGATIVE_INFINITY; + } + + /* Normalize the subnormal number. */ + bits <<= 1; + while ( (bits & 0x0010000000000000L) == 0) { + exp--; + bits <<= 1; + } + } + + + if (exp == -1 || exp == 0) { + if (x < 1.01 && x > 0.99 && hiPrec == null) { + /* The normal method doesn't work well in the range [0.99, 1.01], so call do a straight + polynomial expansion in higer precision. */ + + /* Compute x - 1.0 and split it */ + double xa = x - 1.0; + double xb = xa - x + 1.0; + double tmp = xa * HEX_40000000; + double aa = xa + tmp - tmp; + double ab = xa - aa; + xa = aa; + xb = ab; + + double ya = LN_QUICK_COEF[LN_QUICK_COEF.length-1][0]; + double yb = LN_QUICK_COEF[LN_QUICK_COEF.length-1][1]; + + for (int i = LN_QUICK_COEF.length - 2; i >= 0; i--) { + /* Multiply a = y * x */ + aa = ya * xa; + ab = ya * xb + yb * xa + yb * xb; + /* split, so now y = a */ + tmp = aa * HEX_40000000; + ya = aa + tmp - tmp; + yb = aa - ya + ab; + + /* Add a = y + lnQuickCoef */ + aa = ya + LN_QUICK_COEF[i][0]; + ab = yb + LN_QUICK_COEF[i][1]; + /* Split y = a */ + tmp = aa * HEX_40000000; + ya = aa + tmp - tmp; + yb = aa - ya + ab; + } + + /* Multiply a = y * x */ + aa = ya * xa; + ab = ya * xb + yb * xa + yb * xb; + /* split, so now y = a */ + tmp = aa * HEX_40000000; + ya = aa + tmp - tmp; + yb = aa - ya + ab; + + return ya + yb; + } + } + + // lnm is a log of a number in the range of 1.0 - 2.0, so 0 <= lnm < ln(2) + double lnm[] = LN_MANT[(int)((bits & 0x000ffc0000000000L) >> 42)]; + + /* + double epsilon = x / Double.longBitsToDouble(bits & 0xfffffc0000000000L); + + epsilon -= 1.0; + */ + + // y is the most significant 10 bits of the mantissa + //double y = Double.longBitsToDouble(bits & 0xfffffc0000000000L); + //double epsilon = (x - y) / y; + double epsilon = (bits & 0x3ffffffffffL) / (TWO_POWER_52 + (bits & 0x000ffc0000000000L)); + + double lnza = 0.0; + double lnzb = 0.0; + + if (hiPrec != null) { + /* split epsilon -> x */ + double tmp = epsilon * HEX_40000000; + double aa = epsilon + tmp - tmp; + double ab = epsilon - aa; + double xa = aa; + double xb = ab; + + /* Need a more accurate epsilon, so adjust the division. */ + double numer = bits & 0x3ffffffffffL; + double denom = TWO_POWER_52 + (bits & 0x000ffc0000000000L); + aa = numer - xa*denom - xb * denom; + xb += aa / denom; + + /* Remez polynomial evaluation */ + double ya = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1][0]; + double yb = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1][1]; + + for (int i = LN_HI_PREC_COEF.length - 2; i >= 0; i--) { + /* Multiply a = y * x */ + aa = ya * xa; + ab = ya * xb + yb * xa + yb * xb; + /* split, so now y = a */ + tmp = aa * HEX_40000000; + ya = aa + tmp - tmp; + yb = aa - ya + ab; + + /* Add a = y + lnHiPrecCoef */ + aa = ya + LN_HI_PREC_COEF[i][0]; + ab = yb + LN_HI_PREC_COEF[i][1]; + /* Split y = a */ + tmp = aa * HEX_40000000; + ya = aa + tmp - tmp; + yb = aa - ya + ab; + } + + /* Multiply a = y * x */ + aa = ya * xa; + ab = ya * xb + yb * xa + yb * xb; + + /* split, so now lnz = a */ + /* + tmp = aa * 1073741824.0; + lnza = aa + tmp - tmp; + lnzb = aa - lnza + ab; + */ + lnza = aa + ab; + lnzb = -(lnza - aa - ab); + } else { + /* High precision not required. Eval Remez polynomial + using standard double precision */ + lnza = -0.16624882440418567; + lnza = lnza * epsilon + 0.19999954120254515; + lnza = lnza * epsilon + -0.2499999997677497; + lnza = lnza * epsilon + 0.3333333333332802; + lnza = lnza * epsilon + -0.5; + lnza = lnza * epsilon + 1.0; + lnza = lnza * epsilon; + } + + /* Relative sizes: + * lnzb [0, 2.33E-10] + * lnm[1] [0, 1.17E-7] + * ln2B*exp [0, 1.12E-4] + * lnza [0, 9.7E-4] + * lnm[0] [0, 0.692] + * ln2A*exp [0, 709] + */ + + /* Compute the following sum: + * lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp; + */ + + //return lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp; + double a = LN_2_A*exp; + double b = 0.0; + double c = a+lnm[0]; + double d = -(c-a-lnm[0]); + a = c; + b = b + d; + + c = a + lnza; + d = -(c - a - lnza); + a = c; + b = b + d; + + c = a + LN_2_B*exp; + d = -(c - a - LN_2_B*exp); + a = c; + b = b + d; + + c = a + lnm[1]; + d = -(c - a - lnm[1]); + a = c; + b = b + d; + + c = a + lnzb; + d = -(c - a - lnzb); + a = c; + b = b + d; + + if (hiPrec != null) { + hiPrec[0] = a; + hiPrec[1] = b; + } + + return a + b; + } + + /** Compute log(1 + x). + * @param x a number + * @return log(1 + x) + */ + public static double log1p(final double x) { + double xpa = 1.0 + x; + double xpb = -(xpa - 1.0 - x); + + if (x == -1) { + return x/0.0; // -Infinity + } + + if (x > 0 && 1/x == 0) { // x = Infinity + return x; + } + + if (x>1e-6 || x<-1e-6) { + double hiPrec[] = new double[2]; + + final double lores = log(xpa, hiPrec); + if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN + return lores; + } + + /* Do a taylor series expansion around xpa */ + /* f(x+y) = f(x) + f'(x)*y + f''(x)/2 y^2 */ + double fx1 = xpb/xpa; + + double epsilon = 0.5 * fx1 + 1.0; + epsilon = epsilon * fx1; + + return epsilon + hiPrec[1] + hiPrec[0]; + } + + /* Value is small |x| < 1e6, do a Taylor series centered on 1.0 */ + double y = x * 0.333333333333333 - 0.5; + y = y * x + 1.0; + y = y * x; + + return y; + } + + /** Compute the base 10 logarithm. + * @param x a number + * @return log10(x) + */ + public static double log10(final double x) { + final double hiPrec[] = new double[2]; + + final double lores = log(x, hiPrec); + if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN + return lores; + } + + final double tmp = hiPrec[0] * HEX_40000000; + final double lna = hiPrec[0] + tmp - tmp; + final double lnb = hiPrec[0] - lna + hiPrec[1]; + + final double rln10a = 0.4342944622039795; + final double rln10b = 1.9699272335463627E-8; + + return rln10b * lnb + rln10b * lna + rln10a * lnb + rln10a * lna; + } + + /** + * Power function. Compute x^y. + * + * @param x a double + * @param y a double + * @return double + */ + public static double pow(double x, double y) { + final double lns[] = new double[2]; + + if (y == 0.0) { + return 1.0; + } + + if (x != x) { // X is NaN + return x; + } + + + if (x == 0) { + long bits = Double.doubleToLongBits(x); + if ((bits & 0x8000000000000000L) != 0) { + // -zero + long yi = (long) y; + + if (y < 0 && y == yi && (yi & 1) == 1) { + return Double.NEGATIVE_INFINITY; + } + + if (y < 0 && y == yi && (yi & 1) == 1) { + return -0.0; + } + + if (y > 0 && y == yi && (yi & 1) == 1) { + return -0.0; + } + } + + if (y < 0) { + return Double.POSITIVE_INFINITY; + } + if (y > 0) { + return 0.0; + } + + return Double.NaN; + } + + if (x == Double.POSITIVE_INFINITY) { + if (y != y) { // y is NaN + return y; + } + if (y < 0.0) { + return 0.0; + } else { + return Double.POSITIVE_INFINITY; + } + } + + if (y == Double.POSITIVE_INFINITY) { + if (x * x == 1.0) + return Double.NaN; + + if (x * x > 1.0) { + return Double.POSITIVE_INFINITY; + } else { + return 0.0; + } + } + + if (x == Double.NEGATIVE_INFINITY) { + if (y != y) { // y is NaN + return y; + } + + if (y < 0) { + long yi = (long) y; + if (y == yi && (yi & 1) == 1) { + return -0.0; + } + + return 0.0; + } + + if (y > 0) { + long yi = (long) y; + if (y == yi && (yi & 1) == 1) { + return Double.NEGATIVE_INFINITY; + } + + return Double.POSITIVE_INFINITY; + } + } + + if (y == Double.NEGATIVE_INFINITY) { + + if (x * x == 1.0) { + return Double.NaN; + } + + if (x * x < 1.0) { + return Double.POSITIVE_INFINITY; + } else { + return 0.0; + } + } + + /* Handle special case x<0 */ + if (x < 0) { + // y is an even integer in this case + if (y >= TWO_POWER_52 || y <= -TWO_POWER_52) { + return pow(-x, y); + } + + if (y == (long) y) { + // If y is an integer + return ((long)y & 1) == 0 ? pow(-x, y) : -pow(-x, y); + } else { + return Double.NaN; + } + } + + /* Split y into ya and yb such that y = ya+yb */ + double ya; + double yb; + if (y < 8e298 && y > -8e298) { + double tmp1 = y * HEX_40000000; + ya = y + tmp1 - tmp1; + yb = y - ya; + } else { + double tmp1 = y * 9.31322574615478515625E-10; + double tmp2 = tmp1 * 9.31322574615478515625E-10; + ya = (tmp1 + tmp2 - tmp1) * HEX_40000000 * HEX_40000000; + yb = y - ya; + } + + /* Compute ln(x) */ + final double lores = log(x, lns); + if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN + return lores; + } + + double lna = lns[0]; + double lnb = lns[1]; + + /* resplit lns */ + double tmp1 = lna * HEX_40000000; + double tmp2 = lna + tmp1 - tmp1; + lnb += lna - tmp2; + lna = tmp2; + + // y*ln(x) = (aa+ab) + final double aa = lna * ya; + final double ab = lna * yb + lnb * ya + lnb * yb; + + lna = aa+ab; + lnb = -(lna - aa - ab); + + double z = 1.0 / 120.0; + z = z * lnb + (1.0 / 24.0); + z = z * lnb + (1.0 / 6.0); + z = z * lnb + 0.5; + z = z * lnb + 1.0; + z = z * lnb; + + final double result = exp(lna, z, null); + //result = result + result * z; + return result; + } + + /** xi in the range of [1, 2]. + * 3 5 7 + * x+1 / x x x \ + * ln ----- = 2 * | x + ---- + ---- + ---- + ... | + * 1-x \ 3 5 7 / + * + * So, compute a Remez approximation of the following function + * + * ln ((sqrt(x)+1)/(1-sqrt(x))) / x + * + * This will be an even function with only positive coefficents. + * x is in the range [0 - 1/3]. + * + * Transform xi for input to the above function by setting + * x = (xi-1)/(xi+1). Input to the polynomial is x^2, then + * the result is multiplied by x. + * @param xi number from which log is requested + * @return log(xi) + */ + private static double[] slowLog(double xi) { + double x[] = new double[2]; + double x2[] = new double[2]; + double y[] = new double[2]; + double a[] = new double[2]; + + split(xi, x); + + /* Set X = (x-1)/(x+1) */ + x[0] += 1.0; + resplit(x); + splitReciprocal(x, a); + x[0] -= 2.0; + resplit(x); + splitMult(x, a, y); + x[0] = y[0]; + x[1] = y[1]; + + /* Square X -> X2*/ + splitMult(x, x, x2); + + + //x[0] -= 1.0; + //resplit(x); + + y[0] = LN_SPLIT_COEF[LN_SPLIT_COEF.length-1][0]; + y[1] = LN_SPLIT_COEF[LN_SPLIT_COEF.length-1][1]; + + for (int i = LN_SPLIT_COEF.length-2; i >= 0; i--) { + splitMult(y, x2, a); + y[0] = a[0]; + y[1] = a[1]; + splitAdd(y, LN_SPLIT_COEF[i], a); + y[0] = a[0]; + y[1] = a[1]; + } + + splitMult(y, x, a); + y[0] = a[0]; + y[1] = a[1]; + + return y; + } + + /** + * For x between 0 and pi/4 compute sine. + * @param x number from which sine is requested + * @param result placeholder where to put the result in extended precision + * @return sin(x) + */ + private static double slowSin(final double x, final double result[]) { + final double xs[] = new double[2]; + final double ys[] = new double[2]; + final double facts[] = new double[2]; + final double as[] = new double[2]; + split(x, xs); + ys[0] = ys[1] = 0.0; + + for (int i = 19; i >= 0; i--) { + splitMult(xs, ys, as); + ys[0] = as[0]; ys[1] = as[1]; + + if ( (i & 1) == 0) { + continue; + } + + split(FACT[i], as); + splitReciprocal(as, facts); + + if ( (i & 2) != 0 ) { + facts[0] = -facts[0]; + facts[1] = -facts[1]; + } + + splitAdd(ys, facts, as); + ys[0] = as[0]; ys[1] = as[1]; + } + + if (result != null) { + result[0] = ys[0]; + result[1] = ys[1]; + } + + return ys[0] + ys[1]; + } + + /** + * For x between 0 and pi/4 compute cosine + * @param x number from which cosine is requested + * @param result placeholder where to put the result in extended precision + * @return cos(x) + */ + private static double slowCos(final double x, final double result[]) { + + final double xs[] = new double[2]; + final double ys[] = new double[2]; + final double facts[] = new double[2]; + final double as[] = new double[2]; + split(x, xs); + ys[0] = ys[1] = 0.0; + + for (int i = 19; i >= 0; i--) { + splitMult(xs, ys, as); + ys[0] = as[0]; ys[1] = as[1]; + + if ( (i & 1) != 0) { + continue; + } + + split(FACT[i], as); + splitReciprocal(as, facts); + + if ( (i & 2) != 0 ) { + facts[0] = -facts[0]; + facts[1] = -facts[1]; + } + + splitAdd(ys, facts, as); + ys[0] = as[0]; ys[1] = as[1]; + } + + if (result != null) { + result[0] = ys[0]; + result[1] = ys[1]; + } + + return ys[0] + ys[1]; + } + + /** Build the sine and cosine tables. + */ + private static void buildSinCosTables() { + final double result[] = new double[2]; + + /* Use taylor series for 0 <= x <= 6/8 */ + for (int i = 0; i < 7; i++) { + double x = i / 8.0; + + slowSin(x, result); + SINE_TABLE_A[i] = result[0]; + SINE_TABLE_B[i] = result[1]; + + slowCos(x, result); + COSINE_TABLE_A[i] = result[0]; + COSINE_TABLE_B[i] = result[1]; + } + + /* Use angle addition formula to complete table to 13/8, just beyond pi/2 */ + for (int i = 7; i < 14; i++) { + double xs[] = new double[2]; + double ys[] = new double[2]; + double as[] = new double[2]; + double bs[] = new double[2]; + double temps[] = new double[2]; + + if ( (i & 1) == 0) { + // Even, use double angle + xs[0] = SINE_TABLE_A[i/2]; + xs[1] = SINE_TABLE_B[i/2]; + ys[0] = COSINE_TABLE_A[i/2]; + ys[1] = COSINE_TABLE_B[i/2]; + + /* compute sine */ + splitMult(xs, ys, result); + SINE_TABLE_A[i] = result[0] * 2.0; + SINE_TABLE_B[i] = result[1] * 2.0; + + /* Compute cosine */ + splitMult(ys, ys, as); + splitMult(xs, xs, temps); + temps[0] = -temps[0]; + temps[1] = -temps[1]; + splitAdd(as, temps, result); + COSINE_TABLE_A[i] = result[0]; + COSINE_TABLE_B[i] = result[1]; + } else { + xs[0] = SINE_TABLE_A[i/2]; + xs[1] = SINE_TABLE_B[i/2]; + ys[0] = COSINE_TABLE_A[i/2]; + ys[1] = COSINE_TABLE_B[i/2]; + as[0] = SINE_TABLE_A[i/2+1]; + as[1] = SINE_TABLE_B[i/2+1]; + bs[0] = COSINE_TABLE_A[i/2+1]; + bs[1] = COSINE_TABLE_B[i/2+1]; + + /* compute sine */ + splitMult(xs, bs, temps); + splitMult(ys, as, result); + splitAdd(result, temps, result); + SINE_TABLE_A[i] = result[0]; + SINE_TABLE_B[i] = result[1]; + + /* Compute cosine */ + splitMult(ys, bs, result); + splitMult(xs, as, temps); + temps[0] = -temps[0]; + temps[1] = -temps[1]; + splitAdd(result, temps, result); + COSINE_TABLE_A[i] = result[0]; + COSINE_TABLE_B[i] = result[1]; + } + } + + /* Compute tangent = sine/cosine */ + for (int i = 0; i < 14; i++) { + double xs[] = new double[2]; + double ys[] = new double[2]; + double as[] = new double[2]; + + as[0] = COSINE_TABLE_A[i]; + as[1] = COSINE_TABLE_B[i]; + + splitReciprocal(as, ys); + + xs[0] = SINE_TABLE_A[i]; + xs[1] = SINE_TABLE_B[i]; + + splitMult(xs, ys, as); + + TANGENT_TABLE_A[i] = as[0]; + TANGENT_TABLE_B[i] = as[1]; + } + + } + + /** + * Computes sin(x) - x, where |x| < 1/16. + * Use a Remez polynomial approximation. + * @param x a number smaller than 1/16 + * @return sin(x) - x + */ + private static double polySine(final double x) + { + double x2 = x*x; + + double p = 2.7553817452272217E-6; + p = p * x2 + -1.9841269659586505E-4; + p = p * x2 + 0.008333333333329196; + p = p * x2 + -0.16666666666666666; + //p *= x2; + //p *= x; + p = p * x2 * x; + + return p; + } + + /** + * Computes cos(x) - 1, where |x| < 1/16. + * Use a Remez polynomial approximation. + * @param x a number smaller than 1/16 + * @return cos(x) - 1 + */ + private static double polyCosine(double x) { + double x2 = x*x; + + double p = 2.479773539153719E-5; + p = p * x2 + -0.0013888888689039883; + p = p * x2 + 0.041666666666621166; + p = p * x2 + -0.49999999999999994; + p *= x2; + + return p; + } + + /** + * Compute sine over the first quadrant (0 < x < pi/2). + * Use combination of table lookup and rational polynomial expansion. + * @param xa number from which sine is requested + * @param xb extra bits for x (may be 0.0) + * @return sin(xa + xb) + */ + private static double sinQ(double xa, double xb) { + int idx = (int) ((xa * 8.0) + 0.5); + final double epsilon = xa - EIGHTHS[idx]; //idx*0.125; + + // Table lookups + final double sintA = SINE_TABLE_A[idx]; + final double sintB = SINE_TABLE_B[idx]; + final double costA = COSINE_TABLE_A[idx]; + final double costB = COSINE_TABLE_B[idx]; + + // Polynomial eval of sin(epsilon), cos(epsilon) + double sinEpsA = epsilon; + double sinEpsB = polySine(epsilon); + final double cosEpsA = 1.0; + final double cosEpsB = polyCosine(epsilon); + + // Split epsilon xa + xb = x + final double temp = sinEpsA * HEX_40000000; + double temp2 = (sinEpsA + temp) - temp; + sinEpsB += sinEpsA - temp2; + sinEpsA = temp2; + + /* Compute sin(x) by angle addition formula */ + double result; + + /* Compute the following sum: + * + * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB + + * sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB; + * + * Ranges of elements + * + * xxxtA 0 PI/2 + * xxxtB -1.5e-9 1.5e-9 + * sinEpsA -0.0625 0.0625 + * sinEpsB -6e-11 6e-11 + * cosEpsA 1.0 + * cosEpsB 0 -0.0625 + * + */ + + //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB + + // sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB; + + //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB; + //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB; + double a = 0; + double b = 0; + + double t = sintA; + double c = a + t; + double d = -(c - a - t); + a = c; + b = b + d; + + t = costA * sinEpsA; + c = a + t; + d = -(c - a - t); + a = c; + b = b + d; + + b = b + sintA * cosEpsB + costA * sinEpsB; + /* + t = sintA*cosEpsB; + c = a + t; + d = -(c - a - t); + a = c; + b = b + d; + + t = costA*sinEpsB; + c = a + t; + d = -(c - a - t); + a = c; + b = b + d; + */ + + b = b + sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB; + /* + t = sintB; + c = a + t; + d = -(c - a - t); + a = c; + b = b + d; + + t = costB*sinEpsA; + c = a + t; + d = -(c - a - t); + a = c; + b = b + d; + + t = sintB*cosEpsB; + c = a + t; + d = -(c - a - t); + a = c; + b = b + d; + + t = costB*sinEpsB; + c = a + t; + d = -(c - a - t); + a = c; + b = b + d; + */ + + if (xb != 0.0) { + t = ((costA + costB) * (cosEpsA + cosEpsB) - + (sintA + sintB) * (sinEpsA + sinEpsB)) * xb; // approximate cosine*xb + c = a + t; + d = -(c - a - t); + a = c; + b = b + d; + } + + result = a + b; + + return result; + } + + /** + * Compute cosine in the first quadrant by subtracting input from PI/2 and + * then calling sinQ. This is more accurate as the input approaches PI/2. + * @param xa number from which cosine is requested + * @param xb extra bits for x (may be 0.0) + * @return cos(xa + xb) + */ + private static double cosQ(double xa, double xb) { + final double pi2a = 1.5707963267948966; + final double pi2b = 6.123233995736766E-17; + + final double a = pi2a - xa; + double b = -(a - pi2a + xa); + b += pi2b - xb; + + return sinQ(a, b); + } + + /** + * Compute tangent (or cotangent) over the first quadrant. 0 < x < pi/2 + * Use combination of table lookup and rational polynomial expansion. + * @param xa number from which sine is requested + * @param xb extra bits for x (may be 0.0) + * @param cotanFlag if true, compute the cotangent instead of the tangent + * @return tan(xa+xb) (or cotangent, depending on cotanFlag) + */ + private static double tanQ(double xa, double xb, boolean cotanFlag) { + + int idx = (int) ((xa * 8.0) + 0.5); + final double epsilon = xa - EIGHTHS[idx]; //idx*0.125; + + // Table lookups + final double sintA = SINE_TABLE_A[idx]; + final double sintB = SINE_TABLE_B[idx]; + final double costA = COSINE_TABLE_A[idx]; + final double costB = COSINE_TABLE_B[idx]; + + // Polynomial eval of sin(epsilon), cos(epsilon) + double sinEpsA = epsilon; + double sinEpsB = polySine(epsilon); + final double cosEpsA = 1.0; + final double cosEpsB = polyCosine(epsilon); + + // Split epsilon xa + xb = x + double temp = sinEpsA * HEX_40000000; + double temp2 = (sinEpsA + temp) - temp; + sinEpsB += sinEpsA - temp2; + sinEpsA = temp2; + + /* Compute sin(x) by angle addition formula */ + + /* Compute the following sum: + * + * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB + + * sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB; + * + * Ranges of elements + * + * xxxtA 0 PI/2 + * xxxtB -1.5e-9 1.5e-9 + * sinEpsA -0.0625 0.0625 + * sinEpsB -6e-11 6e-11 + * cosEpsA 1.0 + * cosEpsB 0 -0.0625 + * + */ + + //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB + + // sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB; + + //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB; + //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB; + double a = 0; + double b = 0; + + // Compute sine + double t = sintA; + double c = a + t; + double d = -(c - a - t); + a = c; + b = b + d; + + t = costA*sinEpsA; + c = a + t; + d = -(c - a - t); + a = c; + b = b + d; + + b = b + sintA*cosEpsB + costA*sinEpsB; + b = b + sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB; + + double sina = a + b; + double sinb = -(sina - a - b); + + // Compute cosine + + a = b = c = d = 0.0; + + t = costA*cosEpsA; + c = a + t; + d = -(c - a - t); + a = c; + b = b + d; + + t = -sintA*sinEpsA; + c = a + t; + d = -(c - a - t); + a = c; + b = b + d; + + b = b + costB*cosEpsA + costA*cosEpsB + costB*cosEpsB; + b = b - (sintB*sinEpsA + sintA*sinEpsB + sintB*sinEpsB); + + double cosa = a + b; + double cosb = -(cosa - a - b); + + if (cotanFlag) { + double tmp; + tmp = cosa; cosa = sina; sina = tmp; + tmp = cosb; cosb = sinb; sinb = tmp; + } + + + /* estimate and correct, compute 1.0/(cosa+cosb) */ + /* + double est = (sina+sinb)/(cosa+cosb); + double err = (sina - cosa*est) + (sinb - cosb*est); + est += err/(cosa+cosb); + err = (sina - cosa*est) + (sinb - cosb*est); + */ + + // f(x) = 1/x, f'(x) = -1/x^2 + + double est = sina/cosa; + + /* Split the estimate to get more accurate read on division rounding */ + temp = est * HEX_40000000; + double esta = (est + temp) - temp; + double estb = est - esta; + + temp = cosa * HEX_40000000; + double cosaa = (cosa + temp) - temp; + double cosab = cosa - cosaa; + + //double err = (sina - est*cosa)/cosa; // Correction for division rounding + double err = (sina - esta*cosaa - esta*cosab - estb*cosaa - estb*cosab)/cosa; // Correction for division rounding + err += sinb/cosa; // Change in est due to sinb + err += -sina * cosb / cosa / cosa; // Change in est due to cosb + + if (xb != 0.0) { + // tan' = 1 + tan^2 cot' = -(1 + cot^2) + // Approximate impact of xb + double xbadj = xb + est*est*xb; + if (cotanFlag) { + xbadj = -xbadj; + } + + err += xbadj; + } + + return est+err; + } + + /** Reduce the input argument using the Payne and Hanek method. + * This is good for all inputs 0.0 < x < inf + * Output is remainder after dividing by PI/2 + * The result array should contain 3 numbers. + * result[0] is the integer portion, so mod 4 this gives the quadrant. + * result[1] is the upper bits of the remainder + * result[2] is the lower bits of the remainder + * + * @param x number to reduce + * @param result placeholder where to put the result + */ + private static void reducePayneHanek(double x, double result[]) + { + /* Convert input double to bits */ + long inbits = Double.doubleToLongBits(x); + int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023; + + /* Convert to fixed point representation */ + inbits &= 0x000fffffffffffffL; + inbits |= 0x0010000000000000L; + + /* Normalize input to be between 0.5 and 1.0 */ + exponent++; + inbits <<= 11; + + /* Based on the exponent, get a shifted copy of recip2pi */ + long shpi0; + long shpiA; + long shpiB; + int idx = exponent >> 6; + int shift = exponent - (idx << 6); + + if (shift != 0) { + shpi0 = (idx == 0) ? 0 : (RECIP_2PI[idx-1] << shift); + shpi0 |= RECIP_2PI[idx] >>> (64-shift); + shpiA = (RECIP_2PI[idx] << shift) | (RECIP_2PI[idx+1] >>> (64-shift)); + shpiB = (RECIP_2PI[idx+1] << shift) | (RECIP_2PI[idx+2] >>> (64-shift)); + } else { + shpi0 = (idx == 0) ? 0 : RECIP_2PI[idx-1]; + shpiA = RECIP_2PI[idx]; + shpiB = RECIP_2PI[idx+1]; + } + + /* Multiply input by shpiA */ + long a = inbits >>> 32; + long b = inbits & 0xffffffffL; + + long c = shpiA >>> 32; + long d = shpiA & 0xffffffffL; + + long ac = a * c; + long bd = b * d; + long bc = b * c; + long ad = a * d; + + long prodB = bd + (ad << 32); + long prodA = ac + (ad >>> 32); + + boolean bita = (bd & 0x8000000000000000L) != 0; + boolean bitb = (ad & 0x80000000L ) != 0; + boolean bitsum = (prodB & 0x8000000000000000L) != 0; + + /* Carry */ + if ( (bita && bitb) || + ((bita || bitb) && !bitsum) ) { + prodA++; + } + + bita = (prodB & 0x8000000000000000L) != 0; + bitb = (bc & 0x80000000L ) != 0; + + prodB = prodB + (bc << 32); + prodA = prodA + (bc >>> 32); + + bitsum = (prodB & 0x8000000000000000L) != 0; + + /* Carry */ + if ( (bita && bitb) || + ((bita || bitb) && !bitsum) ) { + prodA++; + } + + /* Multiply input by shpiB */ + c = shpiB >>> 32; + d = shpiB & 0xffffffffL; + ac = a * c; + bc = b * c; + ad = a * d; + + /* Collect terms */ + ac = ac + ((bc + ad) >>> 32); + + bita = (prodB & 0x8000000000000000L) != 0; + bitb = (ac & 0x8000000000000000L ) != 0; + prodB += ac; + bitsum = (prodB & 0x8000000000000000L) != 0; + /* Carry */ + if ( (bita && bitb) || + ((bita || bitb) && !bitsum) ) { + prodA++; + } + + /* Multiply by shpi0 */ + c = shpi0 >>> 32; + d = shpi0 & 0xffffffffL; + + bd = b * d; + bc = b * c; + ad = a * d; + + prodA += bd + ((bc + ad) << 32); + + /* + * prodA, prodB now contain the remainder as a fraction of PI. We want this as a fraction of + * PI/2, so use the following steps: + * 1.) multiply by 4. + * 2.) do a fixed point muliply by PI/4. + * 3.) Convert to floating point. + * 4.) Multiply by 2 + */ + + /* This identifies the quadrant */ + int intPart = (int)(prodA >>> 62); + + /* Multiply by 4 */ + prodA <<= 2; + prodA |= prodB >>> 62; + prodB <<= 2; + + /* Multiply by PI/4 */ + a = prodA >>> 32; + b = prodA & 0xffffffffL; + + c = PI_O_4_BITS[0] >>> 32; + d = PI_O_4_BITS[0] & 0xffffffffL; + + ac = a * c; + bd = b * d; + bc = b * c; + ad = a * d; + + long prod2B = bd + (ad << 32); + long prod2A = ac + (ad >>> 32); + + bita = (bd & 0x8000000000000000L) != 0; + bitb = (ad & 0x80000000L ) != 0; + bitsum = (prod2B & 0x8000000000000000L) != 0; + + /* Carry */ + if ( (bita && bitb) || + ((bita || bitb) && !bitsum) ) { + prod2A++; + } + + bita = (prod2B & 0x8000000000000000L) != 0; + bitb = (bc & 0x80000000L ) != 0; + + prod2B = prod2B + (bc << 32); + prod2A = prod2A + (bc >>> 32); + + bitsum = (prod2B & 0x8000000000000000L) != 0; + + /* Carry */ + if ( (bita && bitb) || + ((bita || bitb) && !bitsum) ) { + prod2A++; + } + + /* Multiply input by pio4bits[1] */ + c = PI_O_4_BITS[1] >>> 32; + d = PI_O_4_BITS[1] & 0xffffffffL; + ac = a * c; + bc = b * c; + ad = a * d; + + /* Collect terms */ + ac = ac + ((bc + ad) >>> 32); + + bita = (prod2B & 0x8000000000000000L) != 0; + bitb = (ac & 0x8000000000000000L ) != 0; + prod2B += ac; + bitsum = (prod2B & 0x8000000000000000L) != 0; + /* Carry */ + if ( (bita && bitb) || + ((bita || bitb) && !bitsum) ) { + prod2A++; + } + + /* Multiply inputB by pio4bits[0] */ + a = prodB >>> 32; + b = prodB & 0xffffffffL; + c = PI_O_4_BITS[0] >>> 32; + d = PI_O_4_BITS[0] & 0xffffffffL; + ac = a * c; + bc = b * c; + ad = a * d; + + /* Collect terms */ + ac = ac + ((bc + ad) >>> 32); + + bita = (prod2B & 0x8000000000000000L) != 0; + bitb = (ac & 0x8000000000000000L ) != 0; + prod2B += ac; + bitsum = (prod2B & 0x8000000000000000L) != 0; + /* Carry */ + if ( (bita && bitb) || + ((bita || bitb) && !bitsum) ) { + prod2A++; + } + + /* Convert to double */ + double tmpA = (prod2A >>> 12) / TWO_POWER_52; // High order 52 bits + double tmpB = (((prod2A & 0xfffL) << 40) + (prod2B >>> 24)) / TWO_POWER_52 / TWO_POWER_52; // Low bits + + double sumA = tmpA + tmpB; + double sumB = -(sumA - tmpA - tmpB); + + /* Multiply by PI/2 and return */ + result[0] = intPart; + result[1] = sumA * 2.0; + result[2] = sumB * 2.0; + } + + /** + * Sine function. + * @param x a number + * @return sin(x) + */ + public static double sin(double x) { + boolean negative = false; + int quadrant = 0; + double xa; + double xb = 0.0; + + /* Take absolute value of the input */ + xa = x; + if (x < 0) { + negative = true; + xa = -xa; + } + + /* Check for zero and negative zero */ + if (xa == 0.0) { + long bits = Double.doubleToLongBits(x); + if (bits < 0) { + return -0.0; + } + return 0.0; + } + + if (xa != xa || xa == Double.POSITIVE_INFINITY) { + return Double.NaN; + } + + /* Perform any argument reduction */ + if (xa > 3294198.0) { + // PI * (2**20) + // Argument too big for CodyWaite reduction. Must use + // PayneHanek. + double reduceResults[] = new double[3]; + reducePayneHanek(xa, reduceResults); + quadrant = ((int) reduceResults[0]) & 3; + xa = reduceResults[1]; + xb = reduceResults[2]; + } else if (xa > 1.5707963267948966) { + /* Inline the Cody/Waite reduction for performance */ + + // Estimate k + //k = (int)(xa / 1.5707963267948966); + int k = (int)(xa * 0.6366197723675814); + + // Compute remainder + double remA; + double remB; + while (true) { + double a = -k * 1.570796251296997; + remA = xa + a; + remB = -(remA - xa - a); + + a = -k * 7.549789948768648E-8; + double b = remA; + remA = a + b; + remB += -(remA - b - a); + + a = -k * 6.123233995736766E-17; + b = remA; + remA = a + b; + remB += -(remA - b - a); + + if (remA > 0.0) + break; + + // Remainder is negative, so decrement k and try again. + // This should only happen if the input is very close + // to an even multiple of pi/2 + k--; + } + quadrant = k & 3; + xa = remA; + xb = remB; + } + + if (negative) { + quadrant ^= 2; // Flip bit 1 + } + + switch (quadrant) { + case 0: + return sinQ(xa, xb); + case 1: + return cosQ(xa, xb); + case 2: + return -sinQ(xa, xb); + case 3: + return -cosQ(xa, xb); + default: + return Double.NaN; + } + } + + /** + * Cosine function + * @param x a number + * @return cos(x) + */ + public static double cos(double x) { + int quadrant = 0; + + /* Take absolute value of the input */ + double xa = x; + if (x < 0) { + xa = -xa; + } + + if (xa != xa || xa == Double.POSITIVE_INFINITY) { + return Double.NaN; + } + + /* Perform any argument reduction */ + double xb = 0; + if (xa > 3294198.0) { + // PI * (2**20) + // Argument too big for CodyWaite reduction. Must use + // PayneHanek. + double reduceResults[] = new double[3]; + reducePayneHanek(xa, reduceResults); + quadrant = ((int) reduceResults[0]) & 3; + xa = reduceResults[1]; + xb = reduceResults[2]; + } else if (xa > 1.5707963267948966) { + /* Inline the Cody/Waite reduction for performance */ + + // Estimate k + //k = (int)(xa / 1.5707963267948966); + int k = (int)(xa * 0.6366197723675814); + + // Compute remainder + double remA; + double remB; + while (true) { + double a = -k * 1.570796251296997; + remA = xa + a; + remB = -(remA - xa - a); + + a = -k * 7.549789948768648E-8; + double b = remA; + remA = a + b; + remB += -(remA - b - a); + + a = -k * 6.123233995736766E-17; + b = remA; + remA = a + b; + remB += -(remA - b - a); + + if (remA > 0.0) + break; + + // Remainder is negative, so decrement k and try again. + // This should only happen if the input is very close + // to an even multiple of pi/2 + k--; + } + quadrant = k & 3; + xa = remA; + xb = remB; + } + + //if (negative) + // quadrant = (quadrant + 2) % 4; + + switch (quadrant) { + case 0: + return cosQ(xa, xb); + case 1: + return -sinQ(xa, xb); + case 2: + return -cosQ(xa, xb); + case 3: + return sinQ(xa, xb); + default: + return Double.NaN; + } + } + + /** + * Tangent function + * @param x a number + * @return tan(x) + */ + public static double tan(double x) { + boolean negative = false; + int quadrant = 0; + + /* Take absolute value of the input */ + double xa = x; + if (x < 0) { + negative = true; + xa = -xa; + } + + /* Check for zero and negative zero */ + if (xa == 0.0) { + long bits = Double.doubleToLongBits(x); + if (bits < 0) { + return -0.0; + } + return 0.0; + } + + if (xa != xa || xa == Double.POSITIVE_INFINITY) { + return Double.NaN; + } + + /* Perform any argument reduction */ + double xb = 0; + if (xa > 3294198.0) { + // PI * (2**20) + // Argument too big for CodyWaite reduction. Must use + // PayneHanek. + double reduceResults[] = new double[3]; + reducePayneHanek(xa, reduceResults); + quadrant = ((int) reduceResults[0]) & 3; + xa = reduceResults[1]; + xb = reduceResults[2]; + } else if (xa > 1.5707963267948966) { + /* Inline the Cody/Waite reduction for performance */ + + // Estimate k + //k = (int)(xa / 1.5707963267948966); + int k = (int)(xa * 0.6366197723675814); + + // Compute remainder + double remA; + double remB; + while (true) { + double a = -k * 1.570796251296997; + remA = xa + a; + remB = -(remA - xa - a); + + a = -k * 7.549789948768648E-8; + double b = remA; + remA = a + b; + remB += -(remA - b - a); + + a = -k * 6.123233995736766E-17; + b = remA; + remA = a + b; + remB += -(remA - b - a); + + if (remA > 0.0) + break; + + // Remainder is negative, so decrement k and try again. + // This should only happen if the input is very close + // to an even multiple of pi/2 + k--; + } + quadrant = k & 3; + xa = remA; + xb = remB; + } + + if (xa > 1.5) { + // Accurracy suffers between 1.5 and PI/2 + final double pi2a = 1.5707963267948966; + final double pi2b = 6.123233995736766E-17; + + final double a = pi2a - xa; + double b = -(a - pi2a + xa); + b += pi2b - xb; + + xa = a + b; + xb = -(xa - a - b); + quadrant ^= 1; + negative ^= true; + } + + double result; + if ((quadrant & 1) == 0) { + result = tanQ(xa, xb, false); + } else { + result = -tanQ(xa, xb, true); + } + + if (negative) { + result = -result; + } + + return result; + } + + /** + * Arctangent function + * @param x a number + * @return atan(x) + */ + public static double atan(double x) { + return atan(x, 0.0, false); + } + + /** Internal helper function to compute arctangent. + * @param xa number from which arctangent is requested + * @param xb extra bits for x (may be 0.0) + * @param leftPlane if true, result angle must be put in the left half plane + * @return atan(xa + xb) (or angle shifted by {@code PI} if leftPlane is true) + */ + private static double atan(double xa, double xb, boolean leftPlane) { + boolean negate = false; + int idx; + + if (xa == 0.0) { // Matches +/- 0.0; return correct sign + return leftPlane ? copySign(Math.PI, xa) : xa; + } + + if (xa < 0) { + // negative + xa = -xa; + xb = -xb; + negate = true; + } + + if (xa > 1.633123935319537E16) { // Very large input + return (negate ^ leftPlane) ? (-Math.PI/2.0) : (Math.PI/2.0); + } + + /* Estimate the closest tabulated arctan value, compute eps = xa-tangentTable */ + if (xa < 1.0) { + idx = (int) (((-1.7168146928204136 * xa * xa + 8.0) * xa) + 0.5); + } else { + double temp = 1.0/xa; + idx = (int) (-((-1.7168146928204136 * temp * temp + 8.0) * temp) + 13.07); + } + double epsA = xa - TANGENT_TABLE_A[idx]; + double epsB = -(epsA - xa + TANGENT_TABLE_A[idx]); + epsB += xb - TANGENT_TABLE_B[idx]; + + double temp = epsA + epsB; + epsB = -(temp - epsA - epsB); + epsA = temp; + + /* Compute eps = eps / (1.0 + xa*tangent) */ + temp = xa * HEX_40000000; + double ya = xa + temp - temp; + double yb = xb + xa - ya; + xa = ya; + xb += yb; + + //if (idx > 8 || idx == 0) + if (idx == 0) { + /* If the slope of the arctan is gentle enough (< 0.45), this approximation will suffice */ + //double denom = 1.0 / (1.0 + xa*tangentTableA[idx] + xb*tangentTableA[idx] + xa*tangentTableB[idx] + xb*tangentTableB[idx]); + double denom = 1.0 / (1.0 + (xa + xb) * (TANGENT_TABLE_A[idx] + TANGENT_TABLE_B[idx])); + //double denom = 1.0 / (1.0 + xa*tangentTableA[idx]); + ya = epsA * denom; + yb = epsB * denom; + } else { + double temp2 = xa * TANGENT_TABLE_A[idx]; + double za = 1.0 + temp2; + double zb = -(za - 1.0 - temp2); + temp2 = xb * TANGENT_TABLE_A[idx] + xa * TANGENT_TABLE_B[idx]; + temp = za + temp2; + zb += -(temp - za - temp2); + za = temp; + + zb += xb * TANGENT_TABLE_B[idx]; + ya = epsA / za; + + temp = ya * HEX_40000000; + final double yaa = (ya + temp) - temp; + final double yab = ya - yaa; + + temp = za * HEX_40000000; + final double zaa = (za + temp) - temp; + final double zab = za - zaa; + + /* Correct for rounding in division */ + yb = (epsA - yaa * zaa - yaa * zab - yab * zaa - yab * zab) / za; + + yb += -epsA * zb / za / za; + yb += epsB / za; + } + + + epsA = ya; + epsB = yb; + + /* Evaluate polynomial */ + double epsA2 = epsA*epsA; + + /* + yb = -0.09001346640161823; + yb = yb * epsA2 + 0.11110718400605211; + yb = yb * epsA2 + -0.1428571349122913; + yb = yb * epsA2 + 0.19999999999273194; + yb = yb * epsA2 + -0.33333333333333093; + yb = yb * epsA2 * epsA; + */ + + yb = 0.07490822288864472; + yb = yb * epsA2 + -0.09088450866185192; + yb = yb * epsA2 + 0.11111095942313305; + yb = yb * epsA2 + -0.1428571423679182; + yb = yb * epsA2 + 0.19999999999923582; + yb = yb * epsA2 + -0.33333333333333287; + yb = yb * epsA2 * epsA; + + + ya = epsA; + + temp = ya + yb; + yb = -(temp - ya - yb); + ya = temp; + + /* Add in effect of epsB. atan'(x) = 1/(1+x^2) */ + yb += epsB / (1.0 + epsA * epsA); + + double result; + double resultb; + + //result = yb + eighths[idx] + ya; + double za = EIGHTHS[idx] + ya; + double zb = -(za - EIGHTHS[idx] - ya); + temp = za + yb; + zb += -(temp - za - yb); + za = temp; + + result = za + zb; + resultb = -(result - za - zb); + + if (leftPlane) { + // Result is in the left plane + final double pia = 1.5707963267948966*2.0; + final double pib = 6.123233995736766E-17*2.0; + + za = pia - result; + zb = -(za - pia + result); + zb += pib - resultb; + + result = za + zb; + resultb = -(result - za - zb); + } + + + if (negate ^ leftPlane) { + result = -result; + } + + return result; + } + + /** + * Two arguments arctangent function + * @param y ordinate + * @param x abscissa + * @return phase angle of point (x,y) between {@code -PI} and {@code PI} + */ + public static double atan2(double y, double x) { + if (x !=x || y != y) { + return Double.NaN; + } + + if (y == 0.0) { + double result = x*y; + double invx = 1.0/x; + double invy = 1.0/y; + + if (invx == 0.0) { // X is infinite + if (x > 0) { + return y; // return +/- 0.0 + } else { + return copySign(Math.PI, y); + } + } + + if (x < 0.0 || invx < 0.0) { + if (y < 0.0 || invy < 0.0) { + return -Math.PI; + } else { + return Math.PI; + } + } else { + return result; + } + } + + // y cannot now be zero + + if (y == Double.POSITIVE_INFINITY) { + if (x == Double.POSITIVE_INFINITY) { + return Math.PI/4.0; + } + + if (x == Double.NEGATIVE_INFINITY) { + return Math.PI*3.0/4.0; + } + + return Math.PI/2.0; + } + + if (y == Double.NEGATIVE_INFINITY) { + if (x == Double.POSITIVE_INFINITY) { + return -Math.PI/4.0; + } + + if (x == Double.NEGATIVE_INFINITY) { + return -Math.PI*3.0/4.0; + } + + return -Math.PI/2.0; + } + + if (x == Double.POSITIVE_INFINITY) { + if (y > 0.0 || 1/y > 0.0) { + return 0.0; + } + + if (y < 0.0 || 1/y < 0.0) { + return -0.0; + } + } + + if (x == Double.NEGATIVE_INFINITY) + { + if (y > 0.0 || 1/y > 0.0) { + return Math.PI; + } + + if (y < 0.0 || 1/y < 0.0) { + return -Math.PI; + } + } + + // Neither y nor x can be infinite or NAN here + + if (x == 0) { + if (y > 0.0 || 1/y > 0.0) { + return Math.PI/2.0; + } + + if (y < 0.0 || 1/y < 0.0) { + return -Math.PI/2.0; + } + } + + // Compute ratio r = y/x + final double r = y/x; + if (Double.isInfinite(r)) { // bypass calculations that can create NaN + return atan(r, 0, x < 0); + } + + double ra = doubleHighPart(r); + double rb = r - ra; + + // Split x + final double xa = doubleHighPart(x); + final double xb = x - xa; + + rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x; + + double temp = ra + rb; + rb = -(temp - ra - rb); + ra = temp; + + if (ra == 0) { // Fix up the sign so atan works correctly + ra = copySign(0.0, y); + } + + // Call atan + double result = atan(ra, rb, x < 0); + + return result; + } + + /** Compute the arc sine of a number. + * @param x number on which evaluation is done + * @return arc sine of x + */ + public static double asin(double x) { + if (x != x) { + return Double.NaN; + } + + if (x > 1.0 || x < -1.0) { + return Double.NaN; + } + + if (x == 1.0) { + return Math.PI/2.0; + } + + if (x == -1.0) { + return -Math.PI/2.0; + } + + if (x == 0.0) { // Matches +/- 0.0; return correct sign + return x; + } + + /* Compute asin(x) = atan(x/sqrt(1-x*x)) */ + + /* Split x */ + double temp = x * HEX_40000000; + final double xa = x + temp - temp; + final double xb = x - xa; + + /* Square it */ + double ya = xa*xa; + double yb = xa*xb*2.0 + xb*xb; + + /* Subtract from 1 */ + ya = -ya; + yb = -yb; + + double za = 1.0 + ya; + double zb = -(za - 1.0 - ya); + + temp = za + yb; + zb += -(temp - za - yb); + za = temp; + + /* Square root */ + double y; + y = sqrt(za); + temp = y * HEX_40000000; + ya = y + temp - temp; + yb = y - ya; + + /* Extend precision of sqrt */ + yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y); + + /* Contribution of zb to sqrt */ + double dx = zb / (2.0*y); + + // Compute ratio r = x/y + double r = x/y; + temp = r * HEX_40000000; + double ra = r + temp - temp; + double rb = r - ra; + + rb += (x - ra*ya - ra*yb - rb*ya - rb*yb) / y; // Correct for rounding in division + rb += -x * dx / y / y; // Add in effect additional bits of sqrt. + + temp = ra + rb; + rb = -(temp - ra - rb); + ra = temp; + + return atan(ra, rb, false); + } + + /** Compute the arc cosine of a number. + * @param x number on which evaluation is done + * @return arc cosine of x + */ + public static double acos(double x) { + if (x != x) { + return Double.NaN; + } + + if (x > 1.0 || x < -1.0) { + return Double.NaN; + } + + if (x == -1.0) { + return Math.PI; + } + + if (x == 1.0) { + return 0.0; + } + + if (x == 0) { + return Math.PI/2.0; + } + + /* Compute acos(x) = atan(sqrt(1-x*x)/x) */ + + /* Split x */ + double temp = x * HEX_40000000; + final double xa = x + temp - temp; + final double xb = x - xa; + + /* Square it */ + double ya = xa*xa; + double yb = xa*xb*2.0 + xb*xb; + + /* Subtract from 1 */ + ya = -ya; + yb = -yb; + + double za = 1.0 + ya; + double zb = -(za - 1.0 - ya); + + temp = za + yb; + zb += -(temp - za - yb); + za = temp; + + /* Square root */ + double y = sqrt(za); + temp = y * HEX_40000000; + ya = y + temp - temp; + yb = y - ya; + + /* Extend precision of sqrt */ + yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y); + + /* Contribution of zb to sqrt */ + yb += zb / (2.0*y); + y = ya+yb; + yb = -(y - ya - yb); + + // Compute ratio r = y/x + double r = y/x; + + // Did r overflow? + if (Double.isInfinite(r)) { // x is effectively zero + return Math.PI/2; // so return the appropriate value + } + + double ra = doubleHighPart(r); + double rb = r - ra; + + rb += (y - ra*xa - ra*xb - rb*xa - rb*xb) / x; // Correct for rounding in division + rb += yb / x; // Add in effect additional bits of sqrt. + + temp = ra + rb; + rb = -(temp - ra - rb); + ra = temp; + + return atan(ra, rb, x<0); + } + + /** Compute the cubic root of a number. + * @param x number on which evaluation is done + * @return cubic root of x + */ + public static double cbrt(double x) { + /* Convert input double to bits */ + long inbits = Double.doubleToLongBits(x); + int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023; + boolean subnormal = false; + + if (exponent == -1023) { + if (x == 0) { + return x; + } + + /* Subnormal, so normalize */ + subnormal = true; + x *= 1.8014398509481984E16; // 2^54 + inbits = Double.doubleToLongBits(x); + exponent = (int) ((inbits >> 52) & 0x7ff) - 1023; + } + + if (exponent == 1024) { + // Nan or infinity. Don't care which. + return x; + } + + /* Divide the exponent by 3 */ + int exp3 = exponent / 3; + + /* p2 will be the nearest power of 2 to x with its exponent divided by 3 */ + double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) | + (long)(((exp3 + 1023) & 0x7ff)) << 52); + + /* This will be a number between 1 and 2 */ + final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L); + + /* Estimate the cube root of mant by polynomial */ + double est = -0.010714690733195933; + est = est * mant + 0.0875862700108075; + est = est * mant + -0.3058015757857271; + est = est * mant + 0.7249995199969751; + est = est * mant + 0.5039018405998233; + + est *= CBRTTWO[exponent % 3 + 2]; + + // est should now be good to about 15 bits of precision. Do 2 rounds of + // Newton's method to get closer, this should get us full double precision + // Scale down x for the purpose of doing newtons method. This avoids over/under flows. + final double xs = x / (p2*p2*p2); + est += (xs - est*est*est) / (3*est*est); + est += (xs - est*est*est) / (3*est*est); + + // Do one round of Newton's method in extended precision to get the last bit right. + double temp = est * HEX_40000000; + double ya = est + temp - temp; + double yb = est - ya; + + double za = ya * ya; + double zb = ya * yb * 2.0 + yb * yb; + temp = za * HEX_40000000; + double temp2 = za + temp - temp; + zb += za - temp2; + za = temp2; + + zb = za * yb + ya * zb + zb * yb; + za = za * ya; + + double na = xs - za; + double nb = -(na - xs + za); + nb -= zb; + + est += (na+nb)/(3*est*est); + + /* Scale by a power of two, so this is exact. */ + est *= p2; + + if (subnormal) { + est *= 3.814697265625E-6; // 2^-18 + } + + return est; + } + + /** + * Convert degrees to radians, with error of less than 0.5 ULP + * @param x angle in degrees + * @return x converted into radians + */ + public static double toRadians(double x) + { + if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign + return x; + } + + // These are PI/180 split into high and low order bits + final double facta = 0.01745329052209854; + final double factb = 1.997844754509471E-9; + + double xa = doubleHighPart(x); + double xb = x - xa; + + double result = xb * factb + xb * facta + xa * factb + xa * facta; + if (result == 0) { + result = result * x; // ensure correct sign if calculation underflows + } + return result; + } + + /** + * Convert radians to degrees, with error of less than 0.5 ULP + * @param x angle in radians + * @return x converted into degrees + */ + public static double toDegrees(double x) + { + if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign + return x; + } + + // These are 180/PI split into high and low order bits + final double facta = 57.2957763671875; + final double factb = 3.145894820876798E-6; + + double xa = doubleHighPart(x); + double xb = x - xa; + + return xb * factb + xb * facta + xa * factb + xa * facta; + } + + /** + * Absolute value. + * @param x number from which absolute value is requested + * @return abs(x) + */ + public static int abs(final int x) { + return (x < 0) ? -x : x; + } + + /** + * Absolute value. + * @param x number from which absolute value is requested + * @return abs(x) + */ + public static long abs(final long x) { + return (x < 0l) ? -x : x; + } + + /** + * Absolute value. + * @param x number from which absolute value is requested + * @return abs(x) + */ + public static float abs(final float x) { + return (x < 0.0f) ? -x : (x == 0.0f) ? 0.0f : x; // -0.0 => +0.0 + } + + /** + * Absolute value. + * @param x number from which absolute value is requested + * @return abs(x) + */ + public static double abs(double x) { + return (x < 0.0) ? -x : (x == 0.0) ? 0.0 : x; // -0.0 => +0.0 + } + + /** + * Compute least significant bit (Unit in Last Position) for a number. + * @param x number from which ulp is requested + * @return ulp(x) + */ + public static double ulp(double x) { + if (Double.isInfinite(x)) { + return Double.POSITIVE_INFINITY; + } + return abs(x - Double.longBitsToDouble(Double.doubleToLongBits(x) ^ 1)); + } + + /** + * Compute least significant bit (Unit in Last Position) for a number. + * @param x number from which ulp is requested + * @return ulp(x) + */ + public static float ulp(float x) { + if (Float.isInfinite(x)) { + return Float.POSITIVE_INFINITY; + } + return abs(x - Float.intBitsToFloat(Float.floatToIntBits(x) ^ 1)); + } + + /** + * Multiply a double number by a power of 2. + * @param d number to multiply + * @param n power of 2 + * @return d × 2<sup>n</sup> + */ + public static double scalb(final double d, final int n) { + + // first simple and fast handling when 2^n can be represented using normal numbers + if ((n > -1023) && (n < 1024)) { + return d * Double.longBitsToDouble(((long) (n + 1023)) << 52); + } + + // handle special cases + if (Double.isNaN(d) || Double.isInfinite(d) || (d == 0)) { + return d; + } + if (n < -2098) { + return (d > 0) ? 0.0 : -0.0; + } + if (n > 2097) { + return (d > 0) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY; + } + + // decompose d + final long bits = Double.doubleToLongBits(d); + final long sign = bits & 0x8000000000000000L; + int exponent = ((int) (bits >>> 52)) & 0x7ff; + long mantissa = bits & 0x000fffffffffffffL; + + // compute scaled exponent + int scaledExponent = exponent + n; + + if (n < 0) { + // we are really in the case n <= -1023 + if (scaledExponent > 0) { + // both the input and the result are normal numbers, we only adjust the exponent + return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa); + } else if (scaledExponent > -53) { + // the input is a normal number and the result is a subnormal number + + // recover the hidden mantissa bit + mantissa = mantissa | (1L << 52); + + // scales down complete mantissa, hence losing least significant bits + final long mostSignificantLostBit = mantissa & (1L << (-scaledExponent)); + mantissa = mantissa >>> (1 - scaledExponent); + if (mostSignificantLostBit != 0) { + // we need to add 1 bit to round up the result + mantissa++; + } + return Double.longBitsToDouble(sign | mantissa); + + } else { + // no need to compute the mantissa, the number scales down to 0 + return (sign == 0L) ? 0.0 : -0.0; + } + } else { + // we are really in the case n >= 1024 + if (exponent == 0) { + + // the input number is subnormal, normalize it + while ((mantissa >>> 52) != 1) { + mantissa = mantissa << 1; + --scaledExponent; + } + ++scaledExponent; + mantissa = mantissa & 0x000fffffffffffffL; + + if (scaledExponent < 2047) { + return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa); + } else { + return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY; + } + + } else if (scaledExponent < 2047) { + return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa); + } else { + return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY; + } + } + + } + + /** + * Multiply a float number by a power of 2. + * @param f number to multiply + * @param n power of 2 + * @return f × 2<sup>n</sup> + */ + public static float scalb(final float f, final int n) { + + // first simple and fast handling when 2^n can be represented using normal numbers + if ((n > -127) && (n < 128)) { + return f * Float.intBitsToFloat((n + 127) << 23); + } + + // handle special cases + if (Float.isNaN(f) || Float.isInfinite(f) || (f == 0f)) { + return f; + } + if (n < -277) { + return (f > 0) ? 0.0f : -0.0f; + } + if (n > 276) { + return (f > 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY; + } + + // decompose f + final int bits = Float.floatToIntBits(f); + final int sign = bits & 0x80000000; + int exponent = (bits >>> 23) & 0xff; + int mantissa = bits & 0x007fffff; + + // compute scaled exponent + int scaledExponent = exponent + n; + + if (n < 0) { + // we are really in the case n <= -127 + if (scaledExponent > 0) { + // both the input and the result are normal numbers, we only adjust the exponent + return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa); + } else if (scaledExponent > -24) { + // the input is a normal number and the result is a subnormal number + + // recover the hidden mantissa bit + mantissa = mantissa | (1 << 23); + + // scales down complete mantissa, hence losing least significant bits + final int mostSignificantLostBit = mantissa & (1 << (-scaledExponent)); + mantissa = mantissa >>> (1 - scaledExponent); + if (mostSignificantLostBit != 0) { + // we need to add 1 bit to round up the result + mantissa++; + } + return Float.intBitsToFloat(sign | mantissa); + + } else { + // no need to compute the mantissa, the number scales down to 0 + return (sign == 0) ? 0.0f : -0.0f; + } + } else { + // we are really in the case n >= 128 + if (exponent == 0) { + + // the input number is subnormal, normalize it + while ((mantissa >>> 23) != 1) { + mantissa = mantissa << 1; + --scaledExponent; + } + ++scaledExponent; + mantissa = mantissa & 0x007fffff; + + if (scaledExponent < 255) { + return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa); + } else { + return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY; + } + + } else if (scaledExponent < 255) { + return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa); + } else { + return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY; + } + } + + } + + /** + * Get the next machine representable number after a number, moving + * in the direction of another number. + * <p> + * The ordering is as follows (increasing): + * <ul> + * <li>-INFINITY</li> + * <li>-MAX_VALUE</li> + * <li>-MIN_VALUE</li> + * <li>-0.0</li> + * <li>+0.0</li> + * <li>+MIN_VALUE</li> + * <li>+MAX_VALUE</li> + * <li>+INFINITY</li> + * <li></li> + * <p> + * If arguments compare equal, then the second argument is returned. + * <p> + * If {@code direction} is greater than {@code d}, + * the smallest machine representable number strictly greater than + * {@code d} is returned; if less, then the largest representable number + * strictly less than {@code d} is returned.</p> + * <p> + * If {@code d} is infinite and direction does not + * bring it back to finite numbers, it is returned unchanged.</p> + * + * @param d base number + * @param direction (the only important thing is whether + * {@code direction} is greater or smaller than {@code d}) + * @return the next machine representable number in the specified direction + */ + public static double nextAfter(double d, double direction) { + + // handling of some important special cases + if (Double.isNaN(d) || Double.isNaN(direction)) { + return Double.NaN; + } else if (d == direction) { + return direction; + } else if (Double.isInfinite(d)) { + return (d < 0) ? -Double.MAX_VALUE : Double.MAX_VALUE; + } else if (d == 0) { + return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE; + } + // special cases MAX_VALUE to infinity and MIN_VALUE to 0 + // are handled just as normal numbers + + final long bits = Double.doubleToLongBits(d); + final long sign = bits & 0x8000000000000000L; + if ((direction < d) ^ (sign == 0L)) { + return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) + 1)); + } else { + return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) - 1)); + } + + } + + /** + * Get the next machine representable number after a number, moving + * in the direction of another number. + * <p> + * The ordering is as follows (increasing): + * <ul> + * <li>-INFINITY</li> + * <li>-MAX_VALUE</li> + * <li>-MIN_VALUE</li> + * <li>-0.0</li> + * <li>+0.0</li> + * <li>+MIN_VALUE</li> + * <li>+MAX_VALUE</li> + * <li>+INFINITY</li> + * <li></li> + * <p> + * If arguments compare equal, then the second argument is returned. + * <p> + * If {@code direction} is greater than {@code f}, + * the smallest machine representable number strictly greater than + * {@code f} is returned; if less, then the largest representable number + * strictly less than {@code f} is returned.</p> + * <p> + * If {@code f} is infinite and direction does not + * bring it back to finite numbers, it is returned unchanged.</p> + * + * @param f base number + * @param direction (the only important thing is whether + * {@code direction} is greater or smaller than {@code f}) + * @return the next machine representable number in the specified direction + */ + public static float nextAfter(final float f, final double direction) { + + // handling of some important special cases + if (Double.isNaN(f) || Double.isNaN(direction)) { + return Float.NaN; + } else if (f == direction) { + return (float) direction; + } else if (Float.isInfinite(f)) { + return (f < 0f) ? -Float.MAX_VALUE : Float.MAX_VALUE; + } else if (f == 0f) { + return (direction < 0) ? -Float.MIN_VALUE : Float.MIN_VALUE; + } + // special cases MAX_VALUE to infinity and MIN_VALUE to 0 + // are handled just as normal numbers + + final int bits = Float.floatToIntBits(f); + final int sign = bits & 0x80000000; + if ((direction < f) ^ (sign == 0)) { + return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) + 1)); + } else { + return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) - 1)); + } + + } + + /** Get the largest whole number smaller than x. + * @param x number from which floor is requested + * @return a double number f such that f is an integer f <= x < f + 1.0 + */ + public static double floor(double x) { + long y; + + if (x != x) { // NaN + return x; + } + + if (x >= TWO_POWER_52 || x <= -TWO_POWER_52) { + return x; + } + + y = (long) x; + if (x < 0 && y != x) { + y--; + } + + if (y == 0) { + return x*y; + } + + return y; + } + + /** Get the smallest whole number larger than x. + * @param x number from which ceil is requested + * @return a double number c such that c is an integer c - 1.0 < x <= c + */ + public static double ceil(double x) { + double y; + + if (x != x) { // NaN + return x; + } + + y = floor(x); + if (y == x) { + return y; + } + + y += 1.0; + + if (y == 0) { + return x*y; + } + + return y; + } + + /** Get the whole number that is the nearest to x, or the even one if x is exactly half way between two integers. + * @param x number from which nearest whole number is requested + * @return a double number r such that r is an integer r - 0.5 <= x <= r + 0.5 + */ + public static double rint(double x) { + double y = floor(x); + double d = x - y; + + if (d > 0.5) { + if (y == -1.0) { + return -0.0; // Preserve sign of operand + } + return y+1.0; + } + if (d < 0.5) { + return y; + } + + /* half way, round to even */ + long z = (long) y; + return (z & 1) == 0 ? y : y + 1.0; + } + + /** Get the closest long to x. + * @param x number from which closest long is requested + * @return closest long to x + */ + public static long round(double x) { + return (long) floor(x + 0.5); + } + + /** Get the closest int to x. + * @param x number from which closest int is requested + * @return closest int to x + */ + public static int round(final float x) { + return (int) floor(x + 0.5f); + } + + /** Compute the minimum of two values + * @param a first value + * @param b second value + * @return a if a is lesser or equal to b, b otherwise + */ + public static int min(final int a, final int b) { + return (a <= b) ? a : b; + } + + /** Compute the minimum of two values + * @param a first value + * @param b second value + * @return a if a is lesser or equal to b, b otherwise + */ + public static long min(final long a, final long b) { + return (a <= b) ? a : b; + } + + /** Compute the minimum of two values + * @param a first value + * @param b second value + * @return a if a is lesser or equal to b, b otherwise + */ + public static float min(final float a, final float b) { + if (a > b) { + return b; + } + if (a < b) { + return a; + } + /* if either arg is NaN, return NaN */ + if (a != b) { + return Float.NaN; + } + /* min(+0.0,-0.0) == -0.0 */ + /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */ + int bits = Float.floatToRawIntBits(a); + if (bits == 0x80000000) { + return a; + } + return b; + } + + /** Compute the minimum of two values + * @param a first value + * @param b second value + * @return a if a is lesser or equal to b, b otherwise + */ + public static double min(final double a, final double b) { + if (a > b) { + return b; + } + if (a < b) { + return a; + } + /* if either arg is NaN, return NaN */ + if (a != b) { + return Double.NaN; + } + /* min(+0.0,-0.0) == -0.0 */ + /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */ + long bits = Double.doubleToRawLongBits(a); + if (bits == 0x8000000000000000L) { + return a; + } + return b; + } + + /** Compute the maximum of two values + * @param a first value + * @param b second value + * @return b if a is lesser or equal to b, a otherwise + */ + public static int max(final int a, final int b) { + return (a <= b) ? b : a; + } + + /** Compute the maximum of two values + * @param a first value + * @param b second value + * @return b if a is lesser or equal to b, a otherwise + */ + public static long max(final long a, final long b) { + return (a <= b) ? b : a; + } + + /** Compute the maximum of two values + * @param a first value + * @param b second value + * @return b if a is lesser or equal to b, a otherwise + */ + public static float max(final float a, final float b) { + if (a > b) { + return a; + } + if (a < b) { + return b; + } + /* if either arg is NaN, return NaN */ + if (a != b) { + return Float.NaN; + } + /* min(+0.0,-0.0) == -0.0 */ + /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */ + int bits = Float.floatToRawIntBits(a); + if (bits == 0x80000000) { + return b; + } + return a; + } + + /** Compute the maximum of two values + * @param a first value + * @param b second value + * @return b if a is lesser or equal to b, a otherwise + */ + public static double max(final double a, final double b) { + if (a > b) { + return a; + } + if (a < b) { + return b; + } + /* if either arg is NaN, return NaN */ + if (a != b) { + return Double.NaN; + } + /* min(+0.0,-0.0) == -0.0 */ + /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */ + long bits = Double.doubleToRawLongBits(a); + if (bits == 0x8000000000000000L) { + return b; + } + return a; + } + + /** + * Returns the hypotenuse of a triangle with sides {@code x} and {@code y} + * - sqrt(<i>x</i><sup>2</sup> +<i>y</i><sup>2</sup>)<br/> + * avoiding intermediate overflow or underflow. + * + * <ul> + * <li> If either argument is infinite, then the result is positive infinity.</li> + * <li> else, if either argument is NaN then the result is NaN.</li> + * </ul> + * + * @param x a value + * @param y a value + * @return sqrt(<i>x</i><sup>2</sup> +<i>y</i><sup>2</sup>) + */ + public static double hypot(final double x, final double y) { + if (Double.isInfinite(x) || Double.isInfinite(y)) { + return Double.POSITIVE_INFINITY; + } else if (Double.isNaN(x) || Double.isNaN(y)) { + return Double.NaN; + } else { + + final int expX = getExponent(x); + final int expY = getExponent(y); + if (expX > expY + 27) { + // y is neglectible with respect to x + return abs(x); + } else if (expY > expX + 27) { + // x is neglectible with respect to y + return abs(y); + } else { + + // find an intermediate scale to avoid both overflow and underflow + final int middleExp = (expX + expY) / 2; + + // scale parameters without losing precision + final double scaledX = scalb(x, -middleExp); + final double scaledY = scalb(y, -middleExp); + + // compute scaled hypotenuse + final double scaledH = sqrt(scaledX * scaledX + scaledY * scaledY); + + // remove scaling + return scalb(scaledH, middleExp); + + } + + } + } + + /** + * Computes the remainder as prescribed by the IEEE 754 standard. + * The remainder value is mathematically equal to {@code x - y*n} + * where {@code n} is the mathematical integer closest to the exact mathematical value + * of the quotient {@code x/y}. + * If two mathematical integers are equally close to {@code x/y} then + * {@code n} is the integer that is even. + * <p> + * <ul> + * <li>If either operand is NaN, the result is NaN.</li> + * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li> + * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li> + * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li> + * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li> + * </ul> + * <p><b>Note:</b> this implementation currently delegates to {@link StrictMath#IEEEremainder} + * @param dividend the number to be divided + * @param divisor the number by which to divide + * @return the remainder, rounded + */ + public static double IEEEremainder(double dividend, double divisor) { + return StrictMath.IEEEremainder(dividend, divisor); // TODO provide our own implementation + } + + /** + * Returns the first argument with the sign of the second argument. + * A NaN {@code sign} argument is treated as positive. + * + * @param magnitude the value to return + * @param sign the sign for the returned value + * @return the magnitude with the same sign as the {@code sign} argument + */ + public static double copySign(double magnitude, double sign){ + long m = Double.doubleToLongBits(magnitude); + long s = Double.doubleToLongBits(sign); + if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK + return magnitude; + } + return -magnitude; // flip sign + } + + /** + * Returns the first argument with the sign of the second argument. + * A NaN {@code sign} argument is treated as positive. + * + * @param magnitude the value to return + * @param sign the sign for the returned value + * @return the magnitude with the same sign as the {@code sign} argument + */ + public static float copySign(float magnitude, float sign){ + int m = Float.floatToIntBits(magnitude); + int s = Float.floatToIntBits(sign); + if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK + return magnitude; + } + return -magnitude; // flip sign + } + + /** + * Return the exponent of a double number, removing the bias. + * <p> + * For double numbers of the form 2<sup>x</sup>, the unbiased + * exponent is exactly x. + * </p> + * @param d number from which exponent is requested + * @return exponent for d in IEEE754 representation, without bias + */ + public static int getExponent(final double d) { + return (int) ((Double.doubleToLongBits(d) >>> 52) & 0x7ff) - 1023; + } + + /** + * Return the exponent of a float number, removing the bias. + * <p> + * For float numbers of the form 2<sup>x</sup>, the unbiased + * exponent is exactly x. + * </p> + * @param f number from which exponent is requested + * @return exponent for d in IEEE754 representation, without bias + */ + public static int getExponent(final float f) { + return ((Float.floatToIntBits(f) >>> 23) & 0xff) - 127; + } + +} |