// Constants used by multiple algorithms

private static final double INFINITY = Double.POSITIVE_INFINITY;

private static final double TWO24 = 0x1.0p24; // 1.67772160000000000000e+07

private static final double TWO54 = 0x1.0p54; // 1.80143985094819840000e+16

private static final double HUGE = 1.0e+300;

/*

private static final int SIGN_BIT = 0x8000_0000;

private static final int EXP_BITS = 0x7ff0_0000;

private static final int EXP_SIGNIF_BITS = 0x7fff_ffff;

private FdLibm() {

throw new UnsupportedOperationException("No FdLibm instances for you.");

}

/**

private static int __LO(double x) {

long transducer = Double.doubleToRawLongBits(x);

return (int)transducer;

}

/**

private static double __LO(double x, int low) {

long transX = Double.doubleToRawLongBits(x);

return Double.longBitsToDouble((transX & 0xFFFF_FFFF_0000_0000L) |

(low & 0x0000_0000_FFFF_FFFFL));

}

/**

private static int __HI(double x) {

long transducer = Double.doubleToRawLongBits(x);

return (int)(transducer >> 32);

}

/**

private static double __HI(double x, int high) {

long transX = Double.doubleToRawLongBits(x);

return Double.longBitsToDouble((transX & 0x0000_0000_FFFF_FFFFL) |

( ((long)high)) << 32 );

}

/**

private static double __HI_LO(int high, int low) {

return Double.longBitsToDouble(((long)high << 32) |

(low & 0xffff_ffffL));

}

/** sin(x)

static final class Sin {

private Sin() {throw new UnsupportedOperationException();}

static double compute(double x) {

double[] y = new double[2];

double z = 0.0;

int n, ix;

// High word of x.

ix = __HI(x);

// |x| ~< pi/4

ix &= EXP_SIGNIF_BITS;

if (ix <= 0x3fe9_21fb) {

return __kernel_sin(x, z, 0);

} else if (ix >= EXP_BITS) { // sin(Inf or NaN) is NaN

return x - x;

} else { // argument reduction needed

n = RemPio2.__ieee754_rem_pio2(x, y);

switch(n & 3) {

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);

default:

return -Cos.__kernel_cos(y[0], y[1]);

}

}

}

/** __kernel_sin( x, y, iy)

private static final double

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

static double __kernel_sin(double x, double y, int iy) {

double z, r, v;

int ix;

ix = __HI(x) & EXP_SIGNIF_BITS; // high word of x

if (ix < 0x3e40_0000) { // |x| < 2**-27

if ((int)x == 0) // generate inexact

return x;

}

z = x*x;

v = z*x;

r = S2 + z*(S3 + z*(S4 + z*(S5 + z*S6)));

if (iy == 0) {

return x + v*(S1 + z*r);

} else {

return x - ((z*(0.5*y - v*r) - y) - v*S1);

}

}

}

/** cos(x)

static final class Cos {

private Cos() {throw new UnsupportedOperationException();}

static double compute(double x) {

double[] y = new double[2];

double z = 0.0;

int n, ix;

// High word of x.

ix = __HI(x);

// |x| ~< pi/4

ix &= EXP_SIGNIF_BITS;

if (ix <= 0x3fe9_21fb) {

return __kernel_cos(x, z);

} else if (ix >= EXP_BITS) { // cos(Inf or NaN) is NaN

return x - x;

} else { // argument reduction needed

n = RemPio2.__ieee754_rem_pio2(x,y);

switch (n & 3) {

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]);

default:

return Sin.__kernel_sin(y[0], y[1], 1);

}

}

}

/**

private static final double

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

static double __kernel_cos(double x, double y) {

double a, hz, z, r, qx = 0.0;

int ix;

ix = __HI(x) & EXP_SIGNIF_BITS; // ix = |x|'s high word

if (ix < 0x3e40_0000) { // if x < 2**27

if (((int)x) == 0) { // generate inexact

return 1.0;

}

}

z = x*x;

r = z*(C1 + z*(C2 + z*(C3 + z*(C4 + z*(C5 + z*C6)))));

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

qx = 0.28125;

} else {

qx = __HI_LO(ix - 0x0020_0000, 0);

}

hz = 0.5*z - qx;

a = 1.0 - qx;

return a - (hz - (z*r - x*y));

}

}

}

/** tan(x)

static final class Tan {

private Tan() {throw new UnsupportedOperationException();}

static double compute(double x) {

double[] y = new double[2];

double z = 0.0;

int n, ix;

// High word of x.

ix = __HI(x);

// |x| ~< pi/4

ix &= EXP_SIGNIF_BITS;

if (ix <= 0x3fe9_21fb) {

return __kernel_tan(x, z, 1);

} else if (ix >= EXP_BITS) { // tan(Inf or NaN) is NaN

return x - x; // NaN

} else { // argument reduction needed

n = RemPio2.__ieee754_rem_pio2(x, y);

return __kernel_tan(y[0], y[1], 1 - ((n & 1) << 1)); // 1 -- n even; -1 -- n odd

}

}

/** __kernel_tan( x, y, k )

private static final double

pio4 = 0x1.921fb54442d18p-1, // 7.85398163397448278999e-01

pio4lo= 0x1.1a62633145c07p-55; // 3.06161699786838301793e-17

@Stable

private static final double[]

T = {

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

};

static double __kernel_tan(double x, double y, int iy) {

double z, r, v, w, s;

int ix, hx;

hx = __HI(x); // high word of x

ix = hx & EXP_SIGNIF_BITS; // 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) {

return x;

} else { // compute -1 / (x+y) carefully

double a, t;

z = w = x + y;

z= __LO(z, 0);

v = y - (z - x);

t = a = -1.0 / w;

t = __LO(t, 0);

s = 1.0 + t * z;

return t + a * (s + t * v);

}

}

}

}

if (ix >= 0x3FE5_9428) { // |x| >= 0.6744

if ( hx < 0) {

x = -x;

y = -y;

}

z = pio4 - x;

w = pio4lo - y;

x = z + w;

y = 0.0;

}

z = x*x;

w = z*z;

/* Break x^5*(T[1]+x^2*T[2]+...) into

r = T[1] + w*(T[3] + w*(T[5] + w*(T[7] + w*(T[9] + w*T[11]))));

v = z*(T[2] + w*(T[4] + w*(T[6] + w*(T[8] + w*(T[10] + w*T[12])))));

s = z*x;

r = y + z*(s*(r + v) + y);

r += T[0]*s;

w = x + r;

if (ix >= 0x3FE5_9428) {

v = (double)iy;

return (double)(1-((hx >> 30) & 2))*(v - 2.0*(x - (w*w/(w + v) - r)));

}

if (iy == 1) {

return w;

} else { /* if were to allow error up to 2 ulp,

// compute -1.0/(x + r) accurately

double a,t;

z = w;

z = __LO(z, 0);

v = r - (z - x); // z + v = r + x

t = a = -1.0/w; // a = -1.0/w

t = __LO(t, 0);

s = 1.0 + t*z;

return t + a*(s + t*v);

}

}

}

/** __ieee754_rem_pio2(x,y)

static final class RemPio2 {

/*

@Stable

private static final int[] two_over_pi = {

0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,

0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,

0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,

0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,

0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,

0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,

0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,

0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,

0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,

0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,

0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,

};

@Stable

private static final int[] npio2_hw = {

0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,

0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,

0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,

0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,

0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,

0x404858EB, 0x404921FB,

};

/*

private static final double

invpio2 = 0x1.45f306dc9c883p-1, // 6.36619772367581382433e-01

pio2_1 = 0x1.921fb544p0, // 1.57079632673412561417e+00

pio2_1t = 0x1.0b4611a626331p-34, // 6.07710050650619224932e-11

pio2_2 = 0x1.0b4611a6p-34, // 6.07710050630396597660e-11

pio2_2t = 0x1.3198a2e037073p-69, // 2.02226624879595063154e-21

pio2_3 = 0x1.3198a2ep-69, // 2.02226624871116645580e-21

pio2_3t = 0x1.b839a252049c1p-104; // 8.47842766036889956997e-32

static int __ieee754_rem_pio2(double x, double[] y) {

double z = 0.0, w, t, r, fn;

double[] tx = new double[3];

int e0, i, j, nx, n, ix, hx;

hx = __HI(x); // high word of x

ix = hx & EXP_SIGNIF_BITS;

if (ix <= 0x3fe9_21fb) { // |x| ~<= pi/4 , no need for reduction

y[0] = x;

y[1] = 0;

return 0;

}

if (ix < 0x4002_d97c) { // |x| < 3pi/4, special case with n=+-1

if (hx > 0) {

z = x - pio2_1;

if (ix != 0x3ff9_21fb) { // 33+53 bit pi is good enough

y[0] = z - pio2_1t;

y[1] = (z - y[0]) - pio2_1t;

} else { // near pi/2, use 33+33+53 bit pi

z -= pio2_2;

y[0] = z - pio2_2t;

y[1] = (z-y[0])-pio2_2t;

}

return 1;

} else { // negative x

z = x + pio2_1;

if (ix != 0x3ff_921fb) { // 33+53 bit pi is good enough

y[0] = z + pio2_1t;

y[1] = (z - y[0]) + pio2_1t;

} else { // near pi/2, use 33+33+53 bit pi

z += pio2_2;

y[0] = z + pio2_2t;

y[1] = (z - y[0]) + pio2_2t;

}

return -1;

}

}

if (ix <= 0x4139_21fb) { // |x| ~<= 2^19*(pi/2), medium size

t = Math.abs(x);

n = (int) (t*invpio2 + 0.5);

fn = (double)n;

r = t - fn*pio2_1;

w = fn*pio2_1t; // 1st round good to 85 bit

if (n < 32 && ix != npio2_hw[n - 1]) {

y[0] = r - w; // quick check no cancellation

} else {

j = ix >> 20;

y[0] = r - w;

i = j - (((__HI(y[0])) >> 20) & 0x7ff);

if (i > 16) { // 2nd iteration needed, good to 118

t = r;

w = fn*pio2_2;

r = t - w;

w = fn*pio2_2t - ((t - r) - w);

y[0] = r - w;

i = j - (((__HI(y[0])) >> 20) & 0x7ff);

if (i > 49) { // 3rd iteration need, 151 bits acc

t = r; // will cover all possible cases

w = fn*pio2_3;

r = t - w;

w = fn*pio2_3t - ((t - r) - w);

y[0] = r - w;

}

}

}

y[1] = (r - y[0]) - w;

if (hx < 0) {

y[0] = -y[0];

y[1] = -y[1];

return -n;

} else {

return n;

}

}

/*

if (ix >= EXP_BITS) { // x is inf or NaN

y[0] = y[1] = x - x;

return 0;

}

// set z = scalbn(|x|, ilogb(x)-23)

z = __LO(z, __LO(x));

e0 = (ix >> 20) - 1046; // e0 = ilogb(z) - 23;

z = __HI(z, ix - (e0 << 20));

for (i=0; i < 2; i++) {

tx[i] = (double)((int)(z));

z = (z - tx[i])*TWO24;

}

tx[2] = z;

nx = 3;

while (tx[nx - 1] == 0.0) { // skip zero term

nx--;

}

n = KernelRemPio2.__kernel_rem_pio2(tx, y, e0, nx, 2, two_over_pi);

if (hx < 0) {

y[0] = -y[0];

y[1] = -y[1];

return -n;

}

return n;

}

}

/**

static final class KernelRemPio2 {

/*

@Stable

private static final int init_jk[] = {2, 3, 4, 6}; // initial value for jk

@Stable

private static final double PIo2[] = {

0x1.921fb4p0, // 1.57079625129699707031e+00

0x1.4442dp-24, // 7.54978941586159635335e-08

0x1.846988p-48, // 5.39030252995776476554e-15

0x1.8cc516p-72, // 3.28200341580791294123e-22

0x1.01b838p-96, // 1.27065575308067607349e-29

0x1.a25204p-120, // 1.22933308981111328932e-36

0x1.382228p-145, // 2.73370053816464559624e-44

0x1.9f31dp-169, // 2.16741683877804819444e-51

};

static final double

twon24 = 0x1.0p-24; // 5.96046447753906250000e-08

static int __kernel_rem_pio2(double[] x, double[] y, int e0, int nx, int prec, final int[] ipio2) {

int jz, jx, jv, jp, jk, carry, n, i, j, k, m, q0, ih;

int[] iq = new int[20];

double z,fw;

double [] f = new double[20];

double [] fq= new double[20];

double [] q = new double[20];

// initialize jk

jk = init_jk[prec];

jp = jk;

// determine jx, jv, q0, note that 3 > q0

jx = nx - 1;

jv = (e0 - 3)/24;

if (jv < 0) {

jv = 0;

}

q0 = e0 - 24*(jv + 1);

// set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk]

j = jv - jx;

m = jx + jk;

for (i = 0; i <= m; i++, j++) {

f[i] = (j < 0) ? 0.0 : (double) ipio2[j];

}

// compute q[0],q[1],...q[jk]

for (i=0; i <= jk; i++) {

for(j = 0, fw = 0.0; j <= jx; j++) {

fw += x[j]*f[jx + i - j];

}

q[i] = fw;

}

jz = jk;

while(true) {

// distill q[] into iq[] reversingly

for(i=0, j=jz, z=q[jz]; j > 0; i++, j--) {

fw = (double)((int)(twon24* z));

iq[i] = (int)(z - TWO24*fw);

z = q[j - 1] + fw;

}

// compute n

z = Math.scalb(z, q0); // actual value of z

z -= 8.0*Math.floor(z*0.125); // trim off integer >= 8

n = (int) z;

z -= (double)n;

ih = 0;

if (q0 > 0) { // need iq[jz - 1] to determine n

i = (iq[jz - 1] >> (24 - q0));

n += i;

iq[jz - 1] -= i << (24 - q0);

ih = iq[jz - 1] >> (23 - q0);

} else if (q0 == 0) {

ih = iq[jz-1]>>23;

} else if (z >= 0.5) {

ih=2;

}

if (ih > 0) { // q > 0.5

n += 1;

carry = 0;

for (i=0; i < jz; i++) { // compute 1-q

j = iq[i];

if (carry == 0) {

if (j != 0) {

carry = 1;

iq[i] = 0x100_0000 - j;

}

} else {

iq[i] = 0xff_ffff - j;

}

}

if (q0 > 0) { // rare case: chance is 1 in 12

switch(q0) {

case 1:

iq[jz-1] &= 0x7f_ffff;

break;

case 2:

iq[jz-1] &= 0x3f_ffff;

break;

}

}

if (ih == 2) {

z = 1.0 - z;

if (carry != 0) {

z -= Math.scalb(1.0, q0);

}

}

}

// check if recomputation is needed

if (z == 0.0) {

j = 0;

for (i = jz - 1; i >= jk; i--) {

j |= iq[i];

}

if (j == 0) { // need recomputation

for(k=1; iq[jk - k] == 0; k++); // k = no. of terms needed

for(i = jz + 1; i <= jz + k; i++) { // add q[jz+1] to q[jz+k]

f[jx + i] = (double) ipio2[jv + i];

for (j=0, fw = 0.0; j <= jx; j++) {

fw += x[j]*f[jx + i - j];

}

q[i] = fw;

}

jz += k;

continue;

} else {

break;

}

} else {

break;

}

}

// chop off zero terms

if (z == 0.0) {

jz -= 1;

q0 -= 24;

while (iq[jz] == 0) {

jz--;

q0-=24;

}

} else { // break z into 24-bit if necessary

z = Math.scalb(z, -q0);

if (z >= TWO24) {

fw = (double)((int)(twon24*z));

iq[jz] = (int)(z - TWO24*fw);

jz += 1;

q0 += 24;

iq[jz] = (int) fw;

} else {

iq[jz] = (int) z;

}

}

// convert integer "bit" chunk to floating-point value

fw = Math.scalb(1.0, q0);

for(i = jz; i>=0; i--) {

q[i] = fw*(double)iq[i];

fw *= twon24;

}

// compute PIo2[0,...,jp]*q[jz,...,0]

for(i = jz; i>=0; i--) {

for (fw = 0.0, k = 0; k <= jp && k <= jz-i; k++) {

fw += PIo2[k]*q[i + k];

}

fq[jz - i] = fw;

}

// compress fq[] into y[]

switch(prec) {

case 0:

fw = 0.0;

for (i = jz; i >=0; i--) {

fw += fq[i];

}

y[0] = (ih == 0)? fw: -fw;

break;

case 1:

case 2:

fw = 0.0;