From 63bb2ce8debdeb7f34bda9a3845fe93cb4d741dd Mon Sep 17 00:00:00 2001 From: Joe Darcy Date: Mon, 30 Jan 2023 20:33:01 +0000 Subject: [PATCH] 8301205: Port fdlibm log10 to Java Reviewed-by: bpb --- .../share/classes/java/lang/FdLibm.java | 79 +++++++++++++++++- .../share/classes/java/lang/StrictMath.java | 6 +- .../java/lang/StrictMath/FdlibmTranslit.java | 82 ++++++++++++++++++- test/jdk/java/lang/StrictMath/Log10Tests.java | 70 ++++++++++++++-- 4 files changed, 222 insertions(+), 15 deletions(-) diff --git a/src/java.base/share/classes/java/lang/FdLibm.java b/src/java.base/share/classes/java/lang/FdLibm.java index f5444c2b05f..7a0b5e624e1 100644 --- a/src/java.base/share/classes/java/lang/FdLibm.java +++ b/src/java.base/share/classes/java/lang/FdLibm.java @@ -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; + } + } } diff --git a/src/java.base/share/classes/java/lang/StrictMath.java b/src/java.base/share/classes/java/lang/StrictMath.java index d891069432e..7820a2e16cf 100644 --- a/src/java.base/share/classes/java/lang/StrictMath.java +++ b/src/java.base/share/classes/java/lang/StrictMath.java @@ -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 diff --git a/test/jdk/java/lang/StrictMath/FdlibmTranslit.java b/test/jdk/java/lang/StrictMath/FdlibmTranslit.java index 5ad48cb4b44..5c04d73fc50 100644 --- a/test/jdk/java/lang/StrictMath/FdlibmTranslit.java +++ b/test/jdk/java/lang/StrictMath/FdlibmTranslit.java @@ -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; + } + } } diff --git a/test/jdk/java/lang/StrictMath/Log10Tests.java b/test/jdk/java/lang/StrictMath/Log10Tests.java index eed6e8a5106..33aec723ec8 100644 --- a/test/jdk/java/lang/StrictMath/Log10Tests.java +++ b/test/jdk/java/lang/StrictMath/Log10Tests.java @@ -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; } }