RFR: JDK-8302027: Port fdlibm trig functions (sin, cos, tan) to Java
Joe Darcy
darcy at openjdk.org
Wed Mar 1 06:09:03 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.
Various diffs to aid the review:
$ diff -w Sin.c Sin.translit.java
1c1
< /* sin(x)
---
> /** sin(x)
31,41c31,34
<
< #include "fdlibm.h"
<
< #ifdef __STDC__
< double sin(double x)
< #else
< double sin(x)
< double x;
< #endif
< {
< double y[2],z=0.0;
---
> static class Sin {
> 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_sin(y[0],y[1],1);
< case 1: return __kernel_cos(y[0],y[1]);
< case 2: return -__kernel_sin(y[0],y[1],1);
---
> case 0: return Sin.__kernel_sin(y[0],y[1],1);
> case 1: return Cos.__kernel_cos(y[0],y[1]);
> case 2: return -Sin.__kernel_sin(y[0],y[1],1);
62c55
< return -__kernel_cos(y[0],y[1]);
---
> return -Cos.__kernel_cos(y[0],y[1]);
66c59,60
< /* __kernel_sin( x, y, iy)
---
>
> /** __kernel_sin( x, y, iy)
93,100c87
<
< #include "fdlibm.h"
<
< #ifdef __STDC__
< static const double
< #else
< static double
< #endif
---
> private static final double
109,115c96
< #ifdef __STDC__
< double __kernel_sin(double x, double y, int iy)
< #else
< double __kernel_sin(x, y, iy)
< double x,y; int iy; /* iy=0 if y is zero */
< #endif
< {
---
> static double __kernel_sin(double x, double y, int iy) {
126a108
> }
$ diff -w Sin.translit.java Sin.fdlibm.java
31a32,33
> private Sin() {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_sin(x,z,0);
<
< /* sin(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_sin(x, z, 0);
> } else if (ix>=0x7ff0_0000) { // sin(Inf or NaN) is NaN
> return x - x;
> } else { // argument reduction needed
88,94c88,93
< half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
< S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
< S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
< S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
< S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
< S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
< S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
---
> S1 = -0x1.5555555555549p-3, // -1.66666666666666324348e-01
> S2 = 0x1.111111110f8a6p-7, // 8.33333333332248946124e-03
> S3 = -0x1.a01a019c161d5p-13, // -1.98412698298579493134e-04
> S4 = 0x1.71de357b1fe7dp-19, // 2.75573137070700676789e-06
> S5 = -0x1.ae5e68a2b9cebp-26, // -2.50507602534068634195e-08
> S6 = 0x1.5d93a5acfd57cp-33; // 1.58969099521155010221e-10
99,101c98,102
< ix = __HI(x)&0x7fffffff; /* high word of x */
< if(ix<0x3e400000) /* |x| < 2**-27 */
< {if((int)x==0) return x;} /* generate inexact */
---
> ix = __HI(x) & 0x7fff_ffff; // high word of x
> if (ix < 0x3e40_0000) { // |x| < 2**-27
> if ((int)x == 0) // generate inexact
> return x;
> }
105,106c106,110
< if(iy==0) return x+v*(S1+z*r);
< else return x-((z*(half*y-v*r)-y)-v*S1);
---
> if (iy == 0) {
> return x + v*(S1 + z*r);
> } else {
> return x - ((z*(0.5*y - v*r) - y) - v*S1);
> }
-------------
PR: https://git.openjdk.org/jdk/pull/12800
More information about the core-libs-dev
mailing list