diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/s_sin.c')
-rw-r--r-- | sysdeps/ieee754/dbl-64/s_sin.c | 883 |
1 files changed, 331 insertions, 552 deletions
diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c index eff120e88d..ca2532fb63 100644 --- a/sysdeps/ieee754/dbl-64/s_sin.c +++ b/sysdeps/ieee754/dbl-64/s_sin.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2015 Free Software Foundation, Inc. + * Copyright (C) 2001-2016 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -32,9 +32,6 @@ /* bsloww1 */ /* bsloww2 */ /* cslow2 */ -/* csloww */ -/* csloww1 */ -/* csloww2 */ /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h usncs.h */ /* branred.c sincos32.c dosincos.c mpa.c */ /* sincos.tbl */ @@ -135,17 +132,14 @@ double __mpcos (double x, double dx, bool reduce_range); static double slow (double x); static double slow1 (double x); static double slow2 (double x); -static double sloww (double x, double dx, double orig); -static double sloww1 (double x, double dx, double orig, int m); +static double sloww (double x, double dx, double orig, int n); +static double sloww1 (double x, double dx, double orig, int m, int n); static double sloww2 (double x, double dx, double orig, int n); static double bsloww (double x, double dx, double orig, int n); static double bsloww1 (double x, double dx, double orig, int n); static double bsloww2 (double x, double dx, double orig, int n); int __branred (double x, double *a, double *aa); static double cslow2 (double x); -static double csloww (double x, double dx, double orig); -static double csloww1 (double x, double dx, double orig, int m); -static double csloww2 (double x, double dx, double orig, int n); /* Given a number partitioned into U and X such that U is an index into the sin/cos table, this macro computes the cosine of the number by combining @@ -276,32 +270,216 @@ reduce_and_compute (double x, unsigned int k) return retval; } +static inline int4 +__always_inline +reduce_sincos_1 (double x, double *a, double *da) +{ + mynumber v; + + double t = (x * hpinv + toint); + double xn = t - toint; + v.x = t; + double y = (x - xn * mp1) - xn * mp2; + int4 n = v.i[LOW_HALF] & 3; + double db = xn * mp3; + double b = y - db; + db = (y - b) - db; + + *a = b; + *da = db; + + return n; +} + +/* Compute sin (A + DA). cos can be computed by shifting the quadrant N + clockwise. */ +static double +__always_inline +do_sincos_1 (double a, double da, double x, int4 n, int4 k) +{ + double xx, retval, res, cor, y; + mynumber u; + int m; + double eps = fabs (x) * 1.2e-30; + + int k1 = (n + k) & 3; + switch (k1) + { /* quarter of unit circle */ + case 2: + a = -a; + da = -da; + case 0: + xx = a * a; + if (xx < 0.01588) + { + /* Taylor series. */ + res = TAYLOR_SIN (xx, a, da, cor); + cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; + retval = (res == res + cor) ? res : sloww (a, da, x, k); + } + else + { + if (a > 0) + m = 1; + else + { + m = 0; + a = -a; + da = -da; + } + u.x = big + a; + y = a - (u.x - big); + res = do_sin (u, y, da, &cor); + cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; + retval = ((res == res + cor) ? ((m) ? res : -res) + : sloww1 (a, da, x, m, k)); + } + break; + + case 1: + case 3: + if (a < 0) + { + a = -a; + da = -da; + } + u.x = big + a; + y = a - (u.x - big) + da; + res = do_cos (u, y, &cor); + cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; + retval = ((res == res + cor) ? ((k1 & 2) ? -res : res) + : sloww2 (a, da, x, n)); + break; + } + + return retval; +} + +static inline int4 +__always_inline +reduce_sincos_2 (double x, double *a, double *da) +{ + mynumber v; + + double t = (x * hpinv + toint); + double xn = t - toint; + v.x = t; + double xn1 = (xn + 8.0e22) - 8.0e22; + double xn2 = xn - xn1; + double y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2); + int4 n = v.i[LOW_HALF] & 3; + double db = xn1 * pp3; + t = y - db; + db = (y - t) - db; + db = (db - xn2 * pp3) - xn * pp4; + double b = t + db; + db = (t - b) + db; + + *a = b; + *da = db; + + return n; +} + +/* Compute sin (A + DA). cos can be computed by shifting the quadrant N + clockwise. */ +static double +__always_inline +do_sincos_2 (double a, double da, double x, int4 n, int4 k) +{ + double res, retval, cor, xx; + mynumber u; + + double eps = 1.0e-24; + + k = (n + k) & 3; + + switch (k) + { + case 2: + a = -a; + da = -da; + /* Fall through. */ + case 0: + xx = a * a; + if (xx < 0.01588) + { + /* Taylor series. */ + res = TAYLOR_SIN (xx, a, da, cor); + cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; + retval = (res == res + cor) ? res : bsloww (a, da, x, n); + } + else + { + double t, db, y; + int m; + if (a > 0) + { + m = 1; + t = a; + db = da; + } + else + { + m = 0; + t = -a; + db = -da; + } + u.x = big + t; + y = t - (u.x - big); + res = do_sin (u, y, db, &cor); + cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; + retval = ((res == res + cor) ? ((m) ? res : -res) + : bsloww1 (a, da, x, n)); + } + break; + + case 1: + case 3: + if (a < 0) + { + a = -a; + da = -da; + } + u.x = big + a; + double y = a - (u.x - big) + da; + res = do_cos (u, y, &cor); + cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; + retval = ((res == res + cor) ? ((n & 2) ? -res : res) + : bsloww2 (a, da, x, n)); + break; + } + + return retval; +} + /*******************************************************************/ /* An ultimate sin routine. Given an IEEE double machine number x */ /* it computes the correctly rounded (to nearest) value of sin(x) */ /*******************************************************************/ +#ifdef IN_SINCOS +static double +#else double SECTION +#endif __sin (double x) { - double xx, res, t, cor, y, s, c, sn, ssn, cs, ccs, xn, a, da, db, eps, xn1, - xn2; - mynumber u, v; - int4 k, m, n; + double xx, res, t, cor, y, s, c, sn, ssn, cs, ccs; + mynumber u; + int4 k, m; double retval = 0; +#ifndef IN_SINCOS SET_RESTORE_ROUND_53BIT (FE_TONEAREST); +#endif u.x = x; m = u.i[HIGH_HALF]; k = 0x7fffffff & m; /* no sign */ if (k < 0x3e500000) /* if x->0 =>sin(x)=x */ { - if (fabs (x) < DBL_MIN) - { - double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); retval = x; } /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/ @@ -353,146 +531,22 @@ __sin (double x) retval = (res == res + 1.020 * cor) ? ((m > 0) ? res : -res) : slow2 (x); } /* else if (k < 0x400368fd) */ +#ifndef IN_SINCOS /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/ else if (k < 0x419921FB) { - t = (x * hpinv + toint); - xn = t - toint; - v.x = t; - y = (x - xn * mp1) - xn * mp2; - n = v.i[LOW_HALF] & 3; - da = xn * mp3; - a = y - da; - da = (y - a) - da; - eps = fabs (x) * 1.2e-30; - - switch (n) - { /* quarter of unit circle */ - case 0: - case 2: - xx = a * a; - if (n) - { - a = -a; - da = -da; - } - if (xx < 0.01588) - { - /* Taylor series. */ - res = TAYLOR_SIN (xx, a, da, cor); - cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; - retval = (res == res + cor) ? res : sloww (a, da, x); - } - else - { - if (a > 0) - m = 1; - else - { - m = 0; - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big); - res = do_sin (u, y, da, &cor); - cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; - retval = ((res == res + cor) ? ((m) ? res : -res) - : sloww1 (a, da, x, m)); - } - break; - - case 1: - case 3: - if (a < 0) - { - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big) + da; - res = do_cos (u, y, &cor); - cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; - retval = ((res == res + cor) ? ((n & 2) ? -res : res) - : sloww2 (a, da, x, n)); - break; - } + double a, da; + int4 n = reduce_sincos_1 (x, &a, &da); + retval = do_sincos_1 (a, da, x, n, 0); } /* else if (k < 0x419921FB ) */ /*---------------------105414350 <|x|< 281474976710656 --------------------*/ else if (k < 0x42F00000) { - t = (x * hpinv + toint); - xn = t - toint; - v.x = t; - xn1 = (xn + 8.0e22) - 8.0e22; - xn2 = xn - xn1; - y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2); - n = v.i[LOW_HALF] & 3; - da = xn1 * pp3; - t = y - da; - da = (y - t) - da; - da = (da - xn2 * pp3) - xn * pp4; - a = t + da; - da = (t - a) + da; - eps = 1.0e-24; - - switch (n) - { - case 0: - case 2: - xx = a * a; - if (n) - { - a = -a; - da = -da; - } - if (xx < 0.01588) - { - /* Taylor series. */ - res = TAYLOR_SIN (xx, a, da, cor); - cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; - retval = (res == res + cor) ? res : bsloww (a, da, x, n); - } - else - { - double t; - if (a > 0) - { - m = 1; - t = a; - db = da; - } - else - { - m = 0; - t = -a; - db = -da; - } - u.x = big + t; - y = t - (u.x - big); - res = do_sin (u, y, db, &cor); - cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; - retval = ((res == res + cor) ? ((m) ? res : -res) - : bsloww1 (a, da, x, n)); - } - break; + double a, da; - case 1: - case 3: - if (a < 0) - { - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big) + da; - res = do_cos (u, y, &cor); - cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; - retval = ((res == res + cor) ? ((n & 2) ? -res : res) - : bsloww2 (a, da, x, n)); - break; - } + int4 n = reduce_sincos_2 (x, &a, &da); + retval = do_sincos_2 (a, da, x, n, 0); } /* else if (k < 0x42F00000 ) */ /* -----------------281474976710656 <|x| <2^1024----------------------------*/ @@ -506,6 +560,7 @@ __sin (double x) __set_errno (EDOM); retval = x / x; } +#endif return retval; } @@ -516,18 +571,23 @@ __sin (double x) /* it computes the correctly rounded (to nearest) value of cos(x) */ /*******************************************************************/ +#ifdef IN_SINCOS +static double +#else double SECTION +#endif __cos (double x) { - double y, xx, res, t, cor, xn, a, da, db, eps, xn1, - xn2; - mynumber u, v; - int4 k, m, n; + double y, xx, res, cor, a, da; + mynumber u; + int4 k, m; double retval = 0; +#ifndef IN_SINCOS SET_RESTORE_ROUND_53BIT (FE_TONEAREST); +#endif u.x = x; m = u.i[HIGH_HALF]; @@ -556,7 +616,7 @@ __cos (double x) { res = TAYLOR_SIN (xx, a, da, cor); cor = (cor > 0) ? 1.02 * cor + 1.0e-31 : 1.02 * cor - 1.0e-31; - retval = (res == res + cor) ? res : csloww (a, da, x); + retval = (res == res + cor) ? res : sloww (a, da, x, 1); } else { @@ -575,150 +635,26 @@ __cos (double x) res = do_sin (u, y, da, &cor); cor = (cor > 0) ? 1.035 * cor + 1.0e-31 : 1.035 * cor - 1.0e-31; retval = ((res == res + cor) ? ((m) ? res : -res) - : csloww1 (a, da, x, m)); + : sloww1 (a, da, x, m, 1)); } } /* else if (k < 0x400368fd) */ +#ifndef IN_SINCOS else if (k < 0x419921FB) { /* 2.426265<|x|< 105414350 */ - t = (x * hpinv + toint); - xn = t - toint; - v.x = t; - y = (x - xn * mp1) - xn * mp2; - n = v.i[LOW_HALF] & 3; - da = xn * mp3; - a = y - da; - da = (y - a) - da; - eps = fabs (x) * 1.2e-30; - - switch (n) - { - case 1: - case 3: - xx = a * a; - if (n == 1) - { - a = -a; - da = -da; - } - if (xx < 0.01588) - { - res = TAYLOR_SIN (xx, a, da, cor); - cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; - retval = (res == res + cor) ? res : csloww (a, da, x); - } - else - { - if (a > 0) - { - m = 1; - } - else - { - m = 0; - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big); - res = do_sin (u, y, da, &cor); - cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; - retval = ((res == res + cor) ? ((m) ? res : -res) - : csloww1 (a, da, x, m)); - } - break; - - case 0: - case 2: - if (a < 0) - { - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big) + da; - res = do_cos (u, y, &cor); - cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; - retval = ((res == res + cor) ? ((n) ? -res : res) - : csloww2 (a, da, x, n)); - break; - } + double a, da; + int4 n = reduce_sincos_1 (x, &a, &da); + retval = do_sincos_1 (a, da, x, n, 1); } /* else if (k < 0x419921FB ) */ else if (k < 0x42F00000) { - t = (x * hpinv + toint); - xn = t - toint; - v.x = t; - xn1 = (xn + 8.0e22) - 8.0e22; - xn2 = xn - xn1; - y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2); - n = v.i[LOW_HALF] & 3; - da = xn1 * pp3; - t = y - da; - da = (y - t) - da; - da = (da - xn2 * pp3) - xn * pp4; - a = t + da; - da = (t - a) + da; - eps = 1.0e-24; - - switch (n) - { - case 1: - case 3: - xx = a * a; - if (n == 1) - { - a = -a; - da = -da; - } - if (xx < 0.01588) - { - res = TAYLOR_SIN (xx, a, da, cor); - cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; - retval = (res == res + cor) ? res : bsloww (a, da, x, n); - } - else - { - double t; - if (a > 0) - { - m = 1; - t = a; - db = da; - } - else - { - m = 0; - t = -a; - db = -da; - } - u.x = big + t; - y = t - (u.x - big); - res = do_sin (u, y, db, &cor); - cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; - retval = ((res == res + cor) ? ((m) ? res : -res) - : bsloww1 (a, da, x, n)); - } - break; + double a, da; - case 0: - case 2: - if (a < 0) - { - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big) + da; - res = do_cos (u, y, &cor); - cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; - retval = ((res == res + cor) ? ((n) ? -res : res) - : bsloww2 (a, da, x, n)); - break; - } + int4 n = reduce_sincos_2 (x, &a, &da); + retval = do_sincos_2 (a, da, x, n, 1); } /* else if (k < 0x42F00000 ) */ /* 281474976710656 <|x| <2^1024 */ @@ -731,6 +667,7 @@ __cos (double x) __set_errno (EDOM); retval = x / x; /* |x| > 2^1024 */ } +#endif return retval; } @@ -748,14 +685,12 @@ slow (double x) res = TAYLOR_SLOW (x, 0, cor); if (res == res + 1.0007 * cor) return res; - else - { - __dubsin (fabs (x), 0, w); - if (w[0] == w[0] + 1.000000001 * w[1]) - return (x > 0) ? w[0] : -w[0]; - else - return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); - } + + __dubsin (fabs (x), 0, w); + if (w[0] == w[0] + 1.000000001 * w[1]) + return (x > 0) ? w[0] : -w[0]; + + return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); } /*******************************************************************************/ @@ -775,14 +710,12 @@ slow1 (double x) res = do_sin_slow (u, y, 0, 0, &cor); if (res == res + cor) return (x > 0) ? res : -res; - else - { - __dubsin (fabs (x), 0, w); - if (w[0] == w[0] + 1.000000005 * w[1]) - return (x > 0) ? w[0] : -w[0]; - else - return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); - } + + __dubsin (fabs (x), 0, w); + if (w[0] == w[0] + 1.000000005 * w[1]) + return (x > 0) ? w[0] : -w[0]; + + return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); } /**************************************************************************/ @@ -813,17 +746,15 @@ slow2 (double x) res = do_cos_slow (u, y, del, 0, &cor); if (res == res + cor) return (x > 0) ? res : -res; - else - { - y = fabs (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]; - else - return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); - } + + y = fabs (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) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); } /***************************************************************************/ @@ -836,12 +767,13 @@ slow2 (double x) static double SECTION -sloww (double x, double dx, double orig) +sloww (double x, double dx, double orig, int k) { double y, t, res, cor, w[2], a, da, xn; mynumber v; int4 n; res = TAYLOR_SLOW (x, dx, cor); + if (cor > 0) cor = 1.0005 * cor + fabs (orig) * 3.1e-30; else @@ -849,46 +781,43 @@ sloww (double x, double dx, double orig) if (res == res + cor) return res; + + (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); + if (w[1] > 0) + cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30; else + cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30; + + if (w[0] == w[0] + cor) + return (x > 0) ? w[0] : -w[0]; + + t = (orig * hpinv + toint); + xn = t - toint; + v.x = t; + y = (orig - xn * mp1) - xn * mp2; + n = (v.i[LOW_HALF] + k) & 3; + da = xn * pp3; + t = y - da; + da = (y - t) - da; + y = xn * pp4; + a = t - y; + da = ((t - a) - y) + da; + + if (n == 2 || n == 1) { - (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); - if (w[1] > 0) - cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30; - else - cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30; + a = -a; + da = -da; + } + (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w); + if (w[1] > 0) + cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40; + else + cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40; - if (w[0] == w[0] + cor) - return (x > 0) ? w[0] : -w[0]; - else - { - t = (orig * hpinv + toint); - xn = t - toint; - v.x = t; - y = (orig - xn * mp1) - xn * mp2; - n = v.i[LOW_HALF] & 3; - da = xn * pp3; - t = y - da; - da = (y - t) - da; - y = xn * pp4; - a = t - y; - da = ((t - a) - y) + da; - if (n & 2) - { - a = -a; - da = -da; - } - (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w); - if (w[1] > 0) - cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40; - else - cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40; + if (w[0] == w[0] + cor) + return (a > 0) ? w[0] : -w[0]; - if (w[0] == w[0] + cor) - return (a > 0) ? w[0] : -w[0]; - else - return __mpsin (orig, 0, true); - } - } + return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); } /***************************************************************************/ @@ -900,7 +829,7 @@ sloww (double x, double dx, double orig) static double SECTION -sloww1 (double x, double dx, double orig, int m) +sloww1 (double x, double dx, double orig, int m, int k) { mynumber u; double w[2], y, cor, res; @@ -911,20 +840,18 @@ sloww1 (double x, double dx, double orig, int m) if (res == res + cor) return (m > 0) ? res : -res; + + __dubsin (x, dx, w); + + if (w[1] > 0) + cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); else - { - __dubsin (x, dx, w); + cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); - if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); - else - cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); + if (w[0] == w[0] + cor) + return (m > 0) ? w[0] : -w[0]; - if (w[0] == w[0] + cor) - return (m > 0) ? w[0] : -w[0]; - else - return __mpsin (orig, 0, true); - } + return (k == 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); } /***************************************************************************/ @@ -947,20 +874,18 @@ sloww2 (double x, double dx, double orig, int n) if (res == res + cor) return (n & 2) ? -res : res; + + __docos (x, dx, w); + + if (w[1] > 0) + cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); else - { - __docos (x, dx, w); + cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); - if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); - else - cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); + if (w[0] == w[0] + cor) + return (n & 2) ? -w[0] : w[0]; - if (w[0] == w[0] + cor) - return (n & 2) ? -w[0] : w[0]; - else - return __mpsin (orig, 0, true); - } + return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true); } /***************************************************************************/ @@ -981,18 +906,17 @@ bsloww (double x, double dx, double orig, int n) cor = (cor > 0) ? 1.0005 * cor + 1.1e-24 : 1.0005 * cor - 1.1e-24; if (res == res + cor) return res; + + (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); + if (w[1] > 0) + cor = 1.000000001 * w[1] + 1.1e-24; else - { - (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); - if (w[1] > 0) - cor = 1.000000001 * w[1] + 1.1e-24; - else - cor = 1.000000001 * w[1] - 1.1e-24; - if (w[0] == w[0] + cor) - return (x > 0) ? w[0] : -w[0]; - else - return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); - } + cor = 1.000000001 * w[1] - 1.1e-24; + + if (w[0] == w[0] + cor) + return (x > 0) ? w[0] : -w[0]; + + return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); } /***************************************************************************/ @@ -1016,20 +940,18 @@ bsloww1 (double x, double dx, double orig, int n) res = do_sin_slow (u, y, dx, 1.1e-24, &cor); if (res == res + cor) return (x > 0) ? res : -res; + + __dubsin (fabs (x), dx, w); + + if (w[1] > 0) + cor = 1.000000005 * w[1] + 1.1e-24; else - { - __dubsin (fabs (x), dx, w); + cor = 1.000000005 * w[1] - 1.1e-24; - if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-24; - else - cor = 1.000000005 * w[1] - 1.1e-24; + if (w[0] == w[0] + cor) + return (x > 0) ? w[0] : -w[0]; - if (w[0] == w[0] + cor) - return (x > 0) ? w[0] : -w[0]; - else - return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); - } + return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); } /***************************************************************************/ @@ -1053,20 +975,18 @@ bsloww2 (double x, double dx, double orig, int n) res = do_cos_slow (u, y, dx, 1.1e-24, &cor); if (res == res + cor) return (n & 2) ? -res : res; + + __docos (fabs (x), dx, w); + + if (w[1] > 0) + cor = 1.000000005 * w[1] + 1.1e-24; else - { - __docos (fabs (x), dx, w); + cor = 1.000000005 * w[1] - 1.1e-24; - if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-24; - else - cor = 1.000000005 * w[1] - 1.1e-24; + if (w[0] == w[0] + cor) + return (n & 2) ? -w[0] : w[0]; - if (w[0] == w[0] + cor) - return (n & 2) ? -w[0] : w[0]; - else - return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true); - } + return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true); } /************************************************************************/ @@ -1087,154 +1007,13 @@ cslow2 (double x) res = do_cos_slow (u, y, 0, 0, &cor); if (res == res + cor) return res; - else - { - y = fabs (x); - __docos (y, 0, w); - if (w[0] == w[0] + 1.000000005 * w[1]) - return w[0]; - else - return __mpcos (x, 0, false); - } -} - -/***************************************************************************/ -/* Routine compute cos(x+dx) (Double-Length number) where x is small enough*/ -/* to use Taylor series around zero and (x+dx) .Routine receive also */ -/* (right argument) the original value of x for computing error of */ -/* result.And if result not accurate enough routine calls other routines */ -/***************************************************************************/ - - -static double -SECTION -csloww (double x, double dx, double orig) -{ - double y, t, res, cor, w[2], a, da, xn; - mynumber v; - int4 n; - - /* Taylor series */ - res = TAYLOR_SLOW (x, dx, cor); - - if (cor > 0) - cor = 1.0005 * cor + fabs (orig) * 3.1e-30; - else - cor = 1.0005 * cor - fabs (orig) * 3.1e-30; - - if (res == res + cor) - return res; - else - { - (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); - - if (w[1] > 0) - cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30; - else - cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30; - - if (w[0] == w[0] + cor) - return (x > 0) ? w[0] : -w[0]; - else - { - t = (orig * hpinv + toint); - xn = t - toint; - v.x = t; - y = (orig - xn * mp1) - xn * mp2; - n = v.i[LOW_HALF] & 3; - da = xn * pp3; - t = y - da; - da = (y - t) - da; - y = xn * pp4; - a = t - y; - da = ((t - a) - y) + da; - if (n == 1) - { - a = -a; - da = -da; - } - (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w); - - if (w[1] > 0) - cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40; - else - cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40; - - if (w[0] == w[0] + cor) - return (a > 0) ? w[0] : -w[0]; - else - return __mpcos (orig, 0, true); - } - } -} - -/***************************************************************************/ -/* Routine compute sin(x+dx) (Double-Length number) where x in first or */ -/* third quarter of unit circle.Routine receive also (right argument) the */ -/* original value of x for computing error of result.And if result not */ -/* accurate enough routine calls other routines */ -/***************************************************************************/ -static double -SECTION -csloww1 (double x, double dx, double orig, int m) -{ - mynumber u; - double w[2], y, cor, res; - - u.x = big + x; - y = x - (u.x - big); - res = do_sin_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); - - if (res == res + cor) - return (m > 0) ? res : -res; - else - { - __dubsin (x, dx, w); - if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); - else - cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); - if (w[0] == w[0] + cor) - return (m > 0) ? w[0] : -w[0]; - else - return __mpcos (orig, 0, true); - } -} - - -/***************************************************************************/ -/* Routine compute sin(x+dx) (Double-Length number) where x in second or */ -/* fourth quarter of unit circle.Routine receive also the original value */ -/* and quarter(n= 1or 3)of x for computing error of result.And if result not*/ -/* accurate enough routine calls other routines */ -/***************************************************************************/ - -static double -SECTION -csloww2 (double x, double dx, double orig, int n) -{ - mynumber u; - double w[2], y, cor, res; - - u.x = big + x; - y = x - (u.x - big); - res = do_cos_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); + y = fabs (x); + __docos (y, 0, w); + if (w[0] == w[0] + 1.000000005 * w[1]) + return w[0]; - if (res == res + cor) - return (n) ? -res : res; - else - { - __docos (x, dx, w); - if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); - else - cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); - if (w[0] == w[0] + cor) - return (n) ? -w[0] : w[0]; - else - return __mpcos (orig, 0, true); - } + return __mpcos (x, 0, false); } #ifndef __cos |