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

8301205: Port fdlibm log10 to Java

Browse Source 
Reviewed-by: bpb
This commit is contained in:
Joe Darcy 2023-01-30 20:33:01 +00:00
parent b84f4c40fd
commit 63bb2ce8de
4 changed files with 222 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

79
src/java.base/share/classes/java/lang/FdLibm.java
View File

@ -1,5 +1,5 @@
/*
* Copyright (c) 1998, 2021, Oracle and/or its affiliates. All rights reserved.
* Copyright (c) 1998, 2023, 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
@ -745,4 +745,81 @@ class FdLibm {
}
}
}
/**
* Return the base 10 logarithm of x
*
* Method :
* Let log10_2hi = leading 40 bits of log10(2) and
* log10_2lo = log10(2) - log10_2hi,
* ivln10 = 1/log(10) rounded.
* Then
* n = ilogb(x),
* if(n<0) n = n+1;
* x = scalbn(x,-n);
* log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))
*
* Note 1:
* To guarantee log10(10**n)=n, where 10**n is normal, the rounding
* mode must set to Round-to-Nearest.
* Note 2:
* [1/log(10)] rounded to 53 bits has error .198 ulps;
* log10 is monotonic at all binary break points.
*
* Special cases:
* log10(x) is NaN with signal if x < 0;
* log10(+INF) is +INF with no signal; log10(0) is -INF with signal;
* log10(NaN) is that NaN with no signal;
* log10(10**N) = N for N=0,1,...,22.
*
* Constants:
* The hexadecimal values are the intended ones for the following constants.
* The decimal values may be used, provided that the compiler will convert
* from decimal to binary accurately enough to produce the hexadecimal values
* shown.
*/
static class Log10 {
private static double two54 = 0x1.0p54; // 1.80143985094819840000e+16;
private static double ivln10 = 0x1.bcb7b1526e50ep-2; // 4.34294481903251816668e-01
private static double log10_2hi = 0x1.34413509f6p-2; // 3.01029995663611771306e-01;
private static double log10_2lo = 0x1.9fef311f12b36p-42; // 3.69423907715893078616e-13;
private Log10() {
throw new UnsupportedOperationException();
}
public static double compute(double x) {
double y, z;
int i, k;
int hx = __HI(x); // high word of x
int lx = __LO(x); // low word of x
k=0;
if (hx < 0x0010_0000) { /* x < 2**-1022 */
if (((hx & 0x7fff_ffff) | lx) == 0) {
return -two54/0.0; /* log(+-0)=-inf */
}
if (hx < 0) {
return (x - x)/0.0; /* log(-#) = NaN */
}
k -= 54;
x *= two54; /* subnormal number, scale up x */
hx = __HI(x);
}
if (hx >= 0x7ff0_0000) {
return x + x;
}
k += (hx >> 20) - 1023;
i = (k & 0x8000_0000) >>> 31; // unsigned shift
hx = (hx & 0x000f_ffff) | ((0x3ff - i) << 20);
y = (double)(k + i);
x = __HI(x, hx); // replace high word of x with hx
z = y * log10_2lo + ivln10 * StrictMath.log(x);
return z + y * log10_2hi;
}
}
}

6
src/java.base/share/classes/java/lang/StrictMath.java
View File

@ -1,5 +1,5 @@
/*
* Copyright (c) 1999, 2022, Oracle and/or its affiliates. All rights reserved.
* Copyright (c) 1999, 2023, 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
@ -281,7 +281,9 @@ public final class StrictMath {
* @return the base 10 logarithm of {@code a}.
* @since 1.5
*/
public static native double log10(double a);
public static double log10(double a) {
return FdLibm.Log10.compute(a);
}
/**
* Returns the correctly rounded positive square root of a

82
test/jdk/java/lang/StrictMath/FdlibmTranslit.java
View File

@ -1,5 +1,5 @@
/*
* Copyright (c) 1998, 2016, Oracle and/or its affiliates. All rights reserved.
* Copyright (c) 1998, 2023, 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
@ -74,6 +74,14 @@ public class FdlibmTranslit {
return Hypot.compute(x, y);
}
public static double cbrt(double x) {
return Cbrt.compute(x);
}
public static double log10(double x) {
return Log10.compute(x);
}
/**
* cbrt(x)
* Return cube root of x
@ -89,7 +97,7 @@ public class FdlibmTranslit {
private static final double F = 1.60714285714285720630e+00; /* 45/28 = 0x3FF9B6DB, 0x6DB6DB6E */
private static final double G = 3.57142857142857150787e-01; /* 5/14 = 0x3FD6DB6D, 0xB6DB6DB7 */
public static strictfp double compute(double x) {
public static double compute(double x) {
int hx;
double r, s, t=0.0, w;
int sign; // unsigned
@ -333,7 +341,7 @@ public class FdlibmTranslit {
private static final double P4 = -1.65339022054652515390e-06; /* 0xBEBBBD41, 0xC5D26BF1 */
private static final double P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
public static strictfp double compute(double x) {
public static double compute(double x) {
double y,hi=0,lo=0,c,t;
int k=0,xsb;
/*unsigned*/ int hx;
@ -384,4 +392,72 @@ public class FdlibmTranslit {
}
}
}
/**
* Return the base 10 logarithm of x
*
* Method :
* Let log10_2hi = leading 40 bits of log10(2) and
* log10_2lo = log10(2) - log10_2hi,
* ivln10 = 1/log(10) rounded.
* Then
* n = ilogb(x),
* if(n<0) n = n+1;
* x = scalbn(x,-n);
* log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))
*
* Note 1:
* To guarantee log10(10**n)=n, where 10**n is normal, the rounding
* mode must set to Round-to-Nearest.
* Note 2:
* [1/log(10)] rounded to 53 bits has error .198 ulps;
* log10 is monotonic at all binary break points.
*
* Special cases:
* log10(x) is NaN with signal if x < 0;
* log10(+INF) is +INF with no signal; log10(0) is -INF with signal;
* log10(NaN) is that NaN with no signal;
* log10(10**N) = N for N=0,1,...,22.
*
* Constants:
* The hexadecimal values are the intended ones for the following constants.
* The decimal values may be used, provided that the compiler will convert
* from decimal to binary accurately enough to produce the hexadecimal values
* shown.
*/
static class Log10 {
private static double two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
private static double ivln10 = 4.34294481903251816668e-01; /* 0x3FDBCB7B, 0x1526E50E */
private static double log10_2hi = 3.01029995663611771306e-01; /* 0x3FD34413, 0x509F6000 */
private static double log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
private static double zero = 0.0;
public static double compute(double x) {
double y,z;
int i,k,hx;
/*unsigned*/ int lx;
hx = __HI(x); /* high word of x */
lx = __LO(x); /* low word of x */
k=0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx)==0)
return -two54/zero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 54; x *= two54; /* subnormal number, scale up x */
hx = __HI(x); /* high word of x */
}
if (hx >= 0x7ff00000) return x+x;
k += (hx>>20)-1023;
i = (k&0x80000000)>>>31; // unsigned shift
hx = (hx&0x000fffff)|((0x3ff-i)<<20);
y = (double)(k+i);
x = __HI(x, hx); //original: __HI(x) = hx;
z = y*log10_2lo + ivln10*StrictMath.log(x); // TOOD: switch to Translit.log when available
return z+y*log10_2hi;
}
}
}

70
test/jdk/java/lang/StrictMath/Log10Tests.java
View File

@ -1,5 +1,5 @@
/*
* Copyright (c) 2003, 2022, Oracle and/or its affiliates. All rights reserved.
* Copyright (c) 2003, 2023, 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
@ -23,10 +23,17 @@
/*
* @test
* @bug 4074599
* @bug 4074599 8301205
* @key randomness
* @library /test/lib
* @build jdk.test.lib.RandomFactory
* @build Tests
* @build FdlibmTranslit
* @build Log10Tests
* @run main Log10Tests
* @summary Tests for StrictMath.log10
*/
import jdk.test.lib.RandomFactory;
/**
* The tests in ../Math/Log10Tests.java test properties that should
@ -42,6 +49,19 @@
public class Log10Tests {
private Log10Tests(){}
public static void main(String... args) {
int failures = 0;
failures += testLog10();
failures += testAgainstTranslit();
if (failures > 0) {
System.err.println("Testing log10 incurred "
+ failures + " failures.");
throw new RuntimeException();
}
}
static int testLog10Case(double input, double expected) {
return Tests.test("StrictMath.log10(double)", input,
StrictMath::log10, expected);
@ -709,15 +729,47 @@ public class Log10Tests {
return failures;
}
public static void main(String [] argv) {
// Initialize shared random number generator
private static java.util.Random random = RandomFactory.getRandom();
/**
* Test StrictMath.log10 against transliteration port of log10.
*/
private static int testAgainstTranslit() {
int failures = 0;
double x;
failures += testLog10();
// Test just above subnormal threshold...
x = Double.MIN_NORMAL;
failures += testRange(x, Math.ulp(x), 1000);
if (failures > 0) {
System.err.println("Testing log10 incurred "
+ failures + " failures.");
throw new RuntimeException();
// ... and just below subnormal threshold ...
x = Math.nextDown(Double.MIN_NORMAL);
failures += testRange(x, -Math.ulp(x), 1000);
// ... and near 1.0 ...
x = 1.0;
failures += testRange(x, Math.ulp(x), 1000);
x = Math.nextDown(1.0);
failures += testRange(x, -Math.ulp(x), 1000);
x = Tests.createRandomDouble(random);
// Make the increment twice the ulp value in case the random
// value is near an exponent threshold. Don't worry about test
// elements overflowing to infinity if the starting value is
// near Double.MAX_VALUE.
failures += testRange(x, 2.0 * Math.ulp(x), 1000);
return failures;
}
private static int testRange(double start, double increment, int count) {
int failures = 0;
double x = start;
for (int i = 0; i < count; i++, x += increment) {
failures += testLog10Case(x, FdlibmTranslit.log10(x));
}
return failures;
}
}