diff options
author | joseph <joseph@7b3dc134-2b1b-0410-93df-9e9f96275f8d> | 2013-05-18 00:51:47 +0000 |
---|---|---|
committer | joseph <joseph@7b3dc134-2b1b-0410-93df-9e9f96275f8d> | 2013-05-18 00:51:47 +0000 |
commit | eab7f6089510455a9b26643c64da331749a15650 (patch) | |
tree | e069c5f33da7c0cffbb68f47ec07b1b10b6789e4 /libc/sysdeps/ieee754 | |
parent | f9b341f7c8c64a0df8707b3cf29b425a25a52d12 (diff) | |
download | eglibc2-eab7f6089510455a9b26643c64da331749a15650.tar.gz |
Merge changes between r22954 and r23097 from /fsf/trunk.
git-svn-id: svn://svn.eglibc.org/trunk@23098 7b3dc134-2b1b-0410-93df-9e9f96275f8d
Diffstat (limited to 'libc/sysdeps/ieee754')
27 files changed, 2273 insertions, 952 deletions
diff --git a/libc/sysdeps/ieee754/dbl-64/e_gamma_r.c b/libc/sysdeps/ieee754/dbl-64/e_gamma_r.c index 987355175..5b17f7b5a 100644 --- a/libc/sysdeps/ieee754/dbl-64/e_gamma_r.c +++ b/libc/sysdeps/ieee754/dbl-64/e_gamma_r.c @@ -19,14 +19,104 @@ #include <math.h> #include <math_private.h> +#include <float.h> +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's + approximation to gamma function. */ + +static const double gamma_coeff[] = + { + 0x1.5555555555555p-4, + -0xb.60b60b60b60b8p-12, + 0x3.4034034034034p-12, + -0x2.7027027027028p-12, + 0x3.72a3c5631fe46p-12, + -0x7.daac36664f1f4p-12, + }; + +#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0])) + +/* Return gamma (X), for positive X less than 184, in the form R * + 2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to + avoid overflow or underflow in intermediate calculations. */ + +static double +gamma_positive (double x, int *exp2_adj) +{ + int local_signgam; + if (x < 0.5) + { + *exp2_adj = 0; + return __ieee754_exp (__ieee754_lgamma_r (x + 1, &local_signgam)) / x; + } + else if (x <= 1.5) + { + *exp2_adj = 0; + return __ieee754_exp (__ieee754_lgamma_r (x, &local_signgam)); + } + else if (x < 6.5) + { + /* Adjust into the range for using exp (lgamma). */ + *exp2_adj = 0; + double n = __ceil (x - 1.5); + double x_adj = x - n; + double eps; + double prod = __gamma_product (x_adj, 0, n, &eps); + return (__ieee754_exp (__ieee754_lgamma_r (x_adj, &local_signgam)) + * prod * (1.0 + eps)); + } + else + { + double eps = 0; + double x_eps = 0; + double x_adj = x; + double prod = 1; + if (x < 12.0) + { + /* Adjust into the range for applying Stirling's + approximation. */ + double n = __ceil (12.0 - x); +#if FLT_EVAL_METHOD != 0 + volatile +#endif + double x_tmp = x + n; + x_adj = x_tmp; + x_eps = (x - (x_adj - n)); + prod = __gamma_product (x_adj - n, x_eps, n, &eps); + } + /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)). + Compute gamma (X_ADJ + X_EPS) using Stirling's approximation, + starting by computing pow (X_ADJ, X_ADJ) with a power of 2 + factored out. */ + double exp_adj = -eps; + double x_adj_int = __round (x_adj); + double x_adj_frac = x_adj - x_adj_int; + int x_adj_log2; + double x_adj_mant = __frexp (x_adj, &x_adj_log2); + if (x_adj_mant < M_SQRT1_2) + { + x_adj_log2--; + x_adj_mant *= 2.0; + } + *exp2_adj = x_adj_log2 * (int) x_adj_int; + double ret = (__ieee754_pow (x_adj_mant, x_adj) + * __ieee754_exp2 (x_adj_log2 * x_adj_frac) + * __ieee754_exp (-x_adj) + * __ieee754_sqrt (2 * M_PI / x_adj) + / prod); + exp_adj += x_eps * __ieee754_log (x); + double bsum = gamma_coeff[NCOEFF - 1]; + double x_adj2 = x_adj * x_adj; + for (size_t i = 1; i <= NCOEFF - 1; i++) + bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i]; + exp_adj += bsum / x_adj; + return ret + ret * __expm1 (exp_adj); + } +} double __ieee754_gamma_r (double x, int *signgamp) { - /* We don't have a real gamma implementation now. We'll use lgamma - and the exp function. But due to the required boundary - conditions we must check some values separately. */ int32_t hx; u_int32_t lx; @@ -51,8 +141,48 @@ __ieee754_gamma_r (double x, int *signgamp) *signgamp = 0; return x - x; } + if (__builtin_expect ((hx & 0x7ff00000) == 0x7ff00000, 0)) + { + /* Positive infinity (return positive infinity) or NaN (return + NaN). */ + *signgamp = 0; + return x + x; + } - /* XXX FIXME. */ - return __ieee754_exp (__ieee754_lgamma_r (x, signgamp)); + if (x >= 172.0) + { + /* Overflow. */ + *signgamp = 0; + return DBL_MAX * DBL_MAX; + } + else if (x > 0.0) + { + *signgamp = 0; + int exp2_adj; + double ret = gamma_positive (x, &exp2_adj); + return __scalbn (ret, exp2_adj); + } + else if (x >= -DBL_EPSILON / 4.0) + { + *signgamp = 0; + return 1.0 / x; + } + else + { + double tx = __trunc (x); + *signgamp = (tx == 2.0 * __trunc (tx / 2.0)) ? -1 : 1; + if (x <= -184.0) + /* Underflow. */ + return DBL_MIN * DBL_MIN; + double frac = tx - x; + if (frac > 0.5) + frac = 1.0 - frac; + double sinpix = (frac <= 0.25 + ? __sin (M_PI * frac) + : __cos (M_PI * (0.5 - frac))); + int exp2_adj; + double ret = M_PI / (-x * sinpix * gamma_positive (-x, &exp2_adj)); + return __scalbn (ret, -exp2_adj); + } } strong_alias (__ieee754_gamma_r, __gamma_r_finite) diff --git a/libc/sysdeps/ieee754/dbl-64/e_remainder.c b/libc/sysdeps/ieee754/dbl-64/e_remainder.c index 39ca0c2d0..2d20bb1df 100644 --- a/libc/sysdeps/ieee754/dbl-64/e_remainder.c +++ b/libc/sysdeps/ieee754/dbl-64/e_remainder.c @@ -51,6 +51,7 @@ double __ieee754_remainder(double x, double y) ky=t.i[HIGH_HALF]; /*------ |x| < 2^1023 and 2^-970 < |y| < 2^1024 ------------------*/ if (kx<0x7fe00000 && ky<0x7ff00000 && ky>=0x03500000) { + SET_RESTORE_ROUND_NOEX (FE_TONEAREST); if (kx+0x00100000<ky) return x; if ((kx-0x01500000)<ky) { z=x/t.x; diff --git a/libc/sysdeps/ieee754/dbl-64/gamma_product.c b/libc/sysdeps/ieee754/dbl-64/gamma_product.c new file mode 100644 index 000000000..2a3fc1aff --- /dev/null +++ b/libc/sysdeps/ieee754/dbl-64/gamma_product.c @@ -0,0 +1,75 @@ +/* Compute a product of X, X+1, ..., with an error estimate. + Copyright (C) 2013 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <float.h> + +/* Calculate X * Y exactly and store the result in *HI + *LO. It is + given that the values are small enough that no overflow occurs and + large enough (or zero) that no underflow occurs. */ + +static void +mul_split (double *hi, double *lo, double x, double y) +{ +#ifdef __FP_FAST_FMA + /* Fast built-in fused multiply-add. */ + *hi = x * y; + *lo = __builtin_fma (x, y, -*hi); +#elif defined FP_FAST_FMA + /* Fast library fused multiply-add, compiler before GCC 4.6. */ + *hi = x * y; + *lo = __fma (x, y, -*hi); +#else + /* Apply Dekker's algorithm. */ + *hi = x * y; +# define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1) + double x1 = x * C; + double y1 = y * C; +# undef C + x1 = (x - x1) + x1; + y1 = (y - y1) + y1; + double x2 = x - x1; + double y2 = y - y1; + *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2; +#endif +} + +/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N + - 1, in the form R * (1 + *EPS) where the return value R is an + approximation to the product and *EPS is set to indicate the + approximate error in the return value. X is such that all the + values X + 1, ..., X + N - 1 are exactly representable, and X_EPS / + X is small enough that factors quadratic in it can be + neglected. */ + +double +__gamma_product (double x, double x_eps, int n, double *eps) +{ + SET_RESTORE_ROUND (FE_TONEAREST); + double ret = x; + *eps = x_eps / x; + for (int i = 1; i < n; i++) + { + *eps += x_eps / (x + i); + double lo; + mul_split (&ret, &lo, ret, x + i); + *eps += lo / ret; + } + return ret; +} diff --git a/libc/sysdeps/ieee754/dbl-64/gamma_productf.c b/libc/sysdeps/ieee754/dbl-64/gamma_productf.c new file mode 100644 index 000000000..46072f16e --- /dev/null +++ b/libc/sysdeps/ieee754/dbl-64/gamma_productf.c @@ -0,0 +1,46 @@ +/* Compute a product of X, X+1, ..., with an error estimate. + Copyright (C) 2013 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <float.h> + +/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N + - 1, in the form R * (1 + *EPS) where the return value R is an + approximation to the product and *EPS is set to indicate the + approximate error in the return value. X is such that all the + values X + 1, ..., X + N - 1 are exactly representable, and X_EPS / + X is small enough that factors quadratic in it can be + neglected. */ + +float +__gamma_productf (float x, float x_eps, int n, float *eps) +{ + double x_full = (double) x + (double) x_eps; + double ret = x_full; + for (int i = 1; i < n; i++) + ret *= x_full + i; + +#if FLT_EVAL_METHOD != 0 + volatile +#endif + float fret = ret; + *eps = (ret - fret) / fret; + + return fret; +} diff --git a/libc/sysdeps/ieee754/dbl-64/s_sin.c b/libc/sysdeps/ieee754/dbl-64/s_sin.c index 5038b7261..5c388c8b9 100644 --- a/libc/sysdeps/ieee754/dbl-64/s_sin.c +++ b/libc/sysdeps/ieee754/dbl-64/s_sin.c @@ -66,299 +66,359 @@ extern const union } __sincostab attribute_hidden; static const double - sn3 = -1.66666666666664880952546298448555E-01, - sn5 = 8.33333214285722277379541354343671E-03, - cs2 = 4.99999999999999999999950396842453E-01, - cs4 = -4.16666666666664434524222570944589E-02, - cs6 = 1.38888874007937613028114285595617E-03; - -void __dubsin(double x, double dx, double w[]); -void __docos(double x, double dx, double w[]); -double __mpsin(double x, double dx); -double __mpcos(double x, double dx); -double __mpsin1(double x); -double __mpcos1(double x); -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); -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); -static double csloww2(double x, double dx, double orig, int n); + sn3 = -1.66666666666664880952546298448555E-01, + sn5 = 8.33333214285722277379541354343671E-03, + cs2 = 4.99999999999999999999950396842453E-01, + cs4 = -4.16666666666664434524222570944589E-02, + cs6 = 1.38888874007937613028114285595617E-03; + +void __dubsin (double x, double dx, double w[]); +void __docos (double x, double dx, double w[]); +double __mpsin (double x, double dx); +double __mpcos (double x, double dx); +double __mpsin1 (double x); +double __mpcos1 (double x); +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); +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); +static double csloww2 (double x, double dx, double orig, int n); + /*******************************************************************/ /* An ultimate sin routine. Given an IEEE double machine number x */ /* it computes the correctly rounded (to nearest) value of sin(x) */ /*******************************************************************/ double SECTION -__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 retval = 0; - - SET_RESTORE_ROUND_53BIT (FE_TONEAREST); - - u.x = x; - m = u.i[HIGH_HALF]; - k = 0x7fffffff&m; /* no sign */ - if (k < 0x3e500000) /* if x->0 =>sin(x)=x */ - { retval = x; goto ret; } +__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 retval = 0; + + SET_RESTORE_ROUND_53BIT (FE_TONEAREST); + + u.x = x; + m = u.i[HIGH_HALF]; + k = 0x7fffffff & m; /* no sign */ + if (k < 0x3e500000) /* if x->0 =>sin(x)=x */ + { + retval = x; + goto ret; + } /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/ - else if (k < 0x3fd00000){ - xx = x*x; - /*Taylor series */ - t = ((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*(xx*x); - res = x+t; - cor = (x-res)+t; - retval = (res == res + 1.07*cor)? res : slow(x); - goto ret; - } /* else if (k < 0x3fd00000) */ + else if (k < 0x3fd00000) + { + xx = x * x; + /*Taylor series. */ + t = (((((s5.x * xx + s4.x) * xx + s3.x) * xx + s2.x) * xx + s1.x) + * (xx * x)); + res = x + t; + cor = (x - res) + t; + retval = (res == res + 1.07 * cor) ? res : slow (x); + goto ret; + } /* else if (k < 0x3fd00000) */ /*---------------------------- 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); - xx=y*y; - s = y + y*xx*(sn3 +xx*sn5); - c = xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=(m>0)?__sincostab.x[k]:-__sincostab.x[k]; - ssn=(m>0)?__sincostab.x[k+1]:-__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - cor=(ssn+s*ccs-sn*c)+cs*s; - res=sn+cor; - cor=(sn-res)+cor; - retval = (res==res+1.096*cor)? res : slow1(x); - goto ret; - } /* else if (k < 0x3feb6000) */ + 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); + xx = y * y; + s = y + y * xx * (sn3 + xx * sn5); + c = xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = (m > 0) ? __sincostab.x[k] : -__sincostab.x[k]; + ssn = (m > 0) ? __sincostab.x[k + 1] : -__sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + cor = (ssn + s * ccs - sn * c) + cs * s; + res = sn + cor; + cor = (sn - res) + cor; + retval = (res == res + 1.096 * cor) ? res : slow1 (x); + goto ret; + } /* else if (k < 0x3feb6000) */ /*----------------------- 0.855469 <|x|<2.426265 ----------------------*/ - else if (k < 0x400368fd ) { - - y = (m>0)? hp0.x-x:hp0.x+x; - if (y>=0) { - u.x = big.x+y; - y = (y-(u.x-big.x))+hp1.x; - } - else { - u.x = big.x-y; - y = (-hp1.x) - (y+(u.x-big.x)); - } - xx=y*y; - s = y + y*xx*(sn3 +xx*sn5); - c = xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - cor=(ccs-s*ssn-cs*c)-sn*s; - res=cs+cor; - cor=(cs-res)+cor; - retval = (res==res+1.020*cor)? ((m>0)?res:-res) : slow2(x); - goto ret; - } /* else if (k < 0x400368fd) */ + else if (k < 0x400368fd) + { + + y = (m > 0) ? hp0.x - x : hp0.x + x; + if (y >= 0) + { + u.x = big.x + y; + y = (y - (u.x - big.x)) + hp1.x; + } + else + { + u.x = big.x - y; + y = (-hp1.x) - (y + (u.x - big.x)); + } + xx = y * y; + s = y + y * xx * (sn3 + xx * sn5); + c = xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + cor = (ccs - s * ssn - cs * c) - sn * s; + res = cs + cor; + cor = (cs - res) + cor; + retval = (res == res + 1.020 * cor) ? ((m > 0) ? res : -res) : slow2 (x); + goto ret; + } /* else if (k < 0x400368fd) */ /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/ - else if (k < 0x419921FB ) { - t = (x*hpinv.x + toint.x); - xn = t - toint.x; - v.x = t; - y = (x - xn*mp1.x) - xn*mp2.x; - n =v.i[LOW_HALF]&3; - da = xn*mp3.x; - a=y-da; - da = (y-a)-da; - eps = ABS(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 */ - t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da; - res = a+t; - cor = (a-res)+t; - cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps; - retval = (res == res + cor)? res : sloww(a,da,x); + else if (k < 0x419921FB) + { + t = (x * hpinv.x + toint.x); + xn = t - toint.x; + v.x = t; + y = (x - xn * mp1.x) - xn * mp2.x; + n = v.i[LOW_HALF] & 3; + da = xn * mp3.x; + a = y - da; + da = (y - a) - da; + eps = ABS (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 */ + t = (((((s5.x * xx + s4.x) * xx + s3.x) * xx + s2.x) * xx + + s1.x) * a - 0.5 * da) * xx + da; + res = a + t; + cor = (a - res) + t; + cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; + retval = (res == res + cor) ? res : sloww (a, da, x); goto ret; } - else { - if (a>0) - {m=1;t=a;db=da;} + else + { + if (a > 0) + { + m = 1; + t = a; + db = da; + } else - {m=0;t=-a;db=-da;} - u.x=big.x+t; - y=t-(u.x-big.x); - xx=y*y; - s = y + (db+y*xx*(sn3 +xx*sn5)); - c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - cor=(ssn+s*ccs-sn*c)+cs*s; - res=sn+cor; - cor=(sn-res)+cor; - cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps; - retval = (res==res+cor)? ((m)?res:-res) : sloww1(a,da,x); + { + m = 0; + t = -a; + db = -da; + } + u.x = big.x + t; + y = t - (u.x - big.x); + xx = y * y; + s = y + (db + y * xx * (sn3 + xx * sn5)); + c = y * db + xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + cor = (ssn + s * ccs - sn * c) + cs * s; + res = sn + cor; + cor = (sn - res) + cor; + cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; + retval = ((res == res + cor) ? ((m) ? res : -res) + : sloww1 (a, da, x)); goto ret; } - break; - - case 1: - case 3: - if (a<0) - {a=-a;da=-da;} - u.x=big.x+a; - y=a-(u.x-big.x)+da; - xx=y*y; - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - s = y + y*xx*(sn3 +xx*sn5); - c = xx*(cs2 +xx*(cs4 + xx*cs6)); - cor=(ccs-s*ssn-cs*c)-sn*s; - res=cs+cor; - cor=(cs-res)+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); - goto ret; - - break; - - } - - } /* else if (k < 0x419921FB ) */ + break; + + case 1: + case 3: + if (a < 0) + { + a = -a; + da = -da; + } + u.x = big.x + a; + y = a - (u.x - big.x) + da; + xx = y * y; + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + s = y + y * xx * (sn3 + xx * sn5); + c = xx * (cs2 + xx * (cs4 + xx * cs6)); + cor = (ccs - s * ssn - cs * c) - sn * s; + res = cs + cor; + cor = (cs - res) + 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)); + goto ret; + + break; + } + + } /* else if (k < 0x419921FB ) */ /*---------------------105414350 <|x|< 281474976710656 --------------------*/ - else if (k < 0x42F00000 ) { - t = (x*hpinv.x + toint.x); - xn = t - toint.x; - 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); - n =v.i[LOW_HALF]&3; - da = xn1*pp3.x; - t=y-da; - da = (y-t)-da; - da = (da - xn2*pp3.x) -xn*pp4.x; - 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) { + else if (k < 0x42F00000) + { + t = (x * hpinv.x + toint.x); + xn = t - toint.x; + 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); + n = v.i[LOW_HALF] & 3; + da = xn1 * pp3.x; + t = y - da; + da = (y - t) - da; + da = (da - xn2 * pp3.x) - xn * pp4.x; + 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 */ - t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da; - res = a+t; - cor = (a-res)+t; - cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps; - retval = (res == res + cor)? res : bsloww(a,da,x,n); + t = (((((s5.x * xx + s4.x) * xx + s3.x) * xx + s2.x) * xx + + s1.x) * a - 0.5 * da) * xx + da; + res = a + t; + cor = (a - res) + t; + cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; + retval = (res == res + cor) ? res : bsloww (a, da, x, n); goto ret; } - else { - if (a>0) {m=1;t=a;db=da;} - else {m=0;t=-a;db=-da;} - u.x=big.x+t; - y=t-(u.x-big.x); - xx=y*y; - s = y + (db+y*xx*(sn3 +xx*sn5)); - c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - cor=(ssn+s*ccs-sn*c)+cs*s; - res=sn+cor; - cor=(sn-res)+cor; - cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps; - retval = (res==res+cor)? ((m)?res:-res) : bsloww1(a,da,x,n); + else + { + if (a > 0) + { + m = 1; + t = a; + db = da; + } + else + { + m = 0; + t = -a; + db = -da; + } + u.x = big.x + t; + y = t - (u.x - big.x); + xx = y * y; + s = y + (db + y * xx * (sn3 + xx * sn5)); + c = y * db + xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + cor = (ssn + s * ccs - sn * c) + cs * s; + res = sn + cor; + cor = (sn - res) + cor; + cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; + retval = ((res == res + cor) ? ((m) ? res : -res) + : bsloww1 (a, da, x, n)); goto ret; - } - break; - - case 1: - case 3: - if (a<0) - {a=-a;da=-da;} - u.x=big.x+a; - y=a-(u.x-big.x)+da; - xx=y*y; - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - s = y + y*xx*(sn3 +xx*sn5); - c = xx*(cs2 +xx*(cs4 + xx*cs6)); - cor=(ccs-s*ssn-cs*c)-sn*s; - res=cs+cor; - cor=(cs-res)+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); - goto ret; - - break; - - } - - } /* else if (k < 0x42F00000 ) */ + } + break; + + case 1: + case 3: + if (a < 0) + { + a = -a; + da = -da; + } + u.x = big.x + a; + y = a - (u.x - big.x) + da; + xx = y * y; + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + s = y + y * xx * (sn3 + xx * sn5); + c = xx * (cs2 + xx * (cs4 + xx * cs6)); + cor = (ccs - s * ssn - cs * c) - sn * s; + res = cs + cor; + cor = (cs - res) + 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)); + goto ret; + + break; + } + } /* else if (k < 0x42F00000 ) */ /* -----------------281474976710656 <|x| <2^1024----------------------------*/ - else if (k < 0x7ff00000) { - - n = __branred(x,&a,&da); - switch (n) { - case 0: - if (a*a < 0.01588) retval = bsloww(a,da,x,n); - else retval = bsloww1(a,da,x,n); - goto ret; - break; - case 2: - if (a*a < 0.01588) retval = bsloww(-a,-da,x,n); - else retval = bsloww1(-a,-da,x,n); - goto ret; - break; - - case 1: - case 3: - retval = bsloww2(a,da,x,n); - goto ret; - break; - } - - } /* else if (k < 0x7ff00000 ) */ + else if (k < 0x7ff00000) + { + n = __branred (x, &a, &da); + switch (n) + { + case 0: + if (a * a < 0.01588) + retval = bsloww (a, da, x, n); + else + retval = bsloww1 (a, da, x, n); + goto ret; + break; + case 2: + if (a * a < 0.01588) + retval = bsloww (-a, -da, x, n); + else + retval = bsloww1 (-a, -da, x, n); + goto ret; + break; -/*--------------------- |x| > 2^1024 ----------------------------------*/ - else { - if (k == 0x7ff00000 && u.i[LOW_HALF] == 0) - __set_errno (EDOM); - retval = x / x; + case 1: + case 3: + retval = bsloww2 (a, da, x, n); goto ret; + break; } + } /* else if (k < 0x7ff00000 ) */ + +/*--------------------- |x| > 2^1024 ----------------------------------*/ + else + { + if (k == 0x7ff00000 && u.i[LOW_HALF] == 0) + __set_errno (EDOM); + retval = x / x; + goto ret; + } - ret: - return retval; +ret: + return retval; } @@ -369,11 +429,12 @@ __sin(double x){ double SECTION -__cos(double x) +__cos (double x) { - double y,xx,res,t,cor,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2; - mynumber u,v; - int4 k,m,n; + double y, xx, res, t, cor, s, c, sn, ssn, cs, ccs, xn, a, da, db, eps, xn1, + xn2; + mynumber u, v; + int4 k, m, n; double retval = 0; @@ -381,251 +442,318 @@ __cos(double x) u.x = x; m = u.i[HIGH_HALF]; - k = 0x7fffffff&m; - - if (k < 0x3e400000 ) { retval = 1.0; goto ret; } /* |x|<2^-27 => cos(x)=1 */ - - else if (k < 0x3feb6000 ) {/* 2^-27 < |x| < 0.855469 */ - y=ABS(x); - u.x = big.x+y; - y = y-(u.x-big.x); - xx=y*y; - s = y + y*xx*(sn3 +xx*sn5); - c = xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - cor=(ccs-s*ssn-cs*c)-sn*s; - res=cs+cor; - cor=(cs-res)+cor; - retval = (res==res+1.020*cor)? res : cslow2(x); - goto ret; - -} /* else if (k < 0x3feb6000) */ - - else if (k < 0x400368fd ) {/* 0.855469 <|x|<2.426265 */; - y=hp0.x-ABS(x); - a=y+hp1.x; - da=(y-a)+hp1.x; - xx=a*a; - if (xx < 0.01588) { - t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da; - res = a+t; - cor = (a-res)+t; - cor = (cor>0)? 1.02*cor+1.0e-31 : 1.02*cor -1.0e-31; - retval = (res == res + cor)? res : csloww(a,da,x); - goto ret; - } - else { - if (a>0) {m=1;t=a;db=da;} - else {m=0;t=-a;db=-da;} - u.x=big.x+t; - y=t-(u.x-big.x); - xx=y*y; - s = y + (db+y*xx*(sn3 +xx*sn5)); - c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - cor=(ssn+s*ccs-sn*c)+cs*s; - res=sn+cor; - cor=(sn-res)+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); - goto ret; -} + k = 0x7fffffff & m; -} /* else if (k < 0x400368fd) */ - - - else if (k < 0x419921FB ) {/* 2.426265<|x|< 105414350 */ - t = (x*hpinv.x + toint.x); - xn = t - toint.x; - v.x = t; - y = (x - xn*mp1.x) - xn*mp2.x; - n =v.i[LOW_HALF]&3; - da = xn*mp3.x; - a=y-da; - da = (y-a)-da; - eps = ABS(x)*1.2e-30; - - switch (n) { - case 1: - case 3: - xx = a*a; - if (n == 1) {a=-a;da=-da;} - if (xx < 0.01588) { - t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da; - res = a+t; - cor = (a-res)+t; - cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps; - retval = (res == res + cor)? res : csloww(a,da,x); - goto ret; - } - else { - if (a>0) {m=1;t=a;db=da;} - else {m=0;t=-a;db=-da;} - u.x=big.x+t; - y=t-(u.x-big.x); - xx=y*y; - s = y + (db+y*xx*(sn3 +xx*sn5)); - c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - cor=(ssn+s*ccs-sn*c)+cs*s; - res=sn+cor; - cor=(sn-res)+cor; - cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps; - retval = (res==res+cor)? ((m)?res:-res) : csloww1(a,da,x); - goto ret; - } - break; - - case 0: - case 2: - if (a<0) {a=-a;da=-da;} - u.x=big.x+a; - y=a-(u.x-big.x)+da; - xx=y*y; - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - s = y + y*xx*(sn3 +xx*sn5); - c = xx*(cs2 +xx*(cs4 + xx*cs6)); - cor=(ccs-s*ssn-cs*c)-sn*s; - res=cs+cor; - cor=(cs-res)+cor; - cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps; - retval = (res==res+cor)? ((n)?-res:res) : csloww2(a,da,x,n); + if (k < 0x3e400000) + { + retval = 1.0; + goto ret; + } /* |x|<2^-27 => cos(x)=1 */ + + else if (k < 0x3feb6000) + { /* 2^-27 < |x| < 0.855469 */ + y = ABS (x); + u.x = big.x + y; + y = y - (u.x - big.x); + xx = y * y; + s = y + y * xx * (sn3 + xx * sn5); + c = xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + cor = (ccs - s * ssn - cs * c) - sn * s; + res = cs + cor; + cor = (cs - res) + cor; + retval = (res == res + 1.020 * cor) ? res : cslow2 (x); goto ret; + } /* else if (k < 0x3feb6000) */ + + else if (k < 0x400368fd) + { /* 0.855469 <|x|<2.426265 */ ; + y = hp0.x - ABS (x); + a = y + hp1.x; + da = (y - a) + hp1.x; + xx = a * a; + if (xx < 0.01588) + { + t = (((((s5.x * xx + s4.x) * xx + s3.x) * xx + s2.x) * xx + s1.x) + * a - 0.5 * da) * xx + da; + res = a + t; + cor = (a - res) + t; + cor = (cor > 0) ? 1.02 * cor + 1.0e-31 : 1.02 * cor - 1.0e-31; + retval = (res == res + cor) ? res : csloww (a, da, x); + goto ret; + } + else + { + if (a > 0) + { + m = 1; + t = a; + db = da; + } + else + { + m = 0; + t = -a; + db = -da; + } + u.x = big.x + t; + y = t - (u.x - big.x); + xx = y * y; + s = y + (db + y * xx * (sn3 + xx * sn5)); + c = y * db + xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + cor = (ssn + s * ccs - sn * c) + cs * s; + res = sn + cor; + cor = (sn - res) + 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)); + goto ret; + } - break; + } /* else if (k < 0x400368fd) */ - } - } /* else if (k < 0x419921FB ) */ - - - else if (k < 0x42F00000 ) { - t = (x*hpinv.x + toint.x); - xn = t - toint.x; - 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); - n =v.i[LOW_HALF]&3; - da = xn1*pp3.x; - t=y-da; - da = (y-t)-da; - da = (da - xn2*pp3.x) -xn*pp4.x; - 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) { - t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da; - res = a+t; - cor = (a-res)+t; - cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps; - retval = (res == res + cor)? res : bsloww(a,da,x,n); - goto ret; - } - else { - if (a>0) {m=1;t=a;db=da;} - else {m=0;t=-a;db=-da;} - u.x=big.x+t; - y=t-(u.x-big.x); - xx=y*y; - s = y + (db+y*xx*(sn3 +xx*sn5)); - c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - cor=(ssn+s*ccs-sn*c)+cs*s; - res=sn+cor; - cor=(sn-res)+cor; - cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps; - retval = (res==res+cor)? ((m)?res:-res) : bsloww1(a,da,x,n); - goto ret; - } - break; - - case 0: - case 2: - if (a<0) {a=-a;da=-da;} - u.x=big.x+a; - y=a-(u.x-big.x)+da; - xx=y*y; - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - s = y + y*xx*(sn3 +xx*sn5); - c = xx*(cs2 +xx*(cs4 + xx*cs6)); - cor=(ccs-s*ssn-cs*c)-sn*s; - res=cs+cor; - cor=(cs-res)+cor; - cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps; - retval = (res==res+cor)? ((n)?-res:res) : bsloww2(a,da,x,n); - goto ret; - break; + else if (k < 0x419921FB) + { /* 2.426265<|x|< 105414350 */ + t = (x * hpinv.x + toint.x); + xn = t - toint.x; + v.x = t; + y = (x - xn * mp1.x) - xn * mp2.x; + n = v.i[LOW_HALF] & 3; + da = xn * mp3.x; + a = y - da; + da = (y - a) - da; + eps = ABS (x) * 1.2e-30; + + switch (n) + { + case 1: + case 3: + xx = a * a; + if (n == 1) + { + a = -a; + da = -da; + } + if (xx < 0.01588) + { + t = (((((s5.x * xx + s4.x) * xx + s3.x) * xx + s2.x) * xx + + s1.x) * a - 0.5 * da) * xx + da; + res = a + t; + cor = (a - res) + t; + cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; + retval = (res == res + cor) ? res : csloww (a, da, x); + goto ret; + } + else + { + if (a > 0) + { + m = 1; + t = a; + db = da; + } + else + { + m = 0; + t = -a; + db = -da; + } + u.x = big.x + t; + y = t - (u.x - big.x); + xx = y * y; + s = y + (db + y * xx * (sn3 + xx * sn5)); + c = y * db + xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + cor = (ssn + s * ccs - sn * c) + cs * s; + res = sn + cor; + cor = (sn - res) + cor; + cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; + retval = ((res == res + cor) ? ((m) ? res : -res) + : csloww1 (a, da, x)); + goto ret; + } + break; + + case 0: + case 2: + if (a < 0) + { + a = -a; + da = -da; + } + u.x = big.x + a; + y = a - (u.x - big.x) + da; + xx = y * y; + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + s = y + y * xx * (sn3 + xx * sn5); + c = xx * (cs2 + xx * (cs4 + xx * cs6)); + cor = (ccs - s * ssn - cs * c) - sn * s; + res = cs + cor; + cor = (cs - res) + cor; + cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; + retval = ((res == res + cor) ? ((n) ? -res : res) + : csloww2 (a, da, x, n)); + goto ret; - } + break; + } + } /* else if (k < 0x419921FB ) */ - } /* else if (k < 0x42F00000 ) */ + else if (k < 0x42F00000) + { + t = (x * hpinv.x + toint.x); + xn = t - toint.x; + 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); + n = v.i[LOW_HALF] & 3; + da = xn1 * pp3.x; + t = y - da; + da = (y - t) - da; + da = (da - xn2 * pp3.x) - xn * pp4.x; + 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) + { + t = (((((s5.x * xx + s4.x) * xx + s3.x) * xx + s2.x) * xx + + s1.x) * a - 0.5 * da) * xx + da; + res = a + t; + cor = (a - res) + t; + cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; + retval = (res == res + cor) ? res : bsloww (a, da, x, n); + goto ret; + } + else + { + if (a > 0) + { + m = 1; + t = a; + db = da; + } + else + { + m = 0; + t = -a; + db = -da; + } + u.x = big.x + t; + y = t - (u.x - big.x); + xx = y * y; + s = y + (db + y * xx * (sn3 + xx * sn5)); + c = y * db + xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + cor = (ssn + s * ccs - sn * c) + cs * s; + res = sn + cor; + cor = (sn - res) + cor; + cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; + retval = ((res == res + cor) ? ((m) ? res : -res) + : bsloww1 (a, da, x, n)); + goto ret; + } + break; + + case 0: + case 2: + if (a < 0) + { + a = -a; + da = -da; + } + u.x = big.x + a; + y = a - (u.x - big.x) + da; + xx = y * y; + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + s = y + y * xx * (sn3 + xx * sn5); + c = xx * (cs2 + xx * (cs4 + xx * cs6)); + cor = (ccs - s * ssn - cs * c) - sn * s; + res = cs + cor; + cor = (cs - res) + cor; + cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; + retval = ((res == res + cor) ? ((n) ? -res : res) + : bsloww2 (a, da, x, n)); + goto ret; + break; + } + } /* else if (k < 0x42F00000 ) */ + + else if (k < 0x7ff00000) + { /* 281474976710656 <|x| <2^1024 */ + + n = __branred (x, &a, &da); + switch (n) + { + case 1: + if (a * a < 0.01588) + retval = bsloww (-a, -da, x, n); + else + retval = bsloww1 (-a, -da, x, n); + goto ret; + break; + case 3: + if (a * a < 0.01588) + retval = bsloww (a, da, x, n); + else + retval = bsloww1 (a, da, x, n); + goto ret; + break; - else if (k < 0x7ff00000) {/* 281474976710656 <|x| <2^1024 */ + case 0: + case 2: + retval = bsloww2 (a, da, x, n); + goto ret; + break; + } + } /* else if (k < 0x7ff00000 ) */ - n = __branred(x,&a,&da); - switch (n) { - case 1: - if (a*a < 0.01588) retval = bsloww(-a,-da,x,n); - else retval = bsloww1(-a,-da,x,n); + else + { + if (k == 0x7ff00000 && u.i[LOW_HALF] == 0) + __set_errno (EDOM); + retval = x / x; /* |x| > 2^1024 */ goto ret; - break; - case 3: - if (a*a < 0.01588) retval = bsloww(a,da,x,n); - else retval = bsloww1(a,da,x,n); - goto ret; - break; - - case 0: - case 2: - retval = bsloww2(a,da,x,n); - goto ret; - break; } - } /* else if (k < 0x7ff00000 ) */ - - - - - else { - if (k == 0x7ff00000 && u.i[LOW_HALF] == 0) - __set_errno (EDOM); - retval = x / x; /* |x| > 2^1024 */ - goto ret; - } - - ret: +ret: return retval; } @@ -636,25 +764,32 @@ __cos(double x) static double SECTION -slow(double x) { -static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ - double y,x1,x2,xx,r,t,res,cor,w[2]; - x1=(x+th2_36)-th2_36; - y = aa.x*x1*x1*x1; - r=x+y; - x2=x-x1; - xx=x*x; - t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + bb.x)*xx + 3.0*aa.x*x1*x2)*x +aa.x*x2*x2*x2; - t=((x-r)+y)+t; - res=r+t; - cor = (r-res)+t; - if (res == res + 1.0007*cor) return res; - else { - __dubsin(ABS(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):-__mpsin(-x,0); - } +slow (double x) +{ + static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ + double y, x1, x2, xx, r, t, res, cor, w[2]; + x1 = (x + th2_36) - th2_36; + y = aa.x * x1 * x1 * x1; + r = x + y; + x2 = x - x1; + xx = x * x; + t = (((((s5.x * xx + s4.x) * xx + s3.x) * xx + s2.x) * xx + bb.x) * xx + + 3.0 * aa.x * x1 * x2) * x + aa.x * x2 * x2 * x2; + t = ((x - r) + y) + t; + res = r + t; + cor = (r - res) + t; + if (res == res + 1.0007 * cor) + return res; + else + { + __dubsin (ABS (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) : -__mpsin (-x, 0); + } } + /*******************************************************************************/ /* 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 */ @@ -662,88 +797,102 @@ static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ static double SECTION -slow1(double x) { +slow1 (double x) +{ mynumber u; - double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res; + double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, c1, c2, xx, cor, res; static const double t22 = 6291456.0; int4 k; - y=ABS(x); - u.x=big.x+y; - y=y-(u.x-big.x); - xx=y*y; - s = y*xx*(sn3 +xx*sn5); - c = xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; /* Data */ - ssn=__sincostab.x[k+1]; /* from */ - cs=__sincostab.x[k+2]; /* tables */ - ccs=__sincostab.x[k+3]; /* __sincostab.tbl */ - y1 = (y+t22)-t22; + y = ABS (x); + u.x = big.x + y; + y = y - (u.x - big.x); + xx = y * y; + s = y * xx * (sn3 + xx * sn5); + c = xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; /* Data */ + ssn = __sincostab.x[k + 1]; /* from */ + cs = __sincostab.x[k + 2]; /* tables */ + ccs = __sincostab.x[k + 3]; /* __sincostab.tbl */ + y1 = (y + t22) - t22; y2 = y - y1; - c1 = (cs+t22)-t22; - c2=(cs-c1)+ccs; - cor=(ssn+s*ccs+cs*s+c2*y+c1*y2)-sn*c; - y=sn+c1*y1; - cor = cor+((sn-y)+c1*y1); - res=y+cor; - cor=(y-res)+cor; - if (res == res+1.0005*cor) return (x>0)?res:-res; - else { - __dubsin(ABS(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):-__mpsin(-x,0); - } + c1 = (cs + t22) - t22; + c2 = (cs - c1) + ccs; + cor = (ssn + s * ccs + cs * s + c2 * y + c1 * y2) - sn * c; + y = sn + c1 * y1; + cor = cor + ((sn - y) + c1 * y1); + res = y + cor; + cor = (y - res) + cor; + if (res == res + 1.0005 * cor) + return (x > 0) ? res : -res; + else + { + __dubsin (ABS (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) : -__mpsin (-x, 0); + } } + /**************************************************************************/ /* 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) { +slow2 (double x) +{ mynumber u; - double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res,del; + double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, e1, e2, xx, cor, res, del; static const double t22 = 6291456.0; int4 k; - y=ABS(x); - y = hp0.x-y; - if (y>=0) { - u.x = big.x+y; - y = y-(u.x-big.x); - del = hp1.x; - } - else { - u.x = big.x-y; - y = -(y+(u.x-big.x)); - del = -hp1.x; - } - xx=y*y; - s = y*xx*(sn3 +xx*sn5); - c = y*del+xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - y1 = (y+t22)-t22; - y2 = (y - y1)+del; - e1 = (sn+t22)-t22; - e2=(sn-e1)+ssn; - cor=(ccs-cs*c-e1*y2-e2*y)-sn*s; - y=cs-e1*y1; - cor = cor+((cs-y)-e1*y1); - res=y+cor; - cor=(y-res)+cor; - if (res == res+1.0005*cor) return (x>0)?res:-res; - else { - y=ABS(x)-hp0.x; - y1=y-hp1.x; - y2=(y-y1)-hp1.x; - __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):-__mpsin(-x,0); - } + y = ABS (x); + y = hp0.x - y; + if (y >= 0) + { + u.x = big.x + y; + y = y - (u.x - big.x); + del = hp1.x; + } + else + { + u.x = big.x - y; + y = -(y + (u.x - big.x)); + del = -hp1.x; + } + xx = y * y; + s = y * xx * (sn3 + xx * sn5); + c = y * del + xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + y1 = (y + t22) - t22; + y2 = (y - y1) + del; + e1 = (sn + t22) - t22; + e2 = (sn - e1) + ssn; + cor = (ccs - cs * c - e1 * y2 - e2 * y) - sn * s; + y = cs - e1 * y1; + cor = cor + ((cs - y) - e1 * y1); + res = y + cor; + cor = (y - res) + cor; + if (res == res + 1.0005 * cor) + return (x > 0) ? res : -res; + else + { + y = ABS (x) - hp0.x; + y1 = y - hp1.x; + y2 = (y - y1) - hp1.x; + __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) : -__mpsin (-x, 0); + } } + /***************************************************************************/ /* Routine compute sin(x+dx) (Double-Length number) where x is small enough*/ /* to use Taylor series around zero and (x+dx) */ @@ -754,46 +903,74 @@ slow2(double x) { static double SECTION -sloww(double x,double dx, double orig) { - static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ - double y,x1,x2,xx,r,t,res,cor,w[2],a,da,xn; - union {int4 i[2]; double x;} v; +sloww (double x, double dx, double orig) +{ + static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ + double y, x1, x2, xx, r, t, res, cor, w[2], a, da, xn; + union + { + int4 i[2]; + double x; + } v; int4 n; - x1=(x+th2_36)-th2_36; - y = aa.x*x1*x1*x1; - r=x+y; - x2=(x-x1)+dx; - xx=x*x; - t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + bb.x)*xx + 3.0*aa.x*x1*x2)*x +aa.x*x2*x2*x2+dx; - t=((x-r)+y)+t; - res=r+t; - cor = (r-res)+t; - cor = (cor>0)? 1.0005*cor+ABS(orig)*3.1e-30 : 1.0005*cor-ABS(orig)*3.1e-30; - if (res == res + cor) return res; - else { - (x>0)? __dubsin(x,dx,w) : __dubsin(-x,-dx,w); - cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-30 : 1.000000001*w[1] - ABS(orig)*1.1e-30; - if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0]; - else { - t = (orig*hpinv.x + toint.x); - xn = t - toint.x; - v.x = t; - y = (orig - xn*mp1.x) - xn*mp2.x; - n =v.i[LOW_HALF]&3; - da = xn*pp3.x; - t=y-da; - da = (y-t)-da; - y = xn*pp4.x; - 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); - cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-40 : 1.000000001*w[1] - ABS(orig)*1.1e-40; - if (w[0] == w[0]+cor) return (a>0)?w[0]:-w[0]; - else return __mpsin1(orig); + x1 = (x + th2_36) - th2_36; + y = aa.x * x1 * x1 * x1; + r = x + y; + x2 = (x - x1) + dx; + xx = x * x; + t = (((((s5.x * xx + s4.x) * xx + s3.x) * xx + s2.x) * xx + bb.x) * xx + + 3.0 * aa.x * x1 * x2) * x + aa.x * x2 * x2 * x2 + dx; + t = ((x - r) + y) + t; + res = r + t; + cor = (r - res) + t; + cor = + (cor > + 0) ? 1.0005 * cor + ABS (orig) * 3.1e-30 : 1.0005 * cor - + ABS (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] + ABS (orig) * 1.1e-30; + else + cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-30; + + if (w[0] == w[0] + cor) + return (x > 0) ? w[0] : -w[0]; + else + { + t = (orig * hpinv.x + toint.x); + xn = t - toint.x; + v.x = t; + y = (orig - xn * mp1.x) - xn * mp2.x; + n = v.i[LOW_HALF] & 3; + da = xn * pp3.x; + t = y - da; + da = (y - t) - da; + y = xn * pp4.x; + 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] + ABS (orig) * 1.1e-40; + else + cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-40; + + if (w[0] == w[0] + cor) + return (a > 0) ? w[0] : -w[0]; + else + return __mpsin1 (orig); + } } - } } + /***************************************************************************/ /* Routine compute sin(x+dx) (Double-Length number) where x in first or */ /* third quarter of unit circle.Routine receive also (right argument) the */ @@ -803,41 +980,58 @@ sloww(double x,double dx, double orig) { static double SECTION -sloww1(double x, double dx, double orig) { +sloww1 (double x, double dx, double orig) +{ mynumber u; - double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res; + double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, c1, c2, xx, cor, res; static const double t22 = 6291456.0; int4 k; - y=ABS(x); - u.x=big.x+y; - y=y-(u.x-big.x); - dx=(x>0)?dx:-dx; - xx=y*y; - s = y*xx*(sn3 +xx*sn5); - c = xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - y1 = (y+t22)-t22; - y2 = (y - y1)+dx; - c1 = (cs+t22)-t22; - c2=(cs-c1)+ccs; - cor=(ssn+s*ccs+cs*s+c2*y+c1*y2-sn*y*dx)-sn*c; - y=sn+c1*y1; - cor = cor+((sn-y)+c1*y1); - res=y+cor; - cor=(y-res)+cor; - cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig); - if (res == res + cor) return (x>0)?res:-res; - else { - __dubsin(ABS(x),dx,w); - cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig); - if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0]; - else return __mpsin1(orig); - } + + y = ABS (x); + u.x = big.x + y; + y = y - (u.x - big.x); + dx = (x > 0) ? dx : -dx; + xx = y * y; + s = y * xx * (sn3 + xx * sn5); + c = xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + y1 = (y + t22) - t22; + y2 = (y - y1) + dx; + c1 = (cs + t22) - t22; + c2 = (cs - c1) + ccs; + cor = (ssn + s * ccs + cs * s + c2 * y + c1 * y2 - sn * y * dx) - sn * c; + y = sn + c1 * y1; + cor = cor + ((sn - y) + c1 * y1); + res = y + cor; + cor = (y - res) + cor; + + if (cor > 0) + cor = 1.0005 * cor + 3.1e-30 * ABS (orig); + else + cor = 1.0005 * cor - 3.1e-30 * ABS (orig); + + if (res == res + cor) + return (x > 0) ? res : -res; + else + { + __dubsin (ABS (x), dx, w); + + if (w[1] > 0) + cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig); + else + cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig); + + if (w[0] == w[0] + cor) + return (x > 0) ? w[0] : -w[0]; + else + return __mpsin1 (orig); + } } + /***************************************************************************/ /* Routine compute sin(x+dx) (Double-Length number) where x in second or */ /* fourth quarter of unit circle.Routine receive also the original value */ @@ -847,42 +1041,59 @@ sloww1(double x, double dx, double orig) { static double SECTION -sloww2(double x, double dx, double orig, int n) { +sloww2 (double x, double dx, double orig, int n) +{ mynumber u; - double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res; + double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, e1, e2, xx, cor, res; static const double t22 = 6291456.0; int4 k; - y=ABS(x); - u.x=big.x+y; - y=y-(u.x-big.x); - dx=(x>0)?dx:-dx; - xx=y*y; - s = y*xx*(sn3 +xx*sn5); - c = y*dx+xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - - y1 = (y+t22)-t22; - y2 = (y - y1)+dx; - e1 = (sn+t22)-t22; - e2=(sn-e1)+ssn; - cor=(ccs-cs*c-e1*y2-e2*y)-sn*s; - y=cs-e1*y1; - cor = cor+((cs-y)-e1*y1); - res=y+cor; - cor=(y-res)+cor; - cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig); - if (res == res + cor) return (n&2)?-res:res; - else { - __docos(ABS(x),dx,w); - cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig); - if (w[0] == w[0]+cor) return (n&2)?-w[0]:w[0]; - else return __mpsin1(orig); - } + + y = ABS (x); + u.x = big.x + y; + y = y - (u.x - big.x); + dx = (x > 0) ? dx : -dx; + xx = y * y; + s = y * xx * (sn3 + xx * sn5); + c = y * dx + xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + + y1 = (y + t22) - t22; + y2 = (y - y1) + dx; + e1 = (sn + t22) - t22; + e2 = (sn - e1) + ssn; + cor = (ccs - cs * c - e1 * y2 - e2 * y) - sn * s; + y = cs - e1 * y1; + cor = cor + ((cs - y) - e1 * y1); + res = y + cor; + cor = (y - res) + cor; + + if (cor > 0) + cor = 1.0005 * cor + 3.1e-30 * ABS (orig); + else + cor = 1.0005 * cor - 3.1e-30 * ABS (orig); + + if (res == res + cor) + return (n & 2) ? -res : res; + else + { + __docos (ABS (x), dx, w); + + if (w[1] > 0) + cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig); + else + cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig); + + if (w[0] == w[0] + cor) + return (n & 2) ? -w[0] : w[0]; + else + return __mpsin1 (orig); + } } + /***************************************************************************/ /* 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) */ @@ -893,26 +1104,36 @@ sloww2(double x, double dx, double orig, int n) { static double SECTION -bsloww(double x,double dx, double orig,int n) { - static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ - double y,x1,x2,xx,r,t,res,cor,w[2]; - x1=(x+th2_36)-th2_36; - y = aa.x*x1*x1*x1; - r=x+y; - x2=(x-x1)+dx; - xx=x*x; - t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + bb.x)*xx + 3.0*aa.x*x1*x2)*x +aa.x*x2*x2*x2+dx; - t=((x-r)+y)+t; - res=r+t; - cor = (r-res)+t; - cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24; - if (res == res + cor) return res; - else { - (x>0)? __dubsin(x,dx,w) : __dubsin(-x,-dx,w); - cor = (w[1]>0)? 1.000000001*w[1] + 1.1e-24 : 1.000000001*w[1] - 1.1e-24; - if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0]; - else return (n&1)?__mpcos1(orig):__mpsin1(orig); - } +bsloww (double x, double dx, double orig, int n) +{ + static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ + double y, x1, x2, xx, r, t, res, cor, w[2]; + + x1 = (x + th2_36) - th2_36; + y = aa.x * x1 * x1 * x1; + r = x + y; + x2 = (x - x1) + dx; + xx = x * x; + t = (((((s5.x * xx + s4.x) * xx + s3.x) * xx + s2.x) * xx + bb.x) * xx + + 3.0 * aa.x * x1 * x2) * x + aa.x * x2 * x2 * x2 + dx; + t = ((x - r) + y) + t; + res = r + t; + cor = (r - res) + t; + cor = (cor > 0) ? 1.0005 * cor + 1.1e-24 : 1.0005 * cor - 1.1e-24; + 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] + 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) ? __mpcos1 (orig) : __mpsin1 (orig); + } } /***************************************************************************/ @@ -924,40 +1145,51 @@ bsloww(double x,double dx, double orig,int n) { static double SECTION -bsloww1(double x, double dx, double orig,int n) { -mynumber u; - double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res; - static const double t22 = 6291456.0; - int4 k; - y=ABS(x); - u.x=big.x+y; - y=y-(u.x-big.x); - dx=(x>0)?dx:-dx; - xx=y*y; - s = y*xx*(sn3 +xx*sn5); - c = xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - y1 = (y+t22)-t22; - y2 = (y - y1)+dx; - c1 = (cs+t22)-t22; - c2=(cs-c1)+ccs; - cor=(ssn+s*ccs+cs*s+c2*y+c1*y2-sn*y*dx)-sn*c; - y=sn+c1*y1; - cor = cor+((sn-y)+c1*y1); - res=y+cor; - cor=(y-res)+cor; - cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24; - if (res == res + cor) return (x>0)?res:-res; - else { - __dubsin(ABS(x),dx,w); - cor = (w[1]>0)? 1.000000005*w[1]+1.1e-24: 1.000000005*w[1]-1.1e-24; - if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0]; - else return (n&1)?__mpcos1(orig):__mpsin1(orig); - } +bsloww1 (double x, double dx, double orig, int n) +{ + mynumber u; + double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, c1, c2, xx, cor, res; + static const double t22 = 6291456.0; + int4 k; + + y = ABS (x); + u.x = big.x + y; + y = y - (u.x - big.x); + dx = (x > 0) ? dx : -dx; + xx = y * y; + s = y * xx * (sn3 + xx * sn5); + c = xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + y1 = (y + t22) - t22; + y2 = (y - y1) + dx; + c1 = (cs + t22) - t22; + c2 = (cs - c1) + ccs; + cor = (ssn + s * ccs + cs * s + c2 * y + c1 * y2 - sn * y * dx) - sn * c; + y = sn + c1 * y1; + cor = cor + ((sn - y) + c1 * y1); + res = y + cor; + cor = (y - res) + cor; + cor = (cor > 0) ? 1.0005 * cor + 1.1e-24 : 1.0005 * cor - 1.1e-24; + if (res == res + cor) + return (x > 0) ? res : -res; + else + { + __dubsin (ABS (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]; + else + return (n & 1) ? __mpcos1 (orig) : __mpsin1 (orig); + } } /***************************************************************************/ @@ -969,41 +1201,52 @@ mynumber u; static double SECTION -bsloww2(double x, double dx, double orig, int n) { -mynumber u; - double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res; - static const double t22 = 6291456.0; - int4 k; - y=ABS(x); - u.x=big.x+y; - y=y-(u.x-big.x); - dx=(x>0)?dx:-dx; - xx=y*y; - s = y*xx*(sn3 +xx*sn5); - c = y*dx+xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - - y1 = (y+t22)-t22; - y2 = (y - y1)+dx; - e1 = (sn+t22)-t22; - e2=(sn-e1)+ssn; - cor=(ccs-cs*c-e1*y2-e2*y)-sn*s; - y=cs-e1*y1; - cor = cor+((cs-y)-e1*y1); - res=y+cor; - cor=(y-res)+cor; - cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24; - if (res == res + cor) return (n&2)?-res:res; - else { - __docos(ABS(x),dx,w); - cor = (w[1]>0)? 1.000000005*w[1]+1.1e-24 : 1.000000005*w[1]-1.1e-24; - if (w[0] == w[0]+cor) return (n&2)?-w[0]:w[0]; - else return (n&1)?__mpsin1(orig):__mpcos1(orig); - } +bsloww2 (double x, double dx, double orig, int n) +{ + mynumber u; + double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, e1, e2, xx, cor, res; + static const double t22 = 6291456.0; + int4 k; + + y = ABS (x); + u.x = big.x + y; + y = y - (u.x - big.x); + dx = (x > 0) ? dx : -dx; + xx = y * y; + s = y * xx * (sn3 + xx * sn5); + c = y * dx + xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + + y1 = (y + t22) - t22; + y2 = (y - y1) + dx; + e1 = (sn + t22) - t22; + e2 = (sn - e1) + ssn; + cor = (ccs - cs * c - e1 * y2 - e2 * y) - sn * s; + y = cs - e1 * y1; + cor = cor + ((cs - y) - e1 * y1); + res = y + cor; + cor = (y - res) + cor; + cor = (cor > 0) ? 1.0005 * cor + 1.1e-24 : 1.0005 * cor - 1.1e-24; + if (res == res + cor) + return (n & 2) ? -res : res; + else + { + __docos (ABS (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]; + else + return (n & 1) ? __mpsin1 (orig) : __mpcos1 (orig); + } } /************************************************************************/ @@ -1013,39 +1256,44 @@ mynumber u; static double SECTION -cslow2(double x) { +cslow2 (double x) +{ mynumber u; - double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res; + double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, e1, e2, xx, cor, res; static const double t22 = 6291456.0; int4 k; - y=ABS(x); - u.x = big.x+y; - y = y-(u.x-big.x); - xx=y*y; - s = y*xx*(sn3 +xx*sn5); - c = xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - y1 = (y+t22)-t22; + + y = ABS (x); + u.x = big.x + y; + y = y - (u.x - big.x); + xx = y * y; + s = y * xx * (sn3 + xx * sn5); + c = xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + y1 = (y + t22) - t22; y2 = y - y1; - e1 = (sn+t22)-t22; - e2=(sn-e1)+ssn; - cor=(ccs-cs*c-e1*y2-e2*y)-sn*s; - y=cs-e1*y1; - cor = cor+((cs-y)-e1*y1); - res=y+cor; - cor=(y-res)+cor; - if (res == res+1.0005*cor) + e1 = (sn + t22) - t22; + e2 = (sn - e1) + ssn; + cor = (ccs - cs * c - e1 * y2 - e2 * y) - sn * s; + y = cs - e1 * y1; + cor = cor + ((cs - y) - e1 * y1); + res = y + cor; + cor = (y - res) + cor; + if (res == res + 1.0005 * cor) return res; - else { - y=ABS(x); - __docos(y,0,w); - if (w[0] == w[0]+1.000000005*w[1]) return w[0]; - else return __mpcos(x,0); - } + else + { + y = ABS (x); + __docos (y, 0, w); + if (w[0] == w[0] + 1.000000005 * w[1]) + return w[0]; + else + return __mpcos (x, 0); + } } /***************************************************************************/ @@ -1058,46 +1306,78 @@ cslow2(double x) { static double SECTION -csloww(double x,double dx, double orig) { - static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ - double y,x1,x2,xx,r,t,res,cor,w[2],a,da,xn; - union {int4 i[2]; double x;} v; +csloww (double x, double dx, double orig) +{ + static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ + double y, x1, x2, xx, r, t, res, cor, w[2], a, da, xn; + union + { + int4 i[2]; + double x; + } v; int4 n; - x1=(x+th2_36)-th2_36; - y = aa.x*x1*x1*x1; - r=x+y; - x2=(x-x1)+dx; - xx=x*x; - /* Taylor series */ - t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + bb.x)*xx + 3.0*aa.x*x1*x2)*x +aa.x*x2*x2*x2+dx; - t=((x-r)+y)+t; - res=r+t; - cor = (r-res)+t; - cor = (cor>0)? 1.0005*cor+ABS(orig)*3.1e-30 : 1.0005*cor-ABS(orig)*3.1e-30; - if (res == res + cor) return res; - else { - (x>0)? __dubsin(x,dx,w) : __dubsin(-x,-dx,w); - cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-30 : 1.000000001*w[1] - ABS(orig)*1.1e-30; - if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0]; - else { - t = (orig*hpinv.x + toint.x); - xn = t - toint.x; - v.x = t; - y = (orig - xn*mp1.x) - xn*mp2.x; - n =v.i[LOW_HALF]&3; - da = xn*pp3.x; - t=y-da; - da = (y-t)-da; - y = xn*pp4.x; - 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); - cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-40 : 1.000000001*w[1] - ABS(orig)*1.1e-40; - if (w[0] == w[0]+cor) return (a>0)?w[0]:-w[0]; - else return __mpcos1(orig); + + x1 = (x + th2_36) - th2_36; + y = aa.x * x1 * x1 * x1; + r = x + y; + x2 = (x - x1) + dx; + xx = x * x; + /* Taylor series */ + t = (((((s5.x * xx + s4.x) * xx + s3.x) * xx + s2.x) * xx + bb.x) * xx + + 3.0 * aa.x * x1 * x2) * x + aa.x * x2 * x2 * x2 + dx; + t = ((x - r) + y) + t; + res = r + t; + cor = (r - res) + t; + + if (cor > 0) + cor = 1.0005 * cor + ABS (orig) * 3.1e-30; + else + cor = 1.0005 * cor - ABS (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] + ABS (orig) * 1.1e-30; + else + cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-30; + + if (w[0] == w[0] + cor) + return (x > 0) ? w[0] : -w[0]; + else + { + t = (orig * hpinv.x + toint.x); + xn = t - toint.x; + v.x = t; + y = (orig - xn * mp1.x) - xn * mp2.x; + n = v.i[LOW_HALF] & 3; + da = xn * pp3.x; + t = y - da; + da = (y - t) - da; + y = xn * pp4.x; + 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] + ABS (orig) * 1.1e-40; + else + cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-40; + + if (w[0] == w[0] + cor) + return (a > 0) ? w[0] : -w[0]; + else + return __mpcos1 (orig); + } } - } } /***************************************************************************/ @@ -1109,40 +1389,54 @@ csloww(double x,double dx, double orig) { static double SECTION -csloww1(double x, double dx, double orig) { +csloww1 (double x, double dx, double orig) +{ mynumber u; - double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res; + double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, c1, c2, xx, cor, res; static const double t22 = 6291456.0; int4 k; - y=ABS(x); - u.x=big.x+y; - y=y-(u.x-big.x); - dx=(x>0)?dx:-dx; - xx=y*y; - s = y*xx*(sn3 +xx*sn5); - c = xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - y1 = (y+t22)-t22; - y2 = (y - y1)+dx; - c1 = (cs+t22)-t22; - c2=(cs-c1)+ccs; - cor=(ssn+s*ccs+cs*s+c2*y+c1*y2-sn*y*dx)-sn*c; - y=sn+c1*y1; - cor = cor+((sn-y)+c1*y1); - res=y+cor; - cor=(y-res)+cor; - cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig); - if (res == res + cor) return (x>0)?res:-res; - else { - __dubsin(ABS(x),dx,w); - cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig); - if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0]; - else return __mpcos1(orig); - } + + y = ABS (x); + u.x = big.x + y; + y = y - (u.x - big.x); + dx = (x > 0) ? dx : -dx; + xx = y * y; + s = y * xx * (sn3 + xx * sn5); + c = xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + y1 = (y + t22) - t22; + y2 = (y - y1) + dx; + c1 = (cs + t22) - t22; + c2 = (cs - c1) + ccs; + cor = (ssn + s * ccs + cs * s + c2 * y + c1 * y2 - sn * y * dx) - sn * c; + y = sn + c1 * y1; + cor = cor + ((sn - y) + c1 * y1); + res = y + cor; + cor = (y - res) + cor; + + if (cor > 0) + cor = 1.0005 * cor + 3.1e-30 * ABS (orig); + else + cor = 1.0005 * cor - 3.1e-30 * ABS (orig); + + if (res == res + cor) + return (x > 0) ? res : -res; + else + { + __dubsin (ABS (x), dx, w); + if (w[1] > 0) + cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig); + else + cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig); + if (w[0] == w[0] + cor) + return (x > 0) ? w[0] : -w[0]; + else + return __mpcos1 (orig); + } } @@ -1155,41 +1449,55 @@ csloww1(double x, double dx, double orig) { static double SECTION -csloww2(double x, double dx, double orig, int n) { +csloww2 (double x, double dx, double orig, int n) +{ mynumber u; - double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res; + double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, e1, e2, xx, cor, res; static const double t22 = 6291456.0; int4 k; - y=ABS(x); - u.x=big.x+y; - y=y-(u.x-big.x); - dx=(x>0)?dx:-dx; - xx=y*y; - s = y*xx*(sn3 +xx*sn5); - c = y*dx+xx*(cs2 +xx*(cs4 + xx*cs6)); - k=u.i[LOW_HALF]<<2; - sn=__sincostab.x[k]; - ssn=__sincostab.x[k+1]; - cs=__sincostab.x[k+2]; - ccs=__sincostab.x[k+3]; - - y1 = (y+t22)-t22; - y2 = (y - y1)+dx; - e1 = (sn+t22)-t22; - e2=(sn-e1)+ssn; - cor=(ccs-cs*c-e1*y2-e2*y)-sn*s; - y=cs-e1*y1; - cor = cor+((cs-y)-e1*y1); - res=y+cor; - cor=(y-res)+cor; - cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig); - if (res == res + cor) return (n)?-res:res; - else { - __docos(ABS(x),dx,w); - cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig); - if (w[0] == w[0]+cor) return (n)?-w[0]:w[0]; - else return __mpcos1(orig); - } + + y = ABS (x); + u.x = big.x + y; + y = y - (u.x - big.x); + dx = (x > 0) ? dx : -dx; + xx = y * y; + s = y * xx * (sn3 + xx * sn5); + c = y * dx + xx * (cs2 + xx * (cs4 + xx * cs6)); + k = u.i[LOW_HALF] << 2; + sn = __sincostab.x[k]; + ssn = __sincostab.x[k + 1]; + cs = __sincostab.x[k + 2]; + ccs = __sincostab.x[k + 3]; + + y1 = (y + t22) - t22; + y2 = (y - y1) + dx; + e1 = (sn + t22) - t22; + e2 = (sn - e1) + ssn; + cor = (ccs - cs * c - e1 * y2 - e2 * y) - sn * s; + y = cs - e1 * y1; + cor = cor + ((cs - y) - e1 * y1); + res = y + cor; + cor = (y - res) + cor; + + if (cor > 0) + cor = 1.0005 * cor + 3.1e-30 * ABS (orig); + else + cor = 1.0005 * cor - 3.1e-30 * ABS (orig); + + if (res == res + cor) + return (n) ? -res : res; + else + { + __docos (ABS (x), dx, w); + if (w[1] > 0) + cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig); + else + cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig); + if (w[0] == w[0] + cor) + return (n) ? -w[0] : w[0]; + else + return __mpcos1 (orig); + } } #ifndef __cos diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/e_fmod.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/e_fmod.c index a630d10fe..f686bb670 100644 --- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/e_fmod.c +++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/e_fmod.c @@ -18,6 +18,7 @@ #include <math.h> #include <math_private.h> +#include <stdint.h> static const double one = 1.0, Zero[] = {0.0, -0.0,}; diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/e_log10.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/e_log10.c index 488a0efae..dcb7b58a1 100644 --- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/e_log10.c +++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/e_log10.c @@ -45,6 +45,7 @@ #include <math.h> #include <math_private.h> +#include <stdint.h> static const double two54 = 1.80143985094819840000e+16; /* 0x4350000000000000 */ static const double ivln10 = 4.34294481903251816668e-01; /* 0x3FDBCB7B1526E50E */ diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/math_private.h b/libc/sysdeps/ieee754/dbl-64/wordsize-64/math_private.h index b66085eb1..4f9219934 100644 --- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/math_private.h +++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/math_private.h @@ -1,6 +1,7 @@ #ifndef _MATH_PRIVATE_H_ #include_next <math_private.h> +#include <stdint.h> #ifndef __isnan extern __always_inline int diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c index f25ede8f9..fcf2e6d5b 100644 --- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c +++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c @@ -16,6 +16,7 @@ #include <math.h> #include <math_private.h> +#include <stdint.h> #undef __finite int diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c index 5beccb0ac..914a3c823 100644 --- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c +++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c @@ -32,6 +32,7 @@ #include <math.h> #include <math_private.h> +#include <stdint.h> /* * floor(x) diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c index 70a620cf6..e80b84ca0 100644 --- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c +++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c @@ -17,6 +17,7 @@ #include <math.h> #include <math_private.h> +#include <stdint.h> #undef __isnan int __isnan(double x) diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_modf.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_modf.c index 89743168c..c309e5627 100644 --- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_modf.c +++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_modf.c @@ -22,6 +22,7 @@ #include <math.h> #include <math_private.h> +#include <stdint.h> static const double one = 1.0; diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c index e9ae82bdb..29e62874b 100644 --- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c +++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c @@ -20,7 +20,7 @@ #include <math.h> #include <math_private.h> - +#include <stdint.h> static const double zero = 0.0; diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c index df674670e..bea796083 100644 --- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c +++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c @@ -20,7 +20,7 @@ #include <math.h> #include <math_private.h> - +#include <stdint.h> static const double huge = 1.0e300; diff --git a/libc/sysdeps/ieee754/flt-32/e_gammaf_r.c b/libc/sysdeps/ieee754/flt-32/e_gammaf_r.c index a312957b0..f58f4c805 100644 --- a/libc/sysdeps/ieee754/flt-32/e_gammaf_r.c +++ b/libc/sysdeps/ieee754/flt-32/e_gammaf_r.c @@ -19,14 +19,97 @@ #include <math.h> #include <math_private.h> +#include <float.h> +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's + approximation to gamma function. */ + +static const float gamma_coeff[] = + { + 0x1.555556p-4f, + -0xb.60b61p-12f, + 0x3.403404p-12f, + }; + +#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0])) + +/* Return gamma (X), for positive X less than 42, in the form R * + 2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to + avoid overflow or underflow in intermediate calculations. */ + +static float +gammaf_positive (float x, int *exp2_adj) +{ + int local_signgam; + if (x < 0.5f) + { + *exp2_adj = 0; + return __ieee754_expf (__ieee754_lgammaf_r (x + 1, &local_signgam)) / x; + } + else if (x <= 1.5f) + { + *exp2_adj = 0; + return __ieee754_expf (__ieee754_lgammaf_r (x, &local_signgam)); + } + else if (x < 2.5f) + { + *exp2_adj = 0; + float x_adj = x - 1; + return (__ieee754_expf (__ieee754_lgammaf_r (x_adj, &local_signgam)) + * x_adj); + } + else + { + float eps = 0; + float x_eps = 0; + float x_adj = x; + float prod = 1; + if (x < 4.0f) + { + /* Adjust into the range for applying Stirling's + approximation. */ + float n = __ceilf (4.0f - x); +#if FLT_EVAL_METHOD != 0 + volatile +#endif + float x_tmp = x + n; + x_adj = x_tmp; + x_eps = (x - (x_adj - n)); + prod = __gamma_productf (x_adj - n, x_eps, n, &eps); + } + /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)). + Compute gamma (X_ADJ + X_EPS) using Stirling's approximation, + starting by computing pow (X_ADJ, X_ADJ) with a power of 2 + factored out. */ + float exp_adj = -eps; + float x_adj_int = __roundf (x_adj); + float x_adj_frac = x_adj - x_adj_int; + int x_adj_log2; + float x_adj_mant = __frexpf (x_adj, &x_adj_log2); + if (x_adj_mant < (float) M_SQRT1_2) + { + x_adj_log2--; + x_adj_mant *= 2.0f; + } + *exp2_adj = x_adj_log2 * (int) x_adj_int; + float ret = (__ieee754_powf (x_adj_mant, x_adj) + * __ieee754_exp2f (x_adj_log2 * x_adj_frac) + * __ieee754_expf (-x_adj) + * __ieee754_sqrtf (2 * (float) M_PI / x_adj) + / prod); + exp_adj += x_eps * __ieee754_logf (x); + float bsum = gamma_coeff[NCOEFF - 1]; + float x_adj2 = x_adj * x_adj; + for (size_t i = 1; i <= NCOEFF - 1; i++) + bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i]; + exp_adj += bsum / x_adj; + return ret + ret * __expm1f (exp_adj); + } +} float __ieee754_gammaf_r (float x, int *signgamp) { - /* We don't have a real gamma implementation now. We'll use lgamma - and the exp function. But due to the required boundary - conditions we must check some values separately. */ int32_t hx; GET_FLOAT_WORD (hx, x); @@ -50,8 +133,49 @@ __ieee754_gammaf_r (float x, int *signgamp) *signgamp = 0; return x - x; } + if (__builtin_expect ((hx & 0x7f800000) == 0x7f800000, 0)) + { + /* Positive infinity (return positive infinity) or NaN (return + NaN). */ + *signgamp = 0; + return x + x; + } - /* XXX FIXME. */ - return __ieee754_expf (__ieee754_lgammaf_r (x, signgamp)); + if (x >= 36.0f) + { + /* Overflow. */ + *signgamp = 0; + return FLT_MAX * FLT_MAX; + } + else if (x > 0.0f) + { + *signgamp = 0; + int exp2_adj; + float ret = gammaf_positive (x, &exp2_adj); + return __scalbnf (ret, exp2_adj); + } + else if (x >= -FLT_EPSILON / 4.0f) + { + *signgamp = 0; + return 1.0f / x; + } + else + { + float tx = __truncf (x); + *signgamp = (tx == 2.0f * __truncf (tx / 2.0f)) ? -1 : 1; + if (x <= -42.0f) + /* Underflow. */ + return FLT_MIN * FLT_MIN; + float frac = tx - x; + if (frac > 0.5f) + frac = 1.0f - frac; + float sinpix = (frac <= 0.25f + ? __sinf ((float) M_PI * frac) + : __cosf ((float) M_PI * (0.5f - frac))); + int exp2_adj; + float ret = (float) M_PI / (-x * sinpix + * gammaf_positive (-x, &exp2_adj)); + return __scalbnf (ret, -exp2_adj); + } } strong_alias (__ieee754_gammaf_r, __gammaf_r_finite) diff --git a/libc/sysdeps/ieee754/k_standard.c b/libc/sysdeps/ieee754/k_standard.c index cd3123046..150921f90 100644 --- a/libc/sysdeps/ieee754/k_standard.c +++ b/libc/sysdeps/ieee754/k_standard.c @@ -837,7 +837,7 @@ __kernel_standard(double x, double y, int type) exc.type = OVERFLOW; exc.name = type < 100 ? "tgamma" : (type < 200 ? "tgammaf" : "tgammal"); - exc.retval = HUGE_VAL; + exc.retval = __copysign (HUGE_VAL, x); if (_LIB_VERSION == _POSIX_) __set_errno (ERANGE); else if (!matherr(&exc)) { diff --git a/libc/sysdeps/ieee754/ldbl-128/e_gammal_r.c b/libc/sysdeps/ieee754/ldbl-128/e_gammal_r.c index b6da31c13..e8d49e987 100644 --- a/libc/sysdeps/ieee754/ldbl-128/e_gammal_r.c +++ b/libc/sysdeps/ieee754/ldbl-128/e_gammal_r.c @@ -20,14 +20,108 @@ #include <math.h> #include <math_private.h> +#include <float.h> +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's + approximation to gamma function. */ + +static const long double gamma_coeff[] = + { + 0x1.5555555555555555555555555555p-4L, + -0xb.60b60b60b60b60b60b60b60b60b8p-12L, + 0x3.4034034034034034034034034034p-12L, + -0x2.7027027027027027027027027028p-12L, + 0x3.72a3c5631fe46ae1d4e700dca8f2p-12L, + -0x7.daac36664f1f207daac36664f1f4p-12L, + 0x1.a41a41a41a41a41a41a41a41a41ap-8L, + -0x7.90a1b2c3d4e5f708192a3b4c5d7p-8L, + 0x2.dfd2c703c0cfff430edfd2c703cp-4L, + -0x1.6476701181f39edbdb9ce625987dp+0L, + 0xd.672219167002d3a7a9c886459cp+0L, + -0x9.cd9292e6660d55b3f712eb9e07c8p+4L, + 0x8.911a740da740da740da740da741p+8L, + -0x8.d0cc570e255bf59ff6eec24b49p+12L, + }; + +#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0])) + +/* Return gamma (X), for positive X less than 1775, in the form R * + 2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to + avoid overflow or underflow in intermediate calculations. */ + +static long double +gammal_positive (long double x, int *exp2_adj) +{ + int local_signgam; + if (x < 0.5L) + { + *exp2_adj = 0; + return __ieee754_expl (__ieee754_lgammal_r (x + 1, &local_signgam)) / x; + } + else if (x <= 1.5L) + { + *exp2_adj = 0; + return __ieee754_expl (__ieee754_lgammal_r (x, &local_signgam)); + } + else if (x < 12.5L) + { + /* Adjust into the range for using exp (lgamma). */ + *exp2_adj = 0; + long double n = __ceill (x - 1.5L); + long double x_adj = x - n; + long double eps; + long double prod = __gamma_productl (x_adj, 0, n, &eps); + return (__ieee754_expl (__ieee754_lgammal_r (x_adj, &local_signgam)) + * prod * (1.0L + eps)); + } + else + { + long double eps = 0; + long double x_eps = 0; + long double x_adj = x; + long double prod = 1; + if (x < 24.0L) + { + /* Adjust into the range for applying Stirling's + approximation. */ + long double n = __ceill (24.0L - x); + x_adj = x + n; + x_eps = (x - (x_adj - n)); + prod = __gamma_productl (x_adj - n, x_eps, n, &eps); + } + /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)). + Compute gamma (X_ADJ + X_EPS) using Stirling's approximation, + starting by computing pow (X_ADJ, X_ADJ) with a power of 2 + factored out. */ + long double exp_adj = -eps; + long double x_adj_int = __roundl (x_adj); + long double x_adj_frac = x_adj - x_adj_int; + int x_adj_log2; + long double x_adj_mant = __frexpl (x_adj, &x_adj_log2); + if (x_adj_mant < M_SQRT1_2l) + { + x_adj_log2--; + x_adj_mant *= 2.0L; + } + *exp2_adj = x_adj_log2 * (int) x_adj_int; + long double ret = (__ieee754_powl (x_adj_mant, x_adj) + * __ieee754_exp2l (x_adj_log2 * x_adj_frac) + * __ieee754_expl (-x_adj) + * __ieee754_sqrtl (2 * M_PIl / x_adj) + / prod); + exp_adj += x_eps * __ieee754_logl (x); + long double bsum = gamma_coeff[NCOEFF - 1]; + long double x_adj2 = x_adj * x_adj; + for (size_t i = 1; i <= NCOEFF - 1; i++) + bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i]; + exp_adj += bsum / x_adj; + return ret + ret * __expm1l (exp_adj); + } +} long double __ieee754_gammal_r (long double x, int *signgamp) { - /* We don't have a real gamma implementation now. We'll use lgamma - and the exp function. But due to the required boundary - conditions we must check some values separately. */ int64_t hx; u_int64_t lx; @@ -51,8 +145,49 @@ __ieee754_gammal_r (long double x, int *signgamp) *signgamp = 0; return x - x; } + if ((hx & 0x7fff000000000000ULL) == 0x7fff000000000000ULL) + { + /* Positive infinity (return positive infinity) or NaN (return + NaN). */ + *signgamp = 0; + return x + x; + } - /* XXX FIXME. */ - return __ieee754_expl (__ieee754_lgammal_r (x, signgamp)); + if (x >= 1756.0L) + { + /* Overflow. */ + *signgamp = 0; + return LDBL_MAX * LDBL_MAX; + } + else if (x > 0.0L) + { + *signgamp = 0; + int exp2_adj; + long double ret = gammal_positive (x, &exp2_adj); + return __scalbnl (ret, exp2_adj); + } + else if (x >= -LDBL_EPSILON / 4.0L) + { + *signgamp = 0; + return 1.0f / x; + } + else + { + long double tx = __truncl (x); + *signgamp = (tx == 2.0L * __truncl (tx / 2.0L)) ? -1 : 1; + if (x <= -1775.0L) + /* Underflow. */ + return LDBL_MIN * LDBL_MIN; + long double frac = tx - x; + if (frac > 0.5L) + frac = 1.0L - frac; + long double sinpix = (frac <= 0.25L + ? __sinl (M_PIl * frac) + : __cosl (M_PIl * (0.5L - frac))); + int exp2_adj; + long double ret = M_PIl / (-x * sinpix + * gammal_positive (-x, &exp2_adj)); + return __scalbnl (ret, -exp2_adj); + } } strong_alias (__ieee754_gammal_r, __gammal_r_finite) diff --git a/libc/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c b/libc/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c index 84846fdc3..ee856acd5 100644 --- a/libc/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c +++ b/libc/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c @@ -184,13 +184,13 @@ static const int32_t two_over_pi[] = { }; static const long double c[] = { -/* 93 bits of pi/2 */ +/* 113 bits of pi/2 */ #define PI_2_1 c[0] - 1.57079632679489661923132169155131424e+00L, /* 3fff921fb54442d18469898cc5100000 */ + 0x1.921fb54442d18469898cc51701b8p+0L, /* pi/2 - PI_2_1 */ #define PI_2_1t c[1] - 8.84372056613570112025531863263659260e-29L, /* 3fa1c06e0e68948127044533e63a0106 */ + 0x3.9a252049c1114cf98e804177d4c8p-116L, }; int32_t __ieee754_rem_pio2l(long double x, long double *y) @@ -213,7 +213,7 @@ int32_t __ieee754_rem_pio2l(long double x, long double *y) { if (hx > 0) { - /* 113 + 93 bit PI is ok */ + /* 113 + 113 bit PI is ok */ z = x - PI_2_1; y[0] = z - PI_2_1t; y[1] = (z - y[0]) - PI_2_1t; @@ -221,7 +221,7 @@ int32_t __ieee754_rem_pio2l(long double x, long double *y) } else { - /* 113 + 93 bit PI is ok */ + /* 113 + 113 bit PI is ok */ z = x + PI_2_1; y[0] = z + PI_2_1t; y[1] = (z - y[0]) + PI_2_1t; diff --git a/libc/sysdeps/ieee754/ldbl-128/gamma_productl.c b/libc/sysdeps/ieee754/ldbl-128/gamma_productl.c new file mode 100644 index 000000000..157dbab9f --- /dev/null +++ b/libc/sysdeps/ieee754/ldbl-128/gamma_productl.c @@ -0,0 +1,75 @@ +/* Compute a product of X, X+1, ..., with an error estimate. + Copyright (C) 2013 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <float.h> + +/* Calculate X * Y exactly and store the result in *HI + *LO. It is + given that the values are small enough that no overflow occurs and + large enough (or zero) that no underflow occurs. */ + +static inline void +mul_split (long double *hi, long double *lo, long double x, long double y) +{ +#ifdef __FP_FAST_FMAL + /* Fast built-in fused multiply-add. */ + *hi = x * y; + *lo = __builtin_fmal (x, y, -*hi); +#elif defined FP_FAST_FMAL + /* Fast library fused multiply-add, compiler before GCC 4.6. */ + *hi = x * y; + *lo = __fmal (x, y, -*hi); +#else + /* Apply Dekker's algorithm. */ + *hi = x * y; +# define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1) + long double x1 = x * C; + long double y1 = y * C; +# undef C + x1 = (x - x1) + x1; + y1 = (y - y1) + y1; + long double x2 = x - x1; + long double y2 = y - y1; + *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2; +#endif +} + +/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N + - 1, in the form R * (1 + *EPS) where the return value R is an + approximation to the product and *EPS is set to indicate the + approximate error in the return value. X is such that all the + values X + 1, ..., X + N - 1 are exactly representable, and X_EPS / + X is small enough that factors quadratic in it can be + neglected. */ + +long double +__gamma_productl (long double x, long double x_eps, int n, long double *eps) +{ + SET_RESTORE_ROUNDL (FE_TONEAREST); + long double ret = x; + *eps = x_eps / x; + for (int i = 1; i < n; i++) + { + *eps += x_eps / (x + i); + long double lo; + mul_split (&ret, &lo, ret, x + i); + *eps += lo / ret; + } + return ret; +} diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c index 52ade9e4a..90d8e3f0d 100644 --- a/libc/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c +++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c @@ -20,14 +20,107 @@ #include <math.h> #include <math_private.h> +#include <float.h> +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's + approximation to gamma function. */ + +static const long double gamma_coeff[] = + { + 0x1.555555555555555555555555558p-4L, + -0xb.60b60b60b60b60b60b60b60b6p-12L, + 0x3.4034034034034034034034034p-12L, + -0x2.7027027027027027027027027p-12L, + 0x3.72a3c5631fe46ae1d4e700dca9p-12L, + -0x7.daac36664f1f207daac36664f2p-12L, + 0x1.a41a41a41a41a41a41a41a41a4p-8L, + -0x7.90a1b2c3d4e5f708192a3b4c5ep-8L, + 0x2.dfd2c703c0cfff430edfd2c704p-4L, + -0x1.6476701181f39edbdb9ce625988p+0L, + 0xd.672219167002d3a7a9c886459cp+0L, + -0x9.cd9292e6660d55b3f712eb9e08p+4L, + 0x8.911a740da740da740da740da74p+8L, + }; + +#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0])) + +/* Return gamma (X), for positive X less than 191, in the form R * + 2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to + avoid overflow or underflow in intermediate calculations. */ + +static long double +gammal_positive (long double x, int *exp2_adj) +{ + int local_signgam; + if (x < 0.5L) + { + *exp2_adj = 0; + return __ieee754_expl (__ieee754_lgammal_r (x + 1, &local_signgam)) / x; + } + else if (x <= 1.5L) + { + *exp2_adj = 0; + return __ieee754_expl (__ieee754_lgammal_r (x, &local_signgam)); + } + else if (x < 11.5L) + { + /* Adjust into the range for using exp (lgamma). */ + *exp2_adj = 0; + long double n = __ceill (x - 1.5L); + long double x_adj = x - n; + long double eps; + long double prod = __gamma_productl (x_adj, 0, n, &eps); + return (__ieee754_expl (__ieee754_lgammal_r (x_adj, &local_signgam)) + * prod * (1.0L + eps)); + } + else + { + long double eps = 0; + long double x_eps = 0; + long double x_adj = x; + long double prod = 1; + if (x < 23.0L) + { + /* Adjust into the range for applying Stirling's + approximation. */ + long double n = __ceill (23.0L - x); + x_adj = x + n; + x_eps = (x - (x_adj - n)); + prod = __gamma_productl (x_adj - n, x_eps, n, &eps); + } + /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)). + Compute gamma (X_ADJ + X_EPS) using Stirling's approximation, + starting by computing pow (X_ADJ, X_ADJ) with a power of 2 + factored out. */ + long double exp_adj = -eps; + long double x_adj_int = __roundl (x_adj); + long double x_adj_frac = x_adj - x_adj_int; + int x_adj_log2; + long double x_adj_mant = __frexpl (x_adj, &x_adj_log2); + if (x_adj_mant < M_SQRT1_2l) + { + x_adj_log2--; + x_adj_mant *= 2.0L; + } + *exp2_adj = x_adj_log2 * (int) x_adj_int; + long double ret = (__ieee754_powl (x_adj_mant, x_adj) + * __ieee754_exp2l (x_adj_log2 * x_adj_frac) + * __ieee754_expl (-x_adj) + * __ieee754_sqrtl (2 * M_PIl / x_adj) + / prod); + exp_adj += x_eps * __ieee754_logl (x); + long double bsum = gamma_coeff[NCOEFF - 1]; + long double x_adj2 = x_adj * x_adj; + for (size_t i = 1; i <= NCOEFF - 1; i++) + bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i]; + exp_adj += bsum / x_adj; + return ret + ret * __expm1l (exp_adj); + } +} long double __ieee754_gammal_r (long double x, int *signgamp) { - /* We don't have a real gamma implementation now. We'll use lgamma - and the exp function. But due to the required boundary - conditions we must check some values separately. */ int64_t hx; u_int64_t lx; @@ -51,8 +144,49 @@ __ieee754_gammal_r (long double x, int *signgamp) *signgamp = 0; return x - x; } + if ((hx & 0x7ff0000000000000ULL) == 0x7ff0000000000000ULL) + { + /* Positive infinity (return positive infinity) or NaN (return + NaN). */ + *signgamp = 0; + return x + x; + } - /* XXX FIXME. */ - return __ieee754_expl (__ieee754_lgammal_r (x, signgamp)); + if (x >= 172.0L) + { + /* Overflow. */ + *signgamp = 0; + return LDBL_MAX * LDBL_MAX; + } + else if (x > 0.0L) + { + *signgamp = 0; + int exp2_adj; + long double ret = gammal_positive (x, &exp2_adj); + return __scalbnl (ret, exp2_adj); + } + else if (x >= -0x1p-110L) + { + *signgamp = 0; + return 1.0f / x; + } + else + { + long double tx = __truncl (x); + *signgamp = (tx == 2.0L * __truncl (tx / 2.0L)) ? -1 : 1; + if (x <= -191.0L) + /* Underflow. */ + return LDBL_MIN * LDBL_MIN; + long double frac = tx - x; + if (frac > 0.5L) + frac = 1.0L - frac; + long double sinpix = (frac <= 0.25L + ? __sinl (M_PIl * frac) + : __cosl (M_PIl * (0.5L - frac))); + int exp2_adj; + long double ret = M_PIl / (-x * sinpix + * gammal_positive (-x, &exp2_adj)); + return __scalbnl (ret, -exp2_adj); + } } strong_alias (__ieee754_gammal_r, __gammal_r_finite) diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c index 692ae2493..6a72d6a85 100644 --- a/libc/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c +++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c @@ -185,13 +185,13 @@ static const int32_t two_over_pi[] = { }; static const long double c[] = { -/* 93 bits of pi/2 */ +/* 106 bits of pi/2 */ #define PI_2_1 c[0] - 1.57079632679489661923132169155131424e+00L, /* 3fff921fb54442d18469898cc5100000 */ + 0x1.921fb54442d18469898cc517018p+0L, /* pi/2 - PI_2_1 */ #define PI_2_1t c[1] - 8.84372056613570112025531863263659260e-29L, /* 3fa1c06e0e68948127044533e63a0106 */ + 0x3.839a252049c1114cf98e804178p-108L, }; int32_t __ieee754_rem_pio2l(long double x, long double *y) @@ -216,7 +216,7 @@ int32_t __ieee754_rem_pio2l(long double x, long double *y) { if (hx > 0) { - /* 113 + 93 bit PI is ok */ + /* 106 + 106 bit PI is ok */ z = x - PI_2_1; y[0] = z - PI_2_1t; y[1] = (z - y[0]) - PI_2_1t; @@ -224,7 +224,7 @@ int32_t __ieee754_rem_pio2l(long double x, long double *y) } else { - /* 113 + 93 bit PI is ok */ + /* 106 + 106 bit PI is ok */ z = x + PI_2_1; y[0] = z + PI_2_1t; y[1] = (z - y[0]) + PI_2_1t; diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/gamma_productl.c b/libc/sysdeps/ieee754/ldbl-128ibm/gamma_productl.c new file mode 100644 index 000000000..7c6186d23 --- /dev/null +++ b/libc/sysdeps/ieee754/ldbl-128ibm/gamma_productl.c @@ -0,0 +1,42 @@ +/* Compute a product of X, X+1, ..., with an error estimate. + Copyright (C) 2013 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> + +/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N + - 1, in the form R * (1 + *EPS) where the return value R is an + approximation to the product and *EPS is set to indicate the + approximate error in the return value. X is such that all the + values X + 1, ..., X + N - 1 are exactly representable, and X_EPS / + X is small enough that factors quadratic in it can be + neglected. */ + +long double +__gamma_productl (long double x, long double x_eps, int n, long double *eps) +{ + long double ret = x; + *eps = x_eps / x; + for (int i = 1; i < n; i++) + { + *eps += x_eps / (x + i); + ret *= x + i; + /* FIXME: no error estimates for the multiplication. */ + } + return ret; +} diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h b/libc/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h index 1cce1fc4d..58eb57cd6 100644 --- a/libc/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h +++ b/libc/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h @@ -4,6 +4,7 @@ #include <sysdeps/ieee754/ldbl-128/math_ldbl.h> #include <ieee754.h> +#include <stdint.h> static inline void ldbl_extract_mantissa (int64_t *hi64, uint64_t *lo64, int *exp, long double x) diff --git a/libc/sysdeps/ieee754/ldbl-96/e_gammal_r.c b/libc/sysdeps/ieee754/ldbl-96/e_gammal_r.c index 0974351a1..7cb3e8563 100644 --- a/libc/sysdeps/ieee754/ldbl-96/e_gammal_r.c +++ b/libc/sysdeps/ieee754/ldbl-96/e_gammal_r.c @@ -19,14 +19,102 @@ #include <math.h> #include <math_private.h> +#include <float.h> +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's + approximation to gamma function. */ + +static const long double gamma_coeff[] = + { + 0x1.5555555555555556p-4L, + -0xb.60b60b60b60b60bp-12L, + 0x3.4034034034034034p-12L, + -0x2.7027027027027028p-12L, + 0x3.72a3c5631fe46aep-12L, + -0x7.daac36664f1f208p-12L, + 0x1.a41a41a41a41a41ap-8L, + -0x7.90a1b2c3d4e5f708p-8L, + }; + +#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0])) + +/* Return gamma (X), for positive X less than 1766, in the form R * + 2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to + avoid overflow or underflow in intermediate calculations. */ + +static long double +gammal_positive (long double x, int *exp2_adj) +{ + int local_signgam; + if (x < 0.5L) + { + *exp2_adj = 0; + return __ieee754_expl (__ieee754_lgammal_r (x + 1, &local_signgam)) / x; + } + else if (x <= 1.5L) + { + *exp2_adj = 0; + return __ieee754_expl (__ieee754_lgammal_r (x, &local_signgam)); + } + else if (x < 7.5L) + { + /* Adjust into the range for using exp (lgamma). */ + *exp2_adj = 0; + long double n = __ceill (x - 1.5L); + long double x_adj = x - n; + long double eps; + long double prod = __gamma_productl (x_adj, 0, n, &eps); + return (__ieee754_expl (__ieee754_lgammal_r (x_adj, &local_signgam)) + * prod * (1.0L + eps)); + } + else + { + long double eps = 0; + long double x_eps = 0; + long double x_adj = x; + long double prod = 1; + if (x < 13.0L) + { + /* Adjust into the range for applying Stirling's + approximation. */ + long double n = __ceill (13.0L - x); + x_adj = x + n; + x_eps = (x - (x_adj - n)); + prod = __gamma_productl (x_adj - n, x_eps, n, &eps); + } + /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)). + Compute gamma (X_ADJ + X_EPS) using Stirling's approximation, + starting by computing pow (X_ADJ, X_ADJ) with a power of 2 + factored out. */ + long double exp_adj = -eps; + long double x_adj_int = __roundl (x_adj); + long double x_adj_frac = x_adj - x_adj_int; + int x_adj_log2; + long double x_adj_mant = __frexpl (x_adj, &x_adj_log2); + if (x_adj_mant < M_SQRT1_2l) + { + x_adj_log2--; + x_adj_mant *= 2.0L; + } + *exp2_adj = x_adj_log2 * (int) x_adj_int; + long double ret = (__ieee754_powl (x_adj_mant, x_adj) + * __ieee754_exp2l (x_adj_log2 * x_adj_frac) + * __ieee754_expl (-x_adj) + * __ieee754_sqrtl (2 * M_PIl / x_adj) + / prod); + exp_adj += x_eps * __ieee754_logl (x); + long double bsum = gamma_coeff[NCOEFF - 1]; + long double x_adj2 = x_adj * x_adj; + for (size_t i = 1; i <= NCOEFF - 1; i++) + bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i]; + exp_adj += bsum / x_adj; + return ret + ret * __expm1l (exp_adj); + } +} long double __ieee754_gammal_r (long double x, int *signgamp) { - /* We don't have a real gamma implementation now. We'll use lgamma - and the exp function. But due to the required boundary - conditions we must check some values separately. */ u_int32_t es, hx, lx; GET_LDOUBLE_WORDS (es, hx, lx, x); @@ -43,22 +131,55 @@ __ieee754_gammal_r (long double x, int *signgamp) *signgamp = 0; return x - x; } - if (__builtin_expect ((es & 0x7fff) == 0x7fff, 0) - && ((hx & 0x7fffffff) | lx) != 0) + if (__builtin_expect ((es & 0x7fff) == 0x7fff, 0)) { - /* NaN, return it. */ + /* Positive infinity (return positive infinity) or NaN (return + NaN). */ *signgamp = 0; - return x; + return x + x; } - if (__builtin_expect ((es & 0x8000) != 0, 0) - && x < 0xffffffff && __rintl (x) == x) + if (__builtin_expect ((es & 0x8000) != 0, 0) && __rintl (x) == x) { /* Return value for integer x < 0 is NaN with invalid exception. */ *signgamp = 0; return (x - x) / (x - x); } - /* XXX FIXME. */ - return __ieee754_expl (__ieee754_lgammal_r (x, signgamp)); + if (x >= 1756.0L) + { + /* Overflow. */ + *signgamp = 0; + return LDBL_MAX * LDBL_MAX; + } + else if (x > 0.0L) + { + *signgamp = 0; + int exp2_adj; + long double ret = gammal_positive (x, &exp2_adj); + return __scalbnl (ret, exp2_adj); + } + else if (x >= -LDBL_EPSILON / 4.0L) + { + *signgamp = 0; + return 1.0f / x; + } + else + { + long double tx = __truncl (x); + *signgamp = (tx == 2.0L * __truncl (tx / 2.0L)) ? -1 : 1; + if (x <= -1766.0L) + /* Underflow. */ + return LDBL_MIN * LDBL_MIN; + long double frac = tx - x; + if (frac > 0.5L) + frac = 1.0L - frac; + long double sinpix = (frac <= 0.25L + ? __sinl (M_PIl * frac) + : __cosl (M_PIl * (0.5L - frac))); + int exp2_adj; + long double ret = M_PIl / (-x * sinpix + * gammal_positive (-x, &exp2_adj)); + return __scalbnl (ret, -exp2_adj); + } } strong_alias (__ieee754_gammal_r, __gammal_r_finite) diff --git a/libc/sysdeps/ieee754/ldbl-96/gamma_product.c b/libc/sysdeps/ieee754/ldbl-96/gamma_product.c new file mode 100644 index 000000000..d464e7084 --- /dev/null +++ b/libc/sysdeps/ieee754/ldbl-96/gamma_product.c @@ -0,0 +1,46 @@ +/* Compute a product of X, X+1, ..., with an error estimate. + Copyright (C) 2013 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <float.h> + +/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N + - 1, in the form R * (1 + *EPS) where the return value R is an + approximation to the product and *EPS is set to indicate the + approximate error in the return value. X is such that all the + values X + 1, ..., X + N - 1 are exactly representable, and X_EPS / + X is small enough that factors quadratic in it can be + neglected. */ + +double +__gamma_product (double x, double x_eps, int n, double *eps) +{ + long double x_full = (long double) x + (long double) x_eps; + long double ret = x_full; + for (int i = 1; i < n; i++) + ret *= x_full + i; + +#if FLT_EVAL_METHOD != 0 + volatile +#endif + double fret = ret; + *eps = (ret - fret) / fret; + + return fret; +} diff --git a/libc/sysdeps/ieee754/ldbl-96/gamma_productl.c b/libc/sysdeps/ieee754/ldbl-96/gamma_productl.c new file mode 100644 index 000000000..157dbab9f --- /dev/null +++ b/libc/sysdeps/ieee754/ldbl-96/gamma_productl.c @@ -0,0 +1,75 @@ +/* Compute a product of X, X+1, ..., with an error estimate. + Copyright (C) 2013 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <float.h> + +/* Calculate X * Y exactly and store the result in *HI + *LO. It is + given that the values are small enough that no overflow occurs and + large enough (or zero) that no underflow occurs. */ + +static inline void +mul_split (long double *hi, long double *lo, long double x, long double y) +{ +#ifdef __FP_FAST_FMAL + /* Fast built-in fused multiply-add. */ + *hi = x * y; + *lo = __builtin_fmal (x, y, -*hi); +#elif defined FP_FAST_FMAL + /* Fast library fused multiply-add, compiler before GCC 4.6. */ + *hi = x * y; + *lo = __fmal (x, y, -*hi); +#else + /* Apply Dekker's algorithm. */ + *hi = x * y; +# define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1) + long double x1 = x * C; + long double y1 = y * C; +# undef C + x1 = (x - x1) + x1; + y1 = (y - y1) + y1; + long double x2 = x - x1; + long double y2 = y - y1; + *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2; +#endif +} + +/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N + - 1, in the form R * (1 + *EPS) where the return value R is an + approximation to the product and *EPS is set to indicate the + approximate error in the return value. X is such that all the + values X + 1, ..., X + N - 1 are exactly representable, and X_EPS / + X is small enough that factors quadratic in it can be + neglected. */ + +long double +__gamma_productl (long double x, long double x_eps, int n, long double *eps) +{ + SET_RESTORE_ROUNDL (FE_TONEAREST); + long double ret = x; + *eps = x_eps / x; + for (int i = 1; i < n; i++) + { + *eps += x_eps / (x + i); + long double lo; + mul_split (&ret, &lo, ret, x + i); + *eps += lo / ret; + } + return ret; +} diff --git a/libc/sysdeps/ieee754/s_lib_version.c b/libc/sysdeps/ieee754/s_lib_version.c index a377ab1f7..7abb3e07a 100644 --- a/libc/sysdeps/ieee754/s_lib_version.c +++ b/libc/sysdeps/ieee754/s_lib_version.c @@ -25,16 +25,17 @@ static char rcsid[] = "$NetBSD: s_lib_version.c,v 1.6 1995/05/10 20:47:44 jtc Ex * define and initialize _LIB_VERSION */ #ifdef _POSIX_MODE -_LIB_VERSION_TYPE _LIB_VERSION = _POSIX_; +_LIB_VERSION_TYPE _LIB_VERSION_INTERNAL = _POSIX_; #else #ifdef _XOPEN_MODE -_LIB_VERSION_TYPE _LIB_VERSION = _XOPEN_; +_LIB_VERSION_TYPE _LIB_VERSION_INTERNAL = _XOPEN_; #else #ifdef _SVID3_MODE -_LIB_VERSION_TYPE _LIB_VERSION = _SVID_; +_LIB_VERSION_TYPE _LIB_VERSION_INTERNAL = _SVID_; #else /* default _IEEE_MODE */ -_LIB_VERSION_TYPE _LIB_VERSION = _IEEE_; +_LIB_VERSION_TYPE _LIB_VERSION_INTERNAL = _IEEE_; #endif #endif #endif +weak_alias (_LIB_VERSION_INTERNAL, _LIB_VERSION) |