| /* |
| * 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; |
| } |
| |
| } |