changeset 12809:8292f92a37b7

7130085: Port fdlibm hypot to Java Reviewed-by: bpb
author darcy
date Wed, 23 Sep 2015 14:14:14 -0700
parents aa7cccf1a672
children c9b4bc199dca
files make/mapfiles/libjava/mapfile-vers src/java.base/share/classes/java/lang/FdLibm.java src/java.base/share/classes/java/lang/StrictMath.java src/java.base/share/native/libfdlibm/e_hypot.c src/java.base/share/native/libfdlibm/w_hypot.c src/java.base/share/native/libjava/StrictMath.c test/java/lang/StrictMath/HypotTests.java
diffstat 7 files changed, 222 insertions(+), 200 deletions(-) [+]
line wrap: on
line diff
--- a/make/mapfiles/libjava/mapfile-vers	Wed Sep 23 15:02:46 2015 -0400
+++ b/make/mapfiles/libjava/mapfile-vers	Wed Sep 23 14:14:14 2015 -0700
@@ -157,7 +157,6 @@
 		Java_java_lang_StrictMath_cosh;
 		Java_java_lang_StrictMath_sinh;
 		Java_java_lang_StrictMath_tanh;
-		Java_java_lang_StrictMath_hypot;
 		Java_java_lang_StrictMath_log1p;
 		Java_java_lang_StrictMath_expm1;
 		Java_java_lang_Object_getClass;
--- a/src/java.base/share/classes/java/lang/FdLibm.java	Wed Sep 23 15:02:46 2015 -0400
+++ b/src/java.base/share/classes/java/lang/FdLibm.java	Wed Sep 23 14:14:14 2015 -0700
@@ -26,7 +26,8 @@
 package java.lang;
 
 /**
- * Port of the "Freely Distributable Math Library", version 5.3, from C to Java.
+ * Port of the "Freely Distributable Math Library", version 5.3, from
+ * C to Java.
  *
  * <p>The C version of fdlibm relied on the idiom of pointer aliasing
  * a 64-bit double floating-point value as a two-element array of
@@ -37,7 +38,7 @@
  * operated on as integer values, the standard library methods for
  * bitwise floating-point to integer conversion,
  * Double.longBitsToDouble and Double.doubleToRawLongBits, are directly
- * or indirectly used .
+ * or indirectly used.
  *
  * <p>The C version of fdlibm also took some pains to signal the
  * correct IEEE 754 exceptional conditions divide by zero, invalid,
@@ -47,13 +48,21 @@
  * handling is not supported natively in the JVM, such coding patterns
  * have been omitted from this port. For example, rather than {@code
  * return huge * huge}, this port will use {@code return INFINITY}.
+ *
+ * <p>Various comparison and arithmetic operations in fdlibm could be
+ * done either based on the integer view of a value or directly on the
+ * floating-point representation. Which idiom is faster may depend on
+ * platform specific factors. However, for code clarity if no other
+ * reason, this port will favor expressing the semantics of those
+ * operations in terms of floating-point operations when convenient to
+ * do so.
  */
 class FdLibm {
     // Constants used by multiple algorithms
     private static final double INFINITY = Double.POSITIVE_INFINITY;
 
     private FdLibm() {
-        throw new UnsupportedOperationException("No instances for you.");
+        throw new UnsupportedOperationException("No FdLibm instances for you.");
     }
 
     /**
@@ -91,13 +100,146 @@
     }
 
     /**
+     * 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 ulp (unit in the last place)
+     */
+    public static class Hypot {
+        public static final double TWO_MINUS_600 = 0x1.0p-600;
+        public static final double TWO_PLUS_600  = 0x1.0p+600;
+
+        public static strictfp double compute(double x, double y) {
+            double a = Math.abs(x);
+            double b = Math.abs(y);
+
+            if (!Double.isFinite(a) || !Double.isFinite(b)) {
+                if (a == INFINITY || b == INFINITY)
+                    return INFINITY;
+                else
+                    return a + b; // Propagate NaN significand bits
+            }
+
+            if (b > a) {
+                double tmp = a;
+                a = b;
+                b = tmp;
+            }
+            assert a >= b;
+
+            // Doing bitwise conversion after screening for NaN allows
+            // the code to not worry about the possibility of
+            // "negative" NaN values.
+
+            // Note: the ha and hb variables are the high-order
+            // 32-bits of a and b stored as integer values. The ha and
+            // hb values are used first for a rough magnitude
+            // comparison of a and b and second for simulating higher
+            // precision by allowing a and b, respectively, to be
+            // decomposed into non-overlapping portions. Both of these
+            // uses could be eliminated. The magnitude comparison
+            // could be eliminated by extracting and comparing the
+            // exponents of a and b or just be performing a
+            // floating-point divide.  Splitting a floating-point
+            // number into non-overlapping portions can be
+            // accomplished by judicious use of multiplies and
+            // additions. For details see T. J. Dekker, A Floating
+            // Point Technique for Extending the Available Precision ,
+            // Numerische Mathematik, vol. 18, 1971, pp.224-242 and
+            // subsequent work.
+
+            int ha = __HI(a);        // high word of a
+            int hb = __HI(b);        // high word of b
+
+            if ((ha - hb) > 0x3c00000) {
+                return a + b;  // x / y > 2**60
+            }
+
+            int k = 0;
+            if (a > 0x1.0p500) {   // a > 2**500
+                // scale a and b by 2**-600
+                ha -= 0x25800000;
+                hb -= 0x25800000;
+                a = a * TWO_MINUS_600;
+                b = b * TWO_MINUS_600;
+                k += 600;
+            }
+            double t1, t2;
+            if (b < 0x1.0p-500) {   // b < 2**-500
+                if (b < Double.MIN_NORMAL) {      // subnormal b or 0 */
+                    if (b == 0.0)
+                        return a;
+                    t1 = 0x1.0p1022;   // 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
+                    a = a * TWO_PLUS_600;
+                    b = b * TWO_PLUS_600;
+                    k -= 600;
+                }
+            }
+            // medium size a and b
+            double 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 {
+                double y1, y2;
+                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) {
+                return Math.powerOfTwoD(k) * w;
+            } else
+                return w;
+        }
+    }
+
+    /**
      * Compute x**y
      *                    n
      * Method:  Let x =  2   * (1+f)
      *      1. Compute and return log2(x) in two pieces:
      *              log2(x) = w1 + w2,
      *         where w1 has 53 - 24 = 29 bit trailing zeros.
-     *      2. Perform y*log2(x) = n+y' by simulating muti-precision
+     *      2. Perform y*log2(x) = n+y' by simulating multi-precision
      *         arithmetic, where |y'| <= 0.5.
      *      3. Return x**y = 2**n*exp(y'*log2)
      *
--- a/src/java.base/share/classes/java/lang/StrictMath.java	Wed Sep 23 15:02:46 2015 -0400
+++ b/src/java.base/share/classes/java/lang/StrictMath.java	Wed Sep 23 14:14:14 2015 -0700
@@ -1329,7 +1329,9 @@
      * without intermediate overflow or underflow
      * @since 1.5
      */
-    public static native double hypot(double x, double y);
+    public static double hypot(double x, double y) {
+        return FdLibm.Hypot.compute(x, y);
+    }
 
     /**
      * Returns <i>e</i><sup>x</sup>&nbsp;-1.  Note that for values of
--- a/src/java.base/share/native/libfdlibm/e_hypot.c	Wed Sep 23 15:02:46 2015 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,128 +0,0 @@
-
-/*
- * Copyright (c) 1998, 2001, 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.
- */
-
-/* __ieee754_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)
- */
-
-#include "fdlibm.h"
-
-#ifdef __STDC__
-        double __ieee754_hypot(double x, double y)
-#else
-        double __ieee754_hypot(x,y)
-        double x, y;
-#endif
-{
-        double a=x,b=y,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;}
-        __HI(a) = ha;   /* a <- |a| */
-        __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;
-           __HI(a) = ha;
-           __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;
-                __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;
-                __HI(a) = ha;
-                __HI(b) = hb;
-            }
-        }
-    /* medium size a and b */
-        w = a-b;
-        if (w>b) {
-            t1 = 0;
-            __HI(t1) = ha;
-            t2 = a-t1;
-            w  = sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
-        } else {
-            a  = a+a;
-            y1 = 0;
-            __HI(y1) = hb;
-            y2 = b - y1;
-            t1 = 0;
-            __HI(t1) = ha+0x00100000;
-            t2 = a - t1;
-            w  = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
-        }
-        if(k!=0) {
-            t1 = 1.0;
-            __HI(t1) += (k<<20);
-            return t1*w;
-        } else return w;
-}
--- a/src/java.base/share/native/libfdlibm/w_hypot.c	Wed Sep 23 15:02:46 2015 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,52 +0,0 @@
-
-/*
- * Copyright (c) 1998, 2001, 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.
- */
-
-/*
- * wrapper hypot(x,y)
- */
-
-#include "fdlibm.h"
-
-
-#ifdef __STDC__
-        double hypot(double x, double y)/* wrapper hypot */
-#else
-        double hypot(x,y)               /* wrapper hypot */
-        double x,y;
-#endif
-{
-#ifdef _IEEE_LIBM
-        return __ieee754_hypot(x,y);
-#else
-        double z;
-        z = __ieee754_hypot(x,y);
-        if(_LIB_VERSION == _IEEE_) return z;
-        if((!finite(z))&&finite(x)&&finite(y))
-            return __kernel_standard(x,y,4); /* hypot overflow */
-        else
-            return z;
-#endif
-}
--- a/src/java.base/share/native/libjava/StrictMath.c	Wed Sep 23 15:02:46 2015 -0400
+++ b/src/java.base/share/native/libjava/StrictMath.c	Wed Sep 23 14:14:14 2015 -0700
@@ -127,14 +127,6 @@
 }
 
 JNIEXPORT jdouble JNICALL
-Java_java_lang_StrictMath_hypot(JNIEnv *env, jclass unused, jdouble x, jdouble y)
-{
-    return (jdouble) jhypot((double)x, (double)y);
-}
-
-
-
-JNIEXPORT jdouble JNICALL
 Java_java_lang_StrictMath_log1p(JNIEnv *env, jclass unused, jdouble d)
 {
     return (jdouble) jlog1p((double)d);
--- a/test/java/lang/StrictMath/HypotTests.java	Wed Sep 23 15:02:46 2015 -0400
+++ b/test/java/lang/StrictMath/HypotTests.java	Wed Sep 23 14:14:14 2015 -0700
@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2003, 2004, 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
@@ -42,9 +42,37 @@
 public class HypotTests {
     private HypotTests(){}
 
+    /**
+     * The hypot implementation is commutative, {@code hypot(a, b) ==
+     * hypot(b, a)}, and independent of sign, {@code hypot(a, b) ==
+     * hypot(-a, b) == hypot(a, -b) == hypot(-a, -b)}.
+     */
     static int testHypotCase(double input1, double input2, double expected) {
-        return Tests.test("StrictMath.hypot(double)", input1, input2,
-                          StrictMath.hypot(input1, input2), expected);
+        int failures = 0;
+        failures += Tests.test("StrictMath.hypot(double)", input1, input2,
+                               StrictMath.hypot(input1, input2), expected);
+
+        failures += Tests.test("StrictMath.hypot(double)", input2, input1,
+                          StrictMath.hypot(input2, input1), expected);
+
+        failures += Tests.test("StrictMath.hypot(double)", -input1, input2,
+                               StrictMath.hypot(-input1, input2), expected);
+
+        failures += Tests.test("StrictMath.hypot(double)", input2, -input1,
+                          StrictMath.hypot(input2, -input1), expected);
+
+        failures += Tests.test("StrictMath.hypot(double)", input1, -input2,
+                               StrictMath.hypot(input1, -input2), expected);
+
+        failures += Tests.test("StrictMath.hypot(double)", -input2, input1,
+                          StrictMath.hypot(-input2, input1), expected);
+
+        failures += Tests.test("StrictMath.hypot(double)", -input1, -input2,
+                               StrictMath.hypot(-input1, -input2), expected);
+
+        failures +=  Tests.test("StrictMath.hypot(double)", -input2, -input1,
+                          StrictMath.hypot(-input2, -input1), expected);
+        return failures;
     }
 
     static int testHypot() {
@@ -611,21 +639,60 @@
             {0x1.8p1,   0x1.8bffffffffff6p6,    0x1.8c2e88e6f44b1p6},
             {0x1.8p1,   0x1.8ffffffffffe8p6,    0x1.902e11d3b5549p6},
             {0x1.8p1,   0x1.8fffffffffffep6,    0x1.902e11d3b556p6},
+
+            // Test near decision points of the fdlibm algorithm
+            {0x1.0000000000001p501,   0x1.000000000000p501,    0x1.6a09e667f3bcdp501},
+            {0x1.0p501,               0x1.0p499,               0x1.07e0f66afed07p501},
+
+            {0x1.0p500,               0x1.0p450,               0x1.0p500},
+            {0x1.0000000000001p500,   0x1.0p450,               0x1.0000000000001p500},
+
+            {0x1.0p500,               0x1.0p440,               0x1.0p500},
+            {0x1.0000000000001p500,   0x1.0p440,               0x1.0000000000001p500},
+            {0x1.0p500,               0x1.0p439,               0x1.0p500},
+            {0x1.0000000000001p500,   0x1.0p439,               0x1.0000000000001p500},
+
+            {0x1.0p-450,              0x1.0p-500,              0x1.0p-450},
+            {0x1.0000000000001p-450,  0x1.0p-500,              0x1.0000000000001p-450},
+            {0x1.0p-450,              0x1.fffffffffffffp-499,  0x1.0p-450},
+            {0x1.0000000000001p-450,  0x1.fffffffffffffp-499,  0x1.0000000000001p-450},
+
+
+            {0x1.0p-450,              0x1.0p-500,              0x1.0p-450},
+            {0x1.0000000000001p-450,  0x1.0p-500,              0x1.0000000000001p-450},
+            {0x1.0p-450,              0x1.fffffffffffffp-499,  0x1.0p-450},
+            {0x1.0000000000001p-450,  0x1.fffffffffffffp-499,  0x1.0000000000001p-450},
+
+            // 0x1.0p-1022 is MIN_NORMAL
+            {0x1.0000000000001p-1022, 0x1.0000000000001p-1022, 0x1.6a09e667f3bcep-1022},
+            {0x1.0000000000001p-1022, 0x1.0p-1022,             0x1.6a09e667f3bcdp-1022},
+            {0x1.0000000000001p-1022, 0x0.fffffffffffffp-1022, 0x1.6a09e667f3bcdp-1022},
+            {0x1.0000000000001p-1022, 0x0.0000000000001P-1022, 0x1.0000000000001p-1022},
+            {0x1.0000000000001p-1022, 0.0,                     0x1.0000000000001p-1022},
+
+            {0x1.0000000000000p-1022, 0x0.fffffffffffffp-1022, 0x1.6a09e667f3bccp-1022},
+            {0x1.0000000000000p-1021, 0x0.fffffffffffffp-1022, 0x1.1e3779b97f4a8p-1021},
+            {0x1.0000000000000p-1020, 0x0.fffffffffffffp-1022, 0x1.07e0f66afed07p-1020},
+
+            // 0x0.0000000000001P-1022 is MIN_VALUE (smallest nonzero number)
+            {0x0.0000000000001p-1022, 0x0.0000000000001p-1022, 0x0.0000000000001p-1022},
+            {0x0.0000000000002p-1022, 0x0.0000000000001p-1022, 0x0.0000000000002p-1022},
+            {0x0.0000000000003p-1022, 0x0.0000000000002p-1022, 0x0.0000000000004p-1022},
         };
 
         for (double[] testCase: testCases)
-            failures+=testHypotCase(testCase[0], testCase[1], testCase[2]);
+            failures += testHypotCase(testCase[0], testCase[1], testCase[2]);
 
         return failures;
     }
 
-    public static void main(String [] argv) {
+    public static void main(String... args) {
         int failures = 0;
 
         failures += testHypot();
 
         if (failures > 0) {
-            System.err.println("Testing log1p incurred "
+            System.err.println("Testing hypot incurred "
                                + failures + " failures.");
             throw new RuntimeException();
         }