--- a/jdk/make/mapfiles/libjava/mapfile-vers Wed Sep 23 15:02:46 2015 -0400
+++ b/jdk/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/jdk/src/java.base/share/classes/java/lang/FdLibm.java Wed Sep 23 15:02:46 2015 -0400
+++ b/jdk/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/jdk/src/java.base/share/classes/java/lang/StrictMath.java Wed Sep 23 15:02:46 2015 -0400
+++ b/jdk/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> -1. Note that for values of
--- a/jdk/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/jdk/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/jdk/src/java.base/share/native/libjava/StrictMath.c Wed Sep 23 15:02:46 2015 -0400
+++ b/jdk/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/jdk/test/java/lang/StrictMath/HypotTests.java Wed Sep 23 15:02:46 2015 -0400
+++ b/jdk/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();
}