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. * *
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}. + * + *
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) * diff -r 1586fc6697da -r b65c2f5d4d01 jdk/src/java.base/share/classes/java/lang/StrictMath.java --- 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 ex -1. Note that for values of diff -r 1586fc6697da -r b65c2f5d4d01 jdk/src/java.base/share/native/libfdlibm/e_hypot.c --- 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; -} diff -r 1586fc6697da -r b65c2f5d4d01 jdk/src/java.base/share/native/libfdlibm/w_hypot.c --- 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 -} diff -r 1586fc6697da -r b65c2f5d4d01 jdk/src/java.base/share/native/libjava/StrictMath.c --- 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); diff -r 1586fc6697da -r b65c2f5d4d01 jdk/test/java/lang/StrictMath/HypotTests.java --- 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(); }