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

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 != x) { // X is NaN
-        if (x < 0) {
+            return x;
-            if (y == (long) y) {
-                // If y is an integer
-                return ((long)y & 1) == 0 ? pow(-x, y) : -pow(-x, y);
-            } else {
-                return Double.NaN;
-            }
         }
 
+
         if (x == 0) {
             long bits = Double.doubleToLongBits(x);
             if ((bits & 0x8000000000000000L) != 0) {
                 // -zero
-                if (y < 0 && y == (long)y)
+                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;
+        double ya;
-        final double ya = y + tmp1 - tmp1;
+        double yb;
-        final double yb = y - ya;
+        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;
+        double tmp1 = lna * 1073741824.0;
-        final double tmp2 = lna + tmp1 - tmp1;
+        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;
+            xa = a + b;
-            xb = 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;
+        //result = yb + eighths[idx] + ya;
-            double zb = -(za - pi2a + ya);
+        double za = EIGHTHES[idx] + ya;
-            temp = za - EIGHTHES[idx];
+        double zb = -(za - EIGHTHES[idx] - ya);
-            zb += -(temp - za + EIGHTHES[idx]);
+        temp = za + yb;
-            za = temp;
+        zb += -(temp - za - yb);
+        za = temp;
 
-            zb += pi2b - yb;
+        result = za + zb;
-            ya = za;
+        resultb = -(result - za - zb);
-            yb = zb;
-
-            result = yb + ya;
-            resultb = -(result - yb - ya);
-        } else {
-            //result = yb + eighths[idx] + ya;
-            double za = EIGHTHES[idx] + ya;
-            double zb = -(za - EIGHTHES[idx] - ya);
-            temp = za + yb;
-            zb += -(temp - za - yb);
-            za = temp;
-
-            result = za + zb;
-            resultb = -(result - za - zb);
-        }
 
         if (leftPlane) {
             // Result is in the left plane
             final double pia = 1.5707963267948966*2.0;
             final double pib = 6.123233995736766E-17*2.0;
 
-            final double za = pia - result;
+            za = pia - result;
-            double zb = -(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.