|
|
|
@ -52,6 +52,7 @@ |
|
|
|
#include "mydefs.h" |
|
|
|
#include "usncs.h" |
|
|
|
#include "MathLib.h" |
|
|
|
#include <math.h> |
|
|
|
#include <math_private.h> |
|
|
|
#include <fenv.h> |
|
|
|
|
|
|
|
@ -355,7 +356,7 @@ __sin (double x) |
|
|
|
da = xn * mp3; |
|
|
|
a = y - da; |
|
|
|
da = (y - a) - da; |
|
|
|
eps = ABS (x) * 1.2e-30; |
|
|
|
eps = fabs (x) * 1.2e-30; |
|
|
|
|
|
|
|
switch (n) |
|
|
|
{ /* quarter of unit circle */ |
|
|
|
@ -530,7 +531,7 @@ __cos (double x) |
|
|
|
|
|
|
|
else if (k < 0x3feb6000) |
|
|
|
{ /* 2^-27 < |x| < 0.855469 */ |
|
|
|
y = ABS (x); |
|
|
|
y = fabs (x); |
|
|
|
u.x = big + y; |
|
|
|
y = y - (u.x - big); |
|
|
|
res = do_cos (u, y, &cor); |
|
|
|
@ -539,7 +540,7 @@ __cos (double x) |
|
|
|
|
|
|
|
else if (k < 0x400368fd) |
|
|
|
{ /* 0.855469 <|x|<2.426265 */ ; |
|
|
|
y = hp0 - ABS (x); |
|
|
|
y = hp0 - fabs (x); |
|
|
|
a = y + hp1; |
|
|
|
da = (y - a) + hp1; |
|
|
|
xx = a * a; |
|
|
|
@ -582,7 +583,7 @@ __cos (double x) |
|
|
|
da = xn * mp3; |
|
|
|
a = y - da; |
|
|
|
da = (y - a) - da; |
|
|
|
eps = ABS (x) * 1.2e-30; |
|
|
|
eps = fabs (x) * 1.2e-30; |
|
|
|
|
|
|
|
switch (n) |
|
|
|
{ |
|
|
|
@ -741,7 +742,7 @@ slow (double x) |
|
|
|
return res; |
|
|
|
else |
|
|
|
{ |
|
|
|
__dubsin (ABS (x), 0, w); |
|
|
|
__dubsin (fabs (x), 0, w); |
|
|
|
if (w[0] == w[0] + 1.000000001 * w[1]) |
|
|
|
return (x > 0) ? w[0] : -w[0]; |
|
|
|
else |
|
|
|
@ -760,7 +761,7 @@ slow1 (double x) |
|
|
|
{ |
|
|
|
mynumber u; |
|
|
|
double w[2], y, cor, res; |
|
|
|
y = ABS (x); |
|
|
|
y = fabs (x); |
|
|
|
u.x = big + y; |
|
|
|
y = y - (u.x - big); |
|
|
|
res = do_sin_slow (u, y, 0, 0, &cor); |
|
|
|
@ -768,7 +769,7 @@ slow1 (double x) |
|
|
|
return (x > 0) ? res : -res; |
|
|
|
else |
|
|
|
{ |
|
|
|
__dubsin (ABS (x), 0, w); |
|
|
|
__dubsin (fabs (x), 0, w); |
|
|
|
if (w[0] == w[0] + 1.000000005 * w[1]) |
|
|
|
return (x > 0) ? w[0] : -w[0]; |
|
|
|
else |
|
|
|
@ -787,7 +788,7 @@ slow2 (double x) |
|
|
|
mynumber u; |
|
|
|
double w[2], y, y1, y2, cor, res, del; |
|
|
|
|
|
|
|
y = ABS (x); |
|
|
|
y = fabs (x); |
|
|
|
y = hp0 - y; |
|
|
|
if (y >= 0) |
|
|
|
{ |
|
|
|
@ -806,7 +807,7 @@ slow2 (double x) |
|
|
|
return (x > 0) ? res : -res; |
|
|
|
else |
|
|
|
{ |
|
|
|
y = ABS (x) - hp0; |
|
|
|
y = fabs (x) - hp0; |
|
|
|
y1 = y - hp1; |
|
|
|
y2 = (y - y1) - hp1; |
|
|
|
__docos (y1, y2, w); |
|
|
|
@ -834,9 +835,9 @@ sloww (double x, double dx, double orig) |
|
|
|
int4 n; |
|
|
|
res = TAYLOR_SLOW (x, dx, cor); |
|
|
|
if (cor > 0) |
|
|
|
cor = 1.0005 * cor + ABS (orig) * 3.1e-30; |
|
|
|
cor = 1.0005 * cor + fabs (orig) * 3.1e-30; |
|
|
|
else |
|
|
|
cor = 1.0005 * cor - ABS (orig) * 3.1e-30; |
|
|
|
cor = 1.0005 * cor - fabs (orig) * 3.1e-30; |
|
|
|
|
|
|
|
if (res == res + cor) |
|
|
|
return res; |
|
|
|
@ -844,9 +845,9 @@ sloww (double x, double dx, double orig) |
|
|
|
{ |
|
|
|
(x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); |
|
|
|
if (w[1] > 0) |
|
|
|
cor = 1.000000001 * w[1] + ABS (orig) * 1.1e-30; |
|
|
|
cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30; |
|
|
|
else |
|
|
|
cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-30; |
|
|
|
cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30; |
|
|
|
|
|
|
|
if (w[0] == w[0] + cor) |
|
|
|
return (x > 0) ? w[0] : -w[0]; |
|
|
|
@ -870,9 +871,9 @@ sloww (double x, double dx, double orig) |
|
|
|
} |
|
|
|
(a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w); |
|
|
|
if (w[1] > 0) |
|
|
|
cor = 1.000000001 * w[1] + ABS (orig) * 1.1e-40; |
|
|
|
cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40; |
|
|
|
else |
|
|
|
cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-40; |
|
|
|
cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40; |
|
|
|
|
|
|
|
if (w[0] == w[0] + cor) |
|
|
|
return (a > 0) ? w[0] : -w[0]; |
|
|
|
@ -898,7 +899,7 @@ sloww1 (double x, double dx, double orig, int m) |
|
|
|
|
|
|
|
u.x = big + x; |
|
|
|
y = x - (u.x - big); |
|
|
|
res = do_sin_slow (u, y, dx, 3.1e-30 * ABS (orig), &cor); |
|
|
|
res = do_sin_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); |
|
|
|
|
|
|
|
if (res == res + cor) |
|
|
|
return (m > 0) ? res : -res; |
|
|
|
@ -907,9 +908,9 @@ sloww1 (double x, double dx, double orig, int m) |
|
|
|
__dubsin (x, dx, w); |
|
|
|
|
|
|
|
if (w[1] > 0) |
|
|
|
cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig); |
|
|
|
cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); |
|
|
|
else |
|
|
|
cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig); |
|
|
|
cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); |
|
|
|
|
|
|
|
if (w[0] == w[0] + cor) |
|
|
|
return (m > 0) ? w[0] : -w[0]; |
|
|
|
@ -934,7 +935,7 @@ sloww2 (double x, double dx, double orig, int n) |
|
|
|
|
|
|
|
u.x = big + x; |
|
|
|
y = x - (u.x - big); |
|
|
|
res = do_cos_slow (u, y, dx, 3.1e-30 * ABS (orig), &cor); |
|
|
|
res = do_cos_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); |
|
|
|
|
|
|
|
if (res == res + cor) |
|
|
|
return (n & 2) ? -res : res; |
|
|
|
@ -943,9 +944,9 @@ sloww2 (double x, double dx, double orig, int n) |
|
|
|
__docos (x, dx, w); |
|
|
|
|
|
|
|
if (w[1] > 0) |
|
|
|
cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig); |
|
|
|
cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); |
|
|
|
else |
|
|
|
cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig); |
|
|
|
cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); |
|
|
|
|
|
|
|
if (w[0] == w[0] + cor) |
|
|
|
return (n & 2) ? -w[0] : w[0]; |
|
|
|
@ -1000,7 +1001,7 @@ bsloww1 (double x, double dx, double orig, int n) |
|
|
|
mynumber u; |
|
|
|
double w[2], y, cor, res; |
|
|
|
|
|
|
|
y = ABS (x); |
|
|
|
y = fabs (x); |
|
|
|
u.x = big + y; |
|
|
|
y = y - (u.x - big); |
|
|
|
dx = (x > 0) ? dx : -dx; |
|
|
|
@ -1009,7 +1010,7 @@ bsloww1 (double x, double dx, double orig, int n) |
|
|
|
return (x > 0) ? res : -res; |
|
|
|
else |
|
|
|
{ |
|
|
|
__dubsin (ABS (x), dx, w); |
|
|
|
__dubsin (fabs (x), dx, w); |
|
|
|
|
|
|
|
if (w[1] > 0) |
|
|
|
cor = 1.000000005 * w[1] + 1.1e-24; |
|
|
|
@ -1037,7 +1038,7 @@ bsloww2 (double x, double dx, double orig, int n) |
|
|
|
mynumber u; |
|
|
|
double w[2], y, cor, res; |
|
|
|
|
|
|
|
y = ABS (x); |
|
|
|
y = fabs (x); |
|
|
|
u.x = big + y; |
|
|
|
y = y - (u.x - big); |
|
|
|
dx = (x > 0) ? dx : -dx; |
|
|
|
@ -1046,7 +1047,7 @@ bsloww2 (double x, double dx, double orig, int n) |
|
|
|
return (n & 2) ? -res : res; |
|
|
|
else |
|
|
|
{ |
|
|
|
__docos (ABS (x), dx, w); |
|
|
|
__docos (fabs (x), dx, w); |
|
|
|
|
|
|
|
if (w[1] > 0) |
|
|
|
cor = 1.000000005 * w[1] + 1.1e-24; |
|
|
|
@ -1072,7 +1073,7 @@ cslow2 (double x) |
|
|
|
mynumber u; |
|
|
|
double w[2], y, cor, res; |
|
|
|
|
|
|
|
y = ABS (x); |
|
|
|
y = fabs (x); |
|
|
|
u.x = big + y; |
|
|
|
y = y - (u.x - big); |
|
|
|
res = do_cos_slow (u, y, 0, 0, &cor); |
|
|
|
@ -1080,7 +1081,7 @@ cslow2 (double x) |
|
|
|
return res; |
|
|
|
else |
|
|
|
{ |
|
|
|
y = ABS (x); |
|
|
|
y = fabs (x); |
|
|
|
__docos (y, 0, w); |
|
|
|
if (w[0] == w[0] + 1.000000005 * w[1]) |
|
|
|
return w[0]; |
|
|
|
@ -1109,9 +1110,9 @@ csloww (double x, double dx, double orig) |
|
|
|
res = TAYLOR_SLOW (x, dx, cor); |
|
|
|
|
|
|
|
if (cor > 0) |
|
|
|
cor = 1.0005 * cor + ABS (orig) * 3.1e-30; |
|
|
|
cor = 1.0005 * cor + fabs (orig) * 3.1e-30; |
|
|
|
else |
|
|
|
cor = 1.0005 * cor - ABS (orig) * 3.1e-30; |
|
|
|
cor = 1.0005 * cor - fabs (orig) * 3.1e-30; |
|
|
|
|
|
|
|
if (res == res + cor) |
|
|
|
return res; |
|
|
|
@ -1120,9 +1121,9 @@ csloww (double x, double dx, double orig) |
|
|
|
(x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); |
|
|
|
|
|
|
|
if (w[1] > 0) |
|
|
|
cor = 1.000000001 * w[1] + ABS (orig) * 1.1e-30; |
|
|
|
cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30; |
|
|
|
else |
|
|
|
cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-30; |
|
|
|
cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30; |
|
|
|
|
|
|
|
if (w[0] == w[0] + cor) |
|
|
|
return (x > 0) ? w[0] : -w[0]; |
|
|
|
@ -1147,9 +1148,9 @@ csloww (double x, double dx, double orig) |
|
|
|
(a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w); |
|
|
|
|
|
|
|
if (w[1] > 0) |
|
|
|
cor = 1.000000001 * w[1] + ABS (orig) * 1.1e-40; |
|
|
|
cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40; |
|
|
|
else |
|
|
|
cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-40; |
|
|
|
cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40; |
|
|
|
|
|
|
|
if (w[0] == w[0] + cor) |
|
|
|
return (a > 0) ? w[0] : -w[0]; |
|
|
|
@ -1175,7 +1176,7 @@ csloww1 (double x, double dx, double orig, int m) |
|
|
|
|
|
|
|
u.x = big + x; |
|
|
|
y = x - (u.x - big); |
|
|
|
res = do_sin_slow (u, y, dx, 3.1e-30 * ABS (orig), &cor); |
|
|
|
res = do_sin_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); |
|
|
|
|
|
|
|
if (res == res + cor) |
|
|
|
return (m > 0) ? res : -res; |
|
|
|
@ -1183,9 +1184,9 @@ csloww1 (double x, double dx, double orig, int m) |
|
|
|
{ |
|
|
|
__dubsin (x, dx, w); |
|
|
|
if (w[1] > 0) |
|
|
|
cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig); |
|
|
|
cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); |
|
|
|
else |
|
|
|
cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig); |
|
|
|
cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); |
|
|
|
if (w[0] == w[0] + cor) |
|
|
|
return (m > 0) ? w[0] : -w[0]; |
|
|
|
else |
|
|
|
@ -1210,7 +1211,7 @@ csloww2 (double x, double dx, double orig, int n) |
|
|
|
|
|
|
|
u.x = big + x; |
|
|
|
y = x - (u.x - big); |
|
|
|
res = do_cos_slow (u, y, dx, 3.1e-30 * ABS (orig), &cor); |
|
|
|
res = do_cos_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); |
|
|
|
|
|
|
|
if (res == res + cor) |
|
|
|
return (n) ? -res : res; |
|
|
|
@ -1218,9 +1219,9 @@ csloww2 (double x, double dx, double orig, int n) |
|
|
|
{ |
|
|
|
__docos (x, dx, w); |
|
|
|
if (w[1] > 0) |
|
|
|
cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig); |
|
|
|
cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); |
|
|
|
else |
|
|
|
cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig); |
|
|
|
cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); |
|
|
|
if (w[0] == w[0] + cor) |
|
|
|
return (n) ? -w[0] : w[0]; |
|
|
|
else |
|
|
|
|