RFR: JDK-8302027: Port fdlibm trig functions (sin, cos, tan) to Java
Joe Darcy
darcy at openjdk.org
Wed Mar 1 06:15:07 UTC 2023
On Wed, 1 Mar 2023 05:28:34 GMT, Joe Darcy <darcy at openjdk.org> wrote:
> Last and certainly not least in the port of FDLIBM to Java, the transcendental methods for sin, cos, and tan.
>
> Some more tests are to be written in the StrictMath directory to verify that the StrictMath algorihtm for sin/cos/tan is being used rather than a different one. However, I wanted to get the rest of the change out for review first.
>
> The sin/cos/tan methods are grouped together since they share the same argument reduction logic. Argument reduction is the process of mapping an argument of a function to an argument in a restricted range (and possibly returning some function of the reduced argument). For sin, cos, and tan, since they are fundamentally periodic with respect to a multiple of pi, argument reduction is done to find the remainder of the original argument with respect to pi/2.
And for cos:
$ diff -w Cos.c Cos.translit.java
1c1
< /* cos(x)
---
> /** cos(x)
31,41c31,34
<
< #include "fdlibm.h"
<
< #ifdef __STDC__
< double cos(double x)
< #else
< double cos(x)
< double x;
< #endif
< {
< double y[2],z=0.0;
---
> static class Cos {
> static double compute(double x) {
> double[] y = new double[2];
> double z=0.0;
56c49
< n = __ieee754_rem_pio2(x,y);
---
> n = RemPio2.__ieee754_rem_pio2(x,y);
58,60c51,53
< case 0: return __kernel_cos(y[0],y[1]);
< case 1: return -__kernel_sin(y[0],y[1],1);
< case 2: return -__kernel_cos(y[0],y[1]);
---
> case 0: return Cos.__kernel_cos(y[0],y[1]);
> case 1: return -Sin.__kernel_sin(y[0],y[1],1);
> case 2: return -Cos.__kernel_cos(y[0],y[1]);
62c55
< return __kernel_sin(y[0],y[1],1);
---
> return Sin.__kernel_sin(y[0],y[1],1);
66c59,60
< /*
---
>
> /**
100,107c94
<
< #include "fdlibm.h"
<
< #ifdef __STDC__
< static const double
< #else
< static double
< #endif
---
> private static final double
116,123c103,104
< #ifdef __STDC__
< double __kernel_cos(double x, double y)
< #else
< double __kernel_cos(x, y)
< double x,y;
< #endif
< {
< double a,hz,z,r,qx;
---
> static double __kernel_cos(double x, double y) {
> double a,hz,z,r,qx = 0.0;
137,138c118,121
< __HI(qx) = ix-0x00200000; /* x/4 */
< __LO(qx) = 0;
---
> //__HI(qx) = ix-0x00200000; /* x/4 */
> qx = __HI(qx, ix-0x00200000);
> // __LO(qx) = 0;
> qx = __LO(qx, 0);
144a128
> }
$ diff -w Cos.translit.java Cos.fdlibm.java
31a32,33
> private Cos() {throw new UnsupportedOperationException();}
>
37c39
< /* High word of x. */
---
> // High word of x.
40,48c42,48
< /* |x| ~< pi/4 */
< ix &= 0x7fffffff;
< if(ix <= 0x3fe921fb) return __kernel_cos(x,z);
<
< /* cos(Inf or NaN) is NaN */
< else if (ix>=0x7ff00000) return x-x;
<
< /* argument reduction needed */
< else {
---
> // |x| ~< pi/4
> ix &= 0x7fff_ffff;
> if (ix <= 0x3fe9_21fb) {
> return __kernel_cos(x, z);
> } else if (ix >= 0x7ff0_0000) { // cos(Inf or NaN) is NaN
> return x-x;
> } else { // argument reduction needed
95,101c95,100
< one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
< C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
< C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
< C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
< C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
< C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
< C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
---
> C1 = 0x1.555555555554cp-5, // 4.16666666666666019037e-02
> C2 = -0x1.6c16c16c15177p-10, // -1.38888888888741095749e-03
> C3 = 0x1.a01a019cb159p-16, // 2.48015872894767294178e-05
> C4 = -0x1.27e4f809c52adp-22, // -2.75573143513906633035e-07
> C5 = 0x1.1ee9ebdb4b1c4p-29, // 2.08757232129817482790e-09
> C6 = -0x1.8fae9be8838d4p-37; // -1.13596475577881948265e-11
106,108c105,109
< ix = __HI(x)&0x7fffffff; /* ix = |x|'s high word*/
< if(ix<0x3e400000) { /* if x < 2**27 */
< if(((int)x)==0) return one; /* generate inexact */
---
> ix = __HI(x) & 0x7fff_ffff; // ix = |x|'s high word
> if (ix < 0x3e40_0000) { // if x < 2**27
> if (((int)x) == 0) { // generate inexact
> return 1.0;
> }
112,115c113,116
< if(ix < 0x3FD33333) /* if |x| < 0.3 */
< return one - (0.5*z - (z*r - x*y));
< else {
< if(ix > 0x3fe90000) { /* x > 0.78125 */
---
> if (ix < 0x3FD3_3333) { // if |x| < 0.3
> return 1.0 - (0.5*z - (z*r - x*y));
> } else {
> if (ix > 0x3fe9_0000) { // x > 0.78125
118,121c119
< //__HI(qx) = ix-0x00200000; /* x/4 */
< qx = __HI(qx, ix-0x00200000);
< // __LO(qx) = 0;
< qx = __LO(qx, 0);
---
> qx = __HI_LO(ix - 0x0020_0000, 0);
124c122
< a = one-qx;
---
> a = 1.0 - qx;
And tan:
$ diff -w Tan.c Tan.translit.java
1c1
< /* tan(x)
---
> /** tan(x)
30,40c30,33
<
< #include "fdlibm.h"
<
< #ifdef __STDC__
< double tan(double x)
< #else
< double tan(x)
< double x;
< #endif
< {
< double y[2],z=0.0;
---
> static class Tan {
> static double compute(double x) {
> double[] y= new double[2];
> double z=0.0;
55c48
< n = __ieee754_rem_pio2(x,y);
---
> n = RemPio2.__ieee754_rem_pio2(x,y);
60c53,54
< /* __kernel_tan( x, y, k )
---
>
> /** __kernel_tan( x, y, k )
93,99c87
<
< #include "fdlibm.h"
< #ifdef __STDC__
< static const double
< #else
< static double
< #endif
---
> private static final double
119,125c107
< #ifdef __STDC__
< double __kernel_tan(double x, double y, int iy)
< #else
< double __kernel_tan(x, y, iy)
< double x,y; int iy;
< #endif
< {
---
> static double __kernel_tan(double x, double y, int iy) {
133c115
< return one / fabs(x);
---
> return one / Math.abs(x);
141c123,124
< __LO(z) = 0;
---
> // __LO(z) = 0;
> z= __LO(z, 0);
144c127,128
< __LO(t) = 0;
---
> //__LO(t) = 0;
> t = __LO(t, 0);
179c163,164
< __LO(z) = 0;
---
> // __LO(z) = 0;
> z = __LO(z, 0);
182c167,168
< __LO(t) = 0;
---
> // __LO(t) = 0;
> t = __LO(t, 0);
186a173
> }
$ diff -w Tan.translit.java Tan.fdlibm.java
30a31,32
> private Tan() {throw new UnsupportedOperationException();}
>
36c38
< /* High word of x. */
---
> // High word of x.
39,47c41,47
< /* |x| ~< pi/4 */
< ix &= 0x7fffffff;
< if(ix <= 0x3fe921fb) return __kernel_tan(x,z,1);
<
< /* tan(Inf or NaN) is NaN */
< else if (ix>=0x7ff00000) return x-x; /* NaN */
<
< /* argument reduction needed */
< else {
---
> // |x| ~< pi/4
> ix &= 0x7fff_ffff;
> if (ix <= 0x3fe9_21fb) {
> return __kernel_tan(x, z, 1);
> } else if (ix >= 0x7ff0_0000) { // tan(Inf or NaN) is NaN
> return x-x; // NaN
> } else { // argument reduction needed
49,50c49
< return __kernel_tan(y[0],y[1],1-((n&1)<<1)); /* 1 -- n even
< -1 -- n odd */
---
> return __kernel_tan(y[0], y[1], 1 - ((n & 1) << 1)); // 1 -- n even; -1 -- n odd
88,90c87,88
< one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
< pio4 = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */
< pio4lo= 3.06161699786838301793e-17, /* 0x3C81A626, 0x33145C07 */
---
> pio4 = 0x1.921fb54442d18p-1, // 7.85398163397448278999e-01
> pio4lo= 0x1.1a62633145c07p-55, // 3.06161699786838301793e-17
92,104c90,102
< 3.33333333333334091986e-01, /* 0x3FD55555, 0x55555563 */
< 1.33333333333201242699e-01, /* 0x3FC11111, 0x1110FE7A */
< 5.39682539762260521377e-02, /* 0x3FABA1BA, 0x1BB341FE */
< 2.18694882948595424599e-02, /* 0x3F9664F4, 0x8406D637 */
< 8.86323982359930005737e-03, /* 0x3F8226E3, 0xE96E8493 */
< 3.59207910759131235356e-03, /* 0x3F6D6D22, 0xC9560328 */
< 1.45620945432529025516e-03, /* 0x3F57DBC8, 0xFEE08315 */
< 5.88041240820264096874e-04, /* 0x3F4344D8, 0xF2F26501 */
< 2.46463134818469906812e-04, /* 0x3F3026F7, 0x1A8D1068 */
< 7.81794442939557092300e-05, /* 0x3F147E88, 0xA03792A6 */
< 7.14072491382608190305e-05, /* 0x3F12B80F, 0x32F0A7E9 */
< -1.85586374855275456654e-05, /* 0xBEF375CB, 0xDB605373 */
< 2.59073051863633712884e-05, /* 0x3EFB2A70, 0x74BF7AD4 */
---
> 0x1.5555555555563p-2, // 3.33333333333334091986e-01
> 0x1.111111110fe7ap-3, // 1.33333333333201242699e-01
> 0x1.ba1ba1bb341fep-5, // 5.39682539762260521377e-02
> 0x1.664f48406d637p-6, // 2.18694882948595424599e-02
> 0x1.226e3e96e8493p-7, // 8.86323982359930005737e-03
> 0x1.d6d22c9560328p-9, // 3.59207910759131235356e-03
> 0x1.7dbc8fee08315p-10, // 1.45620945432529025516e-03
> 0x1.344d8f2f26501p-11, // 5.88041240820264096874e-04
> 0x1.026f71a8d1068p-12, // 2.46463134818469906812e-04
> 0x1.47e88a03792a6p-14, // 7.81794442939557092300e-05
> 0x1.2b80f32f0a7e9p-14, // 7.14072491382608190305e-05
> -0x1.375cbdb605373p-16, // -1.85586374855275456654e-05
> 0x1.b2a7074bf7ad4p-16, // 2.59073051863633712884e-05
110,117c108,115
< hx = __HI(x); /* high word of x */
< ix = hx&0x7fffffff; /* high word of |x| */
< if(ix<0x3e300000) { /* x < 2**-28 */
< if((int)x==0) { /* generate inexact */
< if (((ix | __LO(x)) | (iy + 1)) == 0)
< return one / Math.abs(x);
< else {
< if (iy == 1)
---
> hx = __HI(x); // high word of x
> ix = hx&0x7fff_ffff; // high word of |x|
> if (ix < 0x3e30_0000) { // x < 2**-28
> if ((int)x == 0) { // generate inexact
> if (((ix | __LO(x)) | (iy + 1)) == 0) {
> return 1.0 / Math.abs(x);
> } else {
> if (iy == 1) {
119c117
< else { /* compute -1 / (x+y) carefully */
---
> } else { // compute -1 / (x+y) carefully
123d120
< // __LO(z) = 0;
126,127c123
< t = a = -one / w;
< //__LO(t) = 0;
---
> t = a = -1.0 / w;
129c125
< s = one + t * z;
---
> s = 1.0 + t * z;
135,136c131,135
< if(ix>=0x3FE59428) { /* |x|>=0.6744 */
< if(hx<0) {x = -x; y = -y;}
---
> if (ix >= 0x3FE5_9428) { // |x| >= 0.6744
> if ( hx < 0) {
> x = -x;
> y = -y;
> }
139c138,139
< x = z+w; y = 0.0;
---
> x = z + w;
> y = 0.0;
153c153
< if(ix>=0x3FE59428) {
---
> if (ix >= 0x3FE5_9428) {
157,160c157,161
< if(iy==1) return w;
< else { /* if allow error up to 2 ulp,
< simply return -1.0/(x+r) here */
< /* compute -1.0/(x+r) accurately */
---
> if (iy == 1) {
> return w;
> } else { /* if were to allow error up to 2 ulp,
> could simply return -1.0/(x + r) here */
> // compute -1.0/(x + r) accurately
163d163
< // __LO(z) = 0;
165,167c165,166
< v = r-(z - x); /* z+v = r+x */
< t = a = -1.0/w; /* a = -1.0/w */
< // __LO(t) = 0;
---
> v = r - (z - x); // z + v = r + x
> t = a = -1.0/w; // a = -1.0/w
-------------
PR: https://git.openjdk.org/jdk/pull/12800
More information about the core-libs-dev
mailing list