commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r1004045 - in /commons/proper/math/branches/MATH_2_X/src: main/java/org/apache/commons/math/util/FastMath.java test/java/org/apache/commons/math/util/FastMathTest.java
Date Sun, 03 Oct 2010 20:52:54 GMT
Author: luc
Date: Sun Oct  3 20:52:53 2010
New Revision: 1004045

URL: http://svn.apache.org/viewvc?rev=1004045&view=rev
Log:
applied Bill's patch adding sinh, cosh and tanh functions to FastMath
JIRA: MATH-375

Modified:
    commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java
    commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/util/FastMathTest.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=1004045&r1=1004044&r2=1004045&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 Oct  3 20:52:53 2010
@@ -220,68 +220,286 @@ public class FastMath {
     }
 
     /** Compute the hyperbolic cosine of a number.
-     * @param a number on which evaluation is done
-     * @return hyperbolic cosine of a
+     * @param x number on which evaluation is done
+     * @return hyperbolic cosine of x
      */
-    public static double cosh(final double a) {
-        return 0.5 * (FastMath.exp(a) + FastMath.exp(-a));
+    public static double cosh(double x) {
+      if (x != x) {
+          return x;
+      }
+
+      if (x > 20.0) {
+          return exp(x)/2.0;
+      }
+
+      if (x < -20) {
+          return exp(-x)/2.0;
+      }
+
+      double hiPrec[] = new double[2];
+      if (x < 0.0) {
+          x = -x;
+      }
+      exp(x, 0.0, hiPrec);
+
+      double ya = hiPrec[0] + hiPrec[1];
+      double yb = -(ya - hiPrec[0] - hiPrec[1]);
+
+      double temp = ya * 1073741824.0;
+      double yaa = ya + temp - temp;
+      double yab = ya - yaa;
+
+      // recip = 1/y
+      double recip = 1.0/ya;
+      temp = recip * 1073741824.0;
+      double recipa = recip + temp - temp;
+      double recipb = recip - recipa;
+
+      // Correct for rounding in division
+      recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
+      // Account for yb
+      recipb += -yb * recip * recip;
+
+      // y = y + 1/y
+      temp = ya + recipa;
+      yb += -(temp - ya - recipa);
+      ya = temp;
+      temp = ya + recipb;
+      yb += -(temp - ya - recipb);
+      ya = temp;
+
+      double result = ya + yb;
+      result *= 0.5;
+      return result;
     }
 
     /** Compute the hyperbolic sine of a number.
-     * @param a number on which evaluation is done
-     * @return hyperbolic sine of a
+     * @param x number on which evaluation is done
+     * @return hyperbolic sine of x
      */
-    public static double sinh(double a) {
+    public static double sinh(double x) {
+      boolean negate = false;
+      if (x != x) {
+          return x;
+      }
 
-        boolean negative = false;
-        if (a < 0) {
-            negative = true;
-            a = -a;
-        }
+      if (x > 20.0) {
+          return exp(x)/2.0;
+      }
 
-        double absSinh;
-        if (a > 0.3) {
-            absSinh = 0.5 * (FastMath.exp(a) - FastMath.exp(-a));
-        } else {
-            final double a2 = a * a;
-            if (a > 0.05) {
-                absSinh = a * (1 + a2 * (1 + a2  * (1 + a2 * (1 + a2 * (1 + a2 / 110) / 72)
/ 42) / 20) / 6);
-            } else {
-                absSinh = a * (1 + a2 * (1 + a2  * (1 + a2 / 42) / 20) / 6);
-            }
-        }
+      if (x < -20) {
+          return -exp(-x)/2.0;
+      }
 
-        return negative ? -absSinh : absSinh;
+      if (x == 0) {
+          return x;
+      }
+
+      if (x < 0.0) {
+          x = -x;
+          negate = true;
+      }
 
+      double result;
+
+      if (x > 0.25) {
+          double hiPrec[] = new double[2];
+          exp(x, 0.0, hiPrec);
+
+          double ya = hiPrec[0] + hiPrec[1];
+          double yb = -(ya - hiPrec[0] - hiPrec[1]);
+
+          double temp = ya * 1073741824.0;
+          double yaa = ya + temp - temp;
+          double yab = ya - yaa;
+
+          // recip = 1/y
+          double recip = 1.0/ya;
+          temp = recip * 1073741824.0;
+          double recipa = recip + temp - temp;
+          double recipb = recip - recipa;
+
+          // Correct for rounding in division
+          recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
+          // Account for yb
+          recipb += -yb * recip * recip;
+
+          recipa = -recipa;
+          recipb = -recipb;
+
+          // y = y + 1/y
+          temp = ya + recipa;
+          yb += -(temp - ya - recipa);
+          ya = temp;
+          temp = ya + recipb;
+          yb += -(temp - ya - recipb);
+          ya = temp;
+
+          result = ya + yb;
+          result *= 0.5;
+      }
+      else {
+          double hiPrec[] = new double[2];
+          expm1(x, hiPrec);
+
+          double ya = hiPrec[0] + hiPrec[1];
+          double yb = -(ya - hiPrec[0] - hiPrec[1]);
+
+          /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
+          double denom = 1.0+ya;
+          double denomr = 1.0/denom;
+          double denomb = -(denom - 1.0 - ya) + yb;
+          double ratio = ya*denomr;
+          double ra, rb;
+          double temp = ratio * 1073741824.0;
+          ra = ratio + temp - temp;
+          rb = ratio - ra;
+
+          temp = denom * 1073741824.0;
+          double za = denom + temp - temp;
+          double zb = denom - za;
+
+          rb += (ya - za*ra - za*rb - zb*ra - zb*rb) * denomr;
+
+          // Adjust for yb
+          rb += yb*denomr;                        // numerator
+          rb += -ya * denomb * denomr * denomr;   // denominator
+
+          // y = y - 1/y
+          temp = ya + ra;
+          yb += -(temp - ya - ra);
+          ya = temp;
+          temp = ya + rb;
+          yb += -(temp - ya - rb);
+          ya = temp;
+
+          result = ya + yb;
+          result *= 0.5;
+      }
+
+      if (negate) {
+          result = -result;
+      }
+
+      return result;
     }
 
     /** Compute the hyperbolic tangent of a number.
-     * @param a number on which evaluation is done
-     * @return hyperbolic tangent of a
+     * @param x number on which evaluation is done
+     * @return hyperbolic tangent of x
      */
-    public static double tanh(double a) {
+    public static double tanh(double x) {
+      boolean negate = false;
 
-        boolean negative = false;
-        if (a < 0) {
-            negative = true;
-            a = -a;
-        }
+      if (x != x) {
+          return x;
+      }
 
-        double absTanh;
-        if (a > 0.074) {
-            final double twoA = 2 * a;
-            absTanh = FastMath.expm1(twoA) / (FastMath.exp(twoA) + 1);
-        } else {
-            final double a2 = a * a;
-            if (a > 0.016) {
-                absTanh = a * (1 - a2 * (1 - a2 * (2 - a2 * (17 - a2 * (62 - a2 * 1382 /
55 ) / 9) / 21) / 5) / 3);
-            } else {
-                absTanh = a * (1 - a2 * (1 - a2 * (2 - a2 * 17 / 21) / 5) / 3);
-            }
-        }
+      if (x > 20.0) {
+          return 1.0;
+      }
 
-        return negative ? -absTanh : absTanh;
+      if (x < -20) {
+          return -1.0;
+      }
 
+      if (x == 0) {
+          return x;
+      }
+
+      if (x < 0.0) {
+          x = -x;
+          negate = true;
+      }
+
+      double result;
+      if (x >= 0.5) {
+          double hiPrec[] = new double[2];
+          // tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
+          exp(x*2.0, 0.0, hiPrec);
+
+          double ya = hiPrec[0] + hiPrec[1];
+          double yb = -(ya - hiPrec[0] - hiPrec[1]);
+
+          /* Numerator */
+          double na = -1.0 + ya;
+          double nb = -(na + 1.0 - ya);
+          double temp = na + yb;
+          nb += -(temp - na - yb);
+          na = temp;
+
+          /* Denominator */
+          double da = 1.0 + ya;
+          double db = -(da - 1.0 - ya);
+          temp = da + yb;
+          db += -(temp - da - yb);
+          da = temp;
+
+          temp = da * 1073741824.0;
+          double daa = da + temp - temp;
+          double dab = da - daa;
+
+          // ratio = na/da
+          double ratio = na/da;
+          temp = ratio * 1073741824.0;
+          double ratioa = ratio + temp - temp;
+          double ratiob = ratio - ratioa;
+
+          // Correct for rounding in division
+          ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;
+
+          // Account for nb
+          ratiob += nb / da;
+          // Account for db
+          ratiob += -db * na / da / da;
+
+          result = ratioa + ratiob;
+      }
+      else {
+          double hiPrec[] = new double[2];
+          // tanh(x) = expm1(2x) / (expm1(2x) + 2)
+          expm1(x*2.0, hiPrec);
+
+          double ya = hiPrec[0] + hiPrec[1];
+          double yb = -(ya - hiPrec[0] - hiPrec[1]);
+
+          /* Numerator */
+          double na = ya;
+          double nb = yb;
+
+          /* Denominator */
+          double da = 2.0 + ya;
+          double db = -(da - 2.0 - ya);
+          double temp = da + yb;
+          db += -(temp - da - yb);
+          da = temp;
+
+          temp = da * 1073741824.0;
+          double daa = da + temp - temp;
+          double dab = da - daa;
+
+          // ratio = na/da
+          double ratio = na/da;
+          temp = ratio * 1073741824.0;
+          double ratioa = ratio + temp - temp;
+          double ratiob = ratio - ratioa;
+
+          // Correct for rounding in division
+          ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;
+
+          // Account for nb
+          ratiob += nb / da;
+          // Account for db
+          ratiob += -db * na / da / da;
+
+          result = ratioa + ratiob;
+      }
+      
+      if (negate) {
+          result = -result;
+      }
+
+      return result;
     }
 
     /** Compute the inverse hyperbolic cosine of a number.
@@ -533,6 +751,15 @@ public class FastMath {
      * @return exp(x) - 1
      */
     public static double expm1(double x) {
+      return expm1(x, null);
+    }
+
+    /** Internal helper method for expm1
+     * @param x number to compute shifted exponential
+     * @param hiPrecOut[] receive high precision result for -1.0 < x < 1.0
+     * @return exp(x) - 1
+     */
+    private static double expm1(double x, double hiPrecOut[]) {
         if (x != x || x == 0.0) { // NaN or zero
             return x;
         }
@@ -666,6 +893,11 @@ public class FastMath {
             yb = -rb;
         }
 
+	if (hiPrecOut != null) {
+	  hiPrecOut[0] = ya;
+	  hiPrecOut[1] = yb;
+	}
+
         return ya + yb;
     }
 

Modified: commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/util/FastMathTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/util/FastMathTest.java?rev=1004045&r1=1004044&r2=1004045&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/util/FastMathTest.java
(original)
+++ commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/util/FastMathTest.java
Sun Oct  3 20:52:53 2010
@@ -788,6 +788,84 @@ public class FastMathTest {
         Assert.assertTrue("acos() had errors in excess of " + MAX_ERROR_ULP + " ULP", maxerrulp
< MAX_ERROR_ULP);
     }
 
+    private Dfp cosh(Dfp x) {
+      return DfpMath.exp(x).add(DfpMath.exp(x.negate())).divide(2);
+    }
+
+    private Dfp sinh(Dfp x) {
+      return DfpMath.exp(x).subtract(DfpMath.exp(x.negate())).divide(2);
+    }
+
+    private Dfp tanh(Dfp x) {
+      return sinh(x).divide(cosh(x));
+    }
+
+    @Test
+    public void testSinhAccuracy() {
+        double maxerrulp = 0.0;
+
+        for (int i=0; i<10000; i++) {
+            double x = ((generator.nextDouble() * 16.0) - 8.0) * generator.nextDouble();

+
+            double tst = FastMath.sinh(x);
+            double ref = sinh(field.newDfp(x)).toDouble();
+            double err = (tst - ref) / ref;
+
+            if (err != 0) {
+                double ulp = Math.abs(ref - Double.longBitsToDouble((Double.doubleToLongBits(ref)
^ 1)));
+                double errulp = field.newDfp(tst).subtract(sinh(field.newDfp(x))).divide(field.newDfp(ulp)).toDouble();

+                //System.out.println(x+"\t"+tst+"\t"+ref+"\t"+err+"\t"+errulp); 
+                maxerrulp = Math.max(maxerrulp, Math.abs(errulp));
+            }
+        }
+
+        Assert.assertTrue("sinh() had errors in excess of " + MAX_ERROR_ULP + " ULP", maxerrulp
< MAX_ERROR_ULP);
+    }
+
+    @Test
+    public void testCoshAccuracy() {
+        double maxerrulp = 0.0;
+
+        for (int i=0; i<10000; i++) {
+            double x = ((generator.nextDouble() * 16.0) - 8.0) * generator.nextDouble();

+
+            double tst = FastMath.cosh(x);
+            double ref = cosh(field.newDfp(x)).toDouble();
+            double err = (tst - ref) / ref;
+
+            if (err != 0) {
+                double ulp = Math.abs(ref - Double.longBitsToDouble((Double.doubleToLongBits(ref)
^ 1)));
+                double errulp = field.newDfp(tst).subtract(cosh(field.newDfp(x))).divide(field.newDfp(ulp)).toDouble();

+                //System.out.println(x+"\t"+tst+"\t"+ref+"\t"+err+"\t"+errulp); 
+                maxerrulp = Math.max(maxerrulp, Math.abs(errulp));
+            }
+        }
+
+        Assert.assertTrue("cosh() had errors in excess of " + MAX_ERROR_ULP + " ULP", maxerrulp
< MAX_ERROR_ULP);
+    }
+
+    @Test
+    public void testTanhAccuracy() {
+        double maxerrulp = 0.0;
+
+        for (int i=0; i<10000; i++) {
+            double x = ((generator.nextDouble() * 16.0) - 8.0) * generator.nextDouble();

+
+            double tst = FastMath.tanh(x);
+            double ref = tanh(field.newDfp(x)).toDouble();
+            double err = (tst - ref) / ref;
+
+            if (err != 0) {
+                double ulp = Math.abs(ref - Double.longBitsToDouble((Double.doubleToLongBits(ref)
^ 1)));
+                double errulp = field.newDfp(tst).subtract(tanh(field.newDfp(x))).divide(field.newDfp(ulp)).toDouble();

+                //System.out.println(x+"\t"+tst+"\t"+ref+"\t"+err+"\t"+errulp); 
+                maxerrulp = Math.max(maxerrulp, Math.abs(errulp));
+            }
+        }
+
+        Assert.assertTrue("tanh() had errors in excess of " + MAX_ERROR_ULP + " ULP", maxerrulp
< MAX_ERROR_ULP);
+    }
+
     @Test
     public void testCbrtAccuracy() {
         double maxerrulp = 0.0;
@@ -1017,7 +1095,49 @@ public class FastMathTest {
             time = System.currentTimeMillis() - time;
             System.out.println("FastMath.cbrt " + time + "\t" + x);
 
-           x = 0;
+            x = 0;
+            time = System.currentTimeMillis();
+            for (int i = 0; i < numberOfRuns; i++)
+                x += StrictMath.cosh(i / 1000000.0);
+            time = System.currentTimeMillis() - time;
+            System.out.print("StrictMath.cosh " + time + "\t" + x + "\t");
+
+            x = 0;
+            time = System.currentTimeMillis();
+            for (int i = 0; i < numberOfRuns; i++)
+                x += FastMath.cosh(i / 1000000.0);
+            time = System.currentTimeMillis() - time;
+            System.out.println("FastMath.cosh " + time + "\t" + x);
+
+            x = 0;
+            time = System.currentTimeMillis();
+            for (int i = 0; i < numberOfRuns; i++)
+                x += StrictMath.sinh(i / 1000000.0);
+            time = System.currentTimeMillis() - time;
+            System.out.print("StrictMath.sinh " + time + "\t" + x + "\t");
+
+            x = 0;
+            time = System.currentTimeMillis();
+            for (int i = 0; i < numberOfRuns; i++)
+                x += FastMath.sinh(i / 1000000.0);
+            time = System.currentTimeMillis() - time;
+            System.out.println("FastMath.sinh " + time + "\t" + x);
+
+            x = 0;
+            time = System.currentTimeMillis();
+            for (int i = 0; i < numberOfRuns; i++)
+                x += StrictMath.tanh(i / 1000000.0);
+            time = System.currentTimeMillis() - time;
+            System.out.print("StrictMath.tanh " + time + "\t" + x + "\t");
+
+            x = 0;
+            time = System.currentTimeMillis();
+            for (int i = 0; i < numberOfRuns; i++)
+                x += FastMath.tanh(i / 1000000.0);
+            time = System.currentTimeMillis() - time;
+            System.out.println("FastMath.tanh " + time + "\t" + x);
+
+            x = 0;
             time = System.currentTimeMillis();
             for (int i = 0; i < numberOfRuns; i++)
                 x += StrictMath.expm1(-i / 100000.0);
@@ -1037,6 +1157,7 @@ public class FastMathTest {
                 x += FastMath.expm1(-i / 100000.0);
             time = System.currentTimeMillis() - time;
             System.out.println("FastMath.expm1 " + time + "\t" + x);
+
         }
     }
 



Mime
View raw message