stan/jdk-24
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

8138823: Correct bug in port of fdlibm hypot to Java

Browse Source 
Reviewed-by: bpb
This commit is contained in:
Joe Darcy 2015-10-06 18:39:26 -07:00
parent d12b5c154b
commit f674b491a6
4 changed files with 254 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
jdk/src/java.base/share/classes/java/lang/FdLibm.java
View File

@ -182,7 +182,7 @@ class FdLibm {
}
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;

192
jdk/test/java/lang/StrictMath/FdlibmTranslit.java Normal file
View File

@ -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;
}
}
}

55
jdk/test/java/lang/StrictMath/HypotTests.java
View File

@ -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 @@ public class HypotTests {
{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 @@ public class HypotTests {
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;
}
}

14
jdk/test/java/lang/StrictMath/Tests.java
View File

@ -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 @@ public class Tests {
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);
}
}