From c2134b3a905bae03e87b0fab17f6b5cccf3cbc02 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Sun, 5 Sep 2010 19:27:05 +0000 Subject: [PATCH] fixed errors with infinities added asin/acos git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@992870 13f79535-47bb-0310-9956-ffa450edef68 --- .../apache/commons/math/util/FastMath.java | 348 ++++++++++++++---- 1 file changed, 283 insertions(+), 65 deletions(-) diff --git a/src/main/java/org/apache/commons/math/util/FastMath.java b/src/main/java/org/apache/commons/math/util/FastMath.java index 3d963e54e..09c7cd0bf 100644 --- a/src/main/java/org/apache/commons/math/util/FastMath.java +++ b/src/main/java/org/apache/commons/math/util/FastMath.java @@ -204,22 +204,6 @@ public class FastMath { private FastMath() { } - /** Compute the arc cosine of a number. - * @param a number on which evaluation is done - * @return arc cosine of a - */ - public static double acos(final double a) { - return Math.acos(a); - } - - /** Compute the arc sine of a number. - * @param a number on which evaluation is done - * @return arc sine of a - */ - public static double asin(final double a) { - return Math.asin(a); - } - /** Compute the square root of a number. * @param a number on which evaluation is done * @return square root of a @@ -441,6 +425,10 @@ public class FastMath { intVal = (int) -x; if (intVal > 746) { + if (hiPrec != null) { + hiPrec[0] = 0.0; + hiPrec[1] = 0.0; + } return 0.0; } @@ -474,6 +462,10 @@ public class FastMath { intVal = (int) x; if (intVal > 709) { + if (hiPrec != null) { + hiPrec[0] = Double.POSITIVE_INFINITY; + hiPrec[1] = 0.0; + } return Double.POSITIVE_INFINITY; } @@ -1171,6 +1163,14 @@ public class FastMath { 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]; @@ -1227,22 +1227,28 @@ public class FastMath { return 1.0; } - /* Handle special case x<0 */ - if (x < 0) { - if (y == (long) y) { - // If y is an integer - return ((long)y & 1) == 0 ? pow(-x, y) : -pow(-x, y); - } else { - return Double.NaN; - } + if (x != x) { // X is NaN + return x; } + if (x == 0) { long bits = Double.doubleToLongBits(x); if ((bits & 0x8000000000000000L) != 0) { // -zero - if (y < 0 && y == (long)y) + 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) { @@ -1256,6 +1262,9 @@ public class FastMath { } if (x == Double.POSITIVE_INFINITY) { + if (y != y) { // y is NaN + return y; + } if (y < 0.0) { return 0.0; } else { @@ -1264,6 +1273,9 @@ public class FastMath { } if (y == Double.POSITIVE_INFINITY) { + if (x * x == 1.0) + return Double.NaN; + if (x * x > 1.0) { return Double.POSITIVE_INFINITY; } else { @@ -1271,18 +1283,71 @@ public class FastMath { } } + 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.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 >= 4503599627370496.0 || y <= -4503599627370496.0) { + 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 tmp1 = y * 1073741824.0; - final double ya = y + tmp1 - tmp1; - final double yb = y - ya; + double ya; + double yb; + if (y < 8e298 && y > -8e298) { + double tmp1 = y * 1073741824.0; + ya = y + tmp1 - tmp1; + yb = y - ya; + } else { + double tmp1 = y * 9.31322574615478515625E-10; + double tmp2 = tmp1 * 9.31322574615478515625E-10; + ya = (tmp1 + tmp2 - tmp1) * 1073741824.0 * 1073741824.0; + yb = y - ya; + } /* Compute ln(x) */ log(x, lns); @@ -1290,8 +1355,8 @@ public class FastMath { double lnb = lns[1]; /* resplit lns */ - tmp1 = lna * 1073741824.0; - final double tmp2 = lna + tmp1 - tmp1; + double tmp1 = lna * 1073741824.0; + double tmp2 = lna + tmp1 - tmp1; lnb += lna - tmp2; lna = tmp2; @@ -2375,8 +2440,8 @@ public class FastMath { double b = -(a - pi2a + xa); b += pi2b - xb; - xa = a; - xb = b; + xa = a + b; + xb = -(xa - a - b); quadrant ^= 1; negative ^= true; } @@ -2412,7 +2477,6 @@ public class FastMath { */ private static double atan(double xa, double xb, boolean leftPlane) { boolean negate = false; - boolean recip = false; int idx; if (xa < 0) { @@ -2519,41 +2583,24 @@ public class FastMath { double result; double resultb; - if (recip) { - final double pi2a = 1.5707963267948966; - final double pi2b = 6.123233995736766E-17; - double za = pi2a - ya; - double zb = -(za - pi2a + ya); - temp = za - EIGHTHES[idx]; - zb += -(temp - za + EIGHTHES[idx]); - za = temp; + //result = yb + eighths[idx] + ya; + double za = EIGHTHES[idx] + ya; + double zb = -(za - EIGHTHES[idx] - ya); + temp = za + yb; + zb += -(temp - za - yb); + za = temp; - zb += pi2b - yb; - ya = za; - yb = zb; - - result = yb + ya; - resultb = -(result - yb - ya); - } else { - //result = yb + eighths[idx] + ya; - double za = EIGHTHES[idx] + ya; - double zb = -(za - EIGHTHES[idx] - ya); - temp = za + yb; - zb += -(temp - za - yb); - za = temp; - - result = za + zb; - resultb = -(result - za - zb); - } + result = za + zb; + resultb = -(result - za - zb); if (leftPlane) { // Result is in the left plane final double pia = 1.5707963267948966*2.0; final double pib = 6.123233995736766E-17*2.0; - final double za = pia - result; - double zb = -(za - pia + result); + za = pia - result; + zb = -(za - pia + result); zb += pib - resultb; result = za + zb; @@ -2585,7 +2632,11 @@ public class FastMath { double invy = 1.0/y; if (invx == 0.0) { // X is infinite - return 0.0; + if (x > 0) { + return 0.0; + } else { + return Math.PI; + } } if (result != result) { // y must be infinite @@ -2686,6 +2737,155 @@ public class FastMath { 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; + } + + /* Compute asin(x) = atan(x/sqrt(1-x*x)) */ + + /* Split x */ + double temp = x * 1073741824.0; + 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 * 1073741824.0; + 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 * 1073741824.0; + 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 * 1073741824.0; + 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 * 1073741824.0; + 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; + temp = r * 1073741824.0; + double ra = r + temp - temp; + 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); + } + /** * Convert degrees to radians, with error of less than 0.5 ULP * @param x angle in degrees @@ -2829,15 +3029,23 @@ public class FastMath { public static double floor(double x) { long y; + if (x != x) { // NaN + return x; + } + if (x >= 4503599627370496.0 || x <= -4503599627370496.0) { return x; } y = (long) x; - if (x < 0) { + if (x < 0 && y != x) { y--; } + if (y == 0) { + return x*y; + } + return (double) y; } @@ -2848,12 +3056,22 @@ public class FastMath { public static double ceil(double x) { double y; + if (x != x) { // NaN + return x; + } + y = floor(x); if (y == x) { return y; } - return y + 1.0; + 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.