#include <fenv.h>
/* Helper macros to compute sin of the input values. */
-#define POLYNOMIAL2(xx) ((((s5.x * (xx) + s4.x) * (xx) + s3.x) * (xx) + s2.x) \
- * (xx))
+#define POLYNOMIAL2(xx) ((((s5 * (xx) + s4) * (xx) + s3) * (xx) + s2) * (xx))
-#define POLYNOMIAL(xx) (POLYNOMIAL2 (xx) + s1.x)
+#define POLYNOMIAL(xx) (POLYNOMIAL2 (xx) + s1)
/* The computed polynomial is a variation of the Taylor series expansion for
sin(a):
static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ \
double xx = (x0) * (x0); \
double x1 = ((x0) + th2_36) - th2_36; \
- double y = aa.x * x1 * x1 * x1; \
+ double y = aa * x1 * x1 * x1; \
double r = (x0) + y; \
double x2 = ((x0) - x1) + (dx); \
- double t = (((POLYNOMIAL2 (xx) + bb.x) * xx + 3.0 * aa.x * x1 * x2) \
- * (x0) + aa.x * x2 * x2 * x2 + (dx)); \
+ double t = (((POLYNOMIAL2 (xx) + bb) * xx + 3.0 * aa * x1 * x2) \
+ * (x0) + aa * x2 * x2 * x2 + (dx)); \
t = (((x0) - r) + y) + t; \
double res = r + t; \
(cor) = (r - res) + t; \
/*---------------------------- 0.25<|x|< 0.855469---------------------- */
else if (k < 0x3feb6000)
{
- u.x = (m > 0) ? big.x + x : big.x - x;
- y = (m > 0) ? x - (u.x - big.x) : x + (u.x - big.x);
+ u.x = (m > 0) ? big + x : big - x;
+ y = (m > 0) ? x - (u.x - big) : x + (u.x - big);
xx = y * y;
s = y + y * xx * (sn3 + xx * sn5);
c = xx * (cs2 + xx * (cs4 + xx * cs6));
else if (k < 0x400368fd)
{
- y = (m > 0) ? hp0.x - x : hp0.x + x;
+ y = (m > 0) ? hp0 - x : hp0 + x;
if (y >= 0)
{
- u.x = big.x + y;
- y = (y - (u.x - big.x)) + hp1.x;
+ u.x = big + y;
+ y = (y - (u.x - big)) + hp1;
}
else
{
- u.x = big.x - y;
- y = (-hp1.x) - (y + (u.x - big.x));
+ u.x = big - y;
+ y = (-hp1) - (y + (u.x - big));
}
xx = y * y;
s = y + y * xx * (sn3 + xx * sn5);
/*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
else if (k < 0x419921FB)
{
- t = (x * hpinv.x + toint.x);
- xn = t - toint.x;
+ t = (x * hpinv + toint);
+ xn = t - toint;
v.x = t;
- y = (x - xn * mp1.x) - xn * mp2.x;
+ y = (x - xn * mp1) - xn * mp2;
n = v.i[LOW_HALF] & 3;
- da = xn * mp3.x;
+ da = xn * mp3;
a = y - da;
da = (y - a) - da;
eps = ABS (x) * 1.2e-30;
t = -a;
db = -da;
}
- u.x = big.x + t;
- y = t - (u.x - big.x);
+ u.x = big + t;
+ y = t - (u.x - big);
xx = y * y;
s = y + (db + y * xx * (sn3 + xx * sn5));
c = y * db + xx * (cs2 + xx * (cs4 + xx * cs6));
a = -a;
da = -da;
}
- u.x = big.x + a;
- y = a - (u.x - big.x) + da;
+ u.x = big + a;
+ y = a - (u.x - big) + da;
xx = y * y;
SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
s = y + y * xx * (sn3 + xx * sn5);
/*---------------------105414350 <|x|< 281474976710656 --------------------*/
else if (k < 0x42F00000)
{
- t = (x * hpinv.x + toint.x);
- xn = t - toint.x;
+ t = (x * hpinv + toint);
+ xn = t - toint;
v.x = t;
xn1 = (xn + 8.0e22) - 8.0e22;
xn2 = xn - xn1;
- y = ((((x - xn1 * mp1.x) - xn1 * mp2.x) - xn2 * mp1.x) - xn2 * mp2.x);
+ y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2);
n = v.i[LOW_HALF] & 3;
- da = xn1 * pp3.x;
+ da = xn1 * pp3;
t = y - da;
da = (y - t) - da;
- da = (da - xn2 * pp3.x) - xn * pp4.x;
+ da = (da - xn2 * pp3) - xn * pp4;
a = t + da;
da = (t - a) + da;
eps = 1.0e-24;
t = -a;
db = -da;
}
- u.x = big.x + t;
- y = t - (u.x - big.x);
+ u.x = big + t;
+ y = t - (u.x - big);
xx = y * y;
s = y + (db + y * xx * (sn3 + xx * sn5));
c = y * db + xx * (cs2 + xx * (cs4 + xx * cs6));
a = -a;
da = -da;
}
- u.x = big.x + a;
- y = a - (u.x - big.x) + da;
+ u.x = big + a;
+ y = a - (u.x - big) + da;
xx = y * y;
SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
s = y + y * xx * (sn3 + xx * sn5);
else if (k < 0x3feb6000)
{ /* 2^-27 < |x| < 0.855469 */
y = ABS (x);
- u.x = big.x + y;
- y = y - (u.x - big.x);
+ u.x = big + y;
+ y = y - (u.x - big);
xx = y * y;
s = y + y * xx * (sn3 + xx * sn5);
c = xx * (cs2 + xx * (cs4 + xx * cs6));
else if (k < 0x400368fd)
{ /* 0.855469 <|x|<2.426265 */ ;
- y = hp0.x - ABS (x);
- a = y + hp1.x;
- da = (y - a) + hp1.x;
+ y = hp0 - ABS (x);
+ a = y + hp1;
+ da = (y - a) + hp1;
xx = a * a;
if (xx < 0.01588)
{
t = -a;
db = -da;
}
- u.x = big.x + t;
- y = t - (u.x - big.x);
+ u.x = big + t;
+ y = t - (u.x - big);
xx = y * y;
s = y + (db + y * xx * (sn3 + xx * sn5));
c = y * db + xx * (cs2 + xx * (cs4 + xx * cs6));
else if (k < 0x419921FB)
{ /* 2.426265<|x|< 105414350 */
- t = (x * hpinv.x + toint.x);
- xn = t - toint.x;
+ t = (x * hpinv + toint);
+ xn = t - toint;
v.x = t;
- y = (x - xn * mp1.x) - xn * mp2.x;
+ y = (x - xn * mp1) - xn * mp2;
n = v.i[LOW_HALF] & 3;
- da = xn * mp3.x;
+ da = xn * mp3;
a = y - da;
da = (y - a) - da;
eps = ABS (x) * 1.2e-30;
t = -a;
db = -da;
}
- u.x = big.x + t;
- y = t - (u.x - big.x);
+ u.x = big + t;
+ y = t - (u.x - big);
xx = y * y;
s = y + (db + y * xx * (sn3 + xx * sn5));
c = y * db + xx * (cs2 + xx * (cs4 + xx * cs6));
a = -a;
da = -da;
}
- u.x = big.x + a;
- y = a - (u.x - big.x) + da;
+ u.x = big + a;
+ y = a - (u.x - big) + da;
xx = y * y;
SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
s = y + y * xx * (sn3 + xx * sn5);
else if (k < 0x42F00000)
{
- t = (x * hpinv.x + toint.x);
- xn = t - toint.x;
+ t = (x * hpinv + toint);
+ xn = t - toint;
v.x = t;
xn1 = (xn + 8.0e22) - 8.0e22;
xn2 = xn - xn1;
- y = ((((x - xn1 * mp1.x) - xn1 * mp2.x) - xn2 * mp1.x) - xn2 * mp2.x);
+ y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2);
n = v.i[LOW_HALF] & 3;
- da = xn1 * pp3.x;
+ da = xn1 * pp3;
t = y - da;
da = (y - t) - da;
- da = (da - xn2 * pp3.x) - xn * pp4.x;
+ da = (da - xn2 * pp3) - xn * pp4;
a = t + da;
da = (t - a) + da;
eps = 1.0e-24;
t = -a;
db = -da;
}
- u.x = big.x + t;
- y = t - (u.x - big.x);
+ u.x = big + t;
+ y = t - (u.x - big);
xx = y * y;
s = y + (db + y * xx * (sn3 + xx * sn5));
c = y * db + xx * (cs2 + xx * (cs4 + xx * cs6));
a = -a;
da = -da;
}
- u.x = big.x + a;
- y = a - (u.x - big.x) + da;
+ u.x = big + a;
+ y = a - (u.x - big) + da;
xx = y * y;
SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
s = y + y * xx * (sn3 + xx * sn5);
double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, c1, c2, xx, cor, res;
static const double t22 = 6291456.0;
y = ABS (x);
- u.x = big.x + y;
- y = y - (u.x - big.x);
+ u.x = big + y;
+ y = y - (u.x - big);
xx = y * y;
s = y * xx * (sn3 + xx * sn5);
c = xx * (cs2 + xx * (cs4 + xx * cs6));
double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, e1, e2, xx, cor, res, del;
static const double t22 = 6291456.0;
y = ABS (x);
- y = hp0.x - y;
+ y = hp0 - y;
if (y >= 0)
{
- u.x = big.x + y;
- y = y - (u.x - big.x);
- del = hp1.x;
+ u.x = big + y;
+ y = y - (u.x - big);
+ del = hp1;
}
else
{
- u.x = big.x - y;
- y = -(y + (u.x - big.x));
- del = -hp1.x;
+ u.x = big - y;
+ y = -(y + (u.x - big));
+ del = -hp1;
}
xx = y * y;
s = y * xx * (sn3 + xx * sn5);
return (x > 0) ? res : -res;
else
{
- y = ABS (x) - hp0.x;
- y1 = y - hp1.x;
- y2 = (y - y1) - hp1.x;
+ y = ABS (x) - hp0;
+ y1 = y - hp1;
+ y2 = (y - y1) - hp1;
__docos (y1, y2, w);
if (w[0] == w[0] + 1.000000005 * w[1])
return (x > 0) ? w[0] : -w[0];
return (x > 0) ? w[0] : -w[0];
else
{
- t = (orig * hpinv.x + toint.x);
- xn = t - toint.x;
+ t = (orig * hpinv + toint);
+ xn = t - toint;
v.x = t;
- y = (orig - xn * mp1.x) - xn * mp2.x;
+ y = (orig - xn * mp1) - xn * mp2;
n = v.i[LOW_HALF] & 3;
- da = xn * pp3.x;
+ da = xn * pp3;
t = y - da;
da = (y - t) - da;
- y = xn * pp4.x;
+ y = xn * pp4;
a = t - y;
da = ((t - a) - y) + da;
if (n & 2)
static const double t22 = 6291456.0;
y = ABS (x);
- u.x = big.x + y;
- y = y - (u.x - big.x);
+ u.x = big + y;
+ y = y - (u.x - big);
dx = (x > 0) ? dx : -dx;
xx = y * y;
s = y * xx * (sn3 + xx * sn5);
static const double t22 = 6291456.0;
y = ABS (x);
- u.x = big.x + y;
- y = y - (u.x - big.x);
+ u.x = big + y;
+ y = y - (u.x - big);
dx = (x > 0) ? dx : -dx;
xx = y * y;
s = y * xx * (sn3 + xx * sn5);
static const double t22 = 6291456.0;
y = ABS (x);
- u.x = big.x + y;
- y = y - (u.x - big.x);
+ u.x = big + y;
+ y = y - (u.x - big);
dx = (x > 0) ? dx : -dx;
xx = y * y;
s = y * xx * (sn3 + xx * sn5);
static const double t22 = 6291456.0;
y = ABS (x);
- u.x = big.x + y;
- y = y - (u.x - big.x);
+ u.x = big + y;
+ y = y - (u.x - big);
dx = (x > 0) ? dx : -dx;
xx = y * y;
s = y * xx * (sn3 + xx * sn5);
static const double t22 = 6291456.0;
y = ABS (x);
- u.x = big.x + y;
- y = y - (u.x - big.x);
+ u.x = big + y;
+ y = y - (u.x - big);
xx = y * y;
s = y * xx * (sn3 + xx * sn5);
c = xx * (cs2 + xx * (cs4 + xx * cs6));
return (x > 0) ? w[0] : -w[0];
else
{
- t = (orig * hpinv.x + toint.x);
- xn = t - toint.x;
+ t = (orig * hpinv + toint);
+ xn = t - toint;
v.x = t;
- y = (orig - xn * mp1.x) - xn * mp2.x;
+ y = (orig - xn * mp1) - xn * mp2;
n = v.i[LOW_HALF] & 3;
- da = xn * pp3.x;
+ da = xn * pp3;
t = y - da;
da = (y - t) - da;
- y = xn * pp4.x;
+ y = xn * pp4;
a = t - y;
da = ((t - a) - y) + da;
if (n == 1)
static const double t22 = 6291456.0;
y = ABS (x);
- u.x = big.x + y;
- y = y - (u.x - big.x);
+ u.x = big + y;
+ y = y - (u.x - big);
dx = (x > 0) ? dx : -dx;
xx = y * y;
s = y * xx * (sn3 + xx * sn5);
static const double t22 = 6291456.0;
y = ABS (x);
- u.x = big.x + y;
- y = y - (u.x - big.x);
+ u.x = big + y;
+ y = y - (u.x - big);
dx = (x > 0) ? dx : -dx;
xx = y * y;
s = y * xx * (sn3 + xx * sn5);
#ifndef USNCS_H
#define USNCS_H
-#ifdef BIG_ENDI
-static const mynumber
-
-/**/ s1 = {{0xBFC55555, 0x55555555 }}, /* -0.16666666666666666 */
-/**/ s2 = {{0x3F811111, 0x11110ECE }}, /* 0.0083333333333323288 */
-/**/ s3 = {{0xBF2A01A0, 0x19DB08B8 }}, /* -0.00019841269834414642 */
-/**/ s4 = {{0x3EC71DE2, 0x7B9A7ED9 }}, /* 2.755729806860771e-06 */
-/**/ s5 = {{0xBE5ADDFF, 0xC2FCDF59 }}, /* -2.5022014848318398e-08 */
-/**/ aa = {{0xBFC55580, 0x00000000 }}, /* -0.1666717529296875 */
-/**/ bb = {{0x3ED55555, 0x55556E24 }}, /* 5.0862630208387126e-06 */
-/**/ big = {{0x42c80000, 0x00000000 }}, /* 52776558133248 */
-/**/ hp0 = {{0x3FF921FB, 0x54442D18 }}, /* 1.5707963267948966 */
-/**/ hp1 = {{0x3C91A626, 0x33145C07 }}, /* 6.123233995736766e-17 */
-/**/ mp1 = {{0x3FF921FB, 0x58000000 }}, /* 1.5707963407039642 */
-/**/ mp2 = {{0xBE4DDE97, 0x3C000000 }}, /* -1.3909067564377153e-08 */
-/**/ mp3 = {{0xBC8CB3B3, 0x99D747F2 }}, /* -4.9789962505147994e-17 */
-/**/ pp3 = {{0xBC8CB3B3, 0x98000000 }}, /* -4.9789962314799099e-17 */
-/**/ pp4 = {{0xbacd747f, 0x23e32ed7 }}, /* -1.9034889620193266e-25 */
-/**/ hpinv = {{0x3FE45F30, 0x6DC9C883 }}, /* 0.63661977236758138 */
-/**/ toint = {{0x43380000, 0x00000000 }}; /* 6755399441055744 */
-
-#else
-#ifdef LITTLE_ENDI
-static const mynumber
-
-/**/ s1 = {{0x55555555, 0xBFC55555 }},/* -0.16666666666666666 */
-/**/ s2 = {{0x11110ECE, 0x3F811111 }},/* 0.0083333333333323288 */
-/**/ s3 = {{0x19DB08B8, 0xBF2A01A0 }},/* -0.00019841269834414642 */
-/**/ s4 = {{0x7B9A7ED9, 0x3EC71DE2 }},/* 2.755729806860771e-06 */
-/**/ s5 = {{0xC2FCDF59, 0xBE5ADDFF }},/* -2.5022014848318398e-08 */
-/**/ aa = {{0x00000000, 0xBFC55580 }},/* -0.1666717529296875 */
-/**/ bb = {{0x55556E24, 0x3ED55555 }},/* 5.0862630208387126e-06 */
-/**/ big = {{0x00000000, 0x42c80000 }},/* 52776558133248 */
-/**/ hp0 = {{0x54442D18, 0x3FF921FB }},/* 1.5707963267948966 */
-/**/ hp1 = {{0x33145C07, 0x3C91A626 }},/* 6.123233995736766e-17 */
-/**/ mp1 = {{0x58000000, 0x3FF921FB }},/* 1.5707963407039642 */
-/**/ mp2 = {{0x3C000000, 0xBE4DDE97 }},/* -1.3909067564377153e-08 */
-/**/ mp3 = {{0x99D747F2, 0xBC8CB3B3 }},/* -4.9789962505147994e-17 */
-/**/ pp3 = {{0x98000000, 0xBC8CB3B3 }},/* -4.9789962314799099e-17 */
-/**/ pp4 = {{0x23e32ed7, 0xbacd747f }},/* -1.9034889620193266e-25 */
-/**/ hpinv = {{0x6DC9C883, 0x3FE45F30 }},/* 0.63661977236758138 */
-/**/ toint = {{0x00000000, 0x43380000 }};/* 6755399441055744 */
-
-
-#endif
-#endif
+static const double s1 = -0x1.5555555555555p-3; /* -0.16666666666666666 */
+static const double s2 = 0x1.1111111110ECEp-7; /* 0.0083333333333323288 */
+static const double s3 = -0x1.A01A019DB08B8p-13; /* -0.00019841269834414642 */
+static const double s4 = 0x1.71DE27B9A7ED9p-19; /* 2.755729806860771e-06 */
+static const double s5 = -0x1.ADDFFC2FCDF59p-26; /* -2.5022014848318398e-08 */
+static const double aa = -0x1.5558000000000p-3; /* -0.1666717529296875 */
+static const double bb = 0x1.5555555556E24p-18; /* 5.0862630208387126e-06 */
+static const double big = 0x1.8000000000000p45; /* 52776558133248 */
+static const double hp0 = 0x1.921FB54442D18p0; /* 1.5707963267948966 */
+static const double hp1 = 0x1.1A62633145C07p-54; /* 6.123233995736766e-17 */
+static const double mp1 = 0x1.921FB58000000p0; /* 1.5707963407039642 */
+static const double mp2 = -0x1.DDE973C000000p-27; /* -1.3909067564377153e-08 */
+static const double mp3 = -0x1.CB3B399D747F2p-55; /* -4.9789962505147994e-17 */
+static const double pp3 = -0x1.CB3B398000000p-55; /* -4.9789962314799099e-17 */
+static const double pp4 = -0x1.d747f23e32ed7p-83; /* -1.9034889620193266e-25 */
+static const double hpinv = 0x1.45F306DC9C883p-1; /* 0.63661977236758138 */
+static const double toint = 0x1.8000000000000p52; /* 6755399441055744 */
#endif