--- a/jdk/src/java.base/share/classes/java/lang/FdLibm.java Wed Jul 05 20:52:37 2017 +0200
+++ b/jdk/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/jdk/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/jdk/test/java/lang/StrictMath/HypotTests.java Wed Jul 05 20:52:37 2017 +0200
+++ b/jdk/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;
-
- failures += testHypot();
+ double x = Tests.createRandomDouble(random);
+ double y = Tests.createRandomDouble(random);
- if (failures > 0) {
- System.err.println("Testing hypot incurred "
- + failures + " failures.");
- throw new RuntimeException();
+ // 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);
+
+ // 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/jdk/test/java/lang/StrictMath/Tests.java Wed Jul 05 20:52:37 2017 +0200
+++ b/jdk/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);
+ }
}