diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/s_sin.c')
-rw-r--r-- | sysdeps/ieee754/dbl-64/s_sin.c | 904 |
1 files changed, 102 insertions, 802 deletions
diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c index ca2532fb63..b369ac9f5b 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-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 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 @@ -22,22 +22,11 @@ /* */ /* FUNCTIONS: usin */ /* ucos */ -/* slow */ -/* slow1 */ -/* slow2 */ -/* sloww */ -/* sloww1 */ -/* sloww2 */ -/* bsloww */ -/* bsloww1 */ -/* bsloww2 */ -/* cslow2 */ /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h usncs.h */ -/* branred.c sincos32.c dosincos.c mpa.c */ -/* sincos.tbl */ +/* branred.c sincos.tbl */ /* */ -/* An ultimate sin and routine. Given an IEEE double machine number x */ -/* it computes the correctly rounded (to nearest) value of sin(x) or cos(x) */ +/* An ultimate sin and cos routine. Given an IEEE double machine number x */ +/* it computes sin(x) or cos(x) with ~0.55 ULP. */ /* Assumption: Machine arithmetic operations are performed in */ /* round to nearest mode of IEEE 754 standard. */ /* */ @@ -52,6 +41,8 @@ #include "MathLib.h" #include <math.h> #include <math_private.h> +#include <math-underflow.h> +#include <libm-alias-double.h> #include <fenv.h> /* Helper macros to compute sin of the input values. */ @@ -65,35 +56,11 @@ a - a^3/3! + a^5/5! - a^7/7! + a^9/9! + (1 - a^2) * da / 2 The constants s1, s2, s3, etc. are pre-computed values of 1/3!, 1/5! and so - on. The result is returned to LHS and correction in COR. */ -#define TAYLOR_SIN(xx, a, da, cor) \ + on. The result is returned to LHS. */ +#define TAYLOR_SIN(xx, a, da) \ ({ \ double t = ((POLYNOMIAL (xx) * (a) - 0.5 * (da)) * (xx) + (da)); \ double res = (a) + t; \ - (cor) = ((a) - res) + t; \ - res; \ -}) - -/* This is again a variation of the Taylor series expansion with the term - x^3/3! expanded into the following for better accuracy: - - bb * x ^ 3 + 3 * aa * x * x1 * x2 + aa * x1 ^ 3 + aa * x2 ^ 3 - - The correction term is dx and bb + aa = -1/3! - */ -#define TAYLOR_SLOW(x0, dx, cor) \ -({ \ - 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 * x1 * x1 * x1; \ - double r = (x0) + y; \ - double x2 = ((x0) - x1) + (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; \ res; \ }) @@ -123,156 +90,69 @@ static const double cs4 = -4.16666666666664434524222570944589E-02, cs6 = 1.38888874007937613028114285595617E-03; -static const double t22 = 0x1.8p22; - -void __dubsin (double x, double dx, double w[]); -void __docos (double x, double dx, double w[]); -double __mpsin (double x, double dx, bool reduce_range); -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, 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); -/* 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 - the sin and cos of X (as computed by a variation of the Taylor series) with - the values looked up from the sin/cos table to get the result in RES and a - correction value in COR. */ -static double -do_cos (mynumber u, double x, double *corp) +/* Given a number partitioned into X and DX, this function computes the cosine + of the number by combining the sin and cos of X (as computed by a variation + of the Taylor series) with the values looked up from the sin/cos table to + get the result. */ +static inline double +__always_inline +do_cos (double x, double dx) { - double xx, s, sn, ssn, c, cs, ccs, res, cor; + mynumber u; + + if (x < 0) + dx = -dx; + + u.x = big + fabs (x); + x = fabs (x) - (u.x - big) + dx; + + double xx, s, sn, ssn, c, cs, ccs, cor; xx = x * x; s = x + x * xx * (sn3 + xx * sn5); c = xx * (cs2 + xx * (cs4 + xx * cs6)); SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); cor = (ccs - s * ssn - cs * c) - sn * s; - res = cs + cor; - cor = (cs - res) + cor; - *corp = cor; - return res; + return cs + cor; } -/* A more precise variant of DO_COS where the number is partitioned into U, X - and DX. EPS is the adjustment to the correction COR. */ -static double -do_cos_slow (mynumber u, double x, double dx, double eps, double *corp) +/* Given a number partitioned into X and DX, this function computes the sine of + the number by combining the sin and cos of X (as computed by a variation of + the Taylor series) with the values looked up from the sin/cos table to get + the result. */ +static inline double +__always_inline +do_sin (double x, double dx) { - double xx, y, x1, x2, e1, e2, res, cor; - double s, sn, ssn, c, cs, ccs; - xx = x * x; - s = x * xx * (sn3 + xx * sn5); - c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6)); - SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); - x1 = (x + t22) - t22; - x2 = (x - x1) + dx; - e1 = (sn + t22) - t22; - e2 = (sn - e1) + ssn; - cor = (ccs - cs * c - e1 * x2 - e2 * x) - sn * s; - y = cs - e1 * x1; - cor = cor + ((cs - y) - e1 * x1); - res = y + cor; - cor = (y - res) + cor; - if (cor > 0) - cor = 1.0005 * cor + eps; - else - cor = 1.0005 * cor - eps; - *corp = cor; - return res; -} + double xold = x; + /* Max ULP is 0.501 if |x| < 0.126, otherwise ULP is 0.518. */ + if (fabs (x) < 0.126) + return TAYLOR_SIN (x * x, x, dx); -/* Given a number partitioned into U and X and DX such that U is an index into - the sin/cos table, this macro computes the sine of the number by combining - the sin and cos of X (as computed by a variation of the Taylor series) with - the values looked up from the sin/cos table to get the result in RES and a - correction value in COR. */ -static double -do_sin (mynumber u, double x, double dx, double *corp) -{ - double xx, s, sn, ssn, c, cs, ccs, cor, res; + mynumber u; + + if (x <= 0) + dx = -dx; + u.x = big + fabs (x); + x = fabs (x) - (u.x - big); + + double xx, s, sn, ssn, c, cs, ccs, cor; xx = x * x; s = x + (dx + x * xx * (sn3 + xx * sn5)); c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6)); SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); cor = (ssn + s * ccs - sn * c) + cs * s; - res = sn + cor; - cor = (sn - res) + cor; - *corp = cor; - return res; -} - -/* A more precise variant of res = do_sin where the number is partitioned into U, X - and DX. EPS is the adjustment to the correction COR. */ -static double -do_sin_slow (mynumber u, double x, double dx, double eps, double *corp) -{ - double xx, y, x1, x2, c1, c2, res, cor; - double s, sn, ssn, c, cs, ccs; - xx = x * x; - s = x * xx * (sn3 + xx * sn5); - c = xx * (cs2 + xx * (cs4 + xx * cs6)); - SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); - x1 = (x + t22) - t22; - x2 = (x - x1) + dx; - c1 = (cs + t22) - t22; - c2 = (cs - c1) + ccs; - cor = (ssn + s * ccs + cs * s + c2 * x + c1 * x2 - sn * x * dx) - sn * c; - y = sn + c1 * x1; - cor = cor + ((sn - y) + c1 * x1); - res = y + cor; - cor = (y - res) + cor; - if (cor > 0) - cor = 1.0005 * cor + eps; - else - cor = 1.0005 * cor - eps; - *corp = cor; - return res; -} - -/* Reduce range of X and compute sin of a + da. K is the amount by which to - rotate the quadrants. This allows us to use the same routine to compute cos - by simply rotating the quadrants by 1. */ -static inline double -__always_inline -reduce_and_compute (double x, unsigned int k) -{ - double retval = 0, a, da; - unsigned int n = __branred (x, &a, &da); - k = (n + k) % 4; - switch (k) - { - case 0: - if (a * a < 0.01588) - retval = bsloww (a, da, x, n); - else - retval = bsloww1 (a, da, x, n); - break; - case 2: - if (a * a < 0.01588) - retval = bsloww (-a, -da, x, n); - else - retval = bsloww1 (-a, -da, x, n); - break; - - case 1: - case 3: - retval = bsloww2 (a, da, x, n); - break; - } - return retval; + return __copysign (sn + cor, xold); } +/* Reduce range of x to within PI/2 with abs (x) < 105414350. The high part + is written to *a, the low part to *da. Range reduction is accurate to 136 + bits so that when x is large and *a very close to zero, all 53 bits of *a + are correct. */ static inline int4 __always_inline -reduce_sincos_1 (double x, double *a, double *da) +reduce_sincos (double x, double *a, double *da) { mynumber v; @@ -281,198 +161,54 @@ reduce_sincos_1 (double x, double *a, double *da) 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 b, db, t1, t2; + t1 = xn * pp3; + t2 = y - t1; + db = (y - t2) - t1; - 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; + t1 = xn * pp4; + b = t2 - t1; + db += (t2 - b) - t1; *a = b; *da = db; - return n; } -/* Compute sin (A + DA). cos can be computed by shifting the quadrant N - clockwise. */ +/* Compute sin or cos (A + DA) for the given quadrant N. */ static double __always_inline -do_sincos_2 (double a, double da, double x, int4 n, int4 k) +do_sincos (double a, double da, int4 n) { - 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; + double retval; - 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; - } + if (n & 1) + /* Max ULP is 0.513. */ + retval = do_cos (a, da); + else + /* Max ULP is 0.501 if xx < 0.01588, otherwise ULP is 0.518. */ + retval = do_sin (a, da); - return retval; + return (n & 2) ? -retval : 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 +#ifndef IN_SINCOS double SECTION -#endif __sin (double x) { - double xx, res, t, cor, y, s, c, sn, ssn, cs, ccs; + double t, a, da; mynumber u; - int4 k, m; + int4 k, m, n; double retval = 0; -#ifndef IN_SINCOS SET_RESTORE_ROUND_53BIT (FE_TONEAREST); -#endif u.x = x; m = u.i[HIGH_HALF]; @@ -482,77 +218,34 @@ __sin (double x) math_check_force_underflow (x); retval = x; } - /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/ - else if (k < 0x3fd00000) - { - xx = x * x; - /* Taylor series. */ - t = POLYNOMIAL (xx) * (xx * x); - res = x + t; - cor = (x - res) + t; - retval = (res == res + 1.07 * cor) ? res : slow (x); - } /* else if (k < 0x3fd00000) */ -/*---------------------------- 0.25<|x|< 0.855469---------------------- */ +/*--------------------------- 2^-26<|x|< 0.855469---------------------- */ else if (k < 0x3feb6000) { - 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)); - SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); - if (m <= 0) - { - sn = -sn; - ssn = -ssn; - } - cor = (ssn + s * ccs - sn * c) + cs * s; - res = sn + cor; - cor = (sn - res) + cor; - retval = (res == res + 1.096 * cor) ? res : slow1 (x); + /* Max ULP is 0.548. */ + retval = do_sin (x, 0); } /* else if (k < 0x3feb6000) */ /*----------------------- 0.855469 <|x|<2.426265 ----------------------*/ else if (k < 0x400368fd) { - - y = (m > 0) ? hp0 - x : hp0 + x; - if (y >= 0) - { - u.x = big + y; - y = (y - (u.x - big)) + hp1; - } - else - { - u.x = big - y; - y = (-hp1) - (y + (u.x - big)); - } - res = do_cos (u, y, &cor); - retval = (res == res + 1.020 * cor) ? ((m > 0) ? res : -res) : slow2 (x); + t = hp0 - fabs (x); + /* Max ULP is 0.51. */ + retval = __copysign (do_cos (t, hp1), x); } /* else if (k < 0x400368fd) */ -#ifndef IN_SINCOS /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/ else if (k < 0x419921FB) { - double a, da; - int4 n = reduce_sincos_1 (x, &a, &da); - retval = do_sincos_1 (a, da, x, n, 0); + n = reduce_sincos (x, &a, &da); + retval = do_sincos (a, da, n); } /* else if (k < 0x419921FB ) */ -/*---------------------105414350 <|x|< 281474976710656 --------------------*/ - else if (k < 0x42F00000) - { - double a, da; - - 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----------------------------*/ +/* --------------------105414350 <|x| <2^1024------------------------------*/ else if (k < 0x7ff00000) - retval = reduce_and_compute (x, 0); - + { + n = __branred (x, &a, &da); + retval = do_sincos (a, da, n); + } /*--------------------- |x| > 2^1024 ----------------------------------*/ else { @@ -560,7 +253,6 @@ __sin (double x) __set_errno (EDOM); retval = x / x; } -#endif return retval; } @@ -571,23 +263,17 @@ __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, cor, a, da; + double y, a, da; mynumber u; - int4 k, m; + int4 k, m, n; double retval = 0; -#ifndef IN_SINCOS SET_RESTORE_ROUND_53BIT (FE_TONEAREST); -#endif u.x = x; m = u.i[HIGH_HALF]; @@ -599,11 +285,8 @@ __cos (double x) else if (k < 0x3feb6000) { /* 2^-27 < |x| < 0.855469 */ - y = fabs (x); - u.x = big + y; - y = y - (u.x - big); - res = do_cos (u, y, &cor); - retval = (res == res + 1.020 * cor) ? res : cslow2 (x); + /* Max ULP is 0.51. */ + retval = do_cos (x, 0); } /* else if (k < 0x3feb6000) */ else if (k < 0x400368fd) @@ -611,55 +294,23 @@ __cos (double x) y = hp0 - fabs (x); a = y + hp1; da = (y - a) + hp1; - xx = a * a; - if (xx < 0.01588) - { - 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 : sloww (a, da, x, 1); - } - 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 + 1.0e-31 : 1.035 * cor - 1.0e-31; - retval = ((res == res + cor) ? ((m) ? res : -res) - : sloww1 (a, da, x, m, 1)); - } - + /* Max ULP is 0.501 if xx < 0.01588 or 0.518 otherwise. + Range reduction uses 106 bits here which is sufficient. */ + retval = do_sin (a, da); } /* else if (k < 0x400368fd) */ - -#ifndef IN_SINCOS else if (k < 0x419921FB) { /* 2.426265<|x|< 105414350 */ - double a, da; - int4 n = reduce_sincos_1 (x, &a, &da); - retval = do_sincos_1 (a, da, x, n, 1); + n = reduce_sincos (x, &a, &da); + retval = do_sincos (a, da, n + 1); } /* else if (k < 0x419921FB ) */ - else if (k < 0x42F00000) - { - double a, da; - - 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 */ + /* 105414350 <|x| <2^1024 */ else if (k < 0x7ff00000) - retval = reduce_and_compute (x, 1); + { + n = __branred (x, &a, &da); + retval = do_sincos (a, da, n + 1); + } else { @@ -667,366 +318,15 @@ __cos (double x) __set_errno (EDOM); retval = x / x; /* |x| > 2^1024 */ } -#endif return retval; } -/************************************************************************/ -/* Routine compute sin(x) for 2^-26 < |x|< 0.25 by Taylor with more */ -/* precision and if still doesn't accurate enough by mpsin or dubsin */ -/************************************************************************/ - -static double -SECTION -slow (double x) -{ - double res, cor, w[2]; - res = TAYLOR_SLOW (x, 0, cor); - if (res == res + 1.0007 * cor) - return res; - - __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); -} - -/*******************************************************************************/ -/* Routine compute sin(x) for 0.25<|x|< 0.855469 by __sincostab.tbl and Taylor */ -/* and if result still doesn't accurate enough by mpsin or dubsin */ -/*******************************************************************************/ - -static double -SECTION -slow1 (double x) -{ - mynumber u; - double w[2], y, cor, res; - y = fabs (x); - u.x = big + y; - y = y - (u.x - big); - res = do_sin_slow (u, y, 0, 0, &cor); - if (res == res + cor) - return (x > 0) ? res : -res; - - __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); -} - -/**************************************************************************/ -/* Routine compute sin(x) for 0.855469 <|x|<2.426265 by __sincostab.tbl */ -/* and if result still doesn't accurate enough by mpsin or dubsin */ -/**************************************************************************/ -static double -SECTION -slow2 (double x) -{ - mynumber u; - double w[2], y, y1, y2, cor, res, del; - - y = fabs (x); - y = hp0 - y; - if (y >= 0) - { - u.x = big + y; - y = y - (u.x - big); - del = hp1; - } - else - { - u.x = big - y; - y = -(y + (u.x - big)); - del = -hp1; - } - res = do_cos_slow (u, y, del, 0, &cor); - if (res == res + cor) - return (x > 0) ? res : -res; - - 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); -} - -/***************************************************************************/ -/* Routine compute sin(x+dx) (Double-Length number) where x is small enough*/ -/* to use Taylor series around zero and (x+dx) */ -/* 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 mpsin1 or dubsin */ -/***************************************************************************/ - -static double -SECTION -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 - cor = 1.0005 * cor - fabs (orig) * 3.1e-30; - - 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) - { - 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]; - - return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (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 mpsin1 or dubsin */ -/***************************************************************************/ - -static double -SECTION -sloww1 (double x, double dx, double orig, int m, int k) -{ - 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; - - __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]; - - return (k == 1) ? __mpcos (orig, 0, true) : __mpsin (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 mpsin1 or dubsin */ -/***************************************************************************/ - -static double -SECTION -sloww2 (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); - - 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 - cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); - - if (w[0] == w[0] + cor) - return (n & 2) ? -w[0] : w[0]; - - return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true); -} - -/***************************************************************************/ -/* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */ -/* is small enough to use Taylor series around zero and (x+dx) */ -/* 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 -bsloww (double x, double dx, double orig, int n) -{ - double res, cor, w[2]; - - res = TAYLOR_SLOW (x, dx, cor); - 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 - 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); -} - -/***************************************************************************/ -/* Routine compute sin(x+dx) or cos(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 -bsloww1 (double x, double dx, double orig, int n) -{ - mynumber u; - double w[2], y, cor, res; - - y = fabs (x); - u.x = big + y; - y = y - (u.x - big); - dx = (x > 0) ? dx : -dx; - 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 - cor = 1.000000005 * 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); -} - -/***************************************************************************/ -/* Routine compute sin(x+dx) or cos(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 -bsloww2 (double x, double dx, double orig, int n) -{ - mynumber u; - double w[2], y, cor, res; - - y = fabs (x); - u.x = big + y; - y = y - (u.x - big); - dx = (x > 0) ? dx : -dx; - 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 - cor = 1.000000005 * w[1] - 1.1e-24; - - if (w[0] == w[0] + cor) - return (n & 2) ? -w[0] : w[0]; - - return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true); -} - -/************************************************************************/ -/* Routine compute cos(x) for 2^-27 < |x|< 0.25 by Taylor with more */ -/* precision and if still doesn't accurate enough by mpcos or docos */ -/************************************************************************/ - -static double -SECTION -cslow2 (double x) -{ - mynumber u; - double w[2], y, cor, res; - - y = fabs (x); - u.x = big + y; - y = y - (u.x - big); - res = do_cos_slow (u, y, 0, 0, &cor); - if (res == res + cor) - return res; - - y = fabs (x); - __docos (y, 0, w); - if (w[0] == w[0] + 1.000000005 * w[1]) - return w[0]; - - return __mpcos (x, 0, false); -} - #ifndef __cos -weak_alias (__cos, cos) -# ifdef NO_LONG_DOUBLE -strong_alias (__cos, __cosl) -weak_alias (__cos, cosl) -# endif +libm_alias_double (__cos, cos) #endif #ifndef __sin -weak_alias (__sin, sin) -# ifdef NO_LONG_DOUBLE -strong_alias (__sin, __sinl) -weak_alias (__sin, sinl) -# endif +libm_alias_double (__sin, sin) +#endif + #endif |