diff options
| author | upstream source tree <ports@midipix.org> | 2015-03-15 20:14:05 -0400 |
|---|---|---|
| committer | upstream source tree <ports@midipix.org> | 2015-03-15 20:14:05 -0400 |
| commit | 554fd8c5195424bdbcabf5de30fdc183aba391bd (patch) | |
| tree | 976dc5ab7fddf506dadce60ae936f43f58787092 /libjava/classpath/java/lang/StrictMath.java | |
| download | cbb-gcc-4.6.4-554fd8c5195424bdbcabf5de30fdc183aba391bd.tar.bz2 cbb-gcc-4.6.4-554fd8c5195424bdbcabf5de30fdc183aba391bd.tar.xz | |
obtained gcc-4.6.4.tar.bz2 from upstream website;upstream
verified gcc-4.6.4.tar.bz2.sig;
imported gcc-4.6.4 source tree from verified upstream tarball.
downloading a git-generated archive based on the 'upstream' tag
should provide you with a source tree that is binary identical
to the one extracted from the above tarball.
if you have obtained the source via the command 'git clone',
however, do note that line-endings of files in your working
directory might differ from line-endings of the respective
files in the upstream repository.
Diffstat (limited to 'libjava/classpath/java/lang/StrictMath.java')
| -rw-r--r-- | libjava/classpath/java/lang/StrictMath.java | 2575 |
1 files changed, 2575 insertions, 0 deletions
diff --git a/libjava/classpath/java/lang/StrictMath.java b/libjava/classpath/java/lang/StrictMath.java new file mode 100644 index 000000000..88f5e5750 --- /dev/null +++ b/libjava/classpath/java/lang/StrictMath.java @@ -0,0 +1,2575 @@ +/* java.lang.StrictMath -- common mathematical functions, strict Java + Copyright (C) 1998, 2001, 2002, 2003, 2006 Free Software Foundation, Inc. + +This file is part of GNU Classpath. + +GNU Classpath is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2, or (at your option) +any later version. + +GNU Classpath is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +General Public License for more details. + +You should have received a copy of the GNU General Public License +along with GNU Classpath; see the file COPYING. If not, write to the +Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +02110-1301 USA. + +Linking this library statically or dynamically with other modules is +making a combined work based on this library. Thus, the terms and +conditions of the GNU General Public License cover the whole +combination. + +As a special exception, the copyright holders of this library give you +permission to link this library with independent modules to produce an +executable, regardless of the license terms of these independent +modules, and to copy and distribute the resulting executable under +terms of your choice, provided that you also meet, for each linked +independent module, the terms and conditions of the license of that +module. An independent module is a module which is not derived from +or based on this library. If you modify this library, you may extend +this exception to your version of the library, but you are not +obligated to do so. If you do not wish to do so, delete this +exception statement from your version. */ + +/* + * Some of the algorithms in this class are in the public domain, as part + * of fdlibm (freely-distributable math library), available at + * http://www.netlib.org/fdlibm/, and carry the following copyright: + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +package java.lang; + +import gnu.classpath.Configuration; + +import java.util.Random; + +/** + * Helper class containing useful mathematical functions and constants. + * This class mirrors {@link Math}, but is 100% portable, because it uses + * no native methods whatsoever. Also, these algorithms are all accurate + * to less than 1 ulp, and execute in <code>strictfp</code> mode, while + * Math is allowed to vary in its results for some functions. Unfortunately, + * this usually means StrictMath has less efficiency and speed, as Math can + * use native methods. + * + * <p>The source of the various algorithms used is the fdlibm library, at:<br> + * <a href="http://www.netlib.org/fdlibm/">http://www.netlib.org/fdlibm/</a> + * + * Note that angles are specified in radians. Conversion functions are + * provided for your convenience. + * + * @author Eric Blake (ebb9@email.byu.edu) + * @since 1.3 + */ +public final strictfp class StrictMath +{ + /** + * StrictMath is non-instantiable. + */ + private StrictMath() + { + } + + /** + * A random number generator, initialized on first use. + * + * @see #random() + */ + private static Random rand; + + /** + * The most accurate approximation to the mathematical constant <em>e</em>: + * <code>2.718281828459045</code>. Used in natural log and exp. + * + * @see #log(double) + * @see #exp(double) + */ + public static final double E + = 2.718281828459045; // Long bits 0x4005bf0z8b145769L. + + /** + * The most accurate approximation to the mathematical constant <em>pi</em>: + * <code>3.141592653589793</code>. This is the ratio of a circle's diameter + * to its circumference. + */ + public static final double PI + = 3.141592653589793; // Long bits 0x400921fb54442d18L. + + /** + * Take the absolute value of the argument. (Absolute value means make + * it positive.) + * + * <p>Note that the the largest negative value (Integer.MIN_VALUE) cannot + * be made positive. In this case, because of the rules of negation in + * a computer, MIN_VALUE is what will be returned. + * This is a <em>negative</em> value. You have been warned. + * + * @param i the number to take the absolute value of + * @return the absolute value + * @see Integer#MIN_VALUE + */ + public static int abs(int i) + { + return (i < 0) ? -i : i; + } + + /** + * Take the absolute value of the argument. (Absolute value means make + * it positive.) + * + * <p>Note that the the largest negative value (Long.MIN_VALUE) cannot + * be made positive. In this case, because of the rules of negation in + * a computer, MIN_VALUE is what will be returned. + * This is a <em>negative</em> value. You have been warned. + * + * @param l the number to take the absolute value of + * @return the absolute value + * @see Long#MIN_VALUE + */ + public static long abs(long l) + { + return (l < 0) ? -l : l; + } + + /** + * Take the absolute value of the argument. (Absolute value means make + * it positive.) + * + * @param f the number to take the absolute value of + * @return the absolute value + */ + public static float abs(float f) + { + return (f <= 0) ? 0 - f : f; + } + + /** + * Take the absolute value of the argument. (Absolute value means make + * it positive.) + * + * @param d the number to take the absolute value of + * @return the absolute value + */ + public static double abs(double d) + { + return (d <= 0) ? 0 - d : d; + } + + /** + * Return whichever argument is smaller. + * + * @param a the first number + * @param b a second number + * @return the smaller of the two numbers + */ + public static int min(int a, int b) + { + return (a < b) ? a : b; + } + + /** + * Return whichever argument is smaller. + * + * @param a the first number + * @param b a second number + * @return the smaller of the two numbers + */ + public static long min(long a, long b) + { + return (a < b) ? a : b; + } + + /** + * Return whichever argument is smaller. If either argument is NaN, the + * result is NaN, and when comparing 0 and -0, -0 is always smaller. + * + * @param a the first number + * @param b a second number + * @return the smaller of the two numbers + */ + public static float min(float a, float b) + { + // this check for NaN, from JLS 15.21.1, saves a method call + if (a != a) + return a; + // no need to check if b is NaN; < will work correctly + // recall that -0.0 == 0.0, but [+-]0.0 - [+-]0.0 behaves special + if (a == 0 && b == 0) + return -(-a - b); + return (a < b) ? a : b; + } + + /** + * Return whichever argument is smaller. If either argument is NaN, the + * result is NaN, and when comparing 0 and -0, -0 is always smaller. + * + * @param a the first number + * @param b a second number + * @return the smaller of the two numbers + */ + public static double min(double a, double b) + { + // this check for NaN, from JLS 15.21.1, saves a method call + if (a != a) + return a; + // no need to check if b is NaN; < will work correctly + // recall that -0.0 == 0.0, but [+-]0.0 - [+-]0.0 behaves special + if (a == 0 && b == 0) + return -(-a - b); + return (a < b) ? a : b; + } + + /** + * Return whichever argument is larger. + * + * @param a the first number + * @param b a second number + * @return the larger of the two numbers + */ + public static int max(int a, int b) + { + return (a > b) ? a : b; + } + + /** + * Return whichever argument is larger. + * + * @param a the first number + * @param b a second number + * @return the larger of the two numbers + */ + public static long max(long a, long b) + { + return (a > b) ? a : b; + } + + /** + * Return whichever argument is larger. If either argument is NaN, the + * result is NaN, and when comparing 0 and -0, 0 is always larger. + * + * @param a the first number + * @param b a second number + * @return the larger of the two numbers + */ + public static float max(float a, float b) + { + // this check for NaN, from JLS 15.21.1, saves a method call + if (a != a) + return a; + // no need to check if b is NaN; > will work correctly + // recall that -0.0 == 0.0, but [+-]0.0 - [+-]0.0 behaves special + if (a == 0 && b == 0) + return a - -b; + return (a > b) ? a : b; + } + + /** + * Return whichever argument is larger. If either argument is NaN, the + * result is NaN, and when comparing 0 and -0, 0 is always larger. + * + * @param a the first number + * @param b a second number + * @return the larger of the two numbers + */ + public static double max(double a, double b) + { + // this check for NaN, from JLS 15.21.1, saves a method call + if (a != a) + return a; + // no need to check if b is NaN; > will work correctly + // recall that -0.0 == 0.0, but [+-]0.0 - [+-]0.0 behaves special + if (a == 0 && b == 0) + return a - -b; + return (a > b) ? a : b; + } + + /** + * The trigonometric function <em>sin</em>. The sine of NaN or infinity is + * NaN, and the sine of 0 retains its sign. + * + * @param a the angle (in radians) + * @return sin(a) + */ + public static double sin(double a) + { + if (a == Double.NEGATIVE_INFINITY || ! (a < Double.POSITIVE_INFINITY)) + return Double.NaN; + + if (abs(a) <= PI / 4) + return sin(a, 0); + + // Argument reduction needed. + double[] y = new double[2]; + int n = remPiOver2(a, y); + switch (n & 3) + { + case 0: + return sin(y[0], y[1]); + case 1: + return cos(y[0], y[1]); + case 2: + return -sin(y[0], y[1]); + default: + return -cos(y[0], y[1]); + } + } + + /** + * The trigonometric function <em>cos</em>. The cosine of NaN or infinity is + * NaN. + * + * @param a the angle (in radians). + * @return cos(a). + */ + public static double cos(double a) + { + if (a == Double.NEGATIVE_INFINITY || ! (a < Double.POSITIVE_INFINITY)) + return Double.NaN; + + if (abs(a) <= PI / 4) + return cos(a, 0); + + // Argument reduction needed. + double[] y = new double[2]; + int n = remPiOver2(a, y); + switch (n & 3) + { + case 0: + return cos(y[0], y[1]); + case 1: + return -sin(y[0], y[1]); + case 2: + return -cos(y[0], y[1]); + default: + return sin(y[0], y[1]); + } + } + + /** + * The trigonometric function <em>tan</em>. The tangent of NaN or infinity + * is NaN, and the tangent of 0 retains its sign. + * + * @param a the angle (in radians) + * @return tan(a) + */ + public static double tan(double a) + { + if (a == Double.NEGATIVE_INFINITY || ! (a < Double.POSITIVE_INFINITY)) + return Double.NaN; + + if (abs(a) <= PI / 4) + return tan(a, 0, false); + + // Argument reduction needed. + double[] y = new double[2]; + int n = remPiOver2(a, y); + return tan(y[0], y[1], (n & 1) == 1); + } + + /** + * The trigonometric function <em>arcsin</em>. The range of angles returned + * is -pi/2 to pi/2 radians (-90 to 90 degrees). If the argument is NaN or + * its absolute value is beyond 1, the result is NaN; and the arcsine of + * 0 retains its sign. + * + * @param x the sin to turn back into an angle + * @return arcsin(x) + */ + public static double asin(double x) + { + boolean negative = x < 0; + if (negative) + x = -x; + if (! (x <= 1)) + return Double.NaN; + if (x == 1) + return negative ? -PI / 2 : PI / 2; + if (x < 0.5) + { + if (x < 1 / TWO_27) + return negative ? -x : x; + double t = x * x; + double p = t * (PS0 + t * (PS1 + t * (PS2 + t * (PS3 + t + * (PS4 + t * PS5))))); + double q = 1 + t * (QS1 + t * (QS2 + t * (QS3 + t * QS4))); + return negative ? -x - x * (p / q) : x + x * (p / q); + } + double w = 1 - x; // 1>|x|>=0.5. + double t = w * 0.5; + double p = t * (PS0 + t * (PS1 + t * (PS2 + t * (PS3 + t + * (PS4 + t * PS5))))); + double q = 1 + t * (QS1 + t * (QS2 + t * (QS3 + t * QS4))); + double s = sqrt(t); + if (x >= 0.975) + { + w = p / q; + t = PI / 2 - (2 * (s + s * w) - PI_L / 2); + } + else + { + w = (float) s; + double c = (t - w * w) / (s + w); + p = 2 * s * (p / q) - (PI_L / 2 - 2 * c); + q = PI / 4 - 2 * w; + t = PI / 4 - (p - q); + } + return negative ? -t : t; + } + + /** + * The trigonometric function <em>arccos</em>. The range of angles returned + * is 0 to pi radians (0 to 180 degrees). If the argument is NaN or + * its absolute value is beyond 1, the result is NaN. + * + * @param x the cos to turn back into an angle + * @return arccos(x) + */ + public static double acos(double x) + { + boolean negative = x < 0; + if (negative) + x = -x; + if (! (x <= 1)) + return Double.NaN; + if (x == 1) + return negative ? PI : 0; + if (x < 0.5) + { + if (x < 1 / TWO_57) + return PI / 2; + double z = x * x; + double p = z * (PS0 + z * (PS1 + z * (PS2 + z * (PS3 + z + * (PS4 + z * PS5))))); + double q = 1 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4))); + double r = x - (PI_L / 2 - x * (p / q)); + return negative ? PI / 2 + r : PI / 2 - r; + } + if (negative) // x<=-0.5. + { + double z = (1 + x) * 0.5; + double p = z * (PS0 + z * (PS1 + z * (PS2 + z * (PS3 + z + * (PS4 + z * PS5))))); + double q = 1 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4))); + double s = sqrt(z); + double w = p / q * s - PI_L / 2; + return PI - 2 * (s + w); + } + double z = (1 - x) * 0.5; // x>0.5. + double s = sqrt(z); + double df = (float) s; + double c = (z - df * df) / (s + df); + double p = z * (PS0 + z * (PS1 + z * (PS2 + z * (PS3 + z + * (PS4 + z * PS5))))); + double q = 1 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4))); + double w = p / q * s + c; + return 2 * (df + w); + } + + /** + * The trigonometric function <em>arcsin</em>. The range of angles returned + * is -pi/2 to pi/2 radians (-90 to 90 degrees). If the argument is NaN, the + * result is NaN; and the arctangent of 0 retains its sign. + * + * @param x the tan to turn back into an angle + * @return arcsin(x) + * @see #atan2(double, double) + */ + public static double atan(double x) + { + double lo; + double hi; + boolean negative = x < 0; + if (negative) + x = -x; + if (x >= TWO_66) + return negative ? -PI / 2 : PI / 2; + if (! (x >= 0.4375)) // |x|<7/16, or NaN. + { + if (! (x >= 1 / TWO_29)) // Small, or NaN. + return negative ? -x : x; + lo = hi = 0; + } + else if (x < 1.1875) + { + if (x < 0.6875) // 7/16<=|x|<11/16. + { + x = (2 * x - 1) / (2 + x); + hi = ATAN_0_5H; + lo = ATAN_0_5L; + } + else // 11/16<=|x|<19/16. + { + x = (x - 1) / (x + 1); + hi = PI / 4; + lo = PI_L / 4; + } + } + else if (x < 2.4375) // 19/16<=|x|<39/16. + { + x = (x - 1.5) / (1 + 1.5 * x); + hi = ATAN_1_5H; + lo = ATAN_1_5L; + } + else // 39/16<=|x|<2**66. + { + x = -1 / x; + hi = PI / 2; + lo = PI_L / 2; + } + + // Break sum from i=0 to 10 ATi*z**(i+1) into odd and even poly. + double z = x * x; + double w = z * z; + double s1 = z * (AT0 + w * (AT2 + w * (AT4 + w * (AT6 + w + * (AT8 + w * AT10))))); + double s2 = w * (AT1 + w * (AT3 + w * (AT5 + w * (AT7 + w * AT9)))); + if (hi == 0) + return negative ? x * (s1 + s2) - x : x - x * (s1 + s2); + z = hi - ((x * (s1 + s2) - lo) - x); + return negative ? -z : z; + } + + /** + * A special version of the trigonometric function <em>arctan</em>, for + * converting rectangular coordinates <em>(x, y)</em> to polar + * <em>(r, theta)</em>. This computes the arctangent of x/y in the range + * of -pi to pi radians (-180 to 180 degrees). Special cases:<ul> + * <li>If either argument is NaN, the result is NaN.</li> + * <li>If the first argument is positive zero and the second argument is + * positive, or the first argument is positive and finite and the second + * argument is positive infinity, then the result is positive zero.</li> + * <li>If the first argument is negative zero and the second argument is + * positive, or the first argument is negative and finite and the second + * argument is positive infinity, then the result is negative zero.</li> + * <li>If the first argument is positive zero and the second argument is + * negative, or the first argument is positive and finite and the second + * argument is negative infinity, then the result is the double value + * closest to pi.</li> + * <li>If the first argument is negative zero and the second argument is + * negative, or the first argument is negative and finite and the second + * argument is negative infinity, then the result is the double value + * closest to -pi.</li> + * <li>If the first argument is positive and the second argument is + * positive zero or negative zero, or the first argument is positive + * infinity and the second argument is finite, then the result is the + * double value closest to pi/2.</li> + * <li>If the first argument is negative and the second argument is + * positive zero or negative zero, or the first argument is negative + * infinity and the second argument is finite, then the result is the + * double value closest to -pi/2.</li> + * <li>If both arguments are positive infinity, then the result is the + * double value closest to pi/4.</li> + * <li>If the first argument is positive infinity and the second argument + * is negative infinity, then the result is the double value closest to + * 3*pi/4.</li> + * <li>If the first argument is negative infinity and the second argument + * is positive infinity, then the result is the double value closest to + * -pi/4.</li> + * <li>If both arguments are negative infinity, then the result is the + * double value closest to -3*pi/4.</li> + * + * </ul><p>This returns theta, the angle of the point. To get r, albeit + * slightly inaccurately, use sqrt(x*x+y*y). + * + * @param y the y position + * @param x the x position + * @return <em>theta</em> in the conversion of (x, y) to (r, theta) + * @see #atan(double) + */ + public static double atan2(double y, double x) + { + if (x != x || y != y) + return Double.NaN; + if (x == 1) + return atan(y); + if (x == Double.POSITIVE_INFINITY) + { + if (y == Double.POSITIVE_INFINITY) + return PI / 4; + if (y == Double.NEGATIVE_INFINITY) + return -PI / 4; + return 0 * y; + } + if (x == Double.NEGATIVE_INFINITY) + { + if (y == Double.POSITIVE_INFINITY) + return 3 * PI / 4; + if (y == Double.NEGATIVE_INFINITY) + return -3 * PI / 4; + return (1 / (0 * y) == Double.POSITIVE_INFINITY) ? PI : -PI; + } + if (y == 0) + { + if (1 / (0 * x) == Double.POSITIVE_INFINITY) + return y; + return (1 / y == Double.POSITIVE_INFINITY) ? PI : -PI; + } + if (y == Double.POSITIVE_INFINITY || y == Double.NEGATIVE_INFINITY + || x == 0) + return y < 0 ? -PI / 2 : PI / 2; + + double z = abs(y / x); // Safe to do y/x. + if (z > TWO_60) + z = PI / 2 + 0.5 * PI_L; + else if (x < 0 && z < 1 / TWO_60) + z = 0; + else + z = atan(z); + if (x > 0) + return y > 0 ? z : -z; + return y > 0 ? PI - (z - PI_L) : z - PI_L - PI; + } + + /** + * Returns the hyperbolic sine of <code>x</code> which is defined as + * (exp(x) - exp(-x)) / 2. + * + * Special cases: + * <ul> + * <li>If the argument is NaN, the result is NaN</li> + * <li>If the argument is positive infinity, the result is positive + * infinity.</li> + * <li>If the argument is negative infinity, the result is negative + * infinity.</li> + * <li>If the argument is zero, the result is zero.</li> + * </ul> + * + * @param x the argument to <em>sinh</em> + * @return the hyperbolic sine of <code>x</code> + * + * @since 1.5 + */ + public static double sinh(double x) + { + // Method : + // mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2 + // 1. Replace x by |x| (sinh(-x) = -sinh(x)). + // 2. + // E + E/(E+1) + // 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x) + // 2 + // + // 22 <= x <= lnovft : sinh(x) := exp(x)/2 + // lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2) + // ln2ovft < x : sinh(x) := +inf (overflow) + + double t, w, h; + + long bits; + long h_bits; + long l_bits; + + // handle special cases + if (x != x) + return x; + if (x == Double.POSITIVE_INFINITY) + return Double.POSITIVE_INFINITY; + if (x == Double.NEGATIVE_INFINITY) + return Double.NEGATIVE_INFINITY; + + if (x < 0) + h = - 0.5; + else + h = 0.5; + + bits = Double.doubleToLongBits(x); + h_bits = getHighDWord(bits) & 0x7fffffffL; // ignore sign + l_bits = getLowDWord(bits); + + // |x| in [0, 22], return sign(x) * 0.5 * (E+E/(E+1)) + if (h_bits < 0x40360000L) // |x| < 22 + { + if (h_bits < 0x3e300000L) // |x| < 2^-28 + return x; // for tiny arguments return x + + t = expm1(abs(x)); + + if (h_bits < 0x3ff00000L) + return h * (2.0 * t - t * t / (t + 1.0)); + + return h * (t + t / (t + 1.0)); + } + + // |x| in [22, log(Double.MAX_VALUE)], return 0.5 * exp(|x|) + if (h_bits < 0x40862e42L) + return h * exp(abs(x)); + + // |x| in [log(Double.MAX_VALUE), overflowthreshold] + if ((h_bits < 0x408633ceL) + || ((h_bits == 0x408633ceL) && (l_bits <= 0x8fb9f87dL))) + { + w = exp(0.5 * abs(x)); + t = h * w; + + return t * w; + } + + // |x| > overflowthershold + return h * Double.POSITIVE_INFINITY; + } + + /** + * Returns the hyperbolic cosine of <code>x</code>, which is defined as + * (exp(x) + exp(-x)) / 2. + * + * Special cases: + * <ul> + * <li>If the argument is NaN, the result is NaN</li> + * <li>If the argument is positive infinity, the result is positive + * infinity.</li> + * <li>If the argument is negative infinity, the result is positive + * infinity.</li> + * <li>If the argument is zero, the result is one.</li> + * </ul> + * + * @param x the argument to <em>cosh</em> + * @return the hyperbolic cosine of <code>x</code> + * + * @since 1.5 + */ + public static double cosh(double x) + { + // Method : + // mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2 + // 1. Replace x by |x| (cosh(x) = cosh(-x)). + // 2. + // [ exp(x) - 1 ]^2 + // 0 <= x <= ln2/2 : cosh(x) := 1 + ------------------- + // 2*exp(x) + // + // exp(x) + 1/exp(x) + // ln2/2 <= x <= 22 : cosh(x) := ------------------ + // 2 + // 22 <= x <= lnovft : cosh(x) := exp(x)/2 + // lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2) + // ln2ovft < x : cosh(x) := +inf (overflow) + + double t, w; + long bits; + long hx; + long lx; + + // handle special cases + if (x != x) + return x; + if (x == Double.POSITIVE_INFINITY) + return Double.POSITIVE_INFINITY; + if (x == Double.NEGATIVE_INFINITY) + return Double.POSITIVE_INFINITY; + + bits = Double.doubleToLongBits(x); + hx = getHighDWord(bits) & 0x7fffffffL; // ignore sign + lx = getLowDWord(bits); + + // |x| in [0, 0.5 * ln(2)], return 1 + expm1(|x|)^2 / (2 * exp(|x|)) + if (hx < 0x3fd62e43L) + { + t = expm1(abs(x)); + w = 1.0 + t; + + // for tiny arguments return 1. + if (hx < 0x3c800000L) + return w; + + return 1.0 + (t * t) / (w + w); + } + + // |x| in [0.5 * ln(2), 22], return exp(|x|)/2 + 1 / (2 * exp(|x|)) + if (hx < 0x40360000L) + { + t = exp(abs(x)); + + return 0.5 * t + 0.5 / t; + } + + // |x| in [22, log(Double.MAX_VALUE)], return 0.5 * exp(|x|) + if (hx < 0x40862e42L) + return 0.5 * exp(abs(x)); + + // |x| in [log(Double.MAX_VALUE), overflowthreshold], + // return exp(x/2)/2 * exp(x/2) + if ((hx < 0x408633ceL) + || ((hx == 0x408633ceL) && (lx <= 0x8fb9f87dL))) + { + w = exp(0.5 * abs(x)); + t = 0.5 * w; + + return t * w; + } + + // |x| > overflowthreshold + return Double.POSITIVE_INFINITY; + } + + /** + * Returns the hyperbolic tangent of <code>x</code>, which is defined as + * (exp(x) - exp(-x)) / (exp(x) + exp(-x)), i.e. sinh(x) / cosh(x). + * + Special cases: + * <ul> + * <li>If the argument is NaN, the result is NaN</li> + * <li>If the argument is positive infinity, the result is 1.</li> + * <li>If the argument is negative infinity, the result is -1.</li> + * <li>If the argument is zero, the result is zero.</li> + * </ul> + * + * @param x the argument to <em>tanh</em> + * @return the hyperbolic tagent of <code>x</code> + * + * @since 1.5 + */ + public static double tanh(double x) + { + // Method : + // 0. tanh(x) is defined to be (exp(x) - exp(-x)) / (exp(x) + exp(-x)) + // 1. reduce x to non-negative by tanh(-x) = -tanh(x). + // 2. 0 <= x <= 2^-55 : tanh(x) := x * (1.0 + x) + // -t + // 2^-55 < x <= 1 : tanh(x) := -----; t = expm1(-2x) + // t + 2 + // 2 + // 1 <= x <= 22.0 : tanh(x) := 1 - ----- ; t=expm1(2x) + // t + 2 + // 22.0 < x <= INF : tanh(x) := 1. + + double t, z; + + long bits; + long h_bits; + + // handle special cases + if (x != x) + return x; + if (x == Double.POSITIVE_INFINITY) + return 1.0; + if (x == Double.NEGATIVE_INFINITY) + return -1.0; + + bits = Double.doubleToLongBits(x); + h_bits = getHighDWord(bits) & 0x7fffffffL; // ingnore sign + + if (h_bits < 0x40360000L) // |x| < 22 + { + if (h_bits < 0x3c800000L) // |x| < 2^-55 + return x * (1.0 + x); + + if (h_bits >= 0x3ff00000L) // |x| >= 1 + { + t = expm1(2.0 * abs(x)); + z = 1.0 - 2.0 / (t + 2.0); + } + else // |x| < 1 + { + t = expm1(-2.0 * abs(x)); + z = -t / (t + 2.0); + } + } + else // |x| >= 22 + z = 1.0; + + return (x >= 0) ? z : -z; + } + + /** + * Returns the lower two words of a long. This is intended to be + * used like this: + * <code>getLowDWord(Double.doubleToLongBits(x))</code>. + */ + private static long getLowDWord(long x) + { + return x & 0x00000000ffffffffL; + } + + /** + * Returns the higher two words of a long. This is intended to be + * used like this: + * <code>getHighDWord(Double.doubleToLongBits(x))</code>. + */ + private static long getHighDWord(long x) + { + return (x & 0xffffffff00000000L) >> 32; + } + + /** + * Returns a double with the IEEE754 bit pattern given in the lower + * and higher two words <code>lowDWord</code> and <code>highDWord</code>. + */ + private static double buildDouble(long lowDWord, long highDWord) + { + return Double.longBitsToDouble(((highDWord & 0xffffffffL) << 32) + | (lowDWord & 0xffffffffL)); + } + + /** + * Returns the cube root of <code>x</code>. The sign of the cube root + * is equal to the sign of <code>x</code>. + * + * Special cases: + * <ul> + * <li>If the argument is NaN, the result is NaN</li> + * <li>If the argument is positive infinity, the result is positive + * infinity.</li> + * <li>If the argument is negative infinity, the result is negative + * infinity.</li> + * <li>If the argument is zero, the result is zero with the same + * sign as the argument.</li> + * </ul> + * + * @param x the number to take the cube root of + * @return the cube root of <code>x</code> + * @see #sqrt(double) + * + * @since 1.5 + */ + public static double cbrt(double x) + { + boolean negative = (x < 0); + double r; + double s; + double t; + double w; + + long bits; + long l; + long h; + + // handle the special cases + if (x != x) + return x; + if (x == Double.POSITIVE_INFINITY) + return Double.POSITIVE_INFINITY; + if (x == Double.NEGATIVE_INFINITY) + return Double.NEGATIVE_INFINITY; + if (x == 0) + return x; + + x = abs(x); + bits = Double.doubleToLongBits(x); + + if (bits < 0x0010000000000000L) // subnormal number + { + t = TWO_54; + t *= x; + + // __HI(t)=__HI(t)/3+B2; + bits = Double.doubleToLongBits(t); + h = getHighDWord(bits); + l = getLowDWord(bits); + + h = h / 3 + CBRT_B2; + + t = buildDouble(l, h); + } + else + { + // __HI(t)=__HI(x)/3+B1; + h = getHighDWord(bits); + l = 0; + + h = h / 3 + CBRT_B1; + t = buildDouble(l, h); + } + + // new cbrt to 23 bits + r = t * t / x; + s = CBRT_C + r * t; + t *= CBRT_G + CBRT_F / (s + CBRT_E + CBRT_D / s); + + // chopped to 20 bits and make it larger than cbrt(x) + bits = Double.doubleToLongBits(t); + h = getHighDWord(bits); + + // __LO(t)=0; + // __HI(t)+=0x00000001; + l = 0; + h += 1; + t = buildDouble(l, h); + + // one step newton iteration to 53 bits with error less than 0.667 ulps + s = t * t; // t * t is exact + r = x / s; + w = t + t; + r = (r - t) / (w + r); // r - t is exact + t = t + t * r; + + return negative ? -t : t; + } + + /** + * Take <em>e</em><sup>a</sup>. The opposite of <code>log()</code>. If the + * argument is NaN, the result is NaN; if the argument is positive infinity, + * the result is positive infinity; and if the argument is negative + * infinity, the result is positive zero. + * + * @param x the number to raise to the power + * @return the number raised to the power of <em>e</em> + * @see #log(double) + * @see #pow(double, double) + */ + public static double exp(double x) + { + if (x != x) + return x; + if (x > EXP_LIMIT_H) + return Double.POSITIVE_INFINITY; + if (x < EXP_LIMIT_L) + return 0; + + // Argument reduction. + double hi; + double lo; + int k; + double t = abs(x); + if (t > 0.5 * LN2) + { + if (t < 1.5 * LN2) + { + hi = t - LN2_H; + lo = LN2_L; + k = 1; + } + else + { + k = (int) (INV_LN2 * t + 0.5); + hi = t - k * LN2_H; + lo = k * LN2_L; + } + if (x < 0) + { + hi = -hi; + lo = -lo; + k = -k; + } + x = hi - lo; + } + else if (t < 1 / TWO_28) + return 1; + else + lo = hi = k = 0; + + // Now x is in primary range. + t = x * x; + double c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); + if (k == 0) + return 1 - (x * c / (c - 2) - x); + double y = 1 - (lo - x * c / (2 - c) - hi); + return scale(y, k); + } + + /** + * Returns <em>e</em><sup>x</sup> - 1. + * Special cases: + * <ul> + * <li>If the argument is NaN, the result is NaN.</li> + * <li>If the argument is positive infinity, the result is positive + * infinity</li> + * <li>If the argument is negative infinity, the result is -1.</li> + * <li>If the argument is zero, the result is zero.</li> + * </ul> + * + * @param x the argument to <em>e</em><sup>x</sup> - 1. + * @return <em>e</em> raised to the power <code>x</code> minus one. + * @see #exp(double) + */ + public static double expm1(double x) + { + // Method + // 1. Argument reduction: + // Given x, find r and integer k such that + // + // x = k * ln(2) + r, |r| <= 0.5 * ln(2) + // + // Here a correction term c will be computed to compensate + // the error in r when rounded to a floating-point number. + // + // 2. Approximating expm1(r) by a special rational function on + // the interval [0, 0.5 * ln(2)]: + // Since + // r*(exp(r)+1)/(exp(r)-1) = 2 + r^2/6 - r^4/360 + ... + // we define R1(r*r) by + // r*(exp(r)+1)/(exp(r)-1) = 2 + r^2/6 * R1(r*r) + // That is, + // R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r) + // = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r)) + // = 1 - r^2/60 + r^4/2520 - r^6/100800 + ... + // We use a special Remes algorithm on [0, 0.347] to generate + // a polynomial of degree 5 in r*r to approximate R1. The + // maximum error of this polynomial approximation is bounded + // by 2**-61. In other words, + // R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5 + // where Q1 = -1.6666666666666567384E-2, + // Q2 = 3.9682539681370365873E-4, + // Q3 = -9.9206344733435987357E-6, + // Q4 = 2.5051361420808517002E-7, + // Q5 = -6.2843505682382617102E-9; + // (where z=r*r, and Q1 to Q5 are called EXPM1_Qx in the source) + // with error bounded by + // | 5 | -61 + // | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2 + // | | + // + // expm1(r) = exp(r)-1 is then computed by the following + // specific way which minimize the accumulation rounding error: + // 2 3 + // r r [ 3 - (R1 + R1*r/2) ] + // expm1(r) = r + --- + --- * [--------------------] + // 2 2 [ 6 - r*(3 - R1*r/2) ] + // + // To compensate the error in the argument reduction, we use + // expm1(r+c) = expm1(r) + c + expm1(r)*c + // ~ expm1(r) + c + r*c + // Thus c+r*c will be added in as the correction terms for + // expm1(r+c). Now rearrange the term to avoid optimization + // screw up: + // ( 2 2 ) + // ({ ( r [ R1 - (3 - R1*r/2) ] ) } r ) + // expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- ) + // ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 ) + // ( ) + // + // = r - E + // 3. Scale back to obtain expm1(x): + // From step 1, we have + // expm1(x) = either 2^k*[expm1(r)+1] - 1 + // = or 2^k*[expm1(r) + (1-2^-k)] + // 4. Implementation notes: + // (A). To save one multiplication, we scale the coefficient Qi + // to Qi*2^i, and replace z by (x^2)/2. + // (B). To achieve maximum accuracy, we compute expm1(x) by + // (i) if x < -56*ln2, return -1.0, (raise inexact if x!=inf) + // (ii) if k=0, return r-E + // (iii) if k=-1, return 0.5*(r-E)-0.5 + // (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E) + // else return 1.0+2.0*(r-E); + // (v) if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1) + // (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else + // (vii) return 2^k(1-((E+2^-k)-r)) + + boolean negative = (x < 0); + double y, hi, lo, c, t, e, hxs, hfx, r1; + int k; + + long bits; + long h_bits; + long l_bits; + + c = 0.0; + y = abs(x); + + bits = Double.doubleToLongBits(y); + h_bits = getHighDWord(bits); + l_bits = getLowDWord(bits); + + // handle special cases and large arguments + if (h_bits >= 0x4043687aL) // if |x| >= 56 * ln(2) + { + if (h_bits >= 0x40862e42L) // if |x| >= EXP_LIMIT_H + { + if (h_bits >= 0x7ff00000L) + { + if (((h_bits & 0x000fffffL) | (l_bits & 0xffffffffL)) != 0) + return x; // exp(NaN) = NaN + else + return negative ? -1.0 : x; // exp({+-inf}) = {+inf, -1} + } + + if (x > EXP_LIMIT_H) + return Double.POSITIVE_INFINITY; // overflow + } + + if (negative) // x <= -56 * ln(2) + return -1.0; + } + + // argument reduction + if (h_bits > 0x3fd62e42L) // |x| > 0.5 * ln(2) + { + if (h_bits < 0x3ff0a2b2L) // |x| < 1.5 * ln(2) + { + if (negative) + { + hi = x + LN2_H; + lo = -LN2_L; + k = -1; + } + else + { + hi = x - LN2_H; + lo = LN2_L; + k = 1; + } + } + else + { + k = (int) (INV_LN2 * x + (negative ? - 0.5 : 0.5)); + t = k; + hi = x - t * LN2_H; + lo = t * LN2_L; + } + + x = hi - lo; + c = (hi - x) - lo; + + } + else if (h_bits < 0x3c900000L) // |x| < 2^-54 return x + return x; + else + k = 0; + + // x is now in primary range + hfx = 0.5 * x; + hxs = x * hfx; + r1 = 1.0 + hxs * (EXPM1_Q1 + + hxs * (EXPM1_Q2 + + hxs * (EXPM1_Q3 + + hxs * (EXPM1_Q4 + + hxs * EXPM1_Q5)))); + t = 3.0 - r1 * hfx; + e = hxs * ((r1 - t) / (6.0 - x * t)); + + if (k == 0) + { + return x - (x * e - hxs); // c == 0 + } + else + { + e = x * (e - c) - c; + e -= hxs; + + if (k == -1) + return 0.5 * (x - e) - 0.5; + + if (k == 1) + { + if (x < - 0.25) + return -2.0 * (e - (x + 0.5)); + else + return 1.0 + 2.0 * (x - e); + } + + if (k <= -2 || k > 56) // sufficient to return exp(x) - 1 + { + y = 1.0 - (e - x); + + bits = Double.doubleToLongBits(y); + h_bits = getHighDWord(bits); + l_bits = getLowDWord(bits); + + h_bits += (k << 20); // add k to y's exponent + + y = buildDouble(l_bits, h_bits); + + return y - 1.0; + } + + t = 1.0; + if (k < 20) + { + bits = Double.doubleToLongBits(t); + h_bits = 0x3ff00000L - (0x00200000L >> k); + l_bits = getLowDWord(bits); + + t = buildDouble(l_bits, h_bits); // t = 1 - 2^(-k) + y = t - (e - x); + + bits = Double.doubleToLongBits(y); + h_bits = getHighDWord(bits); + l_bits = getLowDWord(bits); + + h_bits += (k << 20); // add k to y's exponent + + y = buildDouble(l_bits, h_bits); + } + else + { + bits = Double.doubleToLongBits(t); + h_bits = (0x000003ffL - k) << 20; + l_bits = getLowDWord(bits); + + t = buildDouble(l_bits, h_bits); // t = 2^(-k) + + y = x - (e + t); + y += 1.0; + + bits = Double.doubleToLongBits(y); + h_bits = getHighDWord(bits); + l_bits = getLowDWord(bits); + + h_bits += (k << 20); // add k to y's exponent + + y = buildDouble(l_bits, h_bits); + } + } + + return y; + } + + /** + * Take ln(a) (the natural log). The opposite of <code>exp()</code>. If the + * argument is NaN or negative, the result is NaN; if the argument is + * positive infinity, the result is positive infinity; and if the argument + * is either zero, the result is negative infinity. + * + * <p>Note that the way to get log<sub>b</sub>(a) is to do this: + * <code>ln(a) / ln(b)</code>. + * + * @param x the number to take the natural log of + * @return the natural log of <code>a</code> + * @see #exp(double) + */ + public static double log(double x) + { + if (x == 0) + return Double.NEGATIVE_INFINITY; + if (x < 0) + return Double.NaN; + if (! (x < Double.POSITIVE_INFINITY)) + return x; + + // Normalize x. + long bits = Double.doubleToLongBits(x); + int exp = (int) (bits >> 52); + if (exp == 0) // Subnormal x. + { + x *= TWO_54; + bits = Double.doubleToLongBits(x); + exp = (int) (bits >> 52) - 54; + } + exp -= 1023; // Unbias exponent. + bits = (bits & 0x000fffffffffffffL) | 0x3ff0000000000000L; + x = Double.longBitsToDouble(bits); + if (x >= SQRT_2) + { + x *= 0.5; + exp++; + } + x--; + if (abs(x) < 1 / TWO_20) + { + if (x == 0) + return exp * LN2_H + exp * LN2_L; + double r = x * x * (0.5 - 1 / 3.0 * x); + if (exp == 0) + return x - r; + return exp * LN2_H - ((r - exp * LN2_L) - x); + } + double s = x / (2 + x); + double z = s * s; + double w = z * z; + double t1 = w * (LG2 + w * (LG4 + w * LG6)); + double t2 = z * (LG1 + w * (LG3 + w * (LG5 + w * LG7))); + double r = t2 + t1; + if (bits >= 0x3ff6174a00000000L && bits < 0x3ff6b85200000000L) + { + double h = 0.5 * x * x; // Need more accuracy for x near sqrt(2). + if (exp == 0) + return x - (h - s * (h + r)); + return exp * LN2_H - ((h - (s * (h + r) + exp * LN2_L)) - x); + } + if (exp == 0) + return x - s * (x - r); + return exp * LN2_H - ((s * (x - r) - exp * LN2_L) - x); + } + + /** + * Take a square root. If the argument is NaN or negative, the result is + * NaN; if the argument is positive infinity, the result is positive + * infinity; and if the result is either zero, the result is the same. + * + * <p>For other roots, use pow(x, 1/rootNumber). + * + * @param x the numeric argument + * @return the square root of the argument + * @see #pow(double, double) + */ + public static double sqrt(double x) + { + if (x < 0) + return Double.NaN; + if (x == 0 || ! (x < Double.POSITIVE_INFINITY)) + return x; + + // Normalize x. + long bits = Double.doubleToLongBits(x); + int exp = (int) (bits >> 52); + if (exp == 0) // Subnormal x. + { + x *= TWO_54; + bits = Double.doubleToLongBits(x); + exp = (int) (bits >> 52) - 54; + } + exp -= 1023; // Unbias exponent. + bits = (bits & 0x000fffffffffffffL) | 0x0010000000000000L; + if ((exp & 1) == 1) // Odd exp, double x to make it even. + bits <<= 1; + exp >>= 1; + + // Generate sqrt(x) bit by bit. + bits <<= 1; + long q = 0; + long s = 0; + long r = 0x0020000000000000L; // Move r right to left. + while (r != 0) + { + long t = s + r; + if (t <= bits) + { + s = t + r; + bits -= t; + q += r; + } + bits <<= 1; + r >>= 1; + } + + // Use floating add to round correctly. + if (bits != 0) + q += q & 1; + return Double.longBitsToDouble((q >> 1) + ((exp + 1022L) << 52)); + } + + /** + * Raise a number to a power. Special cases:<ul> + * <li>If the second argument is positive or negative zero, then the result + * is 1.0.</li> + * <li>If the second argument is 1.0, then the result is the same as the + * first argument.</li> + * <li>If the second argument is NaN, then the result is NaN.</li> + * <li>If the first argument is NaN and the second argument is nonzero, + * then the result is NaN.</li> + * <li>If the absolute value of the first argument is greater than 1 and + * the second argument is positive infinity, or the absolute value of the + * first argument is less than 1 and the second argument is negative + * infinity, then the result is positive infinity.</li> + * <li>If the absolute value of the first argument is greater than 1 and + * the second argument is negative infinity, or the absolute value of the + * first argument is less than 1 and the second argument is positive + * infinity, then the result is positive zero.</li> + * <li>If the absolute value of the first argument equals 1 and the second + * argument is infinite, then the result is NaN.</li> + * <li>If the first argument is positive zero and the second argument is + * greater than zero, or the first argument is positive infinity and the + * second argument is less than zero, then the result is positive zero.</li> + * <li>If the first argument is positive zero and the second argument is + * less than zero, or the first argument is positive infinity and the + * second argument is greater than zero, then the result is positive + * infinity.</li> + * <li>If the first argument is negative zero and the second argument is + * greater than zero but not a finite odd integer, or the first argument is + * negative infinity and the second argument is less than zero but not a + * finite odd integer, then the result is positive zero.</li> + * <li>If the first argument is negative zero and the second argument is a + * positive finite odd integer, or the first argument is negative infinity + * and the second argument is a negative finite odd integer, then the result + * is negative zero.</li> + * <li>If the first argument is negative zero and the second argument is + * less than zero but not a finite odd integer, or the first argument is + * negative infinity and the second argument is greater than zero but not a + * finite odd integer, then the result is positive infinity.</li> + * <li>If the first argument is negative zero and the second argument is a + * negative finite odd integer, or the first argument is negative infinity + * and the second argument is a positive finite odd integer, then the result + * is negative infinity.</li> + * <li>If the first argument is less than zero and the second argument is a + * finite even integer, then the result is equal to the result of raising + * the absolute value of the first argument to the power of the second + * argument.</li> + * <li>If the first argument is less than zero and the second argument is a + * finite odd integer, then the result is equal to the negative of the + * result of raising the absolute value of the first argument to the power + * of the second argument.</li> + * <li>If the first argument is finite and less than zero and the second + * argument is finite and not an integer, then the result is NaN.</li> + * <li>If both arguments are integers, then the result is exactly equal to + * the mathematical result of raising the first argument to the power of + * the second argument if that result can in fact be represented exactly as + * a double value.</li> + * + * </ul><p>(In the foregoing descriptions, a floating-point value is + * considered to be an integer if and only if it is a fixed point of the + * method {@link #ceil(double)} or, equivalently, a fixed point of the + * method {@link #floor(double)}. A value is a fixed point of a one-argument + * method if and only if the result of applying the method to the value is + * equal to the value.) + * + * @param x the number to raise + * @param y the power to raise it to + * @return x<sup>y</sup> + */ + public static double pow(double x, double y) + { + // Special cases first. + if (y == 0) + return 1; + if (y == 1) + return x; + if (y == -1) + return 1 / x; + if (x != x || y != y) + return Double.NaN; + + // When x < 0, yisint tells if y is not an integer (0), even(1), + // or odd (2). + int yisint = 0; + if (x < 0 && floor(y) == y) + yisint = (y % 2 == 0) ? 2 : 1; + double ax = abs(x); + double ay = abs(y); + + // More special cases, of y. + if (ay == Double.POSITIVE_INFINITY) + { + if (ax == 1) + return Double.NaN; + if (ax > 1) + return y > 0 ? y : 0; + return y < 0 ? -y : 0; + } + if (y == 2) + return x * x; + if (y == 0.5) + return sqrt(x); + + // More special cases, of x. + if (x == 0 || ax == Double.POSITIVE_INFINITY || ax == 1) + { + if (y < 0) + ax = 1 / ax; + if (x < 0) + { + if (x == -1 && yisint == 0) + ax = Double.NaN; + else if (yisint == 1) + ax = -ax; + } + return ax; + } + if (x < 0 && yisint == 0) + return Double.NaN; + + // Now we can start! + double t; + double t1; + double t2; + double u; + double v; + double w; + if (ay > TWO_31) + { + if (ay > TWO_64) // Automatic over/underflow. + return ((ax < 1) ? y < 0 : y > 0) ? Double.POSITIVE_INFINITY : 0; + // Over/underflow if x is not close to one. + if (ax < 0.9999995231628418) + return y < 0 ? Double.POSITIVE_INFINITY : 0; + if (ax >= 1.0000009536743164) + return y > 0 ? Double.POSITIVE_INFINITY : 0; + // Now |1-x| is <= 2**-20, sufficient to compute + // log(x) by x-x^2/2+x^3/3-x^4/4. + t = x - 1; + w = t * t * (0.5 - t * (1 / 3.0 - t * 0.25)); + u = INV_LN2_H * t; + v = t * INV_LN2_L - w * INV_LN2; + t1 = (float) (u + v); + t2 = v - (t1 - u); + } + else + { + long bits = Double.doubleToLongBits(ax); + int exp = (int) (bits >> 52); + if (exp == 0) // Subnormal x. + { + ax *= TWO_54; + bits = Double.doubleToLongBits(ax); + exp = (int) (bits >> 52) - 54; + } + exp -= 1023; // Unbias exponent. + ax = Double.longBitsToDouble((bits & 0x000fffffffffffffL) + | 0x3ff0000000000000L); + boolean k; + if (ax < SQRT_1_5) // |x|<sqrt(3/2). + k = false; + else if (ax < SQRT_3) // |x|<sqrt(3). + k = true; + else + { + k = false; + ax *= 0.5; + exp++; + } + + // Compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5). + u = ax - (k ? 1.5 : 1); + v = 1 / (ax + (k ? 1.5 : 1)); + double s = u * v; + double s_h = (float) s; + double t_h = (float) (ax + (k ? 1.5 : 1)); + double t_l = ax - (t_h - (k ? 1.5 : 1)); + double s_l = v * ((u - s_h * t_h) - s_h * t_l); + // Compute log(ax). + double s2 = s * s; + double r = s_l * (s_h + s) + s2 * s2 + * (L1 + s2 * (L2 + s2 * (L3 + s2 * (L4 + s2 * (L5 + s2 * L6))))); + s2 = s_h * s_h; + t_h = (float) (3.0 + s2 + r); + t_l = r - (t_h - 3.0 - s2); + // u+v = s*(1+...). + u = s_h * t_h; + v = s_l * t_h + t_l * s; + // 2/(3log2)*(s+...). + double p_h = (float) (u + v); + double p_l = v - (p_h - u); + double z_h = CP_H * p_h; + double z_l = CP_L * p_h + p_l * CP + (k ? DP_L : 0); + // log2(ax) = (s+..)*2/(3*log2) = exp + dp_h + z_h + z_l. + t = exp; + t1 = (float) (z_h + z_l + (k ? DP_H : 0) + t); + t2 = z_l - (t1 - t - (k ? DP_H : 0) - z_h); + } + + // Split up y into y1+y2 and compute (y1+y2)*(t1+t2). + boolean negative = x < 0 && yisint == 1; + double y1 = (float) y; + double p_l = (y - y1) * t1 + y * t2; + double p_h = y1 * t1; + double z = p_l + p_h; + if (z >= 1024) // Detect overflow. + { + if (z > 1024 || p_l + OVT > z - p_h) + return negative ? Double.NEGATIVE_INFINITY + : Double.POSITIVE_INFINITY; + } + else if (z <= -1075) // Detect underflow. + { + if (z < -1075 || p_l <= z - p_h) + return negative ? -0.0 : 0; + } + + // Compute 2**(p_h+p_l). + int n = round((float) z); + p_h -= n; + t = (float) (p_l + p_h); + u = t * LN2_H; + v = (p_l - (t - p_h)) * LN2 + t * LN2_L; + z = u + v; + w = v - (z - u); + t = z * z; + t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); + double r = (z * t1) / (t1 - 2) - (w + z * w); + z = scale(1 - (r - z), n); + return negative ? -z : z; + } + + /** + * Get the IEEE 754 floating point remainder on two numbers. This is the + * value of <code>x - y * <em>n</em></code>, where <em>n</em> is the closest + * double to <code>x / y</code> (ties go to the even n); for a zero + * remainder, the sign is that of <code>x</code>. If either argument is NaN, + * the first argument is infinite, or the second argument is zero, the result + * is NaN; if x is finite but y is infinite, the result is x. + * + * @param x the dividend (the top half) + * @param y the divisor (the bottom half) + * @return the IEEE 754-defined floating point remainder of x/y + * @see #rint(double) + */ + public static double IEEEremainder(double x, double y) + { + // Purge off exception values. + if (x == Double.NEGATIVE_INFINITY || ! (x < Double.POSITIVE_INFINITY) + || y == 0 || y != y) + return Double.NaN; + + boolean negative = x < 0; + x = abs(x); + y = abs(y); + if (x == y || x == 0) + return 0 * x; // Get correct sign. + + // Achieve x < 2y, then take first shot at remainder. + if (y < TWO_1023) + x %= y + y; + + // Now adjust x to get correct precision. + if (y < 4 / TWO_1023) + { + if (x + x > y) + { + x -= y; + if (x + x >= y) + x -= y; + } + } + else + { + y *= 0.5; + if (x > y) + { + x -= y; + if (x >= y) + x -= y; + } + } + return negative ? -x : x; + } + + /** + * Take the nearest integer that is that is greater than or equal to the + * argument. If the argument is NaN, infinite, or zero, the result is the + * same; if the argument is between -1 and 0, the result is negative zero. + * Note that <code>Math.ceil(x) == -Math.floor(-x)</code>. + * + * @param a the value to act upon + * @return the nearest integer >= <code>a</code> + */ + public static double ceil(double a) + { + return -floor(-a); + } + + /** + * Take the nearest integer that is that is less than or equal to the + * argument. If the argument is NaN, infinite, or zero, the result is the + * same. Note that <code>Math.ceil(x) == -Math.floor(-x)</code>. + * + * @param a the value to act upon + * @return the nearest integer <= <code>a</code> + */ + public static double floor(double a) + { + double x = abs(a); + if (! (x < TWO_52) || (long) a == a) + return a; // No fraction bits; includes NaN and infinity. + if (x < 1) + return a >= 0 ? 0 * a : -1; // Worry about signed zero. + return a < 0 ? (long) a - 1.0 : (long) a; // Cast to long truncates. + } + + /** + * Take the nearest integer to the argument. If it is exactly between + * two integers, the even integer is taken. If the argument is NaN, + * infinite, or zero, the result is the same. + * + * @param a the value to act upon + * @return the nearest integer to <code>a</code> + */ + public static double rint(double a) + { + double x = abs(a); + if (! (x < TWO_52)) + return a; // No fraction bits; includes NaN and infinity. + if (x <= 0.5) + return 0 * a; // Worry about signed zero. + if (x % 2 <= 0.5) + return (long) a; // Catch round down to even. + return (long) (a + (a < 0 ? -0.5 : 0.5)); // Cast to long truncates. + } + + /** + * Take the nearest integer to the argument. This is equivalent to + * <code>(int) Math.floor(f + 0.5f)</code>. If the argument is NaN, the + * result is 0; otherwise if the argument is outside the range of int, the + * result will be Integer.MIN_VALUE or Integer.MAX_VALUE, as appropriate. + * + * @param f the argument to round + * @return the nearest integer to the argument + * @see Integer#MIN_VALUE + * @see Integer#MAX_VALUE + */ + public static int round(float f) + { + return (int) floor(f + 0.5f); + } + + /** + * Take the nearest long to the argument. This is equivalent to + * <code>(long) Math.floor(d + 0.5)</code>. If the argument is NaN, the + * result is 0; otherwise if the argument is outside the range of long, the + * result will be Long.MIN_VALUE or Long.MAX_VALUE, as appropriate. + * + * @param d the argument to round + * @return the nearest long to the argument + * @see Long#MIN_VALUE + * @see Long#MAX_VALUE + */ + public static long round(double d) + { + return (long) floor(d + 0.5); + } + + /** + * Get a random number. This behaves like Random.nextDouble(), seeded by + * System.currentTimeMillis() when first called. In other words, the number + * is from a pseudorandom sequence, and lies in the range [+0.0, 1.0). + * This random sequence is only used by this method, and is threadsafe, + * although you may want your own random number generator if it is shared + * among threads. + * + * @return a random number + * @see Random#nextDouble() + * @see System#currentTimeMillis() + */ + public static synchronized double random() + { + if (rand == null) + rand = new Random(); + return rand.nextDouble(); + } + + /** + * Convert from degrees to radians. The formula for this is + * radians = degrees * (pi/180); however it is not always exact given the + * limitations of floating point numbers. + * + * @param degrees an angle in degrees + * @return the angle in radians + */ + public static double toRadians(double degrees) + { + return (degrees * PI) / 180; + } + + /** + * Convert from radians to degrees. The formula for this is + * degrees = radians * (180/pi); however it is not always exact given the + * limitations of floating point numbers. + * + * @param rads an angle in radians + * @return the angle in degrees + */ + public static double toDegrees(double rads) + { + return (rads * 180) / PI; + } + + /** + * Constants for scaling and comparing doubles by powers of 2. The compiler + * must automatically inline constructs like (1/TWO_54), so we don't list + * negative powers of two here. + */ + private static final double + TWO_16 = 0x10000, // Long bits 0x40f0000000000000L. + TWO_20 = 0x100000, // Long bits 0x4130000000000000L. + TWO_24 = 0x1000000, // Long bits 0x4170000000000000L. + TWO_27 = 0x8000000, // Long bits 0x41a0000000000000L. + TWO_28 = 0x10000000, // Long bits 0x41b0000000000000L. + TWO_29 = 0x20000000, // Long bits 0x41c0000000000000L. + TWO_31 = 0x80000000L, // Long bits 0x41e0000000000000L. + TWO_49 = 0x2000000000000L, // Long bits 0x4300000000000000L. + TWO_52 = 0x10000000000000L, // Long bits 0x4330000000000000L. + TWO_54 = 0x40000000000000L, // Long bits 0x4350000000000000L. + TWO_57 = 0x200000000000000L, // Long bits 0x4380000000000000L. + TWO_60 = 0x1000000000000000L, // Long bits 0x43b0000000000000L. + TWO_64 = 1.8446744073709552e19, // Long bits 0x43f0000000000000L. + TWO_66 = 7.378697629483821e19, // Long bits 0x4410000000000000L. + TWO_1023 = 8.98846567431158e307; // Long bits 0x7fe0000000000000L. + + /** + * Super precision for 2/pi in 24-bit chunks, for use in + * {@link #remPiOver2(double, double[])}. + */ + private static final int TWO_OVER_PI[] = { + 0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62, + 0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a, + 0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129, + 0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41, + 0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8, + 0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf, + 0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5, + 0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08, + 0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3, + 0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880, + 0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b, + }; + + /** + * Super precision for pi/2 in 24-bit chunks, for use in + * {@link #remPiOver2(double, double[])}. + */ + private static final double PI_OVER_TWO[] = { + 1.570796251296997, // Long bits 0x3ff921fb40000000L. + 7.549789415861596e-8, // Long bits 0x3e74442d00000000L. + 5.390302529957765e-15, // Long bits 0x3cf8469880000000L. + 3.282003415807913e-22, // Long bits 0x3b78cc5160000000L. + 1.270655753080676e-29, // Long bits 0x39f01b8380000000L. + 1.2293330898111133e-36, // Long bits 0x387a252040000000L. + 2.7337005381646456e-44, // Long bits 0x36e3822280000000L. + 2.1674168387780482e-51, // Long bits 0x3569f31d00000000L. + }; + + /** + * More constants related to pi, used in + * {@link #remPiOver2(double, double[])} and elsewhere. + */ + private static final double + PI_L = 1.2246467991473532e-16, // Long bits 0x3ca1a62633145c07L. + PIO2_1 = 1.5707963267341256, // Long bits 0x3ff921fb54400000L. + PIO2_1L = 6.077100506506192e-11, // Long bits 0x3dd0b4611a626331L. + PIO2_2 = 6.077100506303966e-11, // Long bits 0x3dd0b4611a600000L. + PIO2_2L = 2.0222662487959506e-21, // Long bits 0x3ba3198a2e037073L. + PIO2_3 = 2.0222662487111665e-21, // Long bits 0x3ba3198a2e000000L. + PIO2_3L = 8.4784276603689e-32; // Long bits 0x397b839a252049c1L. + + /** + * Natural log and square root constants, for calculation of + * {@link #exp(double)}, {@link #log(double)} and + * {@link #pow(double, double)}. CP is 2/(3*ln(2)). + */ + private static final double + SQRT_1_5 = 1.224744871391589, // Long bits 0x3ff3988e1409212eL. + SQRT_2 = 1.4142135623730951, // Long bits 0x3ff6a09e667f3bcdL. + SQRT_3 = 1.7320508075688772, // Long bits 0x3ffbb67ae8584caaL. + EXP_LIMIT_H = 709.782712893384, // Long bits 0x40862e42fefa39efL. + EXP_LIMIT_L = -745.1332191019411, // Long bits 0xc0874910d52d3051L. + CP = 0.9617966939259756, // Long bits 0x3feec709dc3a03fdL. + CP_H = 0.9617967009544373, // Long bits 0x3feec709e0000000L. + CP_L = -7.028461650952758e-9, // Long bits 0xbe3e2fe0145b01f5L. + LN2 = 0.6931471805599453, // Long bits 0x3fe62e42fefa39efL. + LN2_H = 0.6931471803691238, // Long bits 0x3fe62e42fee00000L. + LN2_L = 1.9082149292705877e-10, // Long bits 0x3dea39ef35793c76L. + INV_LN2 = 1.4426950408889634, // Long bits 0x3ff71547652b82feL. + INV_LN2_H = 1.4426950216293335, // Long bits 0x3ff7154760000000L. + INV_LN2_L = 1.9259629911266175e-8; // Long bits 0x3e54ae0bf85ddf44L. + + /** + * Constants for computing {@link #log(double)}. + */ + private static final double + LG1 = 0.6666666666666735, // Long bits 0x3fe5555555555593L. + LG2 = 0.3999999999940942, // Long bits 0x3fd999999997fa04L. + LG3 = 0.2857142874366239, // Long bits 0x3fd2492494229359L. + LG4 = 0.22222198432149784, // Long bits 0x3fcc71c51d8e78afL. + LG5 = 0.1818357216161805, // Long bits 0x3fc7466496cb03deL. + LG6 = 0.15313837699209373, // Long bits 0x3fc39a09d078c69fL. + LG7 = 0.14798198605116586; // Long bits 0x3fc2f112df3e5244L. + + /** + * Constants for computing {@link #pow(double, double)}. L and P are + * coefficients for series; OVT is -(1024-log2(ovfl+.5ulp)); and DP is ???. + * The P coefficients also calculate {@link #exp(double)}. + */ + private static final double + L1 = 0.5999999999999946, // Long bits 0x3fe3333333333303L. + L2 = 0.4285714285785502, // Long bits 0x3fdb6db6db6fabffL. + L3 = 0.33333332981837743, // Long bits 0x3fd55555518f264dL. + L4 = 0.272728123808534, // Long bits 0x3fd17460a91d4101L. + L5 = 0.23066074577556175, // Long bits 0x3fcd864a93c9db65L. + L6 = 0.20697501780033842, // Long bits 0x3fca7e284a454eefL. + P1 = 0.16666666666666602, // Long bits 0x3fc555555555553eL. + P2 = -2.7777777777015593e-3, // Long bits 0xbf66c16c16bebd93L. + P3 = 6.613756321437934e-5, // Long bits 0x3f11566aaf25de2cL. + P4 = -1.6533902205465252e-6, // Long bits 0xbebbbd41c5d26bf1L. + P5 = 4.1381367970572385e-8, // Long bits 0x3e66376972bea4d0L. + DP_H = 0.5849624872207642, // Long bits 0x3fe2b80340000000L. + DP_L = 1.350039202129749e-8, // Long bits 0x3e4cfdeb43cfd006L. + OVT = 8.008566259537294e-17; // Long bits 0x3c971547652b82feL. + + /** + * Coefficients for computing {@link #sin(double)}. + */ + private static final double + S1 = -0.16666666666666632, // Long bits 0xbfc5555555555549L. + S2 = 8.33333333332249e-3, // Long bits 0x3f8111111110f8a6L. + S3 = -1.984126982985795e-4, // Long bits 0xbf2a01a019c161d5L. + S4 = 2.7557313707070068e-6, // Long bits 0x3ec71de357b1fe7dL. + S5 = -2.5050760253406863e-8, // Long bits 0xbe5ae5e68a2b9cebL. + S6 = 1.58969099521155e-10; // Long bits 0x3de5d93a5acfd57cL. + + /** + * Coefficients for computing {@link #cos(double)}. + */ + private static final double + C1 = 0.0416666666666666, // Long bits 0x3fa555555555554cL. + C2 = -1.388888888887411e-3, // Long bits 0xbf56c16c16c15177L. + C3 = 2.480158728947673e-5, // Long bits 0x3efa01a019cb1590L. + C4 = -2.7557314351390663e-7, // Long bits 0xbe927e4f809c52adL. + C5 = 2.087572321298175e-9, // Long bits 0x3e21ee9ebdb4b1c4L. + C6 = -1.1359647557788195e-11; // Long bits 0xbda8fae9be8838d4L. + + /** + * Coefficients for computing {@link #tan(double)}. + */ + private static final double + T0 = 0.3333333333333341, // Long bits 0x3fd5555555555563L. + T1 = 0.13333333333320124, // Long bits 0x3fc111111110fe7aL. + T2 = 0.05396825397622605, // Long bits 0x3faba1ba1bb341feL. + T3 = 0.021869488294859542, // Long bits 0x3f9664f48406d637L. + T4 = 8.8632398235993e-3, // Long bits 0x3f8226e3e96e8493L. + T5 = 3.5920791075913124e-3, // Long bits 0x3f6d6d22c9560328L. + T6 = 1.4562094543252903e-3, // Long bits 0x3f57dbc8fee08315L. + T7 = 5.880412408202641e-4, // Long bits 0x3f4344d8f2f26501L. + T8 = 2.464631348184699e-4, // Long bits 0x3f3026f71a8d1068L. + T9 = 7.817944429395571e-5, // Long bits 0x3f147e88a03792a6L. + T10 = 7.140724913826082e-5, // Long bits 0x3f12b80f32f0a7e9L. + T11 = -1.8558637485527546e-5, // Long bits 0xbef375cbdb605373L. + T12 = 2.590730518636337e-5; // Long bits 0x3efb2a7074bf7ad4L. + + /** + * Coefficients for computing {@link #asin(double)} and + * {@link #acos(double)}. + */ + private static final double + PS0 = 0.16666666666666666, // Long bits 0x3fc5555555555555L. + PS1 = -0.3255658186224009, // Long bits 0xbfd4d61203eb6f7dL. + PS2 = 0.20121253213486293, // Long bits 0x3fc9c1550e884455L. + PS3 = -0.04005553450067941, // Long bits 0xbfa48228b5688f3bL. + PS4 = 7.915349942898145e-4, // Long bits 0x3f49efe07501b288L. + PS5 = 3.479331075960212e-5, // Long bits 0x3f023de10dfdf709L. + QS1 = -2.403394911734414, // Long bits 0xc0033a271c8a2d4bL. + QS2 = 2.0209457602335057, // Long bits 0x40002ae59c598ac8L. + QS3 = -0.6882839716054533, // Long bits 0xbfe6066c1b8d0159L. + QS4 = 0.07703815055590194; // Long bits 0x3fb3b8c5b12e9282L. + + /** + * Coefficients for computing {@link #atan(double)}. + */ + private static final double + ATAN_0_5H = 0.4636476090008061, // Long bits 0x3fddac670561bb4fL. + ATAN_0_5L = 2.2698777452961687e-17, // Long bits 0x3c7a2b7f222f65e2L. + ATAN_1_5H = 0.982793723247329, // Long bits 0x3fef730bd281f69bL. + ATAN_1_5L = 1.3903311031230998e-17, // Long bits 0x3c7007887af0cbbdL. + AT0 = 0.3333333333333293, // Long bits 0x3fd555555555550dL. + AT1 = -0.19999999999876483, // Long bits 0xbfc999999998ebc4L. + AT2 = 0.14285714272503466, // Long bits 0x3fc24924920083ffL. + AT3 = -0.11111110405462356, // Long bits 0xbfbc71c6fe231671L. + AT4 = 0.09090887133436507, // Long bits 0x3fb745cdc54c206eL. + AT5 = -0.0769187620504483, // Long bits 0xbfb3b0f2af749a6dL. + AT6 = 0.06661073137387531, // Long bits 0x3fb10d66a0d03d51L. + AT7 = -0.058335701337905735, // Long bits 0xbfadde2d52defd9aL. + AT8 = 0.049768779946159324, // Long bits 0x3fa97b4b24760debL. + AT9 = -0.036531572744216916, // Long bits 0xbfa2b4442c6a6c2fL. + AT10 = 0.016285820115365782; // Long bits 0x3f90ad3ae322da11L. + + /** + * Constants for computing {@link #cbrt(double)}. + */ + private static final int + CBRT_B1 = 715094163, // B1 = (682-0.03306235651)*2**20 + CBRT_B2 = 696219795; // B2 = (664-0.03306235651)*2**20 + + /** + * Constants for computing {@link #cbrt(double)}. + */ + private static final double + CBRT_C = 5.42857142857142815906e-01, // Long bits 0x3fe15f15f15f15f1L + CBRT_D = -7.05306122448979611050e-01, // Long bits 0xbfe691de2532c834L + CBRT_E = 1.41428571428571436819e+00, // Long bits 0x3ff6a0ea0ea0ea0fL + CBRT_F = 1.60714285714285720630e+00, // Long bits 0x3ff9b6db6db6db6eL + CBRT_G = 3.57142857142857150787e-01; // Long bits 0x3fd6db6db6db6db7L + + /** + * Constants for computing {@link #expm1(double)} + */ + private static final double + EXPM1_Q1 = -3.33333333333331316428e-02, // Long bits 0xbfa11111111110f4L + EXPM1_Q2 = 1.58730158725481460165e-03, // Long bits 0x3f5a01a019fe5585L + EXPM1_Q3 = -7.93650757867487942473e-05, // Long bits 0xbf14ce199eaadbb7L + EXPM1_Q4 = 4.00821782732936239552e-06, // Long bits 0x3ed0cfca86e65239L + EXPM1_Q5 = -2.01099218183624371326e-07; // Long bits 0xbe8afdb76e09c32dL + + /** + * Helper function for reducing an angle to a multiple of pi/2 within + * [-pi/4, pi/4]. + * + * @param x the angle; not infinity or NaN, and outside pi/4 + * @param y an array of 2 doubles modified to hold the remander x % pi/2 + * @return the quadrant of the result, mod 4: 0: [-pi/4, pi/4], + * 1: [pi/4, 3*pi/4], 2: [3*pi/4, 5*pi/4], 3: [-3*pi/4, -pi/4] + */ + private static int remPiOver2(double x, double[] y) + { + boolean negative = x < 0; + x = abs(x); + double z; + int n; + if (Configuration.DEBUG && (x <= PI / 4 || x != x + || x == Double.POSITIVE_INFINITY)) + throw new InternalError("Assertion failure"); + if (x < 3 * PI / 4) // If |x| is small. + { + z = x - PIO2_1; + if ((float) x != (float) (PI / 2)) // 33+53 bit pi is good enough. + { + y[0] = z - PIO2_1L; + y[1] = z - y[0] - PIO2_1L; + } + else // Near pi/2, use 33+33+53 bit pi. + { + z -= PIO2_2; + y[0] = z - PIO2_2L; + y[1] = z - y[0] - PIO2_2L; + } + n = 1; + } + else if (x <= TWO_20 * PI / 2) // Medium size. + { + n = (int) (2 / PI * x + 0.5); + z = x - n * PIO2_1; + double w = n * PIO2_1L; // First round good to 85 bits. + y[0] = z - w; + if (n >= 32 || (float) x == (float) (w)) + { + if (x / y[0] >= TWO_16) // Second iteration, good to 118 bits. + { + double t = z; + w = n * PIO2_2; + z = t - w; + w = n * PIO2_2L - (t - z - w); + y[0] = z - w; + if (x / y[0] >= TWO_49) // Third iteration, 151 bits accuracy. + { + t = z; + w = n * PIO2_3; + z = t - w; + w = n * PIO2_3L - (t - z - w); + y[0] = z - w; + } + } + } + y[1] = z - y[0] - w; + } + else + { + // All other (large) arguments. + int e0 = (int) (Double.doubleToLongBits(x) >> 52) - 1046; + z = scale(x, -e0); // e0 = ilogb(z) - 23. + double[] tx = new double[3]; + for (int i = 0; i < 2; i++) + { + tx[i] = (int) z; + z = (z - tx[i]) * TWO_24; + } + tx[2] = z; + int nx = 2; + while (tx[nx] == 0) + nx--; + n = remPiOver2(tx, y, e0, nx); + } + if (negative) + { + y[0] = -y[0]; + y[1] = -y[1]; + return -n; + } + return n; + } + + /** + * Helper function for reducing an angle to a multiple of pi/2 within + * [-pi/4, pi/4]. + * + * @param x the positive angle, broken into 24-bit chunks + * @param y an array of 2 doubles modified to hold the remander x % pi/2 + * @param e0 the exponent of x[0] + * @param nx the last index used in x + * @return the quadrant of the result, mod 4: 0: [-pi/4, pi/4], + * 1: [pi/4, 3*pi/4], 2: [3*pi/4, 5*pi/4], 3: [-3*pi/4, -pi/4] + */ + private static int remPiOver2(double[] x, double[] y, int e0, int nx) + { + int i; + int ih; + int n; + double fw; + double z; + int[] iq = new int[20]; + double[] f = new double[20]; + double[] q = new double[20]; + boolean recompute = false; + + // Initialize jk, jz, jv, q0; note that 3>q0. + int jk = 4; + int jz = jk; + int jv = max((e0 - 3) / 24, 0); + int q0 = e0 - 24 * (jv + 1); + + // Set up f[0] to f[nx+jk] where f[nx+jk] = TWO_OVER_PI[jv+jk]. + int j = jv - nx; + int m = nx + jk; + for (i = 0; i <= m; i++, j++) + f[i] = (j < 0) ? 0 : TWO_OVER_PI[j]; + + // Compute q[0],q[1],...q[jk]. + for (i = 0; i <= jk; i++) + { + for (j = 0, fw = 0; j <= nx; j++) + fw += x[j] * f[nx + i - j]; + q[i] = fw; + } + + do + { + // Distill q[] into iq[] reversingly. + for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) + { + fw = (int) (1 / TWO_24 * z); + iq[i] = (int) (z - TWO_24 * fw); + z = q[j - 1] + fw; + } + + // Compute n. + z = scale(z, q0); + z -= 8 * floor(z * 0.125); // Trim off integer >= 8. + n = (int) z; + z -= n; + ih = 0; + if (q0 > 0) // Need iq[jz-1] to determine n. + { + i = iq[jz - 1] >> (24 - q0); + n += i; + iq[jz - 1] -= i << (24 - q0); + ih = iq[jz - 1] >> (23 - q0); + } + else if (q0 == 0) + ih = iq[jz - 1] >> 23; + else if (z >= 0.5) + ih = 2; + + if (ih > 0) // If q > 0.5. + { + n += 1; + int carry = 0; + for (i = 0; i < jz; i++) // Compute 1-q. + { + j = iq[i]; + if (carry == 0) + { + if (j != 0) + { + carry = 1; + iq[i] = 0x1000000 - j; + } + } + else + iq[i] = 0xffffff - j; + } + switch (q0) + { + case 1: // Rare case: chance is 1 in 12 for non-default. + iq[jz - 1] &= 0x7fffff; + break; + case 2: + iq[jz - 1] &= 0x3fffff; + } + if (ih == 2) + { + z = 1 - z; + if (carry != 0) + z -= scale(1, q0); + } + } + + // Check if recomputation is needed. + if (z == 0) + { + j = 0; + for (i = jz - 1; i >= jk; i--) + j |= iq[i]; + if (j == 0) // Need recomputation. + { + int k; // k = no. of terms needed. + for (k = 1; iq[jk - k] == 0; k++) + ; + + for (i = jz + 1; i <= jz + k; i++) // Add q[jz+1] to q[jz+k]. + { + f[nx + i] = TWO_OVER_PI[jv + i]; + for (j = 0, fw = 0; j <= nx; j++) + fw += x[j] * f[nx + i - j]; + q[i] = fw; + } + jz += k; + recompute = true; + } + } + } + while (recompute); + + // Chop off zero terms. + if (z == 0) + { + jz--; + q0 -= 24; + while (iq[jz] == 0) + { + jz--; + q0 -= 24; + } + } + else // Break z into 24-bit if necessary. + { + z = scale(z, -q0); + if (z >= TWO_24) + { + fw = (int) (1 / TWO_24 * z); + iq[jz] = (int) (z - TWO_24 * fw); + jz++; + q0 += 24; + iq[jz] = (int) fw; + } + else + iq[jz] = (int) z; + } + + // Convert integer "bit" chunk to floating-point value. + fw = scale(1, q0); + for (i = jz; i >= 0; i--) + { + q[i] = fw * iq[i]; + fw *= 1 / TWO_24; + } + + // Compute PI_OVER_TWO[0,...,jk]*q[jz,...,0]. + double[] fq = new double[20]; + for (i = jz; i >= 0; i--) + { + fw = 0; + for (int k = 0; k <= jk && k <= jz - i; k++) + fw += PI_OVER_TWO[k] * q[i + k]; + fq[jz - i] = fw; + } + + // Compress fq[] into y[]. + fw = 0; + for (i = jz; i >= 0; i--) + fw += fq[i]; + y[0] = (ih == 0) ? fw : -fw; + fw = fq[0] - fw; + for (i = 1; i <= jz; i++) + fw += fq[i]; + y[1] = (ih == 0) ? fw : -fw; + return n; + } + + /** + * Helper method for scaling a double by a power of 2. + * + * @param x the double + * @param n the scale; |n| < 2048 + * @return x * 2**n + */ + private static double scale(double x, int n) + { + if (Configuration.DEBUG && abs(n) >= 2048) + throw new InternalError("Assertion failure"); + if (x == 0 || x == Double.NEGATIVE_INFINITY + || ! (x < Double.POSITIVE_INFINITY) || n == 0) + return x; + long bits = Double.doubleToLongBits(x); + int exp = (int) (bits >> 52) & 0x7ff; + if (exp == 0) // Subnormal x. + { + x *= TWO_54; + exp = ((int) (Double.doubleToLongBits(x) >> 52) & 0x7ff) - 54; + } + exp += n; + if (exp > 0x7fe) // Overflow. + return Double.POSITIVE_INFINITY * x; + if (exp > 0) // Normal. + return Double.longBitsToDouble((bits & 0x800fffffffffffffL) + | ((long) exp << 52)); + if (exp <= -54) + return 0 * x; // Underflow. + exp += 54; // Subnormal result. + x = Double.longBitsToDouble((bits & 0x800fffffffffffffL) + | ((long) exp << 52)); + return x * (1 / TWO_54); + } + + /** + * Helper trig function; computes sin in range [-pi/4, pi/4]. + * + * @param x angle within about pi/4 + * @param y tail of x, created by remPiOver2 + * @return sin(x+y) + */ + private static double sin(double x, double y) + { + if (Configuration.DEBUG && abs(x + y) > 0.7854) + throw new InternalError("Assertion failure"); + if (abs(x) < 1 / TWO_27) + return x; // If |x| ~< 2**-27, already know answer. + + double z = x * x; + double v = z * x; + double r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6))); + if (y == 0) + return x + v * (S1 + z * r); + return x - ((z * (0.5 * y - v * r) - y) - v * S1); + } + + /** + * Helper trig function; computes cos in range [-pi/4, pi/4]. + * + * @param x angle within about pi/4 + * @param y tail of x, created by remPiOver2 + * @return cos(x+y) + */ + private static double cos(double x, double y) + { + if (Configuration.DEBUG && abs(x + y) > 0.7854) + throw new InternalError("Assertion failure"); + x = abs(x); + if (x < 1 / TWO_27) + return 1; // If |x| ~< 2**-27, already know answer. + + double z = x * x; + double r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6))))); + + if (x < 0.3) + return 1 - (0.5 * z - (z * r - x * y)); + + double qx = (x > 0.78125) ? 0.28125 : (x * 0.25); + return 1 - qx - ((0.5 * z - qx) - (z * r - x * y)); + } + + /** + * Helper trig function; computes tan in range [-pi/4, pi/4]. + * + * @param x angle within about pi/4 + * @param y tail of x, created by remPiOver2 + * @param invert true iff -1/tan should be returned instead + * @return tan(x+y) + */ + private static double tan(double x, double y, boolean invert) + { + // PI/2 is irrational, so no double is a perfect multiple of it. + if (Configuration.DEBUG && (abs(x + y) > 0.7854 || (x == 0 && invert))) + throw new InternalError("Assertion failure"); + boolean negative = x < 0; + if (negative) + { + x = -x; + y = -y; + } + if (x < 1 / TWO_28) // If |x| ~< 2**-28, already know answer. + return (negative ? -1 : 1) * (invert ? -1 / x : x); + + double z; + double w; + boolean large = x >= 0.6744; + if (large) + { + z = PI / 4 - x; + w = PI_L / 4 - y; + x = z + w; + y = 0; + } + z = x * x; + w = z * z; + // Break x**5*(T1+x**2*T2+...) into + // x**5(T1+x**4*T3+...+x**20*T11) + // + x**5(x**2*(T2+x**4*T4+...+x**22*T12)). + double r = T1 + w * (T3 + w * (T5 + w * (T7 + w * (T9 + w * T11)))); + double v = z * (T2 + w * (T4 + w * (T6 + w * (T8 + w * (T10 + w * T12))))); + double s = z * x; + r = y + z * (s * (r + v) + y); + r += T0 * s; + w = x + r; + if (large) + { + v = invert ? -1 : 1; + return (negative ? -1 : 1) * (v - 2 * (x - (w * w / (w + v) - r))); + } + if (! invert) + return w; + + // Compute -1.0/(x+r) accurately. + z = (float) w; + v = r - (z - x); + double a = -1 / w; + double t = (float) a; + return t + a * (1 + t * z + t * v); + } + + /** + * <p> + * Returns the sign of the argument as follows: + * </p> + * <ul> + * <li>If <code>a</code> is greater than zero, the result is 1.0.</li> + * <li>If <code>a</code> is less than zero, the result is -1.0.</li> + * <li>If <code>a</code> is <code>NaN</code>, the result is <code>NaN</code>. + * <li>If <code>a</code> is positive or negative zero, the result is the + * same.</li> + * </ul> + * + * @param a the numeric argument. + * @return the sign of the argument. + * @since 1.5. + */ + public static double signum(double a) + { + // There's no difference. + return Math.signum(a); + } + + /** + * <p> + * Returns the sign of the argument as follows: + * </p> + * <ul> + * <li>If <code>a</code> is greater than zero, the result is 1.0f.</li> + * <li>If <code>a</code> is less than zero, the result is -1.0f.</li> + * <li>If <code>a</code> is <code>NaN</code>, the result is <code>NaN</code>. + * <li>If <code>a</code> is positive or negative zero, the result is the + * same.</li> + * </ul> + * + * @param a the numeric argument. + * @return the sign of the argument. + * @since 1.5. + */ + public static float signum(float a) + { + // There's no difference. + return Math.signum(a); + } + + /** + * Return the ulp for the given double argument. The ulp is the + * difference between the argument and the next larger double. Note + * that the sign of the double argument is ignored, that is, + * ulp(x) == ulp(-x). If the argument is a NaN, then NaN is returned. + * If the argument is an infinity, then +Inf is returned. If the + * argument is zero (either positive or negative), then + * {@link Double#MIN_VALUE} is returned. + * @param d the double whose ulp should be returned + * @return the difference between the argument and the next larger double + * @since 1.5 + */ + public static double ulp(double d) + { + // There's no difference. + return Math.ulp(d); + } + + /** + * Return the ulp for the given float argument. The ulp is the + * difference between the argument and the next larger float. Note + * that the sign of the float argument is ignored, that is, + * ulp(x) == ulp(-x). If the argument is a NaN, then NaN is returned. + * If the argument is an infinity, then +Inf is returned. If the + * argument is zero (either positive or negative), then + * {@link Float#MIN_VALUE} is returned. + * @param f the float whose ulp should be returned + * @return the difference between the argument and the next larger float + * @since 1.5 + */ + public static float ulp(float f) + { + // There's no difference. + return Math.ulp(f); + } +} |