changeset 12856:f1799292ea4b

8138823: Correct bug in port of fdlibm hypot to Java Reviewed-by: bpb
author darcy
date Tue, 06 Oct 2015 18:39:26 -0700
parents fe40b31c0e52
children ecfc81a3741c
files src/java.base/share/classes/java/lang/FdLibm.java test/java/lang/StrictMath/FdlibmTranslit.java test/java/lang/StrictMath/HypotTests.java test/java/lang/StrictMath/Tests.java
diffstat 4 files changed, 254 insertions(+), 9 deletions(-) [+]
line wrap: on
line diff
--- a/src/java.base/share/classes/java/lang/FdLibm.java	Tue Oct 06 13:20:21 2015 -0700
+++ b/src/java.base/share/classes/java/lang/FdLibm.java	Tue Oct 06 18:39:26 2015 -0700
@@ -182,7 +182,7 @@
             }
 
             int k = 0;
-            if (a > 0x1.0p500) {   // a > 2**500
+            if (a > 0x1.00000_ffff_ffffp500) {   // a > ~2**500
                 // scale a and b by 2**-600
                 ha -= 0x25800000;
                 hb -= 0x25800000;
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/java/lang/StrictMath/FdlibmTranslit.java	Tue Oct 06 18:39:26 2015 -0700
@@ -0,0 +1,192 @@
+/*
+ * Copyright (c) 1998, 2015, Oracle and/or its affiliates. All rights reserved.
+ * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
+ *
+ * This code is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License version 2 only, as
+ * published by the Free Software Foundation.  Oracle designates this
+ * particular file as subject to the "Classpath" exception as provided
+ * by Oracle in the LICENSE file that accompanied this code.
+ *
+ * This code 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
+ * version 2 for more details (a copy is included in the LICENSE file that
+ * accompanied this code).
+ *
+ * You should have received a copy of the GNU General Public License version
+ * 2 along with this work; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
+ * or visit www.oracle.com if you need additional information or have any
+ * questions.
+ */
+
+/**
+ * A transliteration of the "Freely Distributable Math Library"
+ * algorithms from C into Java. That is, this port of the algorithms
+ * is as close to the C originals as possible while still being
+ * readable legal Java.
+ */
+public class FdlibmTranslit {
+    private FdlibmTranslit() {
+        throw new UnsupportedOperationException("No FdLibmTranslit instances for you.");
+    }
+
+    /**
+     * Return the low-order 32 bits of the double argument as an int.
+     */
+    private static int __LO(double x) {
+        long transducer = Double.doubleToRawLongBits(x);
+        return (int)transducer;
+    }
+
+    /**
+     * Return a double with its low-order bits of the second argument
+     * and the high-order bits of the first argument..
+     */
+    private static double __LO(double x, int low) {
+        long transX = Double.doubleToRawLongBits(x);
+        return Double.longBitsToDouble((transX & 0xFFFF_FFFF_0000_0000L)|low );
+    }
+
+    /**
+     * Return the high-order 32 bits of the double argument as an int.
+     */
+    private static int __HI(double x) {
+        long transducer = Double.doubleToRawLongBits(x);
+        return (int)(transducer >> 32);
+    }
+
+    /**
+     * Return a double with its high-order bits of the second argument
+     * and the low-order bits of the first argument..
+     */
+    private static double __HI(double x, int high) {
+        long transX = Double.doubleToRawLongBits(x);
+        return Double.longBitsToDouble((transX & 0x0000_0000_FFFF_FFFFL)|( ((long)high)) << 32 );
+    }
+
+    public static double hypot(double x, double y) {
+        return Hypot.compute(x, y);
+    }
+
+    /**
+     * hypot(x,y)
+     *
+     * Method :
+     *      If (assume round-to-nearest) z = x*x + y*y
+     *      has error less than sqrt(2)/2 ulp, than
+     *      sqrt(z) has error less than 1 ulp (exercise).
+     *
+     *      So, compute sqrt(x*x + y*y) with some care as
+     *      follows to get the error below 1 ulp:
+     *
+     *      Assume x > y > 0;
+     *      (if possible, set rounding to round-to-nearest)
+     *      1. if x > 2y  use
+     *              x1*x1 + (y*y + (x2*(x + x1))) for x*x + y*y
+     *      where x1 = x with lower 32 bits cleared, x2 = x - x1; else
+     *      2. if x <= 2y use
+     *              t1*y1 + ((x-y) * (x-y) + (t1*y2 + t2*y))
+     *      where t1 = 2x with lower 32 bits cleared, t2 = 2x - t1,
+     *      y1= y with lower 32 bits chopped, y2 = y - y1.
+     *
+     *      NOTE: scaling may be necessary if some argument is too
+     *            large or too tiny
+     *
+     * Special cases:
+     *      hypot(x,y) is INF if x or y is +INF or -INF; else
+     *      hypot(x,y) is NAN if x or y is NAN.
+     *
+     * Accuracy:
+     *      hypot(x,y) returns sqrt(x^2 + y^2) with error less
+     *      than 1 ulps (units in the last place)
+     */
+    static class Hypot {
+        public static double compute(double x, double y) {
+            double a = x;
+            double b = y;
+            double t1, t2, y1, y2, w;
+            int j, k, ha, hb;
+
+            ha = __HI(x) & 0x7fffffff;        // high word of  x
+            hb = __HI(y) & 0x7fffffff;        // high word of  y
+            if(hb > ha) {
+                a = y;
+                b = x;
+                j = ha;
+                ha = hb;
+                hb = j;
+            } else {
+                a = x;
+                b = y;
+            }
+            a = __HI(a, ha);   // a <- |a|
+            b = __HI(b, hb);   // b <- |b|
+            if ((ha - hb) > 0x3c00000) {
+                return a + b;  // x / y > 2**60
+            }
+            k=0;
+            if (ha > 0x5f300000) {   // a>2**500
+                if (ha >= 0x7ff00000) {       // Inf or NaN
+                    w = a + b;                // for sNaN
+                    if (((ha & 0xfffff) | __LO(a)) == 0)
+                        w = a;
+                    if (((hb ^ 0x7ff00000) | __LO(b)) == 0)
+                        w = b;
+                    return w;
+                }
+                // scale a and b by 2**-600
+                ha -= 0x25800000;
+                hb -= 0x25800000;
+                k += 600;
+                a = __HI(a, ha);
+                b = __HI(b, hb);
+            }
+            if (hb < 0x20b00000) {   // b < 2**-500
+                if (hb <= 0x000fffff) {      // subnormal b or 0 */
+                    if ((hb | (__LO(b))) == 0)
+                        return a;
+                    t1 = 0;
+                    t1 = __HI(t1, 0x7fd00000);  // t1=2^1022
+                    b *= t1;
+                    a *= t1;
+                    k -= 1022;
+                } else {            // scale a and b by 2^600
+                    ha += 0x25800000;       // a *= 2^600
+                    hb += 0x25800000;       // b *= 2^600
+                    k -= 600;
+                    a = __HI(a, ha);
+                    b = __HI(b, hb);
+                }
+            }
+            // medium size a and b
+            w = a - b;
+            if (w > b) {
+                t1 = 0;
+                t1 = __HI(t1, ha);
+                t2 = a - t1;
+                w  = Math.sqrt(t1*t1 - (b*(-b) - t2 * (a + t1)));
+            } else {
+                a  = a + a;
+                y1 = 0;
+                y1 = __HI(y1, hb);
+                y2 = b - y1;
+                t1 = 0;
+                t1 = __HI(t1, ha + 0x00100000);
+                t2 = a - t1;
+                w  = Math.sqrt(t1*y1 - (w*(-w) - (t1*y2 + t2*b)));
+            }
+            if (k != 0) {
+                t1 = 1.0;
+                int t1_hi = __HI(t1);
+                t1_hi += (k << 20);
+                t1 = __HI(t1, t1_hi);
+                return t1 * w;
+            } else
+                return w;
+        }
+    }
+}
--- a/test/java/lang/StrictMath/HypotTests.java	Tue Oct 06 13:20:21 2015 -0700
+++ b/test/java/lang/StrictMath/HypotTests.java	Tue Oct 06 18:39:26 2015 -0700
@@ -24,10 +24,19 @@
 /*
  * @test
  * @bug 4851638
+ * @key randomness
  * @summary Tests for StrictMath.hypot
+ * @library /lib/testlibrary/
+ * @build jdk.testlibrary.*
+ * @build Tests
+ * @build FdlibmTranslit
+ * @build HypotTests
+ * @run main HypotTests
  * @author Joseph D. Darcy
  */
 
+import jdk.testlibrary.RandomFactory;
+
 /**
  * The tests in ../Math/HypotTests.java test properties that should
  * hold for any hypot implementation, including the FDLIBM-based one
@@ -42,6 +51,19 @@
 public class HypotTests {
     private HypotTests(){}
 
+    public static void main(String... args) {
+        int failures = 0;
+
+        failures += testHypot();
+        failures += testAgainstTranslit();
+
+        if (failures > 0) {
+            System.err.println("Testing hypot incurred "
+                               + failures + " failures.");
+            throw new RuntimeException();
+        }
+    }
+
     /**
      * The hypot implementation is commutative, {@code hypot(a, b) ==
      * hypot(b, a)}, and independent of sign, {@code hypot(a, b) ==
@@ -663,6 +685,12 @@
             {0x1.0p-450,              0x1.fffffffffffffp-499,  0x1.0p-450},
             {0x1.0000000000001p-450,  0x1.fffffffffffffp-499,  0x1.0000000000001p-450},
 
+            {0x1.00000_ffff_0000p500,  0x1.fffffffffffffp499,  0x1.6a09f1b837ccfp500},
+            {0x1.00000_0000_0001p500,  0x1.fffffffffffffp499,  0x1.6a09e667f3bcdp500},
+            {0x1.00000_ffff_ffffp500,  0x1.fffffffffffffp499,  0x1.6a09f1b8431d3p500},
+            {0x1.00001_0000_0000p500,  0x1.fffffffffffffp499,  0x1.6a09f1b8431d5p500},
+
+
             // 0x1.0p-1022 is MIN_NORMAL
             {0x1.0000000000001p-1022, 0x1.0000000000001p-1022, 0x1.6a09e667f3bcep-1022},
             {0x1.0000000000001p-1022, 0x1.0p-1022,             0x1.6a09e667f3bcdp-1022},
@@ -686,15 +714,30 @@
         return failures;
     }
 
-    public static void main(String... args) {
+    // Initialize shared random number generator
+    private static java.util.Random random = RandomFactory.getRandom();
+
+    /**
+     * Test StrictMath.hypot against transliteration port of hypot.
+     */
+    private static int testAgainstTranslit() {
         int failures = 0;
+        double x = Tests.createRandomDouble(random);
+        double y = Tests.createRandomDouble(random);
 
-        failures += testHypot();
+        // Make the increment twice the ulp value in case the random
+        // value is near an exponent threshold.
+        double increment_x = 2.0 * Math.ulp(x);
+        double increment_y = 2.0 * Math.ulp(y);
 
-        if (failures > 0) {
-            System.err.println("Testing hypot incurred "
-                               + failures + " failures.");
-            throw new RuntimeException();
+        // Don't worry about x or y overflowing to infinity if their
+        // exponent is MAX_EXPONENT.
+        for (int i = 0; i < 200; i++, x += increment_x) {
+            for (int j = 0; j < 200; j++, y += increment_y) {
+                failures += testHypotCase(x, y, FdlibmTranslit.hypot(x, y));
+            }
         }
+
+        return failures;
     }
 }
--- a/test/java/lang/StrictMath/Tests.java	Tue Oct 06 13:20:21 2015 -0700
+++ b/test/java/lang/StrictMath/Tests.java	Tue Oct 06 18:39:26 2015 -0700
@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2003, Oracle and/or its affiliates. All rights reserved.
+ * Copyright (c) 2003, 2015, Oracle and/or its affiliates. All rights reserved.
  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
  *
  * This code is free software; you can redistribute it and/or modify it
@@ -60,6 +60,16 @@
             return 0;
     }
 
+    /**
+     * Returns a double over the normalized range of floating-point values.
+     * @return a double over the normalized range of floating-point values
+     */
+    static double createRandomDouble(java.util.Random random) {
+        final int EXPONENT_RANGE = Double.MAX_EXPONENT - Double.MIN_EXPONENT + 1;
 
-
+        int targetExponent = Double.MIN_EXPONENT + random.nextInt(EXPONENT_RANGE + 1);
+        double tmp = random.nextDouble(); // Double in the range of [0.0, 1.0)
+        int tmpExponent = Math.getExponent(tmp);
+        return Math.scalb(tmp, targetExponent - tmpExponent);
+    }
 }