diff --git a/src/java.base/share/classes/java/lang/FdLibm.java b/src/java.base/share/classes/java/lang/FdLibm.java index 460353d4156..9526e7b695f 100644 --- a/src/java.base/share/classes/java/lang/FdLibm.java +++ b/src/java.base/share/classes/java/lang/FdLibm.java @@ -60,7 +60,9 @@ package java.lang; class FdLibm { // Constants used by multiple algorithms private static final double INFINITY = Double.POSITIVE_INFINITY; - private static final double TWO54 = 0x1.0p54; // 1.80143985094819840000e+16 + private static final double TWO54 = 0x1.0p54; // 1.80143985094819840000e+16 + private static final double HUGE = 1.0e+300; + private FdLibm() { throw new UnsupportedOperationException("No FdLibm instances for you."); @@ -102,6 +104,296 @@ class FdLibm { ( ((long)high)) << 32 ); } + /** Returns the arcsine of x. + * + * Method : + * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ... + * we approximate asin(x) on [0,0.5] by + * asin(x) = x + x*x^2*R(x^2) + * where + * R(x^2) is a rational approximation of (asin(x)-x)/x^3 + * and its remez error is bounded by + * |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75) + * + * For x in [0.5,1] + * asin(x) = pi/2-2*asin(sqrt((1-x)/2)) + * Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2; + * then for x>0.98 + * asin(x) = pi/2 - 2*(s+s*z*R(z)) + * = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo) + * For x<=0.98, let pio4_hi = pio2_hi/2, then + * f = hi part of s; + * c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z) + * and + * asin(x) = pi/2 - 2*(s+s*z*R(z)) + * = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo) + * = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c)) + * + * Special cases: + * if x is NaN, return x itself; + * if |x|>1, return NaN with invalid signal. + * + */ + static class Asin { + private Asin() {throw new UnsupportedOperationException();} + + private static final double + pio2_hi = 0x1.921fb54442d18p0, // 1.57079632679489655800e+00 + pio2_lo = 0x1.1a62633145c07p-54, // 6.12323399573676603587e-17 + pio4_hi = 0x1.921fb54442d18p-1, // 7.85398163397448278999e-01 + // coefficient for R(x^2) + pS0 = 0x1.5555555555555p-3, // 1.66666666666666657415e-01 + pS1 = -0x1.4d61203eb6f7dp-2, // -3.25565818622400915405e-01 + pS2 = 0x1.9c1550e884455p-3, // 2.01212532134862925881e-01 + pS3 = -0x1.48228b5688f3bp-5, // -4.00555345006794114027e-02 + pS4 = 0x1.9efe07501b288p-11, // 7.91534994289814532176e-04 + pS5 = 0x1.23de10dfdf709p-15, // 3.47933107596021167570e-05 + qS1 = -0x1.33a271c8a2d4bp1, // -2.40339491173441421878e+00 + qS2 = 0x1.02ae59c598ac8p1, // 2.02094576023350569471e+00 + qS3 = -0x1.6066c1b8d0159p-1, // -6.88283971605453293030e-01 + qS4 = 0x1.3b8c5b12e9282p-4; // 7.70381505559019352791e-02 + + static double compute(double x) { + double t = 0, w, p, q, c, r, s; + int hx, ix; + hx = __HI(x); + ix = hx & 0x7fff_ffff; + if (ix >= 0x3ff0_0000) { // |x| >= 1 + if(((ix - 0x3ff0_0000) | __LO(x)) == 0) { + // asin(1) = +-pi/2 with inexact + return x*pio2_hi + x*pio2_lo; + } + return (x - x)/(x - x); // asin(|x| > 1) is NaN + } else if (ix < 0x3fe0_0000) { // |x| < 0.5 + if (ix < 0x3e40_0000) { // if |x| < 2**-27 + if (HUGE + x > 1.0) {// return x with inexact if x != 0 + return x; + } + } else { + t = x*x; + } + p = t*(pS0 + t*(pS1 + t*(pS2 + t*(pS3 + t*(pS4 + t*pS5))))); + q = 1.0 + t*(qS1 + t*(qS2 + t*(qS3 + t*qS4))); + w = p/q; + return x + x*w; + } + // 1 > |x| >= 0.5 + w = 1.0 - Math.abs(x); + t = w*0.5; + p = t*(pS0 + t*(pS1 + t*(pS2 + t*(pS3 + t*(pS4 + t*pS5))))); + q = 1.0 + t*(qS1 + t*(qS2 + t*(qS3 + t*qS4))); + s = Math.sqrt(t); + if (ix >= 0x3FEF_3333) { // if |x| > 0.975 + w = p/q; + t = pio2_hi - (2.0*(s + s*w) - pio2_lo); + } else { + w = s; + w = __LO(w, 0); + c = (t - w*w)/(s + w); + r = p/q; + p = 2.0*s*r - (pio2_lo - 2.0*c); + q = pio4_hi - 2.0*w; + t = pio4_hi - (p - q); + } + return (hx > 0) ? t : -t; + } + } + + /** Returns the arccosine of x. + * Method : + * acos(x) = pi/2 - asin(x) + * acos(-x) = pi/2 + asin(x) + * For |x| <= 0.5 + * acos(x) = pi/2 - (x + x*x^2*R(x^2)) (see asin.c) + * For x > 0.5 + * acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2))) + * = 2asin(sqrt((1-x)/2)) + * = 2s + 2s*z*R(z) ...z=(1-x)/2, s=sqrt(z) + * = 2f + (2c + 2s*z*R(z)) + * where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term + * for f so that f+c ~ sqrt(z). + * For x <- 0.5 + * acos(x) = pi - 2asin(sqrt((1-|x|)/2)) + * = pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z) + * + * Special cases: + * if x is NaN, return x itself; + * if |x|>1, return NaN with invalid signal. + * + * Function needed: sqrt + */ + static class Acos { + private Acos() {throw new UnsupportedOperationException();} + + private static final double + pio2_hi = 0x1.921fb54442d18p0, // 1.57079632679489655800e+00 + pio2_lo = 0x1.1a62633145c07p-54, // 6.12323399573676603587e-17 + pS0 = 0x1.5555555555555p-3, // 1.66666666666666657415e-01 + pS1 = -0x1.4d61203eb6f7dp-2, // -3.25565818622400915405e-01 + pS2 = 0x1.9c1550e884455p-3, // 2.01212532134862925881e-01 + pS3 = -0x1.48228b5688f3bp-5, // -4.00555345006794114027e-02 + pS4 = 0x1.9efe07501b288p-11, // 7.91534994289814532176e-04 + pS5 = 0x1.23de10dfdf709p-15, // 3.47933107596021167570e-05 + qS1 = -0x1.33a271c8a2d4bp1, // -2.40339491173441421878e+00 + qS2 = 0x1.02ae59c598ac8p1, // 2.02094576023350569471e+00 + qS3 = -0x1.6066c1b8d0159p-1, // -6.88283971605453293030e-01 + qS4 = 0x1.3b8c5b12e9282p-4; // 7.70381505559019352791e-02 + + static double compute(double x) { + double z, p, q, r, w, s, c, df; + int hx, ix; + hx = __HI(x); + ix = hx & 0x7fff_ffff; + if (ix >= 0x3ff0_0000) { // |x| >= 1 + if (((ix - 0x3ff0_0000) | __LO(x)) == 0) { // |x| == 1 + if (hx > 0) {// acos(1) = 0 + return 0.0; + }else { // acos(-1)= pi + return Math.PI + 2.0*pio2_lo; + } + } + return (x-x)/(x-x); // acos(|x| > 1) is NaN + } + if (ix < 0x3fe0_0000) { // |x| < 0.5 + if (ix <= 0x3c60_0000) { // if |x| < 2**-57 + return pio2_hi + pio2_lo; + } + z = x*x; + p = z*(pS0 + z*(pS1 + z*(pS2 + z*(pS3 + z*(pS4 + z*pS5))))); + q = 1.0 + z*(qS1 + z*(qS2 + z*(qS3 + z*qS4))); + r = p/q; + return pio2_hi - (x - (pio2_lo - x*r)); + } else if (hx < 0) { // x < -0.5 + z = (1.0 + x)*0.5; + p = z*(pS0 + z*(pS1 + z*(pS2 + z*(pS3 + z*(pS4 + z*pS5))))); + q = 1.0 + z*(qS1 + z*(qS2 + z*(qS3 + z*qS4))); + s = Math.sqrt(z); + r = p/q; + w = r*s - pio2_lo; + return Math.PI - 2.0*(s+w); + } else { // x > 0.5 + z = (1.0 - x)*0.5; + s = Math.sqrt(z); + df = s; + df = __LO(df, 0); + c = (z - df*df)/(s + df); + p = z*(pS0 + z*(pS1 + z*(pS2 + z*(pS3 + z*(pS4 + z*pS5))))); + q = 1.0 + z*(qS1 + z*(qS2 + z*(qS3 + z*qS4))); + r = p/q; + w = r*s + c; + return 2.0*(df + w); + } + } + } + + /* Returns the arctangent of x. + * Method + * 1. Reduce x to positive by atan(x) = -atan(-x). + * 2. According to the integer k=4t+0.25 chopped, t=x, the argument + * is further reduced to one of the following intervals and the + * arctangent of t is evaluated by the corresponding formula: + * + * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) + * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) ) + * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) ) + * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) ) + * [39/16,INF] atan(x) = atan(INF) + atan( -1/t ) + * + * 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 Atan { + private Atan() {throw new UnsupportedOperationException();} + + private static final double atanhi[] = { + 0x1.dac670561bb4fp-2, // atan(0.5)hi 4.63647609000806093515e-01 + 0x1.921fb54442d18p-1, // atan(1.0)hi 7.85398163397448278999e-01 + 0x1.f730bd281f69bp-1, // atan(1.5)hi 9.82793723247329054082e-01 + 0x1.921fb54442d18p0, // atan(inf)hi 1.57079632679489655800e+00 + }; + + private static final double atanlo[] = { + 0x1.a2b7f222f65e2p-56, // atan(0.5)lo 2.26987774529616870924e-17 + 0x1.1a62633145c07p-55, // atan(1.0)lo 3.06161699786838301793e-17 + 0x1.007887af0cbbdp-56, // atan(1.5)lo 1.39033110312309984516e-17 + 0x1.1a62633145c07p-54, // atan(inf)lo 6.12323399573676603587e-17 + }; + + private static final double aT[] = { + 0x1.555555555550dp-2, // 3.33333333333329318027e-01 + -0x1.999999998ebc4p-3, // -1.99999999998764832476e-01 + 0x1.24924920083ffp-3, // 1.42857142725034663711e-01 + -0x1.c71c6fe231671p-4, // -1.11111104054623557880e-01 + 0x1.745cdc54c206ep-4, // 9.09088713343650656196e-02 + -0x1.3b0f2af749a6dp-4, // -7.69187620504482999495e-02 + 0x1.10d66a0d03d51p-4, // 6.66107313738753120669e-02 + -0x1.dde2d52defd9ap-5, // -5.83357013379057348645e-02 + 0x1.97b4b24760debp-5, // 4.97687799461593236017e-02 + -0x1.2b4442c6a6c2fp-5, // -3.65315727442169155270e-02 + 0x1.0ad3ae322da11p-6, // 1.62858201153657823623e-02 + }; + + static double compute(double x) { + double w, s1, s2, z; + int ix, hx, id; + + hx = __HI(x); + ix = hx & 0x7fff_ffff; + if (ix >= 0x4410_0000) { // if |x| >= 2^66 + if (ix > 0x7ff0_0000 || + (ix == 0x7ff0_0000 && (__LO(x) != 0))) { + return x+x; // NaN + } + if (hx > 0) { + return atanhi[3] + atanlo[3]; + } else { + return -atanhi[3] - atanlo[3]; + } + } if (ix < 0x3fdc_0000) { // |x| < 0.4375 + if (ix < 0x3e20_0000) { // |x| < 2^-29 + if (HUGE + x > 1.0) { // raise inexact + return x; + } + } + id = -1; + } else { + x = Math.abs(x); + if (ix < 0x3ff3_0000) { // |x| < 1.1875 + if (ix < 0x3fe60000) { // 7/16 <= |x| < 11/16 + id = 0; + x = (2.0*x - 1.0)/(2.0 + x); + } else { // 11/16 <= |x| < 19/16 + id = 1; + x = (x - 1.0)/(x + 1.0); + } + } else { + if (ix < 0x4003_8000) { // |x| < 2.4375 + id = 2; + x = (x - 1.5)/(1.0 + 1.5*x); + } else { // 2.4375 <= |x| < 2^66 + id = 3; + x = -1.0/x; + } + } + } + // end of argument reduction + z = x*x; + w = z*z; + // break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly + s1 = z*(aT[0] + w*(aT[2] + w*(aT[4] + w*(aT[6] + w*(aT[8] + w*aT[10]))))); + s2 = w*(aT[1] + w*(aT[3] + w*(aT[5] + w*(aT[7] + w*aT[9])))); + if (id < 0) { + return x - x*(s1 + s2); + } else { + z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x); + return (hx < 0) ? -z: z; + } + } + } + /** * cbrt(x) * Return cube root 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 f08b016b629..a190c92539f 100644 --- a/src/java.base/share/classes/java/lang/StrictMath.java +++ b/src/java.base/share/classes/java/lang/StrictMath.java @@ -164,7 +164,9 @@ public final class StrictMath { * @param a the value whose arc sine is to be returned. * @return the arc sine of the argument. */ - public static native double asin(double a); + public static double asin(double a) { + return FdLibm.Asin.compute(a); + } /** * Returns the arc cosine of a value; the returned angle is in the @@ -177,7 +179,9 @@ public final class StrictMath { * @param a the value whose arc cosine is to be returned. * @return the arc cosine of the argument. */ - public static native double acos(double a); + public static double acos(double a) { + return FdLibm.Acos.compute(a); + } /** * Returns the arc tangent of a value; the returned angle is in the @@ -193,7 +197,9 @@ public final class StrictMath { * @param a the value whose arc tangent is to be returned. * @return the arc tangent of the argument. */ - public static native double atan(double a); + public static double atan(double a) { + return FdLibm.Atan.compute(a); + } /** * Converts an angle measured in degrees to an approximately diff --git a/test/jdk/java/lang/Math/InverseTrigTests.java b/test/jdk/java/lang/Math/InverseTrigTests.java new file mode 100644 index 00000000000..e779ce20c7e --- /dev/null +++ b/test/jdk/java/lang/Math/InverseTrigTests.java @@ -0,0 +1,177 @@ +/* + * 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 8302026 + * @build Tests + * @build InverseTrigTests + * @run main InverseTrigTests + * @summary Tests for {Math, StrictMath}.{asin, acos, atan} + */ + +public class InverseTrigTests { + private InverseTrigTests(){} + + public static void main(String... args) { + int failures = 0; + + failures += testAsinSpecialCases(); + failures += testAcosSpecialCases(); + failures += testAtanSpecialCases(); + + if (failures > 0) { + System.err.println("Testing inverse trig mthods 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.asin: + * + * "Special cases: + * + * If the argument is NaN or its absolute value is greater than 1, + * then the result is NaN. + * + * If the argument is zero, then the result is a zero with the + * same sign as the argument." + */ + private static int testAsinSpecialCases() { + int failures = 0; + + double [][] testCases = { + {NaNd, NaNd}, + {Math.nextUp(1.0), NaNd}, + {Math.nextDown(-1.0), NaNd}, + { InfinityD, NaNd}, + {-InfinityD, NaNd}, + + {-0.0, -0.0}, + {+0.0, +0.0}, + }; + + for(int i = 0; i < testCases.length; i++) { + failures += testAsinCase(testCases[i][0], + testCases[i][1]); + } + + return failures; + } + + private static int testAsinCase(double input, double expected) { + int failures=0; + + failures+=Tests.test("Math.asin", input, Math::asin, expected); + failures+=Tests.test("StrictMath.asin", input, StrictMath::asin, expected); + + return failures; + } + + /** + * From the spec for Math.acos: + * + * "Special case: + * + * If the argument is NaN or its absolute value is greater than 1, + * then the result is NaN. + * + * If the argument is 1.0, the result is positive zero." + */ + private static int testAcosSpecialCases() { + int failures = 0; + + double [][] testCases = { + {NaNd, NaNd}, + {Math.nextUp(1.0), NaNd}, + {Math.nextDown(-1.0), NaNd}, + {InfinityD, NaNd}, + {-InfinityD, NaNd}, + + {1.0, +0.0}, + }; + + for(int i = 0; i < testCases.length; i++) { + failures += testAcosCase(testCases[i][0], + testCases[i][1]); + } + + return failures; + } + + private static int testAcosCase(double input, double expected) { + int failures=0; + + failures+=Tests.test("Math.acos", input, Math::acos, expected); + failures+=Tests.test("StrictMath.acos", input, StrictMath::acos, expected); + + return failures; + } + + /** + * From the spec for Math.atan: + * + * "Special cases: + * + * If the argument is NaN, then the result is NaN. + * + * If the argument is zero, then the result is a zero with the + * same sign as the argument. + * + * If the argument is infinite, then the result is the closest + * value to pi/2 with the same sign as the input." + */ + private static int testAtanSpecialCases() { + int failures = 0; + + double [][] testCases = { + {NaNd, NaNd}, + + {-0.0, -0.0}, + {+0.0, +0.0}, + + { InfinityD, +Math.PI/2.0}, + {-InfinityD, -Math.PI/2.0}, + }; + + for(int i = 0; i < testCases.length; i++) { + failures += testAtanCase(testCases[i][0], + testCases[i][1]); + } + + return failures; + } + + private static int testAtanCase(double input, double expected) { + int failures=0; + + failures+=Tests.test("Math.atan", input, Math::atan, expected); + failures+=Tests.test("StrictMath.atan", input, StrictMath::atan, expected); + + return failures; + } +} diff --git a/test/jdk/java/lang/StrictMath/ExhaustingTests.java b/test/jdk/java/lang/StrictMath/ExhaustingTests.java index 09e182d4e16..457e91d2c99 100644 --- a/test/jdk/java/lang/StrictMath/ExhaustingTests.java +++ b/test/jdk/java/lang/StrictMath/ExhaustingTests.java @@ -23,7 +23,7 @@ /* * @test - * @bug 8301833 + * @bug 8301833 8302026 * @build Tests * @build FdlibmTranslit * @build ExhaustingTests @@ -84,9 +84,9 @@ public class ExhaustingTests { // new UnaryTestCase("cos", FdlibmTranslit::cos, StrictMath::cos, DEFAULT_SHIFT), // new UnaryTestCase("tan", FdlibmTranslit::tan, StrictMath::tan, DEFAULT_SHIFT), - // new UnaryTestCase("asin", FdlibmTranslit::asin, StrictMath::asin, DEFAULT_SHIFT), - // new UnaryTestCase("acos", FdlibmTranslit::acos, StrictMath::acos, DEFAULT_SHIFT), - // new UnaryTestCase("atan", FdlibmTranslit::atan, StrictMath::atan, DEFAULT_SHIFT), + new UnaryTestCase("asin", FdlibmTranslit::asin, StrictMath::asin, DEFAULT_SHIFT), + new UnaryTestCase("acos", FdlibmTranslit::acos, StrictMath::acos, DEFAULT_SHIFT), + new UnaryTestCase("atan", FdlibmTranslit::atan, StrictMath::atan, DEFAULT_SHIFT), }; for (var testCase : testCases) { diff --git a/test/jdk/java/lang/StrictMath/FdlibmTranslit.java b/test/jdk/java/lang/StrictMath/FdlibmTranslit.java index d20dbb16bb1..4dcd98ad4d6 100644 --- a/test/jdk/java/lang/StrictMath/FdlibmTranslit.java +++ b/test/jdk/java/lang/StrictMath/FdlibmTranslit.java @@ -70,6 +70,18 @@ public class FdlibmTranslit { ( ((long)high)) << 32 ); } + public static double asin(double x) { + return Asin.compute(x); + } + + public static double acos(double x) { + return Acos.compute(x); + } + + public static double atan(double x) { + return Atan.compute(x); + } + public static double hypot(double x, double y) { return Hypot.compute(x, y); } @@ -94,6 +106,279 @@ public class FdlibmTranslit { return Expm1.compute(x); } + /** Returns the arcsine of x. + * + * Method : + * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ... + * we approximate asin(x) on [0,0.5] by + * asin(x) = x + x*x^2*R(x^2) + * where + * R(x^2) is a rational approximation of (asin(x)-x)/x^3 + * and its remez error is bounded by + * |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75) + * + * For x in [0.5,1] + * asin(x) = pi/2-2*asin(sqrt((1-x)/2)) + * Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2; + * then for x>0.98 + * asin(x) = pi/2 - 2*(s+s*z*R(z)) + * = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo) + * For x<=0.98, let pio4_hi = pio2_hi/2, then + * f = hi part of s; + * c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z) + * and + * asin(x) = pi/2 - 2*(s+s*z*R(z)) + * = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo) + * = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c)) + * + * Special cases: + * if x is NaN, return x itself; + * if |x|>1, return NaN with invalid signal. + * + */ + static class Asin { + private static final double + one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ + huge = 1.000e+300, + pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ + pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */ + pio4_hi = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */ + /* coefficient for R(x^2) */ + pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */ + pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */ + pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */ + pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */ + pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */ + pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */ + qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */ + qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */ + qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ + qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ + + static double compute(double x) { + double t=0,w,p,q,c,r,s; + int hx,ix; + hx = __HI(x); + ix = hx&0x7fffffff; + if(ix>= 0x3ff00000) { /* |x|>= 1 */ + if(((ix-0x3ff00000)|__LO(x))==0) + /* asin(1)=+-pi/2 with inexact */ + return x*pio2_hi+x*pio2_lo; + return (x-x)/(x-x); /* asin(|x|>1) is NaN */ + } else if (ix<0x3fe00000) { /* |x|<0.5 */ + if(ix<0x3e400000) { /* if |x| < 2**-27 */ + if(huge+x>one) return x;/* return x with inexact if x!=0*/ + } else + t = x*x; + p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); + q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); + w = p/q; + return x+x*w; + } + /* 1> |x|>= 0.5 */ + w = one-Math.abs(x); + t = w*0.5; + p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); + q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); + s = Math.sqrt(t); + if(ix>=0x3FEF3333) { /* if |x| > 0.975 */ + w = p/q; + t = pio2_hi-(2.0*(s+s*w)-pio2_lo); + } else { + w = s; + // __LO(w) = 0; + w = __LO(w, 0); + c = (t-w*w)/(s+w); + r = p/q; + p = 2.0*s*r-(pio2_lo-2.0*c); + q = pio4_hi-2.0*w; + t = pio4_hi-(p-q); + } + if(hx>0) return t; else return -t; + } + } + + /** Returns the arccosine of x. + * Method : + * acos(x) = pi/2 - asin(x) + * acos(-x) = pi/2 + asin(x) + * For |x|<=0.5 + * acos(x) = pi/2 - (x + x*x^2*R(x^2)) (see asin.c) + * For x>0.5 + * acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2))) + * = 2asin(sqrt((1-x)/2)) + * = 2s + 2s*z*R(z) ...z=(1-x)/2, s=sqrt(z) + * = 2f + (2c + 2s*z*R(z)) + * where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term + * for f so that f+c ~ sqrt(z). + * For x<-0.5 + * acos(x) = pi - 2asin(sqrt((1-|x|)/2)) + * = pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z) + * + * Special cases: + * if x is NaN, return x itself; + * if |x|>1, return NaN with invalid signal. + * + * Function needed: sqrt + */ + static class Acos { + private static final double + one= 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ + pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ + pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ + pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */ + pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */ + pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */ + pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */ + pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */ + pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */ + pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */ + qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */ + qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */ + qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ + qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ + + static double compute(double x) { + double z,p,q,r,w,s,c,df; + int hx,ix; + hx = __HI(x); + ix = hx&0x7fffffff; + if(ix>=0x3ff00000) { /* |x| >= 1 */ + if(((ix-0x3ff00000)|__LO(x))==0) { /* |x|==1 */ + if(hx>0) return 0.0; /* acos(1) = 0 */ + else return pi+2.0*pio2_lo; /* acos(-1)= pi */ + } + return (x-x)/(x-x); /* acos(|x|>1) is NaN */ + } + if(ix<0x3fe00000) { /* |x| < 0.5 */ + if(ix<=0x3c600000) return pio2_hi+pio2_lo;/*if|x|<2**-57*/ + z = x*x; + p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); + q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); + r = p/q; + return pio2_hi - (x - (pio2_lo-x*r)); + } else if (hx<0) { /* x < -0.5 */ + z = (one+x)*0.5; + p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); + q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); + s = Math.sqrt(z); + r = p/q; + w = r*s-pio2_lo; + return pi - 2.0*(s+w); + } else { /* x > 0.5 */ + z = (one-x)*0.5; + s = Math.sqrt(z); + df = s; + // __LO(df) = 0; + df = __LO(df, 0); + c = (z-df*df)/(s+df); + p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); + q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); + r = p/q; + w = r*s+c; + return 2.0*(df+w); + } + } + } + + /* Returns the arctangent of x. + * Method + * 1. Reduce x to positive by atan(x) = -atan(-x). + * 2. According to the integer k=4t+0.25 chopped, t=x, the argument + * is further reduced to one of the following intervals and the + * arctangent of t is evaluated by the corresponding formula: + * + * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) + * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) ) + * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) ) + * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) ) + * [39/16,INF] atan(x) = atan(INF) + atan( -1/t ) + * + * 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 Atan { + private static final double atanhi[] = { + 4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */ + 7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */ + 9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */ + 1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */ + }; + + private static final double atanlo[] = { + 2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */ + 3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */ + 1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */ + 6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */ + }; + + private static final double aT[] = { + 3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */ + -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */ + 1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */ + -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */ + 9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */ + -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */ + 6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */ + -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */ + 4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */ + -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */ + 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */ + }; + + private static final double + one = 1.0, + huge = 1.0e300; + + static double compute(double x) { + double w,s1,s2,z; + int ix,hx,id; + + hx = __HI(x); + ix = hx&0x7fffffff; + if(ix>=0x44100000) { /* if |x| >= 2^66 */ + if(ix>0x7ff00000|| + (ix==0x7ff00000&&(__LO(x)!=0))) + return x+x; /* NaN */ + if(hx>0) return atanhi[3]+atanlo[3]; + else return -atanhi[3]-atanlo[3]; + } if (ix < 0x3fdc0000) { /* |x| < 0.4375 */ + if (ix < 0x3e200000) { /* |x| < 2^-29 */ + if(huge+x>one) return x; /* raise inexact */ + } + id = -1; + } else { + x = Math.abs(x); + if (ix < 0x3ff30000) { /* |x| < 1.1875 */ + if (ix < 0x3fe60000) { /* 7/16 <=|x|<11/16 */ + id = 0; x = (2.0*x-one)/(2.0+x); + } else { /* 11/16<=|x|< 19/16 */ + id = 1; x = (x-one)/(x+one); + } + } else { + if (ix < 0x40038000) { /* |x| < 2.4375 */ + id = 2; x = (x-1.5)/(one+1.5*x); + } else { /* 2.4375 <= |x| < 2^66 */ + id = 3; x = -1.0/x; + } + }} + /* end of argument reduction */ + z = x*x; + w = z*z; + /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */ + s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10]))))); + s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9])))); + if (id<0) return x - x*(s1+s2); + else { + z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x); + return (hx<0)? -z:z; + } + } + } + /** * cbrt(x) * Return cube root of x diff --git a/test/jdk/java/lang/StrictMath/InverseTrigTests.java b/test/jdk/java/lang/StrictMath/InverseTrigTests.java new file mode 100644 index 00000000000..ca5bc59d2c9 --- /dev/null +++ b/test/jdk/java/lang/StrictMath/InverseTrigTests.java @@ -0,0 +1,226 @@ +/* + * 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. + */ +import jdk.test.lib.RandomFactory; +import java.util.function.DoubleUnaryOperator; + +/* + * @test + * @bug 8302026 + * @key randomness + * @library /test/lib + * @build jdk.test.lib.RandomFactory + * @build Tests + * @build FdlibmTranslit + * @build InverseTrigTests + * @run main InverseTrigTests + * @summary Tests for StrictMath.{asin, acos, atan} + */ + +/** + * The tests in ../Math/InverseTrigTests.java test properties that + * should hold for any implementation of the inverse trig functions + * ason, acos, and atan, including the FDLIBM-based ones required by + * the StrictMath class. Therefore, the test cases in + * ../Math/InverseTrig.java are run against both the Math and + * StrictMath versions of the inverse trig methods. The role of this + * test is to verify that the FDLIBM algorithms are being used by + * running golden file tests on values that may vary from one + * conforming implementation of the hyperbolics to another. + */ + +public class InverseTrigTests { + private InverseTrigTests(){} + + public static void main(String... args) { + int failures = 0; + + failures += testAgainstTranslitCommon(); + + failures += testAgainstTranslitAsin(); + failures += testAgainstTranslitAcos(); + failures += testAgainstTranslitAtan(); + + if (failures > 0) { + System.err.println("Testing the inverse trig functions incurred " + + failures + " failures."); + throw new RuntimeException(); + } + } + + /** + * Bundle together groups of testing methods. + */ + private static enum InverseTrigTest { + ASIN(InverseTrigTests::testAsinCase, FdlibmTranslit::asin), + ACOS(InverseTrigTests::testAcosCase, FdlibmTranslit::acos), + ATAN(InverseTrigTests::testAtanCase, FdlibmTranslit::atan); + + private DoubleDoubleToInt testCase; + private DoubleUnaryOperator transliteration; + + InverseTrigTest(DoubleDoubleToInt testCase, DoubleUnaryOperator transliteration) { + this.testCase = testCase; + this.transliteration = transliteration; + } + + public DoubleDoubleToInt testCase() {return testCase;} + public DoubleUnaryOperator transliteration() {return transliteration;} + } + + // Initialize shared random number generator + private static java.util.Random random = RandomFactory.getRandom(); + + /** + * Test against shared points of interest. + */ + private static int testAgainstTranslitCommon() { + int failures = 0; + double[] pointsOfInterest = { + Double.MIN_NORMAL, + 1.0, + Tests.createRandomDouble(random), + }; + + for (var testMethods : InverseTrigTest.values()) { + for (double testPoint : pointsOfInterest) { + failures += testRangeMidpoint(testPoint, Math.ulp(testPoint), 1000, testMethods); + } + } + + return failures; + } + + /** + * Test StrictMath.asin against transliteration port of asin. + */ + private static int testAgainstTranslitAsin() { + int failures = 0; + + // Probe near decision points in the FDLIBM algorithm. + double[] decisionPoints = { + 0x1p-27, + -0x1p-27, + + 0.5, + -0.5, + + 0.975, + -0.975, + }; + + for (double testPoint : decisionPoints) { + failures += testRangeMidpoint(testPoint, Math.ulp(testPoint), 1000, InverseTrigTest.ASIN); + } + + return failures; + } + + /** + * Test StrictMath.acos against transliteration port of acos. + */ + private static int testAgainstTranslitAcos() { + int failures = 0; + + // Probe near decision points in the FDLIBM algorithm. + double[] decisionPoints = { + 0.5, + -0.5, + + 0x1.0p-57, + -0x1.0p-57, + }; + + for (double testPoint : decisionPoints) { + failures += testRangeMidpoint(testPoint, Math.ulp(testPoint), 1000, InverseTrigTest.ACOS); + } + + return failures; + } + + /** + * Test StrictMath.atan against transliteration port of atan + */ + private static int testAgainstTranslitAtan() { + int failures = 0; + + // Probe near decision points in the FDLIBM algorithm. + double[] decisionPoints = { + 0.0, + + 7.0/16.0, + 11.0/16.0, + 19.0/16.0, + 39.0/16.0, + + 0x1.0p66, + 0x1.0p-29, + }; + + for (double testPoint : decisionPoints) { + failures += testRangeMidpoint(testPoint, Math.ulp(testPoint), 1000, InverseTrigTest.ATAN); + } + + return failures; + } + + private interface DoubleDoubleToInt { + int apply(double x, double y); + } + + private static int testRange(double start, double increment, int count, + InverseTrigTest testMethods) { + int failures = 0; + double x = start; + for (int i = 0; i < count; i++, x += increment) { + failures += + testMethods.testCase().apply(x, testMethods.transliteration().applyAsDouble(x)); + } + return failures; + } + + private static int testRangeMidpoint(double midpoint, double increment, int count, + InverseTrigTest testMethods) { + int failures = 0; + double x = midpoint - increment*(count / 2) ; + for (int i = 0; i < count; i++, x += increment) { + failures += + testMethods.testCase().apply(x, testMethods.transliteration().applyAsDouble(x)); + } + return failures; + } + + private static int testAsinCase(double input, double expected) { + return Tests.test("StrictMath.asin(double)", input, + StrictMath::asin, expected); + } + + private static int testAcosCase(double input, double expected) { + return Tests.test("StrictMath.acos(double)", input, + StrictMath::acos, expected); + } + + private static int testAtanCase(double input, double expected) { + return Tests.test("StrictMath.atan(double)", input, + StrictMath::atan, expected); + } +}