diff --git a/src/java.base/share/classes/java/lang/FdLibm.java b/src/java.base/share/classes/java/lang/FdLibm.java index 94bda83df6c..460353d4156 100644 --- a/src/java.base/share/classes/java/lang/FdLibm.java +++ b/src/java.base/share/classes/java/lang/FdLibm.java @@ -747,6 +747,146 @@ class FdLibm { } } + /** + * Return the (natural) logarithm of x + * + * Method : + * 1. Argument Reduction: find k and f such that + * x = 2^k * (1+f), + * where sqrt(2)/2 < 1+f < sqrt(2) . + * + * 2. Approximation of log(1+f). + * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) + * = 2s + 2/3 s**3 + 2/5 s**5 + ....., + * = 2s + s*R + * We use a special Reme algorithm on [0,0.1716] to generate + * a polynomial of degree 14 to approximate R The maximum error + * of this polynomial approximation is bounded by 2**-58.45. In + * other words, + * 2 4 6 8 10 12 14 + * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s + * (the values of Lg1 to Lg7 are listed in the program) + * and + * | 2 14 | -58.45 + * | Lg1*s +...+Lg7*s - R(z) | <= 2 + * | | + * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. + * In order to guarantee error in log below 1ulp, we compute log + * by + * log(1+f) = f - s*(f - R) (if f is not too large) + * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) + * + * 3. Finally, log(x) = k*ln2 + log(1+f). + * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) + * Here ln2 is split into two floating point number: + * ln2_hi + ln2_lo, + * where n*ln2_hi is always exact for |n| < 2000. + * + * Special cases: + * log(x) is NaN with signal if x < 0 (including -INF) ; + * log(+INF) is +INF; log(0) is -INF with signal; + * log(NaN) is that NaN with no signal. + * + * Accuracy: + * according to an error analysis, the error is always less than + * 1 ulp (unit in the last place). + * + * 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 final class Log { + private Log() {throw new UnsupportedOperationException();} + + private static final double + ln2_hi = 0x1.62e42feep-1, // 6.93147180369123816490e-01 + ln2_lo = 0x1.a39ef35793c76p-33, // 1.90821492927058770002e-10 + + Lg1 = 0x1.5555555555593p-1, // 6.666666666666735130e-01 + Lg2 = 0x1.999999997fa04p-2, // 3.999999999940941908e-01 + Lg3 = 0x1.2492494229359p-2, // 2.857142874366239149e-01 + Lg4 = 0x1.c71c51d8e78afp-3, // 2.222219843214978396e-01 + Lg5 = 0x1.7466496cb03dep-3, // 1.818357216161805012e-01 + Lg6 = 0x1.39a09d078c69fp-3, // 1.531383769920937332e-01 + Lg7 = 0x1.2f112df3e5244p-3; // 1.479819860511658591e-01 + + private static final double zero = 0.0; + + static double compute(double x) { + double hfsq, f, s, z, R, w, t1, t2, dk; + int k, hx, i, j; + /*unsigned*/ int lx; + + hx = __HI(x); // high word of x + lx = __LO(x); // low word of x + + k=0; + if (hx < 0x0010_0000) { // x < 2**-1022 + if (((hx & 0x7fff_ffff) | lx) == 0) { // log(+-0) = -inf + return -TWO54/zero; + } + if (hx < 0) { // log(-#) = NaN + return (x - x)/zero; + } + k -= 54; + x *= TWO54; // subnormal number, scale up x + hx = __HI(x); // high word of x + } + if (hx >= 0x7ff0_0000) { + return x + x; + } + k += (hx >> 20) - 1023; + hx &= 0x000f_ffff; + i = (hx + 0x9_5f64) & 0x10_0000; + x =__HI(x, hx | (i ^ 0x3ff0_0000)); // normalize x or x/2 + k += (i >> 20); + f = x - 1.0; + if ((0x000f_ffff & (2 + hx)) < 3) {// |f| < 2**-20 + if (f == zero) { + if (k == 0) { + return zero; + } else { + dk = (double)k; + return dk*ln2_hi + dk*ln2_lo; + } + } + R = f*f*(0.5 - 0.33333333333333333*f); + if (k == 0) { + return f - R; + } else { + dk = (double)k; + return dk*ln2_hi - ((R - dk*ln2_lo) - f); + } + } + s = f/(2.0 + f); + dk = (double)k; + z = s*s; + i = hx - 0x6_147a; + w = z*z; + j = 0x6b851 - hx; + t1= w*(Lg2 + w*(Lg4 + w*Lg6)); + t2= z*(Lg1 + w*(Lg3 + w*(Lg5 + w*Lg7))); + i |= j; + R = t2 + t1; + if (i > 0) { + hfsq = 0.5*f*f; + if (k == 0) { + return f-(hfsq - s*(hfsq + R)); + } else { + return dk*ln2_hi - ((hfsq - (s*(hfsq + R) + dk*ln2_lo)) - f); + } + } else { + if (k == 0) { + return f - s*(f - R); + } else { + return dk*ln2_hi - ((s*(f - R) - dk*ln2_lo) - f); + } + } + } + } + /** * Return the base 10 logarithm of x * diff --git a/src/java.base/share/classes/java/lang/StrictMath.java b/src/java.base/share/classes/java/lang/StrictMath.java index 71b79a26891..6cec10e2f85 100644 --- a/src/java.base/share/classes/java/lang/StrictMath.java +++ b/src/java.base/share/classes/java/lang/StrictMath.java @@ -259,7 +259,9 @@ public final class StrictMath { * @return the value ln {@code a}, the natural logarithm of * {@code a}. */ - public static native double log(double a); + public static double log(double a) { + return FdLibm.Log.compute(a); + } /** * Returns the base 10 logarithm of a {@code double} value. diff --git a/test/jdk/java/lang/Math/LogTests.java b/test/jdk/java/lang/Math/LogTests.java new file mode 100644 index 00000000000..eaaa58d80ed --- /dev/null +++ b/test/jdk/java/lang/Math/LogTests.java @@ -0,0 +1,95 @@ +/* + * 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 + * under the terms of the GNU General Public License version 2 only, as + * published by the Free Software Foundation. + * + * 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. + */ + +/* + * @test + * @bug 8301202 + * @build Tests + * @build LogTests + * @run main LogTests + * @summary Tests for {Math, StrictMath}.log + */ + +public class LogTests { + private LogTests(){} + + public static void main(String... args) { + int failures = 0; + + failures += testLogSpecialCases(); + + if (failures > 0) { + System.err.println("Testing log incurred " + + failures + " failures."); + throw new RuntimeException(); + } + } + + private static final double infinityD = Double.POSITIVE_INFINITY; + private static final double NaNd = Double.NaN; + + /** + * From the spec for Math.log: + * "Special cases: + * + * If the argument is NaN or less than zero, then the result is NaN. + * If the argument is positive infinity, then the result is positive infinity. + * If the argument is positive zero or negative zero, then the result is negative infinity. + * If the argument is 1.0, then the result is positive zero. + */ + private static int testLogSpecialCases() { + int failures = 0; + + double [][] testCases = { + {Double.NaN, NaNd}, + {Double.NEGATIVE_INFINITY, NaNd}, + {-Double.MAX_VALUE, NaNd}, + {-1.0, NaNd}, + {-Double.MIN_NORMAL, NaNd}, + {-Double.MIN_VALUE, NaNd}, + + {Double.POSITIVE_INFINITY, infinityD}, + + {-0.0, -infinityD}, + {+0.0, -infinityD}, + + {+1.0, 0.0}, + }; + + for(int i = 0; i < testCases.length; i++) { + failures += testLogCase(testCases[i][0], + testCases[i][1]); + } + + return failures; + } + + private static int testLogCase(double input, double expected) { + int failures=0; + + failures+=Tests.test("Math.log", input, Math::log, expected); + failures+=Tests.test("StrictMath.log", input, StrictMath::log, expected); + + return failures; + } +} diff --git a/test/jdk/java/lang/StrictMath/FdlibmTranslit.java b/test/jdk/java/lang/StrictMath/FdlibmTranslit.java index 565bd69dc2f..d20dbb16bb1 100644 --- a/test/jdk/java/lang/StrictMath/FdlibmTranslit.java +++ b/test/jdk/java/lang/StrictMath/FdlibmTranslit.java @@ -78,6 +78,10 @@ public class FdlibmTranslit { return Cbrt.compute(x); } + public static double log(double x) { + return Log.compute(x); + } + public static double log10(double x) { return Log10.compute(x); } @@ -401,6 +405,125 @@ public class FdlibmTranslit { } } + /** + * Return the logarithm of x + * + * Method : + * 1. Argument Reduction: find k and f such that + * x = 2^k * (1+f), + * where sqrt(2)/2 < 1+f < sqrt(2) . + * + * 2. Approximation of log(1+f). + * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) + * = 2s + 2/3 s**3 + 2/5 s**5 + ....., + * = 2s + s*R + * We use a special Reme algorithm on [0,0.1716] to generate + * a polynomial of degree 14 to approximate R The maximum error + * of this polynomial approximation is bounded by 2**-58.45. In + * other words, + * 2 4 6 8 10 12 14 + * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s + * (the values of Lg1 to Lg7 are listed in the program) + * and + * | 2 14 | -58.45 + * | Lg1*s +...+Lg7*s - R(z) | <= 2 + * | | + * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. + * In order to guarantee error in log below 1ulp, we compute log + * by + * log(1+f) = f - s*(f - R) (if f is not too large) + * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) + * + * 3. Finally, log(x) = k*ln2 + log(1+f). + * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) + * Here ln2 is split into two floating point number: + * ln2_hi + ln2_lo, + * where n*ln2_hi is always exact for |n| < 2000. + * + * Special cases: + * log(x) is NaN with signal if x < 0 (including -INF) ; + * log(+INF) is +INF; log(0) is -INF with signal; + * log(NaN) is that NaN with no signal. + * + * Accuracy: + * according to an error analysis, the error is always less than + * 1 ulp (unit in the last place). + * + * 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. + */ + private static final class Log { + private static final double + ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ + ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ + two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ + Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ + Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ + Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ + Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ + Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ + Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ + Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ + + private static double zero = 0.0; + + static double compute(double x) { + double hfsq,f,s,z,R,w,t1,t2,dk; + int k,hx,i,j; + /*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; + hx &= 0x000fffff; + i = (hx+0x95f64)&0x100000; + // __HI(x) = hx|(i^0x3ff00000); /* normalize x or x/2 */ + x =__HI(x, hx|(i^0x3ff00000)); /* normalize x or x/2 */ + k += (i>>20); + f = x-1.0; + if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */ + if(f==zero) { + if (k==0) return zero; + else {dk=(double)k; return dk*ln2_hi+dk*ln2_lo;} + } + R = f*f*(0.5-0.33333333333333333*f); + if(k==0) return f-R; else {dk=(double)k; + return dk*ln2_hi-((R-dk*ln2_lo)-f);} + } + s = f/(2.0+f); + dk = (double)k; + z = s*s; + i = hx-0x6147a; + w = z*z; + j = 0x6b851-hx; + t1= w*(Lg2+w*(Lg4+w*Lg6)); + t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); + i |= j; + R = t2+t1; + if(i>0) { + hfsq=0.5*f*f; + if(k==0) return f-(hfsq-s*(hfsq+R)); else + return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f); + } else { + if(k==0) return f-s*(f-R); else + return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f); + } + } + } + /** * Return the base 10 logarithm of x * @@ -464,7 +587,7 @@ public class FdlibmTranslit { 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 + z = y*log10_2lo + ivln10*log(x); return z+y*log10_2hi; } } diff --git a/test/jdk/java/lang/StrictMath/LogTests.java b/test/jdk/java/lang/StrictMath/LogTests.java new file mode 100644 index 00000000000..a3e29b5cd08 --- /dev/null +++ b/test/jdk/java/lang/StrictMath/LogTests.java @@ -0,0 +1,136 @@ +/* + * Copyright (c) 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 + * under the terms of the GNU General Public License version 2 only, as + * published by the Free Software Foundation. + * + * 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. + */ + +/* + * @test + * @bug 8301202 + * @key randomness + * @library /test/lib + * @build jdk.test.lib.RandomFactory + * @build Tests + * @build FdlibmTranslit + * @build LogTests + * @run main LogTests + * @summary Tests for StrictMath.log + */ + +import jdk.test.lib.RandomFactory; + +public class LogTests { + private LogTests(){} + + public static void main(String... args) { + int failures = 0; + + failures += testLog(); + failures += testAgainstTranslit(); + + if (failures > 0) { + System.err.println("Testing log incurred " + + failures + " failures."); + throw new RuntimeException(); + } + } + + static int testLogCase(double input, double expected) { + return Tests.test("StrictMath.log(double)", input, + StrictMath::log, expected); + } + + // Inputs where Math.log and StrictMath.log differ for at least + // one Math.log implementation. + static int testLog() { + int failures = 0; + + double [][] testCases = { + {0x1.000000089cd6fp-43, -0x1.dce2a0697a102p4}, + {0x1.0000000830698p182, 0x1.f89c7428dd67ap6}, + {0x1.0000000744b3ap632, 0x1.b611ab2bd53cep8}, + {0x1.000000037d81fp766, 0x1.0979b1dbc4a42p9}, + {0x1.000000024028p991, 0x1.577455642bb92p9}, + }; + + for (double[] testCase: testCases) + failures+=testLogCase(testCase[0], testCase[1]); + + return failures; + } + + // Initialize shared random number generator + private static java.util.Random random = RandomFactory.getRandom(); + + /** + * Test StrictMath.log against transliteration port of log. + */ + private static int testAgainstTranslit() { + int failures = 0; + double x; + + // Test just above subnormal threshold... + x = Double.MIN_NORMAL; + failures += testRange(x, Math.ulp(x), 1000); + + // ... and just below subnormal threshold ... + x = Math.nextDown(Double.MIN_NORMAL); + failures += testRange(x, -Math.ulp(x), 1000); + + // Probe near decision points in the FDLIBM algorithm. + double[] decisionPoints = { + 0x1.0p-1022, + + 0x1.0p-20, + }; + + for (double testPoint : decisionPoints) { + failures += testRangeMidpoint(testPoint, Math.ulp(testPoint), 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 += testLogCase(x, FdlibmTranslit.log(x)); + } + return failures; + } + + private static int testRangeMidpoint(double midpoint, double increment, int count) { + int failures = 0; + double x = midpoint - increment*(count / 2) ; + for (int i = 0; i < count; i++, x += increment) { + failures += testLogCase(x, FdlibmTranslit.log(x)); + } + return failures; + } +}