You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2010/09/05 21:26:35 UTC
svn commit: r992869 -
/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java
Author: luc
Date: Sun Sep 5 19:26:35 2010
New Revision: 992869
URL: http://svn.apache.org/viewvc?rev=992869&view=rev
Log:
fixed errors with infinities
added asin/acos
Modified:
commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java
Modified: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java?rev=992869&r1=992868&r2=992869&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java Sun Sep 5 19:26:35 2010
@@ -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;
- 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 = 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.