You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2010/08/30 00:04:15 UTC
svn commit: r990658 [6/10] - in /commons/proper/math/trunk: ./
src/main/java/org/apache/commons/math/analysis/
src/main/java/org/apache/commons/math/analysis/integration/
src/main/java/org/apache/commons/math/analysis/interpolation/
src/main/java/org/a...
Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMath.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMath.java?rev=990658&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMath.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMath.java Sun Aug 29 22:04:09 2010
@@ -0,0 +1,2967 @@
+/*
+ * 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 StrictMath.
+ * @version $Revision$ $Date$
+ * @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 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 long PI_O_4_BITS[] = new long[] {
+ (0xc90fdaa2L << 32) | 0x2168c234L,
+ (0xc4c6628bL << 32) | 0x80dc1cd1L };
+
+ /** Eighthes.
+ * 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 EIGHTHES[] = {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};
+
+ // Initialize tables
+ static {
+ int i;
+
+ // Generate an array of factorials
+ FACT[0] = 1.0;
+ for (i = 1; i < 20; 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 < 1025; 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 < 1024; 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() {
+ }
+
+ /** Compute the arc cosine of a number.
+ * @param a number on which evaluation is done
+ * @return arc cosine of a
+ */
+ public static double acos(final double a) {
+ return Math.acos(a);
+ }
+
+ /** Compute the arc sine of a number.
+ * @param a number on which evaluation is done
+ * @return arc sine of a
+ */
+ public static double asin(final double a) {
+ return Math.asin(a);
+ }
+
+ /** Compute the square root of a number.
+ * @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 cubic root of a number.
+ * @param a number on which evaluation is done
+ * @return cubic root of a
+ */
+ public static double cbrt(final double a) {
+ return Math.cbrt(a);
+ }
+
+ /** Compute the hyperbolic cosine of a number.
+ * @param a number on which evaluation is done
+ * @return hyperbolic cosine of a
+ */
+ public static double cosh(final double a) {
+ return 0.5 * (FastMath.exp(a) + FastMath.exp(-a));
+ }
+
+ /** Compute the hyperbolic sine of a number.
+ * @param a number on which evaluation is done
+ * @return hyperbolic sine of a
+ */
+ public static double sinh(double a) {
+
+ boolean negative = false;
+ if (a < 0) {
+ negative = true;
+ a = -a;
+ }
+
+ double absSinh;
+ if (a > 0.3) {
+ absSinh = 0.5 * (FastMath.exp(a) - FastMath.exp(-a));
+ } else {
+ final double a2 = a * a;
+ if (a > 0.05) {
+ absSinh = a * (1 + a2 * (1 + a2 * (1 + a2 * (1 + a2 * (1 + a2 / 110) / 72) / 42) / 20) / 6);
+ } else {
+ absSinh = a * (1 + a2 * (1 + a2 * (1 + a2 / 42) / 20) / 6);
+ }
+ }
+
+ return negative ? -absSinh : absSinh;
+
+ }
+
+ /** Compute the hyperbolic tangent of a number.
+ * @param a number on which evaluation is done
+ * @return hyperbolic tangent of a
+ */
+ public static double tanh(double a) {
+
+ boolean negative = false;
+ if (a < 0) {
+ negative = true;
+ a = -a;
+ }
+
+ double absTanh;
+ if (a > 0.074) {
+ final double twoA = 2 * a;
+ absTanh = FastMath.expm1(twoA) / (FastMath.exp(twoA) + 1);
+ } else {
+ final double a2 = a * a;
+ if (a > 0.016) {
+ absTanh = a * (1 - a2 * (1 - a2 * (2 - a2 * (17 - a2 * (62 - a2 * 1382 / 55 ) / 9) / 21) / 5) / 3);
+ } else {
+ absTanh = a * (1 - a2 * (1 - a2 * (2 - a2 * 17 / 21) / 5) / 3);
+ }
+ }
+
+ return negative ? -absTanh : absTanh;
+
+ }
+
+ /** 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, +1 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 : (Double.isNaN(a) ? Double.NaN : 0.0));
+ }
+
+ /** 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);
+ }
+
+ /** Returns a pseudo-random number between 0.0 and 1.0.
+ * @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) {
+ 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) {
+ 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) {
+ 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 * 1073741824.0;
+ 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 * 1073741824.0;
+ 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 * 1073741824.0;
+ final double ra = ratio + temp - temp;
+ double rb = ratio - ra;
+
+ temp = denom * 1073741824.0;
+ 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;
+ }
+
+ 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 * 1073741824.0;
+ split[0] = (d + a) - a;
+ split[1] = d - split[0];
+ } else {
+ final double a = d * 9.31322574615478515625E-10;
+ split[0] = (d + a - d) * 1073741824.0;
+ 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 * 1073741824.0;
+ a[0] = (c + z) - z;
+ a[1] = c - a[0] + d;
+ } else {
+ double z = c * 9.31322574615478515625E-10;
+ a[0] = (c + z - c) * 1073741824.0;
+ 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) {
+ 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 * 1073741824.0;
+ 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 * 1073741824.0;
+ 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 * 1073741824.0;
+ 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 * 1073741824.0;
+ ya = aa + tmp - tmp;
+ yb = aa - ya + ab;
+
+ if (hiPrec != null) {
+ hiPrec[0] = ya;
+ hiPrec[1] = yb;
+ }
+
+ 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 = (double)(bits & 0x3ffffffffffL) / (4503599627370496.0 + (bits & 0x000ffc0000000000L));
+
+ double lnza = 0.0;
+ double lnzb = 0.0;
+
+ if (hiPrec != null) {
+ /* split epsilon -> x */
+ double tmp = epsilon * 1073741824.0;
+ 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 = (double)(bits & 0x3ffffffffffL);
+ double denom = 4503599627370496.0 + (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 * 1073741824.0;
+ 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 * 1073741824.0;
+ 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>1e-6 || x<-1e-6) {
+ double hiPrec[] = new double[2];
+
+ log(xpa, hiPrec);
+
+ /* 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];
+
+ log(x, hiPrec);
+
+ final double tmp = hiPrec[0] * 1073741824.0;
+ 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;
+ }
+
+ /* Handle special case x<0 */
+ if (x < 0) {
+ if (y == (long) y) {
+ // If y is an integer
+ return ((long)y & 1) == 0 ? pow(-x, y) : -pow(-x, y);
+ } else {
+ return Double.NaN;
+ }
+ }
+
+ if (x == 0) {
+ long bits = Double.doubleToLongBits(x);
+ if ((bits & 0x8000000000000000L) != 0) {
+ // -zero
+ if (y < 0 && y == (long)y)
+ return Double.NEGATIVE_INFINITY;
+ }
+
+ if (y < 0) {
+ return Double.POSITIVE_INFINITY;
+ }
+ if (y > 0) {
+ return 0.0;
+ }
+
+ return Double.NaN;
+ }
+
+ if (x == Double.POSITIVE_INFINITY) {
+ if (y < 0.0) {
+ return 0.0;
+ } else {
+ return Double.POSITIVE_INFINITY;
+ }
+ }
+
+ if (y == Double.POSITIVE_INFINITY) {
+ if (x * x > 1.0) {
+ return Double.POSITIVE_INFINITY;
+ } else {
+ return 0.0;
+ }
+ }
+
+ if (y == Double.NEGATIVE_INFINITY) {
+ if (x*x < 1.0) {
+ return Double.NEGATIVE_INFINITY;
+ } else {
+ return 0.0;
+ }
+ }
+
+ /* Split y into ya and yb such that y = ya+yb */
+ double tmp1 = y * 1073741824.0;
+ final double ya = y + tmp1 - tmp1;
+ final double yb = y - ya;
+
+ /* Compute ln(x) */
+ log(x, lns);
+ double lna = lns[0];
+ double lnb = lns[1];
+
+ /* resplit lns */
+ tmp1 = lna * 1073741824.0;
+ final 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 - EIGHTHES[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 * 1073741824.0;
+ 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 - EIGHTHES[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 * 1073741824.0;
+ 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 * 1073741824.0;
+ double esta = (est + temp) - temp;
+ double estb = est - esta;
+
+ temp = cosa * 1073741824.0;
+ 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) / 4503599627370496.0; // High order 52 bits
+ double tmpB = (((prod2A & 0xfffL) << 40) + (prod2B >>> 24)) / 4503599627370496.0 / 4503599627370496.0; // 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;
+ xb = 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 π if leftPlane is true)
+ */
+ private static double atan(double xa, double xb, boolean leftPlane) {
+ boolean negate = false;
+ boolean recip = false;
+ int idx;
+
+ 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 * 1073741824.0;
+ 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 * 1073741824.0;
+ final double yaa = (ya + temp) - temp;
+ final double yab = ya - yaa;
+
+ temp = za * 1073741824.0;
+ 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;
+ if (recip) {
+ final double pi2a = 1.5707963267948966;
+ final double pi2b = 6.123233995736766E-17;
+
+ double za = pi2a - ya;
+ double zb = -(za - pi2a + ya);
+ temp = za - EIGHTHES[idx];
+ zb += -(temp - za + EIGHTHES[idx]);
+ za = temp;
+
+ zb += pi2b - yb;
+ ya = za;
+ yb = zb;
+
+ result = yb + ya;
+ resultb = -(result - yb - ya);
+ } else {
+ //result = yb + eighths[idx] + ya;
+ double za = EIGHTHES[idx] + ya;
+ double zb = -(za - EIGHTHES[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;
+
+ final double za = pia - result;
+ double 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 -π and π
+ */
+ 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
+ return 0.0;
+ }
+
+ if (result != result) { // y must be infinite
+ return x/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;
+ }
+ }
+
+ 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;
+ }
+ }
+
+ 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;
+ }
+ }
+
+ if (x > 8e298 || x < -8e298) { // This would cause split of x to fail
+ x *= 9.31322574615478515625E-10;
+ y *= 9.31322574615478515625E-10;
+ }
+
+ // Split y
+ double temp = x * 1073741824.0;
+ final double xa = x + temp - temp;
+ final double xb = x - xa;
+
+ // Compute ratio r = x/y
+ final double r = y/x;
+ temp = r * 1073741824.0;
+ double ra = r + temp - temp;
+ double rb = r - ra;
+
+ rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x;
+
+ temp = ra + rb;
+ rb = -(temp - ra - rb);
+ ra = temp;
+
+ // Call atan
+ double result = atan(ra, rb, x < 0);
+
+ return result;
+ }
+
+ /**
+ * 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)
+ {
+ final double facta = 0.01745329052209854;
+ final double factb = 1.997844754509471E-9;
+
+ double temp = x * 1073741824.0;
+ double xa = x + temp - temp;
+ double xb = x - xa;
+
+ return xb * factb + xb * facta + xa * factb + xa * facta;
+ }
+
+ /**
+ * 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)
+ {
+ final double facta = 57.2957763671875;
+ final double factb = 3.145894820876798E-6;
+
+ double temp = x * 1073741824.0;
+ double xa = x + temp - temp;
+ 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;
+ }
+
+ /**
+ * 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;
+ }
+
+ /**
+ * 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) {
+ return abs(x - Double.longBitsToDouble(Double.doubleToLongBits(x) ^ 1));
+ }
+
+ /**
+ * Get the next machine representable number after a number, moving
+ * in the direction of another number.
+ * <p>
+ * If <code>direction</code> is greater than or equal to<code>d</code>,
+ * the smallest machine representable number strictly greater than
+ * <code>d</code> is returned; otherwise the largest representable number
+ * strictly less than <code>d</code> is returned.</p>
+ * <p>
+ * If <code>d</code> is NaN or Infinite, it is returned unchanged.</p>
+ *
+ * @param d base number
+ * @param direction (the only important thing is whether
+ * direction is greater or smaller than 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.isInfinite(d)) {
+ return d;
+ } 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
+
+ // split the double in raw components
+ long bits = Double.doubleToLongBits(d);
+ long sign = bits & 0x8000000000000000L;
+ long exponent = bits & 0x7ff0000000000000L;
+ long mantissa = bits & 0x000fffffffffffffL;
+
+ if (d * (direction - d) >= 0) {
+ // we should increase the mantissa
+ if (mantissa == 0x000fffffffffffffL) {
+ return Double.longBitsToDouble(sign |
+ (exponent + 0x0010000000000000L));
+ } else {
+ return Double.longBitsToDouble(sign |
+ exponent | (mantissa + 1));
+ }
+ } else {
+ // we should decrease the mantissa
+ if (mantissa == 0L) {
+ return Double.longBitsToDouble(sign |
+ (exponent - 0x0010000000000000L) |
+ 0x000fffffffffffffL);
+ } else {
+ return Double.longBitsToDouble(sign |
+ exponent | (mantissa - 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
[... 141 lines stripped ...]