diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_acosh.c | 29 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_asin.c | 30 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_atan2.c | 374 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_atanh.c | 122 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_cosh.c | 60 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_exp2.c | 4 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_fmod.c | 43 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_gamma_r.c | 10 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_hypot.c | 25 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_j0.c | 212 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_j1.c | 255 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_jn.c | 82 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_lgamma_r.c | 85 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_log.c | 20 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_log10.c | 42 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_log2.c | 34 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_pow.c | 9 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_remainder.c | 11 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_sinh.c | 36 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_sqrt.c | 3 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/halfulp.c | 18 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_asinh.c | 40 |
22 files changed, 616 insertions, 928 deletions
diff --git a/sysdeps/ieee754/dbl-64/e_acosh.c b/sysdeps/ieee754/dbl-64/e_acosh.c index 27c29cd8c9..f474e9ab5c 100644 --- a/sysdeps/ieee754/dbl-64/e_acosh.c +++ b/sysdeps/ieee754/dbl-64/e_acosh.c @@ -5,18 +5,14 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_acosh.c,v 1.9 1995/05/12 04:57:18 jtc Exp $"; -#endif - /* __ieee754_acosh(x) * Method : - * Based on + * Based on * acosh(x) = log [ x + sqrt(x*x-1) ] * we have * acosh(x) := log(x)+ln2, if x is large; else @@ -31,21 +27,13 @@ static char rcsid[] = "$NetBSD: e_acosh.c,v 1.9 1995/05/12 04:57:18 jtc Exp $"; #include "math.h" #include "math_private.h" -#ifdef __STDC__ -static const double -#else -static double -#endif +static const double one = 1.0, ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */ -#ifdef __STDC__ - double __ieee754_acosh(double x) -#else - double __ieee754_acosh(x) - double x; -#endif -{ +double +__ieee754_acosh(double x) +{ double t; int32_t hx; u_int32_t lx; @@ -54,8 +42,8 @@ ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */ return (x-x)/(x-x); } else if(hx >=0x41b00000) { /* x > 2**28 */ if(hx >=0x7ff00000) { /* x is inf of NaN */ - return x+x; - } else + return x+x; + } else return __ieee754_log(x)+ln2; /* acosh(huge)=log(2x) */ } else if(((hx-0x3ff00000)|lx)==0) { return 0.0; /* acosh(1) = 0 */ @@ -67,3 +55,4 @@ ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */ return __log1p(t+__sqrt(2.0*t+t*t)); } } +strong_alias (__ieee754_acosh, __acosh_finite) diff --git a/sysdeps/ieee754/dbl-64/e_asin.c b/sysdeps/ieee754/dbl-64/e_asin.c index ce5d227b71..02efb7ad2e 100644 --- a/sysdeps/ieee754/dbl-64/e_asin.c +++ b/sysdeps/ieee754/dbl-64/e_asin.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -209,9 +209,9 @@ double __ieee754_asin(double x){ else xx = -x - asncs.x[n]; t = asncs.x[n+1]*xx; p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+ - xx*(asncs.x[n+5]+xx*(asncs.x[n+6] + xx*(asncs.x[n+5]+xx*(asncs.x[n+6] +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+ - xx*asncs.x[n+9])))))))+asncs.x[n+10]; + xx*asncs.x[n+9])))))))+asncs.x[n+10]; t+=p; res =asncs.x[n+11] +t; cor = (asncs.x[n+11]-res)+t; @@ -248,9 +248,9 @@ double __ieee754_asin(double x){ else xx = -x - asncs.x[n]; t = asncs.x[n+1]*xx; p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+ - xx*(asncs.x[n+5]+xx*(asncs.x[n+6] + xx*(asncs.x[n+5]+xx*(asncs.x[n+6] +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+ - xx*(asncs.x[n+9]+xx*asncs.x[n+10]))))))))+asncs.x[n+11]; + xx*(asncs.x[n+9]+xx*asncs.x[n+10]))))))))+asncs.x[n+11]; t+=p; res =asncs.x[n+12] +t; cor = (asncs.x[n+12]-res)+t; @@ -324,6 +324,7 @@ double __ieee754_asin(double x){ return u.x/v.x; /* NaN */ } } +strong_alias (__ieee754_asin, __asin_finite) /*******************************************************************/ /* */ @@ -397,7 +398,7 @@ double __ieee754_acos(double x) else xx = -x - asncs.x[n]; t = asncs.x[n+1]*xx; p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+ - xx*(asncs.x[n+5]+xx*asncs.x[n+6]))))+asncs.x[n+7]; + xx*(asncs.x[n+5]+xx*asncs.x[n+6]))))+asncs.x[n+7]; t+=p; y = (m>0)?(hp0.x-asncs.x[n+8]):(hp0.x+asncs.x[n+8]); t = (m>0)?(hp1.x-t):(hp1.x+t); @@ -433,8 +434,8 @@ double __ieee754_acos(double x) else {xx = -x - asncs.x[n]; eps=1.02; } t = asncs.x[n+1]*xx; p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+ - xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+ - xx*asncs.x[n+7])))))+asncs.x[n+8]; + xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+ + xx*asncs.x[n+7])))))+asncs.x[n+8]; t+=p; y = (m>0)?(hp0.x-asncs.x[n+9]):(hp0.x+asncs.x[n+9]); t = (m>0)?(hp1.x-t):(hp1.x+t); @@ -468,8 +469,8 @@ double __ieee754_acos(double x) else {xx = -x - asncs.x[n]; eps = 1.01; } t = asncs.x[n+1]*xx; p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+ - xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+ - xx*asncs.x[n+8]))))))+asncs.x[n+9]; + xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+ + xx*asncs.x[n+8]))))))+asncs.x[n+9]; t+=p; y = (m>0)?(hp0.x-asncs.x[n+10]):(hp0.x+asncs.x[n+10]); t = (m>0)?(hp1.x-t):(hp1.x+t); @@ -503,9 +504,9 @@ double __ieee754_acos(double x) else {xx = -x - asncs.x[n]; eps =1.005; } t = asncs.x[n+1]*xx; p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+ - xx*(asncs.x[n+5]+xx*(asncs.x[n+6] + xx*(asncs.x[n+5]+xx*(asncs.x[n+6] +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+ - xx*asncs.x[n+9])))))))+asncs.x[n+10]; + xx*asncs.x[n+9])))))))+asncs.x[n+10]; t+=p; y = (m>0)?(hp0.x-asncs.x[n+11]):(hp0.x+asncs.x[n+11]); t = (m>0)?(hp1.x-t):(hp1.x+t); @@ -539,9 +540,9 @@ double __ieee754_acos(double x) else {xx = -x - asncs.x[n]; eps=1.005;} t = asncs.x[n+1]*xx; p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+ - xx*(asncs.x[n+5]+xx*(asncs.x[n+6] + xx*(asncs.x[n+5]+xx*(asncs.x[n+6] +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+xx*(asncs.x[n+9]+ - xx*asncs.x[n+10]))))))))+asncs.x[n+11]; + xx*asncs.x[n+10]))))))))+asncs.x[n+11]; t+=p; y = (m>0)?(hp0.x-asncs.x[n+12]):(hp0.x+asncs.x[n+12]); t = (m>0)?(hp1.x-t):(hp1.x+t); @@ -635,3 +636,4 @@ double __ieee754_acos(double x) return u.x/v.x; } } +strong_alias (__ieee754_acos, __acos_finite) diff --git a/sysdeps/ieee754/dbl-64/e_atan2.c b/sysdeps/ieee754/dbl-64/e_atan2.c index 9e1a794ec8..784fc5a0c3 100644 --- a/sysdeps/ieee754/dbl-64/e_atan2.c +++ b/sysdeps/ieee754/dbl-64/e_atan2.c @@ -63,7 +63,7 @@ double __ieee754_atan2(double y,double x) { #endif static const int pr[MM]={6,8,10,20,32}; double ax,ay,u,du,u9,ua,v,vv,dv,t1,t2,t3,t4,t5,t6,t7,t8, - z,zz,cor,s1,ss1,s2,ss2; + z,zz,cor,s1,ss1,s2,ss2; #if 0 double z1,z2; #endif @@ -73,7 +73,7 @@ double __ieee754_atan2(double y,double x) { #endif static const int ep= 59768832, /* 57*16**5 */ - em=-59768832; /* -57*16**5 */ + em=-59768832; /* -57*16**5 */ /* x=NaN or y=NaN */ num.d = x; ux = num.i[HIGH_HALF]; dx = num.i[LOW_HALF]; @@ -102,23 +102,23 @@ double __ieee754_atan2(double y,double x) { if (ux==0x7ff00000) { if (dx==0x00000000) { if (uy==0x7ff00000) { - if (dy==0x00000000) return qpi.d; } + if (dy==0x00000000) return qpi.d; } else if (uy==0xfff00000) { - if (dy==0x00000000) return mqpi.d; } + if (dy==0x00000000) return mqpi.d; } else { - if ((uy&0x80000000)==0x00000000) return ZERO; - else return MZERO; } + if ((uy&0x80000000)==0x00000000) return ZERO; + else return MZERO; } } } else if (ux==0xfff00000) { if (dx==0x00000000) { if (uy==0x7ff00000) { - if (dy==0x00000000) return tqpi.d; } + if (dy==0x00000000) return tqpi.d; } else if (uy==0xfff00000) { - if (dy==0x00000000) return mtqpi.d; } + if (dy==0x00000000) return mtqpi.d; } else { - if ((uy&0x80000000)==0x00000000) return opi.d; - else return mopi.d; } + if ((uy&0x80000000)==0x00000000) return opi.d; + else return mopi.d; } } } @@ -156,108 +156,108 @@ double __ieee754_atan2(double y,double x) { /* (i) x>0, abs(y)< abs(x): atan(ay/ax) */ if (ay<ax) { if (u<inv16.d) { - v=u*u; zz=du+u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d))))); - if ((z=u+(zz-u1.d*u)) == u+(zz+u1.d*u)) return signArctan2(y,z); + v=u*u; zz=du+u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d))))); + if ((z=u+(zz-u1.d*u)) == u+(zz+u1.d*u)) return signArctan2(y,z); - MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8) - s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d)))); - ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(u,du,s2,ss2,s1,ss1,t1,t2) - if ((z=s1+(ss1-u5.d*s1)) == s1+(ss1+u5.d*s1)) return signArctan2(y,z); - return atan2Mp(x,y,pr); + MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8) + s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d)))); + ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(u,du,s2,ss2,s1,ss1,t1,t2) + if ((z=s1+(ss1-u5.d*s1)) == s1+(ss1+u5.d*s1)) return signArctan2(y,z); + return atan2Mp(x,y,pr); } else { - i=(TWO52+TWO8*u)-TWO52; i-=16; - t3=u-cij[i][0].d; - EADD(t3,du,v,dv) - t1=cij[i][1].d; t2=cij[i][2].d; - zz=v*t2+(dv*t2+v*v*(cij[i][3].d+v*(cij[i][4].d+ - v*(cij[i][5].d+v* cij[i][6].d)))); - if (i<112) { - if (i<48) u9=u91.d; /* u < 1/4 */ - else u9=u92.d; } /* 1/4 <= u < 1/2 */ - else { - if (i<176) u9=u93.d; /* 1/2 <= u < 3/4 */ - else u9=u94.d; } /* 3/4 <= u <= 1 */ - if ((z=t1+(zz-u9*t1)) == t1+(zz+u9*t1)) return signArctan2(y,z); + i=(TWO52+TWO8*u)-TWO52; i-=16; + t3=u-cij[i][0].d; + EADD(t3,du,v,dv) + t1=cij[i][1].d; t2=cij[i][2].d; + zz=v*t2+(dv*t2+v*v*(cij[i][3].d+v*(cij[i][4].d+ + v*(cij[i][5].d+v* cij[i][6].d)))); + if (i<112) { + if (i<48) u9=u91.d; /* u < 1/4 */ + else u9=u92.d; } /* 1/4 <= u < 1/2 */ + else { + if (i<176) u9=u93.d; /* 1/2 <= u < 3/4 */ + else u9=u94.d; } /* 3/4 <= u <= 1 */ + if ((z=t1+(zz-u9*t1)) == t1+(zz+u9*t1)) return signArctan2(y,z); - t1=u-hij[i][0].d; - EADD(t1,du,v,vv) - s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+ - v*(hij[i][14].d+v* hij[i][15].d)))); - ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2) - if ((z=s2+(ss2-ub.d*s2)) == s2+(ss2+ub.d*s2)) return signArctan2(y,z); - return atan2Mp(x,y,pr); + t1=u-hij[i][0].d; + EADD(t1,du,v,vv) + s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+ + v*(hij[i][14].d+v* hij[i][15].d)))); + ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2) + if ((z=s2+(ss2-ub.d*s2)) == s2+(ss2+ub.d*s2)) return signArctan2(y,z); + return atan2Mp(x,y,pr); } } /* (ii) x>0, abs(x)<=abs(y): pi/2-atan(ax/ay) */ else { if (u<inv16.d) { - v=u*u; - zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d))))); - ESUB(hpi.d,u,t2,cor) - t3=((hpi1.d+cor)-du)-zz; - if ((z=t2+(t3-u2.d)) == t2+(t3+u2.d)) return signArctan2(y,z); + v=u*u; + zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d))))); + ESUB(hpi.d,u,t2,cor) + t3=((hpi1.d+cor)-du)-zz; + if ((z=t2+(t3-u2.d)) == t2+(t3+u2.d)) return signArctan2(y,z); - MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8) - s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d)))); - ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(u,du,s2,ss2,s1,ss1,t1,t2) - SUB2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2) - if ((z=s2+(ss2-u6.d)) == s2+(ss2+u6.d)) return signArctan2(y,z); - return atan2Mp(x,y,pr); + MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8) + s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d)))); + ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(u,du,s2,ss2,s1,ss1,t1,t2) + SUB2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2) + if ((z=s2+(ss2-u6.d)) == s2+(ss2+u6.d)) return signArctan2(y,z); + return atan2Mp(x,y,pr); } else { - i=(TWO52+TWO8*u)-TWO52; i-=16; - v=(u-cij[i][0].d)+du; - zz=hpi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+ - v*(cij[i][5].d+v* cij[i][6].d)))); - t1=hpi.d-cij[i][1].d; - if (i<112) ua=ua1.d; /* w < 1/2 */ - else ua=ua2.d; /* w >= 1/2 */ - if ((z=t1+(zz-ua)) == t1+(zz+ua)) return signArctan2(y,z); + i=(TWO52+TWO8*u)-TWO52; i-=16; + v=(u-cij[i][0].d)+du; + zz=hpi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+ + v*(cij[i][5].d+v* cij[i][6].d)))); + t1=hpi.d-cij[i][1].d; + if (i<112) ua=ua1.d; /* w < 1/2 */ + else ua=ua2.d; /* w >= 1/2 */ + if ((z=t1+(zz-ua)) == t1+(zz+ua)) return signArctan2(y,z); - t1=u-hij[i][0].d; - EADD(t1,du,v,vv) - s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+ - v*(hij[i][14].d+v* hij[i][15].d)))); - ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2) - SUB2(hpi.d,hpi1.d,s2,ss2,s1,ss1,t1,t2) - if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d)) return signArctan2(y,z); - return atan2Mp(x,y,pr); + t1=u-hij[i][0].d; + EADD(t1,du,v,vv) + s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+ + v*(hij[i][14].d+v* hij[i][15].d)))); + ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2) + SUB2(hpi.d,hpi1.d,s2,ss2,s1,ss1,t1,t2) + if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d)) return signArctan2(y,z); + return atan2Mp(x,y,pr); } } } @@ -266,112 +266,114 @@ double __ieee754_atan2(double y,double x) { /* (iii) x<0, abs(x)< abs(y): pi/2+atan(ax/ay) */ if (ax<ay) { if (u<inv16.d) { - v=u*u; - zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d))))); - EADD(hpi.d,u,t2,cor) - t3=((hpi1.d+cor)+du)+zz; - if ((z=t2+(t3-u3.d)) == t2+(t3+u3.d)) return signArctan2(y,z); + v=u*u; + zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d))))); + EADD(hpi.d,u,t2,cor) + t3=((hpi1.d+cor)+du)+zz; + if ((z=t2+(t3-u3.d)) == t2+(t3+u3.d)) return signArctan2(y,z); - MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8) - s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d)))); - ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(u,du,s2,ss2,s1,ss1,t1,t2) - ADD2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2) - if ((z=s2+(ss2-u7.d)) == s2+(ss2+u7.d)) return signArctan2(y,z); - return atan2Mp(x,y,pr); + MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8) + s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d)))); + ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(u,du,s2,ss2,s1,ss1,t1,t2) + ADD2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2) + if ((z=s2+(ss2-u7.d)) == s2+(ss2+u7.d)) return signArctan2(y,z); + return atan2Mp(x,y,pr); } else { - i=(TWO52+TWO8*u)-TWO52; i-=16; - v=(u-cij[i][0].d)+du; - zz=hpi1.d+v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+ - v*(cij[i][5].d+v* cij[i][6].d)))); - t1=hpi.d+cij[i][1].d; - if (i<112) ua=ua1.d; /* w < 1/2 */ - else ua=ua2.d; /* w >= 1/2 */ - if ((z=t1+(zz-ua)) == t1+(zz+ua)) return signArctan2(y,z); + i=(TWO52+TWO8*u)-TWO52; i-=16; + v=(u-cij[i][0].d)+du; + zz=hpi1.d+v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+ + v*(cij[i][5].d+v* cij[i][6].d)))); + t1=hpi.d+cij[i][1].d; + if (i<112) ua=ua1.d; /* w < 1/2 */ + else ua=ua2.d; /* w >= 1/2 */ + if ((z=t1+(zz-ua)) == t1+(zz+ua)) return signArctan2(y,z); - t1=u-hij[i][0].d; - EADD(t1,du,v,vv) - s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+ - v*(hij[i][14].d+v* hij[i][15].d)))); - ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2) - ADD2(hpi.d,hpi1.d,s2,ss2,s1,ss1,t1,t2) - if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d)) return signArctan2(y,z); - return atan2Mp(x,y,pr); + t1=u-hij[i][0].d; + EADD(t1,du,v,vv) + s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+ + v*(hij[i][14].d+v* hij[i][15].d)))); + ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2) + ADD2(hpi.d,hpi1.d,s2,ss2,s1,ss1,t1,t2) + if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d)) return signArctan2(y,z); + return atan2Mp(x,y,pr); } } /* (iv) x<0, abs(y)<=abs(x): pi-atan(ax/ay) */ else { if (u<inv16.d) { - v=u*u; - zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d))))); - ESUB(opi.d,u,t2,cor) - t3=((opi1.d+cor)-du)-zz; - if ((z=t2+(t3-u4.d)) == t2+(t3+u4.d)) return signArctan2(y,z); + v=u*u; + zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d))))); + ESUB(opi.d,u,t2,cor) + t3=((opi1.d+cor)-du)-zz; + if ((z=t2+(t3-u4.d)) == t2+(t3+u4.d)) return signArctan2(y,z); - MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8) - s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d)))); - ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(u,du,s2,ss2,s1,ss1,t1,t2) - SUB2(opi.d,opi1.d,s1,ss1,s2,ss2,t1,t2) - if ((z=s2+(ss2-u8.d)) == s2+(ss2+u8.d)) return signArctan2(y,z); - return atan2Mp(x,y,pr); + MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8) + s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d)))); + ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(u,du,s2,ss2,s1,ss1,t1,t2) + SUB2(opi.d,opi1.d,s1,ss1,s2,ss2,t1,t2) + if ((z=s2+(ss2-u8.d)) == s2+(ss2+u8.d)) return signArctan2(y,z); + return atan2Mp(x,y,pr); } else { - i=(TWO52+TWO8*u)-TWO52; i-=16; - v=(u-cij[i][0].d)+du; - zz=opi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+ - v*(cij[i][5].d+v* cij[i][6].d)))); - t1=opi.d-cij[i][1].d; - if (i<112) ua=ua1.d; /* w < 1/2 */ - else ua=ua2.d; /* w >= 1/2 */ - if ((z=t1+(zz-ua)) == t1+(zz+ua)) return signArctan2(y,z); + i=(TWO52+TWO8*u)-TWO52; i-=16; + v=(u-cij[i][0].d)+du; + zz=opi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+ + v*(cij[i][5].d+v* cij[i][6].d)))); + t1=opi.d-cij[i][1].d; + if (i<112) ua=ua1.d; /* w < 1/2 */ + else ua=ua2.d; /* w >= 1/2 */ + if ((z=t1+(zz-ua)) == t1+(zz+ua)) return signArctan2(y,z); - t1=u-hij[i][0].d; - EADD(t1,du,v,vv) - s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+ - v*(hij[i][14].d+v* hij[i][15].d)))); - ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2) - SUB2(opi.d,opi1.d,s2,ss2,s1,ss1,t1,t2) - if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d)) return signArctan2(y,z); - return atan2Mp(x,y,pr); + t1=u-hij[i][0].d; + EADD(t1,du,v,vv) + s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+ + v*(hij[i][14].d+v* hij[i][15].d)))); + ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2) + SUB2(opi.d,opi1.d,s2,ss2,s1,ss1,t1,t2) + if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d)) return signArctan2(y,z); + return atan2Mp(x,y,pr); } } } } +strong_alias (__ieee754_atan2, __atan2_finite) + /* Treat the Denormalized case */ static double normalized(double ax,double ay,double y, double z) { int p; diff --git a/sysdeps/ieee754/dbl-64/e_atanh.c b/sysdeps/ieee754/dbl-64/e_atanh.c index fa4fe675c9..de3bc8f144 100644 --- a/sysdeps/ieee754/dbl-64/e_atanh.c +++ b/sysdeps/ieee754/dbl-64/e_atanh.c @@ -1,74 +1,70 @@ -/* @(#)e_atanh.c 5.1 93/09/24 */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ +/* Copyright (C) 2011 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@gmail.com>, 2011. + + 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, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_atanh.c,v 1.8 1995/05/10 20:44:55 jtc Exp $"; -#endif /* __ieee754_atanh(x) - * Method : - * 1.Reduced x to positive by atanh(-x) = -atanh(x) - * 2.For x>=0.5 - * 1 2x x - * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) - * 2 1 - x 1 - x - * - * For x<0.5 - * atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) - * - * Special cases: - * atanh(x) is NaN if |x| > 1 with signal; - * atanh(NaN) is that NaN with no signal; - * atanh(+-1) is +-INF with signal. - * + Method : + 1.Reduced x to positive by atanh(-x) = -atanh(x) + 2.For x>=0.5 + 1 2x x + atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) + 2 1 - x 1 - x + + For x<0.5 + atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) + + Special cases: + atanh(x) is NaN if |x| > 1 with signal; + atanh(NaN) is that NaN with no signal; + atanh(+-1) is +-INF with signal. + */ +#include <inttypes.h> #include "math.h" #include "math_private.h" -#ifdef __STDC__ -static const double one = 1.0, huge = 1e300; -#else -static double one = 1.0, huge = 1e300; -#endif - -#ifdef __STDC__ -static const double zero = 0.0; -#else -static double zero = 0.0; -#endif +static const double huge = 1e300; -#ifdef __STDC__ - double __ieee754_atanh(double x) -#else - double __ieee754_atanh(x) - double x; -#endif +double +__ieee754_atanh (double x) { - double t; - int32_t hx,ix; - u_int32_t lx; - EXTRACT_WORDS(hx,lx,x); - ix = hx&0x7fffffff; - if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */ - return (x-x)/(x-x); - if(ix==0x3ff00000) - return x/zero; - if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */ - SET_HIGH_WORD(x,ix); - if(ix<0x3fe00000) { /* x < 0.5 */ - t = x+x; - t = 0.5*__log1p(t+t*x/(one-x)); - } else - t = 0.5*__log1p((x+x)/(one-x)); - if(hx>=0) return t; else return -t; + double xa = fabs (x); + double t; + if (xa < 0.5) + { + if (__builtin_expect (xa < 0x1.0p-28, 0) && (huge + x) > 0.0) + return x; + + t = xa + xa; + t = 0.5 * __log1p (t + t * xa / (1.0 - xa)); + } + else if (__builtin_expect (xa < 1.0, 1)) + t = 0.5 * __log1p ((xa + xa) / (1.0 - xa)); + else + { + if (xa > 1.0) + return (x - x) / (x - x); + + return x / 0.0; + } + + return __copysign (t, x); } +strong_alias (__ieee754_atanh, __atanh_finite) diff --git a/sysdeps/ieee754/dbl-64/e_cosh.c b/sysdeps/ieee754/dbl-64/e_cosh.c index 65106b9989..180ca42881 100644 --- a/sysdeps/ieee754/dbl-64/e_cosh.c +++ b/sysdeps/ieee754/dbl-64/e_cosh.c @@ -1,11 +1,11 @@ -/* @(#)e_cosh.c 5.1 93/09/24 */ +/* Optimized by Ulrich Drepper <drepper@gmail.com>, 2011 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -15,18 +15,18 @@ static char rcsid[] = "$NetBSD: e_cosh.c,v 1.7 1995/05/10 20:44:58 jtc Exp $"; #endif /* __ieee754_cosh(x) - * Method : + * Method : * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2 - * 1. Replace x by |x| (cosh(x) = cosh(-x)). - * 2. - * [ exp(x) - 1 ]^2 + * 1. Replace x by |x| (cosh(x) = cosh(-x)). + * 2. + * [ exp(x) - 1 ]^2 * 0 <= x <= ln2/2 : cosh(x) := 1 + ------------------- - * 2*exp(x) + * 2*exp(x) * - * exp(x) + 1/exp(x) + * exp(x) + 1/exp(x) * ln2/2 <= x <= 22 : cosh(x) := ------------------- - * 2 - * 22 <= x <= lnovft : cosh(x) := exp(x)/2 + * 2 + * 22 <= x <= lnovft : cosh(x) := exp(x)/2 * lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2) * ln2ovft < x : cosh(x) := huge*huge (overflow) * @@ -38,19 +38,11 @@ static char rcsid[] = "$NetBSD: e_cosh.c,v 1.7 1995/05/10 20:44:58 jtc Exp $"; #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const double one = 1.0, half=0.5, huge = 1.0e300; -#else -static double one = 1.0, half=0.5, huge = 1.0e300; -#endif -#ifdef __STDC__ - double __ieee754_cosh(double x) -#else - double __ieee754_cosh(x) - double x; -#endif -{ +double +__ieee754_cosh (double x) +{ double t,w; int32_t ix; u_int32_t lx; @@ -59,19 +51,17 @@ static double one = 1.0, half=0.5, huge = 1.0e300; GET_HIGH_WORD(ix,x); ix &= 0x7fffffff; - /* x is INF or NaN */ - if(ix>=0x7ff00000) return x*x; - - /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ - if(ix<0x3fd62e43) { - t = __expm1(fabs(x)); - w = one+t; - if (ix<0x3c800000) return w; /* cosh(tiny) = 1 */ - return one+(t*t)/(w+w); - } - - /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */ + /* |x| in [0,22] */ if (ix < 0x40360000) { + /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ + if(ix<0x3fd62e43) { + t = __expm1(fabs(x)); + w = one+t; + if (ix<0x3c800000) return w; /* cosh(tiny) = 1 */ + return one+(t*t)/(w+w); + } + + /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */ t = __ieee754_exp(fabs(x)); return half*t+half/t; } @@ -87,6 +77,10 @@ static double one = 1.0, half=0.5, huge = 1.0e300; return t*w; } + /* x is INF or NaN */ + if(ix>=0x7ff00000) return x*x; + /* |x| > overflowthresold, cosh(x) overflow */ return huge*huge; } +strong_alias (__ieee754_cosh, __cosh_finite) diff --git a/sysdeps/ieee754/dbl-64/e_exp2.c b/sysdeps/ieee754/dbl-64/e_exp2.c index ce6368be43..674cdb058c 100644 --- a/sysdeps/ieee754/dbl-64/e_exp2.c +++ b/sysdeps/ieee754/dbl-64/e_exp2.c @@ -1,5 +1,6 @@ /* Double-precision floating point 2^x. - Copyright (C) 1997,1998,2000,2001,2005,2006 Free Software Foundation, Inc. + Copyright (C) 1997,1998,2000,2001,2005,2006,2011 + Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Geoffrey Keating <geoffk@ozemail.com.au> @@ -133,3 +134,4 @@ __ieee754_exp2 (double x) /* Return x, if x is a NaN or Inf; or overflow, otherwise. */ return TWO1023*x; } +strong_alias (__ieee754_exp2, __exp2_finite) diff --git a/sysdeps/ieee754/dbl-64/e_fmod.c b/sysdeps/ieee754/dbl-64/e_fmod.c index 2ce613574a..a575f616bc 100644 --- a/sysdeps/ieee754/dbl-64/e_fmod.c +++ b/sysdeps/ieee754/dbl-64/e_fmod.c @@ -5,16 +5,12 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_fmod.c,v 1.8 1995/05/10 20:45:07 jtc Exp $"; -#endif - -/* +/* * __ieee754_fmod(x,y) * Return x mod y in exact arithmetic * Method: shift and subtract @@ -23,18 +19,10 @@ static char rcsid[] = "$NetBSD: e_fmod.c,v 1.8 1995/05/10 20:45:07 jtc Exp $"; #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const double one = 1.0, Zero[] = {0.0, -0.0,}; -#else -static double one = 1.0, Zero[] = {0.0, -0.0,}; -#endif -#ifdef __STDC__ - double __ieee754_fmod(double x, double y) -#else - double __ieee754_fmod(x,y) - double x,y ; -#endif +double +__ieee754_fmod (double x, double y) { int32_t n,hx,hy,hz,ix,iy,sx,i; u_int32_t lx,ly,lz; @@ -51,7 +39,7 @@ static double one = 1.0, Zero[] = {0.0, -0.0,}; return (x*y)/(x*y); if(hx<=hy) { if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */ - if(lx==ly) + if(lx==ly) return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/ } @@ -74,25 +62,25 @@ static double one = 1.0, Zero[] = {0.0, -0.0,}; } else iy = (hy>>20)-1023; /* set up {hx,lx}, {hy,ly} and align y to x */ - if(ix >= -1022) + if(ix >= -1022) hx = 0x00100000|(0x000fffff&hx); else { /* subnormal x, shift x to normal */ n = -1022-ix; if(n<=31) { - hx = (hx<<n)|(lx>>(32-n)); - lx <<= n; + hx = (hx<<n)|(lx>>(32-n)); + lx <<= n; } else { hx = lx<<(n-32); lx = 0; } } - if(iy >= -1022) + if(iy >= -1022) hy = 0x00100000|(0x000fffff&hy); else { /* subnormal y, shift y to normal */ n = -1022-iy; if(n<=31) { - hy = (hy<<n)|(ly>>(32-n)); - ly <<= n; + hy = (hy<<n)|(ly>>(32-n)); + ly <<= n; } else { hy = ly<<(n-32); ly = 0; @@ -105,17 +93,17 @@ static double one = 1.0, Zero[] = {0.0, -0.0,}; hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;} else { - if((hz|lz)==0) /* return sign(x)*0 */ + if((hz|lz)==0) /* return sign(x)*0 */ return Zero[(u_int32_t)sx>>31]; - hx = hz+hz+(lz>>31); lx = lz+lz; + hx = hz+hz+(lz>>31); lx = lz+lz; } } hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; if(hz>=0) {hx=hz;lx=lz;} /* convert back to floating value and restore the sign */ - if((hx|lx)==0) /* return sign(x)*0 */ - return Zero[(u_int32_t)sx>>31]; + if((hx|lx)==0) /* return sign(x)*0 */ + return Zero[(u_int32_t)sx>>31]; while(hx<0x00100000) { /* normalize x */ hx = hx+hx+(lx>>31); lx = lx+lx; iy -= 1; @@ -138,3 +126,4 @@ static double one = 1.0, Zero[] = {0.0, -0.0,}; } return x; /* exact output */ } +strong_alias (__ieee754_fmod, __fmod_finite) diff --git a/sysdeps/ieee754/dbl-64/e_gamma_r.c b/sysdeps/ieee754/dbl-64/e_gamma_r.c index f32309350c..c4b7470e5b 100644 --- a/sysdeps/ieee754/dbl-64/e_gamma_r.c +++ b/sysdeps/ieee754/dbl-64/e_gamma_r.c @@ -1,5 +1,5 @@ /* Implementation of gamma function according to ISO C. - Copyright (C) 1997, 1999, 2001, 2004 Free Software Foundation, Inc. + Copyright (C) 1997, 1999, 2001, 2004, 2011 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -33,19 +33,20 @@ __ieee754_gamma_r (double x, int *signgamp) EXTRACT_WORDS (hx, lx, x); - if (((hx & 0x7fffffff) | lx) == 0) + if (__builtin_expect (((hx & 0x7fffffff) | lx) == 0, 0)) { /* Return value for x == 0 is Inf with divide by zero exception. */ *signgamp = 0; return 1.0 / x; } - if (hx < 0 && (u_int32_t) hx < 0xfff00000 && __rint (x) == x) + if (__builtin_expect (hx < 0, 0) + && (u_int32_t) hx < 0xfff00000 && __rint (x) == x) { /* Return value for integer x < 0 is NaN with invalid exception. */ *signgamp = 0; return (x - x) / (x - x); } - if ((unsigned int) hx == 0xfff00000 && lx==0) + if (__builtin_expect ((unsigned int) hx == 0xfff00000 && lx==0, 0)) { /* x == -Inf. According to ISO this is NaN. */ *signgamp = 0; @@ -55,3 +56,4 @@ __ieee754_gamma_r (double x, int *signgamp) /* XXX FIXME. */ return __ieee754_exp (__ieee754_lgamma_r (x, signgamp)); } +strong_alias (__ieee754_gamma_r, __gamma_r_finite) diff --git a/sysdeps/ieee754/dbl-64/e_hypot.c b/sysdeps/ieee754/dbl-64/e_hypot.c index 76a77ec33a..a89ccaa35e 100644 --- a/sysdeps/ieee754/dbl-64/e_hypot.c +++ b/sysdeps/ieee754/dbl-64/e_hypot.c @@ -10,10 +10,6 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; -#endif - /* __ieee754_hypot(x,y) * * Method : @@ -42,19 +38,15 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; * hypot(x,y) is NAN if x or y is NAN. * * Accuracy: - * hypot(x,y) returns sqrt(x^2+y^2) with error less - * than 1 ulps (units in the last place) + * hypot(x,y) returns sqrt(x^2+y^2) with error less + * than 1 ulps (units in the last place) */ #include "math.h" #include "math_private.h" -#ifdef __STDC__ - double __ieee754_hypot(double x, double y) -#else - double __ieee754_hypot(x,y) - double x, y; -#endif +double +__ieee754_hypot(double x, double y) { double a,b,t1,t2,y1,y2,w; int32_t j,k,ha,hb; @@ -68,7 +60,7 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; SET_HIGH_WORD(b,hb); /* b <- |b| */ if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */ k=0; - if(ha > 0x5f300000) { /* a>2**500 */ + if(__builtin_expect(ha > 0x5f300000, 0)) { /* a>2**500 */ if(ha >= 0x7ff00000) { /* Inf or NaN */ u_int32_t low; w = a+b; /* for sNaN */ @@ -83,9 +75,9 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; SET_HIGH_WORD(a,ha); SET_HIGH_WORD(b,hb); } - if(hb < 0x20b00000) { /* b < 2**-500 */ + if(__builtin_expect(hb < 0x20b00000, 0)) { /* b < 2**-500 */ if(hb <= 0x000fffff) { /* subnormal b or 0 */ - u_int32_t low; + u_int32_t low; GET_LOW_WORD(low,b); if((hb|low)==0) return a; t1=0; @@ -94,7 +86,7 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; a *= t1; k -= 1022; } else { /* scale a and b by 2^600 */ - ha += 0x25800000; /* a *= 2^600 */ + ha += 0x25800000; /* a *= 2^600 */ hb += 0x25800000; /* b *= 2^600 */ k -= 600; SET_HIGH_WORD(a,ha); @@ -126,3 +118,4 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; return t1*w; } else return w; } +strong_alias (__ieee754_hypot, __hypot_finite) diff --git a/sysdeps/ieee754/dbl-64/e_j0.c b/sysdeps/ieee754/dbl-64/e_j0.c index 302df49d62..5ebf2056bf 100644 --- a/sysdeps/ieee754/dbl-64/e_j0.c +++ b/sysdeps/ieee754/dbl-64/e_j0.c @@ -13,10 +13,6 @@ for performance improvement on pipelined processors. */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_j0.c,v 1.8 1995/05/10 20:45:23 jtc Exp $"; -#endif - /* __ieee754_j0(x), __ieee754_y0(x) * Bessel function of the first and second kinds of order zero. * Method -- j0(x): @@ -26,16 +22,16 @@ static char rcsid[] = "$NetBSD: e_j0.c,v 1.8 1995/05/10 20:45:23 jtc Exp $"; * j0(x) = 1-z/4+ z^2*R0/S0, where z = x*x; * (precision: |j0-1+z/4-z^2R0/S0 |<2**-63.67 ) * for x in (2,inf) - * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0)) - * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) + * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0)) + * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) * as follow: * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) * = 1/sqrt(2) * (cos(x) + sin(x)) * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4) * = 1/sqrt(2) * (sin(x) - cos(x)) - * (To avoid cancellation, use + * (To avoid cancellation, use * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one.) + * to compute the worse one.) * * 3 Special cases * j0(nan)= nan @@ -56,8 +52,8 @@ static char rcsid[] = "$NetBSD: e_j0.c,v 1.8 1995/05/10 20:45:23 jtc Exp $"; * Note: For tiny x, U/V = u0 and j0(x)~1, hence * y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27) * 2. For x>=2. - * y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0)) - * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) + * y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0)) + * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) * by the method mentioned above. * 3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0. */ @@ -65,22 +61,14 @@ static char rcsid[] = "$NetBSD: e_j0.c,v 1.8 1995/05/10 20:45:23 jtc Exp $"; #include "math.h" #include "math_private.h" -#ifdef __STDC__ static double pzero(double), qzero(double); -#else -static double pzero(), qzero(); -#endif -#ifdef __STDC__ static const double -#else -static double -#endif -huge = 1e300, +huge = 1e300, one = 1.0, invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */ tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ - /* R0/S0 on [0, 2.00] */ + /* R0/S0 on [0, 2.00] */ R[] = {0.0, 0.0, 1.56249999999999947958e-02, /* 0x3F8FFFFF, 0xFFFFFFFD */ -1.89979294238854721751e-04, /* 0xBF28E6A5, 0xB61AC6E9 */ 1.82954049532700665670e-06, /* 0x3EBEB1D1, 0x0C503919 */ @@ -90,18 +78,10 @@ S[] = {0.0, 1.56191029464890010492e-02, /* 0x3F8FFCE8, 0x82C8C2A4 */ 5.13546550207318111446e-07, /* 0x3EA13B54, 0xCE84D5A9 */ 1.16614003333790000205e-09}; /* 0x3E1408BC, 0xF4745D8F */ -#ifdef __STDC__ static const double zero = 0.0; -#else -static double zero = 0.0; -#endif -#ifdef __STDC__ - double __ieee754_j0(double x) -#else - double __ieee754_j0(x) - double x; -#endif +double +__ieee754_j0(double x) { double z, s,c,ss,cc,r,u,v,r1,r2,s1,s2,z2,z4; int32_t hx,ix; @@ -117,7 +97,7 @@ static double zero = 0.0; if(ix<0x7fe00000) { /* make sure x+x not overflow */ z = -__cos(x+x); if ((s*c)<zero) cc = z/ss; - else ss = z/cc; + else ss = z/cc; } /* * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) @@ -132,8 +112,8 @@ static double zero = 0.0; } if(ix<0x3f200000) { /* |x| < 2**-13 */ if(huge+x>one) { /* raise inexact if x != 0 */ - if(ix<0x3e400000) return one; /* |x|<2**-27 */ - else return one - 0.25*x*x; + if(ix<0x3e400000) return one; /* |x|<2**-27 */ + else return one - 0.25*x*x; } } z = x*x; @@ -155,12 +135,9 @@ static double zero = 0.0; return((one+u)*(one-u)+z*(r/s)); } } +strong_alias (__ieee754_j0, __j0_finite) -#ifdef __STDC__ static const double -#else -static double -#endif U[] = {-7.38042951086872317523e-02, /* 0xBFB2E4D6, 0x99CBD01F */ 1.76666452509181115538e-01, /* 0x3FC69D01, 0x9DE9E3FC */ -1.38185671945596898896e-02, /* 0xBF8C4CE8, 0xB16CFA97 */ @@ -173,52 +150,48 @@ V[] = {1.27304834834123699328e-02, /* 0x3F8A1270, 0x91C9C71A */ 2.59150851840457805467e-07, /* 0x3E91642D, 0x7FF202FD */ 4.41110311332675467403e-10}; /* 0x3DFE5018, 0x3BD6D9EF */ -#ifdef __STDC__ - double __ieee754_y0(double x) -#else - double __ieee754_y0(x) - double x; -#endif +double +__ieee754_y0(double x) { double z, s,c,ss,cc,u,v,z2,z4,z6,u1,u2,u3,v1,v2; int32_t hx,ix,lx; EXTRACT_WORDS(hx,lx,x); - ix = 0x7fffffff&hx; + ix = 0x7fffffff&hx; /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf. */ if(ix>=0x7ff00000) return one/(x+x*x); - if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception. */ - if(hx<0) return zero/(zero*x); - if(ix >= 0x40000000) { /* |x| >= 2.0 */ - /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) - * where x0 = x-pi/4 - * Better formula: - * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) - * = 1/sqrt(2) * (sin(x) + cos(x)) - * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one. - */ + if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception. */ + if(hx<0) return zero/(zero*x); + if(ix >= 0x40000000) { /* |x| >= 2.0 */ + /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) + * where x0 = x-pi/4 + * Better formula: + * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) + * = 1/sqrt(2) * (sin(x) + cos(x)) + * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) + * = 1/sqrt(2) * (sin(x) - cos(x)) + * To avoid cancellation, use + * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) + * to compute the worse one. + */ __sincos (x, &s, &c); - ss = s-c; - cc = s+c; + ss = s-c; + cc = s+c; /* * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) */ - if(ix<0x7fe00000) { /* make sure x+x not overflow */ - z = -__cos(x+x); - if ((s*c)<zero) cc = z/ss; - else ss = z/cc; - } - if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x); - else { - u = pzero(x); v = qzero(x); - z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x); - } - return z; + if(ix<0x7fe00000) { /* make sure x+x not overflow */ + z = -__cos(x+x); + if ((s*c)<zero) cc = z/ss; + else ss = z/cc; + } + if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x); + else { + u = pzero(x); v = qzero(x); + z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x); + } + return z; } if(ix<=0x3e400000) { /* x < 2**-27 */ return(U[0] + tpi*__ieee754_log(x)); @@ -238,21 +211,18 @@ V[] = {1.27304834834123699328e-02, /* 0x3F8A1270, 0x91C9C71A */ #endif return(u/v + tpi*(__ieee754_j0(x)*__ieee754_log(x))); } +strong_alias (__ieee754_y0, __y0_finite) /* The asymptotic expansions of pzero is * 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x. * For x >= 2, We approximate pzero by - * pzero(x) = 1 + (R/S) + * pzero(x) = 1 + (R/S) * where R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10 - * S = 1 + pS0*s^2 + ... + pS4*s^10 + * S = 1 + pS0*s^2 + ... + pS4*s^10 * and * | pzero(x)-1-R/S | <= 2 ** ( -60.26) */ -#ifdef __STDC__ static const double pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#else -static double pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#endif 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ -7.03124999999900357484e-02, /* 0xBFB1FFFF, 0xFFFFFD32 */ -8.08167041275349795626e+00, /* 0xC02029D0, 0xB44FA779 */ @@ -260,11 +230,7 @@ static double pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -2.48521641009428822144e+03, /* 0xC0A36A6E, 0xCD4DCAFC */ -5.25304380490729545272e+03, /* 0xC0B4850B, 0x36CC643D */ }; -#ifdef __STDC__ static const double pS8[5] = { -#else -static double pS8[5] = { -#endif 1.16534364619668181717e+02, /* 0x405D2233, 0x07A96751 */ 3.83374475364121826715e+03, /* 0x40ADF37D, 0x50596938 */ 4.05978572648472545552e+04, /* 0x40E3D2BB, 0x6EB6B05F */ @@ -272,11 +238,7 @@ static double pS8[5] = { 4.76277284146730962675e+04, /* 0x40E74177, 0x4F2C49DC */ }; -#ifdef __STDC__ static const double pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#else -static double pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#endif -1.14125464691894502584e-11, /* 0xBDA918B1, 0x47E495CC */ -7.03124940873599280078e-02, /* 0xBFB1FFFF, 0xE69AFBC6 */ -4.15961064470587782438e+00, /* 0xC010A370, 0xF90C6BBF */ @@ -284,11 +246,7 @@ static double pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -3.31231299649172967747e+02, /* 0xC074B3B3, 0x6742CC63 */ -3.46433388365604912451e+02, /* 0xC075A6EF, 0x28A38BD7 */ }; -#ifdef __STDC__ static const double pS5[5] = { -#else -static double pS5[5] = { -#endif 6.07539382692300335975e+01, /* 0x404E6081, 0x0C98C5DE */ 1.05125230595704579173e+03, /* 0x40906D02, 0x5C7E2864 */ 5.97897094333855784498e+03, /* 0x40B75AF8, 0x8FBE1D60 */ @@ -296,11 +254,7 @@ static double pS5[5] = { 2.40605815922939109441e+03, /* 0x40A2CC1D, 0xC70BE864 */ }; -#ifdef __STDC__ static const double pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -#else -static double pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -#endif -2.54704601771951915620e-09, /* 0xBE25E103, 0x6FE1AA86 */ -7.03119616381481654654e-02, /* 0xBFB1FFF6, 0xF7C0E24B */ -2.40903221549529611423e+00, /* 0xC00345B2, 0xAEA48074 */ @@ -308,11 +262,7 @@ static double pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -5.80791704701737572236e+01, /* 0xC04D0A22, 0x420A1A45 */ -3.14479470594888503854e+01, /* 0xC03F72AC, 0xA892D80F */ }; -#ifdef __STDC__ static const double pS3[5] = { -#else -static double pS3[5] = { -#endif 3.58560338055209726349e+01, /* 0x4041ED92, 0x84077DD3 */ 3.61513983050303863820e+02, /* 0x40769839, 0x464A7C0E */ 1.19360783792111533330e+03, /* 0x4092A66E, 0x6D1061D6 */ @@ -320,11 +270,7 @@ static double pS3[5] = { 1.73580930813335754692e+02, /* 0x4065B296, 0xFC379081 */ }; -#ifdef __STDC__ static const double pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#else -static double pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#endif -8.87534333032526411254e-08, /* 0xBE77D316, 0xE927026D */ -7.03030995483624743247e-02, /* 0xBFB1FF62, 0x495E1E42 */ -1.45073846780952986357e+00, /* 0xBFF73639, 0x8A24A843 */ @@ -332,11 +278,7 @@ static double pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -1.11931668860356747786e+01, /* 0xC02662E6, 0xC5246303 */ -3.23364579351335335033e+00, /* 0xC009DE81, 0xAF8FE70F */ }; -#ifdef __STDC__ static const double pS2[5] = { -#else -static double pS2[5] = { -#endif 2.22202997532088808441e+01, /* 0x40363865, 0x908B5959 */ 1.36206794218215208048e+02, /* 0x4061069E, 0x0EE8878F */ 2.70470278658083486789e+02, /* 0x4070E786, 0x42EA079B */ @@ -344,18 +286,10 @@ static double pS2[5] = { 1.46576176948256193810e+01, /* 0x402D50B3, 0x44391809 */ }; -#ifdef __STDC__ - static double pzero(double x) -#else - static double pzero(x) - double x; -#endif +static double +pzero(double x) { -#ifdef __STDC__ const double *p,*q; -#else - double *p,*q; -#endif double z,r,s,z2,z4,r1,r2,r3,s1,s2,s3; int32_t ix; GET_HIGH_WORD(ix,x); @@ -385,17 +319,13 @@ static double pS2[5] = { /* For x >= 8, the asymptotic expansions of qzero is * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. * We approximate pzero by - * qzero(x) = s*(-1.25 + (R/S)) + * qzero(x) = s*(-1.25 + (R/S)) * where R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10 - * S = 1 + qS0*s^2 + ... + qS5*s^12 + * S = 1 + qS0*s^2 + ... + qS5*s^12 * and * | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22) */ -#ifdef __STDC__ static const double qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#else -static double qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#endif 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ 7.32421874999935051953e-02, /* 0x3FB2BFFF, 0xFFFFFE2C */ 1.17682064682252693899e+01, /* 0x40278952, 0x5BB334D6 */ @@ -403,11 +333,7 @@ static double qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 8.85919720756468632317e+03, /* 0x40C14D99, 0x3E18F46D */ 3.70146267776887834771e+04, /* 0x40E212D4, 0x0E901566 */ }; -#ifdef __STDC__ static const double qS8[6] = { -#else -static double qS8[6] = { -#endif 1.63776026895689824414e+02, /* 0x406478D5, 0x365B39BC */ 8.09834494656449805916e+03, /* 0x40BFA258, 0x4E6B0563 */ 1.42538291419120476348e+05, /* 0x41016652, 0x54D38C3F */ @@ -416,11 +342,7 @@ static double qS8[6] = { -3.43899293537866615225e+05, /* 0xC114FD6D, 0x2C9530C5 */ }; -#ifdef __STDC__ static const double qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#else -static double qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#endif 1.84085963594515531381e-11, /* 0x3DB43D8F, 0x29CC8CD9 */ 7.32421766612684765896e-02, /* 0x3FB2BFFF, 0xD172B04C */ 5.83563508962056953777e+00, /* 0x401757B0, 0xB9953DD3 */ @@ -428,11 +350,7 @@ static double qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 1.02724376596164097464e+03, /* 0x40900CF9, 0x9DC8C481 */ 1.98997785864605384631e+03, /* 0x409F17E9, 0x53C6E3A6 */ }; -#ifdef __STDC__ static const double qS5[6] = { -#else -static double qS5[6] = { -#endif 8.27766102236537761883e+01, /* 0x4054B1B3, 0xFB5E1543 */ 2.07781416421392987104e+03, /* 0x40A03BA0, 0xDA21C0CE */ 1.88472887785718085070e+04, /* 0x40D267D2, 0x7B591E6D */ @@ -441,11 +359,7 @@ static double qS5[6] = { -5.35434275601944773371e+03, /* 0xC0B4EA57, 0xBEDBC609 */ }; -#ifdef __STDC__ static const double qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -#else -static double qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -#endif 4.37741014089738620906e-09, /* 0x3E32CD03, 0x6ADECB82 */ 7.32411180042911447163e-02, /* 0x3FB2BFEE, 0x0E8D0842 */ 3.34423137516170720929e+00, /* 0x400AC0FC, 0x61149CF5 */ @@ -453,11 +367,7 @@ static double qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ 1.70808091340565596283e+02, /* 0x406559DB, 0xE25EFD1F */ 1.66733948696651168575e+02, /* 0x4064D77C, 0x81FA21E0 */ }; -#ifdef __STDC__ static const double qS3[6] = { -#else -static double qS3[6] = { -#endif 4.87588729724587182091e+01, /* 0x40486122, 0xBFE343A6 */ 7.09689221056606015736e+02, /* 0x40862D83, 0x86544EB3 */ 3.70414822620111362994e+03, /* 0x40ACF04B, 0xE44DFC63 */ @@ -466,11 +376,7 @@ static double qS3[6] = { -1.49247451836156386662e+02, /* 0xC062A7EB, 0x201CF40F */ }; -#ifdef __STDC__ static const double qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#else -static double qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#endif 1.50444444886983272379e-07, /* 0x3E84313B, 0x54F76BDB */ 7.32234265963079278272e-02, /* 0x3FB2BEC5, 0x3E883E34 */ 1.99819174093815998816e+00, /* 0x3FFFF897, 0xE727779C */ @@ -478,11 +384,7 @@ static double qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 3.16662317504781540833e+01, /* 0x403FAA8E, 0x29FBDC4A */ 1.62527075710929267416e+01, /* 0x403040B1, 0x71814BB4 */ }; -#ifdef __STDC__ static const double qS2[6] = { -#else -static double qS2[6] = { -#endif 3.03655848355219184498e+01, /* 0x403E5D96, 0xF7C07AED */ 2.69348118608049844624e+02, /* 0x4070D591, 0xE4D14B40 */ 8.44783757595320139444e+02, /* 0x408A6645, 0x22B3BF22 */ @@ -491,18 +393,10 @@ static double qS2[6] = { -5.31095493882666946917e+00, /* 0xC0153E6A, 0xF8B32931 */ }; -#ifdef __STDC__ - static double qzero(double x) -#else - static double qzero(x) - double x; -#endif +static double +qzero(double x) { -#ifdef __STDC__ const double *p,*q; -#else - double *p,*q; -#endif double s,r,z,z2,z4,z6,r1,r2,r3,s1,s2,s3; int32_t ix; GET_HIGH_WORD(ix,x); diff --git a/sysdeps/ieee754/dbl-64/e_j1.c b/sysdeps/ieee754/dbl-64/e_j1.c index 8a3b2ffd19..fdc6b5b896 100644 --- a/sysdeps/ieee754/dbl-64/e_j1.c +++ b/sysdeps/ieee754/dbl-64/e_j1.c @@ -13,10 +13,6 @@ for performance improvement on pipelined processors. */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_j1.c,v 1.8 1995/05/10 20:45:27 jtc Exp $"; -#endif - /* __ieee754_j1(x), __ieee754_y1(x) * Bessel function of the first and second kinds of order zero. * Method -- j1(x): @@ -26,17 +22,17 @@ static char rcsid[] = "$NetBSD: e_j1.c,v 1.8 1995/05/10 20:45:27 jtc Exp $"; * j1(x) = x/2 + x*z*R0/S0, where z = x*x; * (precision: |j1/x - 1/2 - R0/S0 |<2**-61.51 ) * for x in (2,inf) - * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1)) - * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) - * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) + * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1)) + * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) + * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) * as follow: * cos(x1) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) * = 1/sqrt(2) * (sin(x) - cos(x)) * sin(x1) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) * = -1/sqrt(2) * (sin(x) + cos(x)) - * (To avoid cancellation, use + * (To avoid cancellation, use * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one.) + * to compute the worse one.) * * 3 Special cases * j1(nan)= nan @@ -57,25 +53,17 @@ static char rcsid[] = "$NetBSD: e_j1.c,v 1.8 1995/05/10 20:45:27 jtc Exp $"; * Note: For tiny x, 1/x dominate y1 and hence * y1(tiny) = -2/pi/tiny, (choose tiny<2**-54) * 3. For x>=2. - * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) - * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) + * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) + * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) * by method mentioned above. */ #include "math.h" #include "math_private.h" -#ifdef __STDC__ static double pone(double), qone(double); -#else -static double pone(), qone(); -#endif -#ifdef __STDC__ static const double -#else -static double -#endif huge = 1e300, one = 1.0, invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */ @@ -91,25 +79,17 @@ S[] = {0.0, 1.91537599538363460805e-02, /* 0x3F939D0B, 0x12637E53 */ 5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */ 1.23542274426137913908e-11}; /* 0x3DAB2ACF, 0xCFB97ED8 */ -#ifdef __STDC__ static const double zero = 0.0; -#else -static double zero = 0.0; -#endif -#ifdef __STDC__ - double __ieee754_j1(double x) -#else - double __ieee754_j1(x) - double x; -#endif +double +__ieee754_j1(double x) { double z, s,c,ss,cc,r,u,v,y,r1,r2,s1,s2,s3,z2,z4; int32_t hx,ix; GET_HIGH_WORD(hx,x); ix = hx&0x7fffffff; - if(ix>=0x7ff00000) return one/x; + if(__builtin_expect(ix>=0x7ff00000, 0)) return one/x; y = fabs(x); if(ix >= 0x40000000) { /* |x| >= 2.0 */ __sincos (y, &s, &c); @@ -118,7 +98,7 @@ static double zero = 0.0; if(ix<0x7fe00000) { /* make sure y+y not overflow */ z = __cos(y+y); if ((s*c)>zero) cc = z/ss; - else ss = z/cc; + else ss = z/cc; } /* * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x) @@ -130,9 +110,9 @@ static double zero = 0.0; z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrt(y); } if(hx<0) return -z; - else return z; + else return z; } - if(ix<0x3e400000) { /* |x|<2**-27 */ + if(__builtin_expect(ix<0x3e400000, 0)) { /* |x|<2**-27 */ if(huge+x>one) return 0.5*x;/* inexact if x!=0 necessary */ } z = x*x; @@ -144,7 +124,7 @@ static double zero = 0.0; r1 = z*R[0]; z2=z*z; r2 = R[1]+z*R[2]; z4=z2*z2; r = r1 + z2*r2 + z4*R[3]; - r *= x; + r *= x; s1 = one+z*S[1]; s2 = S[2]+z*S[3]; s3 = S[4]+z*S[5]; @@ -152,23 +132,16 @@ static double zero = 0.0; #endif return(x*0.5+r/s); } +strong_alias (__ieee754_j1, __j1_finite) -#ifdef __STDC__ static const double U0[5] = { -#else -static double U0[5] = { -#endif -1.96057090646238940668e-01, /* 0xBFC91866, 0x143CBC8A */ 5.04438716639811282616e-02, /* 0x3FA9D3C7, 0x76292CD1 */ -1.91256895875763547298e-03, /* 0xBF5F55E5, 0x4844F50F */ 2.35252600561610495928e-05, /* 0x3EF8AB03, 0x8FA6B88E */ -9.19099158039878874504e-08, /* 0xBE78AC00, 0x569105B8 */ }; -#ifdef __STDC__ static const double V0[5] = { -#else -static double V0[5] = { -#endif 1.99167318236649903973e-02, /* 0x3F94650D, 0x3F4DA9F0 */ 2.02552581025135171496e-04, /* 0x3F2A8C89, 0x6C257764 */ 1.35608801097516229404e-06, /* 0x3EB6C05A, 0x894E8CA6 */ @@ -176,56 +149,53 @@ static double V0[5] = { 1.66559246207992079114e-11, /* 0x3DB25039, 0xDACA772A */ }; -#ifdef __STDC__ - double __ieee754_y1(double x) -#else - double __ieee754_y1(x) - double x; -#endif +double +__ieee754_y1(double x) { double z, s,c,ss,cc,u,v,u1,u2,v1,v2,v3,z2,z4; int32_t hx,ix,lx; EXTRACT_WORDS(hx,lx,x); - ix = 0x7fffffff&hx; + ix = 0x7fffffff&hx; /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ - if(ix>=0x7ff00000) return one/(x+x*x); - if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception. */; - if(hx<0) return zero/(zero*x); - if(ix >= 0x40000000) { /* |x| >= 2.0 */ + if(__builtin_expect(ix>=0x7ff00000, 0)) return one/(x+x*x); + if(__builtin_expect((ix|lx)==0, 0)) + return -HUGE_VAL+x; /* -inf and overflow exception. */; + if(__builtin_expect(hx<0, 0)) return zero/(zero*x); + if(ix >= 0x40000000) { /* |x| >= 2.0 */ __sincos (x, &s, &c); - ss = -s-c; - cc = s-c; - if(ix<0x7fe00000) { /* make sure x+x not overflow */ - z = __cos(x+x); - if ((s*c)>zero) cc = z/ss; - else ss = z/cc; - } - /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0)) - * where x0 = x-3pi/4 - * Better formula: - * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) - * = -1/sqrt(2) * (cos(x) + sin(x)) - * To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one. - */ - if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x); - else { - u = pone(x); v = qone(x); - z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x); - } - return z; - } - if(ix<=0x3c900000) { /* x < 2**-54 */ - return(-tpi/x); - } - z = x*x; + ss = -s-c; + cc = s-c; + if(ix<0x7fe00000) { /* make sure x+x not overflow */ + z = __cos(x+x); + if ((s*c)>zero) cc = z/ss; + else ss = z/cc; + } + /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0)) + * where x0 = x-3pi/4 + * Better formula: + * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) + * = 1/sqrt(2) * (sin(x) - cos(x)) + * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) + * = -1/sqrt(2) * (cos(x) + sin(x)) + * To avoid cancellation, use + * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) + * to compute the worse one. + */ + if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x); + else { + u = pone(x); v = qone(x); + z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x); + } + return z; + } + if(__builtin_expect(ix<=0x3c900000, 0)) { /* x < 2**-54 */ + return(-tpi/x); + } + z = x*x; #ifdef DO_NOT_USE_THIS - u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); - v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); + u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); + v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); #else u1 = U0[0]+z*U0[1];z2=z*z; u2 = U0[2]+z*U0[3];z4=z2*z2; @@ -235,24 +205,21 @@ static double V0[5] = { v3 = V0[3]+z*V0[4]; v = v1 + z2*v2 + z4*v3; #endif - return(x*(u/v) + tpi*(__ieee754_j1(x)*__ieee754_log(x)-one/x)); + return(x*(u/v) + tpi*(__ieee754_j1(x)*__ieee754_log(x)-one/x)); } +strong_alias (__ieee754_y1, __y1_finite) /* For x >= 8, the asymptotic expansions of pone is * 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x. * We approximate pone by - * pone(x) = 1 + (R/S) + * pone(x) = 1 + (R/S) * where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10 - * S = 1 + ps0*s^2 + ... + ps4*s^10 + * S = 1 + ps0*s^2 + ... + ps4*s^10 * and * | pone(x)-1-R/S | <= 2 ** ( -60.06) */ -#ifdef __STDC__ static const double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#else -static double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#endif 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ 1.17187499999988647970e-01, /* 0x3FBDFFFF, 0xFFFFFCCE */ 1.32394806593073575129e+01, /* 0x402A7A9D, 0x357F7FCE */ @@ -260,11 +227,7 @@ static double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 3.87474538913960532227e+03, /* 0x40AE457D, 0xA3A532CC */ 7.91447954031891731574e+03, /* 0x40BEEA7A, 0xC32782DD */ }; -#ifdef __STDC__ static const double ps8[5] = { -#else -static double ps8[5] = { -#endif 1.14207370375678408436e+02, /* 0x405C8D45, 0x8E656CAC */ 3.65093083420853463394e+03, /* 0x40AC85DC, 0x964D274F */ 3.69562060269033463555e+04, /* 0x40E20B86, 0x97C5BB7F */ @@ -272,11 +235,7 @@ static double ps8[5] = { 3.08042720627888811578e+04, /* 0x40DE1511, 0x697A0B2D */ }; -#ifdef __STDC__ static const double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#else -static double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#endif 1.31990519556243522749e-11, /* 0x3DAD0667, 0xDAE1CA7D */ 1.17187493190614097638e-01, /* 0x3FBDFFFF, 0xE2C10043 */ 6.80275127868432871736e+00, /* 0x401B3604, 0x6E6315E3 */ @@ -284,11 +243,7 @@ static double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 5.17636139533199752805e+02, /* 0x40802D16, 0xD052D649 */ 5.28715201363337541807e+02, /* 0x408085B8, 0xBB7E0CB7 */ }; -#ifdef __STDC__ static const double ps5[5] = { -#else -static double ps5[5] = { -#endif 5.92805987221131331921e+01, /* 0x404DA3EA, 0xA8AF633D */ 9.91401418733614377743e+02, /* 0x408EFB36, 0x1B066701 */ 5.35326695291487976647e+03, /* 0x40B4E944, 0x5706B6FB */ @@ -296,11 +251,7 @@ static double ps5[5] = { 1.50404688810361062679e+03, /* 0x40978030, 0x036F5E51 */ }; -#ifdef __STDC__ static const double pr3[6] = { -#else -static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -#endif 3.02503916137373618024e-09, /* 0x3E29FC21, 0xA7AD9EDD */ 1.17186865567253592491e-01, /* 0x3FBDFFF5, 0x5B21D17B */ 3.93297750033315640650e+00, /* 0x400F76BC, 0xE85EAD8A */ @@ -308,11 +259,7 @@ static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ 9.10550110750781271918e+01, /* 0x4056C385, 0x4D2C1837 */ 4.85590685197364919645e+01, /* 0x4048478F, 0x8EA83EE5 */ }; -#ifdef __STDC__ static const double ps3[5] = { -#else -static double ps3[5] = { -#endif 3.47913095001251519989e+01, /* 0x40416549, 0xA134069C */ 3.36762458747825746741e+02, /* 0x40750C33, 0x07F1A75F */ 1.04687139975775130551e+03, /* 0x40905B7C, 0x5037D523 */ @@ -320,11 +267,7 @@ static double ps3[5] = { 1.03787932439639277504e+02, /* 0x4059F26D, 0x7C2EED53 */ }; -#ifdef __STDC__ static const double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#else -static double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#endif 1.07710830106873743082e-07, /* 0x3E7CE9D4, 0xF65544F4 */ 1.17176219462683348094e-01, /* 0x3FBDFF42, 0xBE760D83 */ 2.36851496667608785174e+00, /* 0x4002F2B7, 0xF98FAEC0 */ @@ -332,11 +275,7 @@ static double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 1.76939711271687727390e+01, /* 0x4031B1A8, 0x177F8EE2 */ 5.07352312588818499250e+00, /* 0x40144B49, 0xA574C1FE */ }; -#ifdef __STDC__ static const double ps2[5] = { -#else -static double ps2[5] = { -#endif 2.14364859363821409488e+01, /* 0x40356FBD, 0x8AD5ECDC */ 1.25290227168402751090e+02, /* 0x405F5293, 0x14F92CD5 */ 2.32276469057162813669e+02, /* 0x406D08D8, 0xD5A2DBD9 */ @@ -344,30 +283,22 @@ static double ps2[5] = { 8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */ }; -#ifdef __STDC__ - static double pone(double x) -#else - static double pone(x) - double x; -#endif +static double +pone(double x) { -#ifdef __STDC__ const double *p,*q; -#else - double *p,*q; -#endif double z,r,s,r1,r2,r3,s1,s2,s3,z2,z4; - int32_t ix; + int32_t ix; GET_HIGH_WORD(ix,x); ix &= 0x7fffffff; - if(ix>=0x40200000) {p = pr8; q= ps8;} - else if(ix>=0x40122E8B){p = pr5; q= ps5;} - else if(ix>=0x4006DB6D){p = pr3; q= ps3;} - else if(ix>=0x40000000){p = pr2; q= ps2;} - z = one/(x*x); + if(ix>=0x40200000) {p = pr8; q= ps8;} + else if(ix>=0x40122E8B){p = pr5; q= ps5;} + else if(ix>=0x4006DB6D){p = pr3; q= ps3;} + else if(ix>=0x40000000){p = pr2; q= ps2;} + z = one/(x*x); #ifdef DO_NOT_USE_THIS - r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); + r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); + s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); #else r1 = p[0]+z*p[1]; z2=z*z; r2 = p[2]+z*p[3]; z4=z2*z2; @@ -378,25 +309,21 @@ static double ps2[5] = { s3 = q[3]+z*q[4]; s = s1 + z2*s2 + z4*s3; #endif - return one+ r/s; + return one+ r/s; } /* For x >= 8, the asymptotic expansions of qone is * 3/8 s - 105/1024 s^3 - ..., where s = 1/x. * We approximate pone by - * qone(x) = s*(0.375 + (R/S)) + * qone(x) = s*(0.375 + (R/S)) * where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10 - * S = 1 + qs1*s^2 + ... + qs6*s^12 + * S = 1 + qs1*s^2 + ... + qs6*s^12 * and * | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13) */ -#ifdef __STDC__ static const double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#else -static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#endif 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ -1.02539062499992714161e-01, /* 0xBFBA3FFF, 0xFFFFFDF3 */ -1.62717534544589987888e+01, /* 0xC0304591, 0xA26779F7 */ @@ -404,11 +331,7 @@ static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -1.18498066702429587167e+04, /* 0xC0C724E7, 0x40F87415 */ -4.84385124285750353010e+04, /* 0xC0E7A6D0, 0x65D09C6A */ }; -#ifdef __STDC__ static const double qs8[6] = { -#else -static double qs8[6] = { -#endif 1.61395369700722909556e+02, /* 0x40642CA6, 0xDE5BCDE5 */ 7.82538599923348465381e+03, /* 0x40BE9162, 0xD0D88419 */ 1.33875336287249578163e+05, /* 0x4100579A, 0xB0B75E98 */ @@ -417,11 +340,7 @@ static double qs8[6] = { -2.94490264303834643215e+05, /* 0xC111F969, 0x0EA5AA18 */ }; -#ifdef __STDC__ static const double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#else -static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#endif -2.08979931141764104297e-11, /* 0xBDB6FA43, 0x1AA1A098 */ -1.02539050241375426231e-01, /* 0xBFBA3FFF, 0xCB597FEF */ -8.05644828123936029840e+00, /* 0xC0201CE6, 0xCA03AD4B */ @@ -429,11 +348,7 @@ static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -1.37319376065508163265e+03, /* 0xC09574C6, 0x6931734F */ -2.61244440453215656817e+03, /* 0xC0A468E3, 0x88FDA79D */ }; -#ifdef __STDC__ static const double qs5[6] = { -#else -static double qs5[6] = { -#endif 8.12765501384335777857e+01, /* 0x405451B2, 0xFF5A11B2 */ 1.99179873460485964642e+03, /* 0x409F1F31, 0xE77BF839 */ 1.74684851924908907677e+04, /* 0x40D10F1F, 0x0D64CE29 */ @@ -442,11 +357,7 @@ static double qs5[6] = { -4.71918354795128470869e+03, /* 0xC0B26F2E, 0xFCFFA004 */ }; -#ifdef __STDC__ static const double qr3[6] = { -#else -static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -#endif -5.07831226461766561369e-09, /* 0xBE35CFA9, 0xD38FC84F */ -1.02537829820837089745e-01, /* 0xBFBA3FEB, 0x51AEED54 */ -4.61011581139473403113e+00, /* 0xC01270C2, 0x3302D9FF */ @@ -454,11 +365,7 @@ static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -2.28244540737631695038e+02, /* 0xC06C87D3, 0x4718D55F */ -2.19210128478909325622e+02, /* 0xC06B66B9, 0x5F5C1BF6 */ }; -#ifdef __STDC__ static const double qs3[6] = { -#else -static double qs3[6] = { -#endif 4.76651550323729509273e+01, /* 0x4047D523, 0xCCD367E4 */ 6.73865112676699709482e+02, /* 0x40850EEB, 0xC031EE3E */ 3.38015286679526343505e+03, /* 0x40AA684E, 0x448E7C9A */ @@ -467,11 +374,7 @@ static double qs3[6] = { -1.35201191444307340817e+02, /* 0xC060E670, 0x290A311F */ }; -#ifdef __STDC__ static const double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#else -static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#endif -1.78381727510958865572e-07, /* 0xBE87F126, 0x44C626D2 */ -1.02517042607985553460e-01, /* 0xBFBA3E8E, 0x9148B010 */ -2.75220568278187460720e+00, /* 0xC0060484, 0x69BB4EDA */ @@ -479,11 +382,7 @@ static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -4.23253133372830490089e+01, /* 0xC04529A3, 0xDE104AAA */ -2.13719211703704061733e+01, /* 0xC0355F36, 0x39CF6E52 */ }; -#ifdef __STDC__ static const double qs2[6] = { -#else -static double qs2[6] = { -#endif 2.95333629060523854548e+01, /* 0x403D888A, 0x78AE64FF */ 2.52981549982190529136e+02, /* 0x406F9F68, 0xDB821CBA */ 7.57502834868645436472e+02, /* 0x4087AC05, 0xCE49A0F7 */ @@ -492,18 +391,10 @@ static double qs2[6] = { -4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */ }; -#ifdef __STDC__ - static double qone(double x) -#else - static double qone(x) - double x; -#endif +static double +qone(double x) { -#ifdef __STDC__ const double *p,*q; -#else - double *p,*q; -#endif double s,r,z,r1,r2,r3,s1,s2,s3,z2,z4,z6; int32_t ix; GET_HIGH_WORD(ix,x); diff --git a/sysdeps/ieee754/dbl-64/e_jn.c b/sysdeps/ieee754/dbl-64/e_jn.c index d9d6f91762..f8b8e2ee6a 100644 --- a/sysdeps/ieee754/dbl-64/e_jn.c +++ b/sysdeps/ieee754/dbl-64/e_jn.c @@ -10,10 +10,6 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_jn.c,v 1.9 1995/05/10 20:45:34 jtc Exp $"; -#endif - /* * __ieee754_jn(n, x), __ieee754_yn(n, x) * floating point Bessel's function of the 1st and 2nd kind @@ -43,27 +39,15 @@ static char rcsid[] = "$NetBSD: e_jn.c,v 1.9 1995/05/10 20:45:34 jtc Exp $"; #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const double -#else -static double -#endif invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */ two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */ one = 1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */ -#ifdef __STDC__ static const double zero = 0.00000000000000000000e+00; -#else -static double zero = 0.00000000000000000000e+00; -#endif -#ifdef __STDC__ - double __ieee754_jn(int n, double x) -#else - double __ieee754_jn(n,x) - int n; double x; -#endif +double +__ieee754_jn(int n, double x) { int32_t i,hx,ix,lx, sgn; double a, b, temp, di; @@ -75,7 +59,8 @@ static double zero = 0.00000000000000000000e+00; EXTRACT_WORDS(hx,lx,x); ix = 0x7fffffff&hx; /* if J(n,NaN) is NaN */ - if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x; + if(__builtin_expect((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000, 0)) + return x+x; if(n<0){ n = -n; x = -x; @@ -85,8 +70,9 @@ static double zero = 0.00000000000000000000e+00; if(n==1) return(__ieee754_j1(x)); sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */ x = fabs(x); - if((ix|lx)==0||ix>=0x7ff00000) /* if x is 0 or inf */ - b = zero; + if(__builtin_expect((ix|lx)==0||ix>=0x7ff00000,0)) + /* if x is 0 or inf */ + b = zero; else if((double)n<=x) { /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ if(ix>=0x52D00000) { /* x > 2**302 */ @@ -99,7 +85,7 @@ static double zero = 0.00000000000000000000e+00; * n sin(xn)*sqt2 cos(xn)*sqt2 * ---------------------------------- * 0 s-c c+s - * 1 -s-c -c+s + * 1 -s-c -c+s * 2 -s+c -c-s * 3 s+c c-s */ @@ -114,13 +100,13 @@ static double zero = 0.00000000000000000000e+00; } b = invsqrtpi*temp/__ieee754_sqrt(x); } else { - a = __ieee754_j0(x); - b = __ieee754_j1(x); - for(i=1;i<n;i++){ + a = __ieee754_j0(x); + b = __ieee754_j1(x); + for(i=1;i<n;i++){ temp = b; b = b*((double)(i+i)/x) - a; /* avoid underflow */ a = temp; - } + } } } else { if(ix<0x3e100000) { /* x < 2**-29 */ @@ -139,11 +125,11 @@ static double zero = 0.00000000000000000000e+00; } } else { /* use backward recurrence */ - /* x x^2 x^2 + /* x x^2 x^2 * J(n,x)/J(n-1,x) = ---- ------ ------ ..... * 2n - 2(n+1) - 2(n+2) * - * 1 1 1 + * 1 1 1 * (for large x) = ---- ------ ------ ..... * 2n 2(n+1) 2(n+2) * -- - ------ - ------ - @@ -156,7 +142,7 @@ static double zero = 0.00000000000000000000e+00; * 1 * w - ----------------- * 1 - * w+h - --------- + * w+h - --------- * w+2h - ... * * To determine how many terms needed, let @@ -193,19 +179,19 @@ static double zero = 0.00000000000000000000e+00; v = two/x; tmp = tmp*__ieee754_log(fabs(v*tmp)); if(tmp<7.09782712893383973096e+02) { - for(i=n-1,di=(double)(i+i);i>0;i--){ - temp = b; + for(i=n-1,di=(double)(i+i);i>0;i--){ + temp = b; b *= di; b = b/x - a; - a = temp; + a = temp; di -= two; - } + } } else { - for(i=n-1,di=(double)(i+i);i>0;i--){ - temp = b; + for(i=n-1,di=(double)(i+i);i>0;i--){ + temp = b; b *= di; b = b/x - a; - a = temp; + a = temp; di -= two; /* scale b to avoid spurious overflow */ if(b>1e100) { @@ -213,7 +199,7 @@ static double zero = 0.00000000000000000000e+00; t /= b; b = one; } - } + } } /* j0() and j1() suffer enormous loss of precision at and * near zero; however, we know that their zero points never @@ -229,13 +215,10 @@ static double zero = 0.00000000000000000000e+00; } if(sgn==1) return -b; else return b; } +strong_alias (__ieee754_jn, __jn_finite) -#ifdef __STDC__ - double __ieee754_yn(int n, double x) -#else - double __ieee754_yn(n,x) - int n; double x; -#endif +double +__ieee754_yn(int n, double x) { int32_t i,hx,ix,lx; int32_t sign; @@ -244,9 +227,11 @@ static double zero = 0.00000000000000000000e+00; EXTRACT_WORDS(hx,lx,x); ix = 0x7fffffff&hx; /* if Y(n,NaN) is NaN */ - if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x; - if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception. */; - if(hx<0) return zero/(zero*x); + if(__builtin_expect((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000,0)) + return x+x; + if(__builtin_expect((ix|lx)==0, 0)) + return -HUGE_VAL+x; /* -inf and overflow exception. */; + if(__builtin_expect(hx<0, 0)) return zero/(zero*x); sign = 1; if(n<0){ n = -n; @@ -254,7 +239,7 @@ static double zero = 0.00000000000000000000e+00; } if(n==0) return(__ieee754_y0(x)); if(n==1) return(sign*__ieee754_y1(x)); - if(ix==0x7ff00000) return zero; + if(__builtin_expect(ix==0x7ff00000, 0)) return zero; if(ix>=0x52D00000) { /* x > 2**302 */ /* (x >> n**2) * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) @@ -265,7 +250,7 @@ static double zero = 0.00000000000000000000e+00; * n sin(xn)*sqt2 cos(xn)*sqt2 * ---------------------------------- * 0 s-c c+s - * 1 -s-c -c+s + * 1 -s-c -c+s * 2 -s+c -c-s * 3 s+c c-s */ @@ -294,3 +279,4 @@ static double zero = 0.00000000000000000000e+00; } if(sign>0) return b; else return -b; } +strong_alias (__ieee754_yn, __yn_finite) diff --git a/sysdeps/ieee754/dbl-64/e_lgamma_r.c b/sysdeps/ieee754/dbl-64/e_lgamma_r.c index a298a5a2a4..e26ce8a247 100644 --- a/sysdeps/ieee754/dbl-64/e_lgamma_r.c +++ b/sysdeps/ieee754/dbl-64/e_lgamma_r.c @@ -10,19 +10,15 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_lgamma_r.c,v 1.7 1995/05/10 20:45:42 jtc Exp $"; -#endif - /* __ieee754_lgamma_r(x, signgamp) * Reentrant version of the logarithm of the Gamma function * with user provide pointer for the sign of Gamma(x). * * Method: * 1. Argument Reduction for 0 < x <= 8 - * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may - * reduce x to a number in [1.5,2.5] by - * lgamma(1+s) = log(s) + lgamma(s) + * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may + * reduce x to a number in [1.5,2.5] by + * lgamma(1+s) = log(s) + lgamma(s) * for example, * lgamma(7.3) = log(6.3) + lgamma(6.3) * = log(6.3*5.3) + lgamma(5.3) @@ -56,15 +52,15 @@ static char rcsid[] = "$NetBSD: e_lgamma_r.c,v 1.7 1995/05/10 20:45:42 jtc Exp $ * Let z = 1/x, then we approximation * f(z) = lgamma(x) - (x-0.5)(log(x)-1) * by - * 3 5 11 + * 3 5 11 * w = w0 + w1*z + w2*z + w3*z + ... + w6*z * where * |w - f(z)| < 2**-58.74 * * 4. For negative x, since (G is gamma function) * -x*G(-x)*G(x) = pi/sin(pi*x), - * we have - * G(x) = pi/(sin(pi*x)*(-x)*G(-x)) + * we have + * G(x) = pi/(sin(pi*x)*(-x)*G(-x)) * since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0 * Hence, for x<0, signgam = sign(sin(pi*x)) and * lgamma(x) = log(|Gamma(x)|) @@ -77,18 +73,14 @@ static char rcsid[] = "$NetBSD: e_lgamma_r.c,v 1.7 1995/05/10 20:45:42 jtc Exp $ * lgamma(1)=lgamma(2)=0 * lgamma(x) ~ -log(x) for tiny x * lgamma(0) = lgamma(inf) = inf - * lgamma(-integer) = +-inf + * lgamma(-integer) = +-inf * */ #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const double -#else -static double -#endif two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ @@ -156,18 +148,10 @@ w4 = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */ w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */ w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */ -#ifdef __STDC__ static const double zero= 0.00000000000000000000e+00; -#else -static double zero= 0.00000000000000000000e+00; -#endif -#ifdef __STDC__ - static double sin_pi(double x) -#else - static double sin_pi(x) - double x; -#endif +static double +sin_pi(double x) { double y,z; int n,ix; @@ -188,16 +172,16 @@ static double zero= 0.00000000000000000000e+00; y = 2.0*(y - __floor(y)); /* y = |x| mod 2.0 */ n = (int) (y*4.0); } else { - if(ix>=0x43400000) { - y = zero; n = 0; /* y must be even */ - } else { - if(ix<0x43300000) z = y+two52; /* exact */ + if(ix>=0x43400000) { + y = zero; n = 0; /* y must be even */ + } else { + if(ix<0x43300000) z = y+two52; /* exact */ GET_LOW_WORD(n,z); n &= 1; - y = n; - n<<= 2; - } - } + y = n; + n<<= 2; + } + } switch (n) { case 0: y = __sin(pi*y); break; case 1: @@ -212,12 +196,8 @@ static double zero= 0.00000000000000000000e+00; } -#ifdef __STDC__ - double __ieee754_lgamma_r(double x, int *signgamp) -#else - double __ieee754_lgamma_r(x,signgamp) - double x; int *signgamp; -#endif +double +__ieee754_lgamma_r(double x, int *signgamp) { double t,y,z,nadj,p,p1,p2,p3,q,r,w; int i,hx,lx,ix; @@ -227,21 +207,23 @@ static double zero= 0.00000000000000000000e+00; /* purge off +-inf, NaN, +-0, and negative arguments */ *signgamp = 1; ix = hx&0x7fffffff; - if(ix>=0x7ff00000) return x*x; - if((ix|lx)==0) + if(__builtin_expect(ix>=0x7ff00000, 0)) return x*x; + if(__builtin_expect((ix|lx)==0, 0)) { if (hx < 0) *signgamp = -1; return one/fabs(x); } - if(ix<0x3b900000) { /* |x|<2**-70, return -log(|x|) */ + if(__builtin_expect(ix<0x3b900000, 0)) { + /* |x|<2**-70, return -log(|x|) */ if(hx<0) { - *signgamp = -1; - return -__ieee754_log(-x); + *signgamp = -1; + return -__ieee754_log(-x); } else return -__ieee754_log(x); } if(hx<0) { - if(ix>=0x43300000) /* |x|>=2**52, must be -integer */ + if(__builtin_expect(ix>=0x43300000, 0)) + /* |x|>=2**52, must be -integer */ return x/zero; t = sin_pi(x); if(t==zero) return one/fabsf(t); /* -integer */ @@ -254,15 +236,15 @@ static double zero= 0.00000000000000000000e+00; if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0; /* for x < 2.0 */ else if(ix<0x40000000) { - if(ix<=0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */ + if(ix<=0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */ r = -__ieee754_log(x); if(ix>=0x3FE76944) {y = one-x; i= 0;} else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;} - else {y = x; i=2;} + else {y = x; i=2;} } else { - r = zero; - if(ix>=0x3FFBB4C3) {y=2.0-x;i=0;} /* [1.7316,2] */ - else if(ix>=0x3FF3B4C4) {y=x-tc;i=1;} /* [1.23,1.73] */ + r = zero; + if(ix>=0x3FFBB4C3) {y=2.0-x;i=0;} /* [1.7316,2] */ + else if(ix>=0x3FF3B4C4) {y=x-tc;i=1;} /* [1.23,1.73] */ else {y=x-one;i=2;} } switch(i) { @@ -286,7 +268,7 @@ static double zero= 0.00000000000000000000e+00; r += (-0.5*y + p1/p2); } } - else if(ix<0x40200000) { /* x < 8.0 */ + else if(ix<0x40200000) { /* x < 8.0 */ i = (int)x; t = zero; y = x-(double)i; @@ -315,3 +297,4 @@ static double zero= 0.00000000000000000000e+00; if(hx<0) r = nadj - r; return r; } +strong_alias (__ieee754_lgamma_r, __lgamma_r_finite) diff --git a/sysdeps/ieee754/dbl-64/e_log.c b/sysdeps/ieee754/dbl-64/e_log.c index 1a9967b546..5d320db994 100644 --- a/sysdeps/ieee754/dbl-64/e_log.c +++ b/sysdeps/ieee754/dbl-64/e_log.c @@ -55,9 +55,9 @@ double __ieee754_log(double x) { int k; #endif double dbl_n,u,p0,q,r0,w,nln2a,luai,lubi,lvaj,lvbj, - sij,ssij,ttij,A,B,B0,y,y1,y2,polI,polII,sa,sb, - t1,t2,t3,t4,t5,t6,t7,t8,t,ra,rb,ww, - a0,aa0,s1,s2,ss2,s3,ss3,a1,aa1,a,aa,b,bb,c; + sij,ssij,ttij,A,B,B0,y,y1,y2,polI,polII,sa,sb, + t1,t2,t3,t4,t5,t6,t7,t8,t,ra,rb,ww, + a0,aa0,s1,s2,ss2,s3,ss3,a1,aa1,a,aa,b,bb,c; number num; mp_no mpx,mpy,mpy1,mpy2,mperr; @@ -69,12 +69,15 @@ double __ieee754_log(double x) { num.d = x; ux = num.i[HIGH_HALF]; dx = num.i[LOW_HALF]; n=0; if (ux < 0x00100000) { - if (((ux & 0x7fffffff) | dx) == 0) return MHALF/ZERO; /* return -INF */ - if (ux < 0) return (x-x)/ZERO; /* return NaN */ + if (__builtin_expect(((ux & 0x7fffffff) | dx) == 0, 0)) + return MHALF/ZERO; /* return -INF */ + if (__builtin_expect(ux < 0, 0)) + return (x-x)/ZERO; /* return NaN */ n -= 54; x *= two54.d; /* scale x */ num.d = x; } - if (ux >= 0x7ff00000) return x+x; /* INF or NaN */ + if (__builtin_expect(ux >= 0x7ff00000, 0)) + return x+x; /* INF or NaN */ /* Regular values of x */ @@ -90,7 +93,7 @@ double __ieee754_log(double x) { /* Evaluate polynomial II */ polII = (b0.d+w*(b1.d+w*(b2.d+w*(b3.d+w*(b4.d+ - w*(b5.d+w*(b6.d+w*(b7.d+w*b8.d))))))))*w*w*w; + w*(b5.d+w*(b6.d+w*(b7.d+w*b8.d))))))))*w*w*w; c = (aa+bb)+polII; /* End stage I, case abs(x-1) < 0.03 */ @@ -99,7 +102,7 @@ double __ieee754_log(double x) { /*--- Stage II, the case abs(x-1) < 0.03 */ a = d11.d+w*(d12.d+w*(d13.d+w*(d14.d+w*(d15.d+w*(d16.d+ - w*(d17.d+w*(d18.d+w*(d19.d+w*d20.d)))))))); + w*(d17.d+w*(d18.d+w*(d19.d+w*d20.d)))))))); EMULV(w,a,s2,ss2,t1,t2,t3,t4,t5) ADD2(d10.d,dd10.d,s2,ss2,s3,ss3,t1,t2) MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) @@ -201,3 +204,4 @@ double __ieee754_log(double x) { } return y1; } +strong_alias (__ieee754_log, __log_finite) diff --git a/sysdeps/ieee754/dbl-64/e_log10.c b/sysdeps/ieee754/dbl-64/e_log10.c index e8a3278eaf..6a630bcef7 100644 --- a/sysdeps/ieee754/dbl-64/e_log10.c +++ b/sysdeps/ieee754/dbl-64/e_log10.c @@ -10,10 +10,6 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_log10.c,v 1.9 1995/05/10 20:45:51 jtc Exp $"; -#endif - /* __ieee754_log10(x) * Return the base 10 logarithm of x * @@ -50,28 +46,16 @@ static char rcsid[] = "$NetBSD: e_log10.c,v 1.9 1995/05/10 20:45:51 jtc Exp $"; #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const double -#else -static double -#endif two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ ivln10 = 4.34294481903251816668e-01, /* 0x3FDBCB7B, 0x1526E50E */ log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */ log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */ -#ifdef __STDC__ static const double zero = 0.0; -#else -static double zero = 0.0; -#endif -#ifdef __STDC__ - double __ieee754_log10(double x) -#else - double __ieee754_log10(x) - double x; -#endif +double +__ieee754_log10(double x) { double y,z; int32_t i,k,hx; @@ -79,20 +63,22 @@ static double zero = 0.0; EXTRACT_WORDS(hx,lx,x); - k=0; - if (hx < 0x00100000) { /* x < 2**-1022 */ - if (((hx&0x7fffffff)|lx)==0) - return -two54/(x-x); /* log(+-0)=-inf */ - if (hx<0) return (x-x)/(x-x); /* log(-#) = NaN */ - k -= 54; x *= two54; /* subnormal number, scale up x */ + k=0; + if (hx < 0x00100000) { /* x < 2**-1022 */ + if (__builtin_expect(((hx&0x7fffffff)|lx)==0, 0)) + return -two54/(x-x); /* log(+-0)=-inf */ + if (__builtin_expect(hx<0, 0)) + return (x-x)/(x-x); /* log(-#) = NaN */ + k -= 54; x *= two54; /* subnormal number, scale up x */ GET_HIGH_WORD(hx,x); - } - if (hx >= 0x7ff00000) return x+x; + } + if (__builtin_expect(hx >= 0x7ff00000, 0)) return x+x; k += (hx>>20)-1023; i = ((u_int32_t)k&0x80000000)>>31; - hx = (hx&0x000fffff)|((0x3ff-i)<<20); - y = (double)(k+i); + hx = (hx&0x000fffff)|((0x3ff-i)<<20); + y = (double)(k+i); SET_HIGH_WORD(x,hx); z = y*log10_2lo + ivln10*__ieee754_log(x); return z+y*log10_2hi; } +strong_alias (__ieee754_log10, __log10_finite) diff --git a/sysdeps/ieee754/dbl-64/e_log2.c b/sysdeps/ieee754/dbl-64/e_log2.c index f05d0ce966..be41cb4e68 100644 --- a/sysdeps/ieee754/dbl-64/e_log2.c +++ b/sysdeps/ieee754/dbl-64/e_log2.c @@ -21,14 +21,14 @@ * 2. Approximation of log(1+f). * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) * = 2s + 2/3 s**3 + 2/5 s**5 + ....., - * = 2s + s*R + * = 2s + s*R * We use a special Reme algorithm on [0,0.1716] to generate - * a polynomial of degree 14 to approximate R The maximum error + * a polynomial of degree 14 to approximate R The maximum error * of this polynomial approximation is bounded by 2**-58.45. In * other words, - * 2 4 6 8 10 12 14 + * 2 4 6 8 10 12 14 * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s - * (the values of Lg1 to Lg7 are listed in the program) + * (the values of Lg1 to Lg7 are listed in the program) * and * | 2 14 | -58.45 * | Lg1*s +...+Lg7*s - R(z) | <= 2 @@ -57,11 +57,7 @@ #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const double -#else -static double -#endif ln2 = 0.69314718055994530942, two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ @@ -72,18 +68,10 @@ Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ -#ifdef __STDC__ static const double zero = 0.0; -#else -static double zero = 0.0; -#endif -#ifdef __STDC__ - double __ieee754_log2(double x) -#else - double __ieee754_log2(x) - double x; -#endif +double +__ieee754_log2(double x) { double hfsq,f,s,z,R,w,t1,t2,dk; int32_t k,hx,i,j; @@ -93,13 +81,14 @@ static double zero = 0.0; k=0; if (hx < 0x00100000) { /* x < 2**-1022 */ - if (((hx&0x7fffffff)|lx)==0) + if (__builtin_expect(((hx&0x7fffffff)|lx)==0, 0)) return -two54/(x-x); /* log(+-0)=-inf */ - if (hx<0) return (x-x)/(x-x); /* log(-#) = NaN */ + if (__builtin_expect(hx<0, 0)) + return (x-x)/(x-x); /* log(-#) = NaN */ k -= 54; x *= two54; /* subnormal number, scale up x */ GET_HIGH_WORD(hx,x); } - if (hx >= 0x7ff00000) return x+x; + if (__builtin_expect(hx >= 0x7ff00000, 0)) return x+x; k += (hx>>20)-1023; hx &= 0x000fffff; i = (hx+0x95f64)&0x100000; @@ -112,7 +101,7 @@ static double zero = 0.0; R = f*f*(0.5-0.33333333333333333*f); return dk-(R-f)/ln2; } - s = f/(2.0+f); + s = f/(2.0+f); z = s*s; i = hx-0x6147a; w = z*z; @@ -128,3 +117,4 @@ static double zero = 0.0; return dk-((s*(f-R))-f)/ln2; } } +strong_alias (__ieee754_log2, __log2_finite) diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c index 1e159f2c0b..83a5eff5c2 100644 --- a/sysdeps/ieee754/dbl-64/e_pow.c +++ b/sysdeps/ieee754/dbl-64/e_pow.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001, 2002, 2004 Free Software Foundation + * Copyright (C) 2001, 2002, 2004, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -22,12 +22,12 @@ /* */ /* FUNCTIONS: upow */ /* power1 */ -/* my_log2 */ +/* my_log2 */ /* log1 */ /* checkint */ /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h */ /* halfulp.c mpexp.c mplog.c slowexp.c slowpow.c mpa.c */ -/* uexp.c upow.c */ +/* uexp.c upow.c */ /* root.tbl uexp.tbl upow.tbl */ /* An ultimate power routine. Given two IEEE double machine numbers y,x */ /* it computes the correctly rounded (to nearest) value of x^y. */ @@ -77,7 +77,7 @@ double __ieee754_pow(double x, double y) { /* else */ if(((u.i[HIGH_HALF]>0 && u.i[HIGH_HALF]<0x7ff00000)|| /* x>0 and not x->0 */ (u.i[HIGH_HALF]==0 && u.i[LOW_HALF]!=0)) && - /* 2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */ + /* 2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */ (v.i[HIGH_HALF]&0x7fffffff) < 0x4ff00000) { /* if y<-1 or y>1 */ z = log1(x,&aa,&error); /* x^y =e^(y log (X)) */ t = y*134217729.0; @@ -153,6 +153,7 @@ double __ieee754_pow(double x, double y) { if (y<0) return (x<1.0)?INF.x:0; return 0; /* unreachable, to make the compiler happy */ } +strong_alias (__ieee754_pow, __pow_finite) /**************************************************************************/ /* Computing x^y using more accurate but more slow log routine */ diff --git a/sysdeps/ieee754/dbl-64/e_remainder.c b/sysdeps/ieee754/dbl-64/e_remainder.c index cc06e18ce8..d1782a15cf 100644 --- a/sysdeps/ieee754/dbl-64/e_remainder.c +++ b/sysdeps/ieee754/dbl-64/e_remainder.c @@ -1,8 +1,8 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation - * + * Copyright (C) 2001, 2011 Free Software Foundation + * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by * the Free Software Foundation; either version 2.1 of the License, or @@ -96,9 +96,9 @@ double __ieee754_remainder(double x, double y) u.x=(u.x-d*w.x)-d*ww.x; if (ABS(u.x)<0.5*t.x) return (u.x!=0)?u.x:((x>0)?ZERO.x:nZERO.x); else - if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x; - else - {z=u.x/t.x; d=(z+big.x)-big.x; return ((u.x-d*w.x)-d*ww.x);} + if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x; + else + {z=u.x/t.x; d=(z+big.x)-big.x; return ((u.x-d*w.x)-d*ww.x);} } } /* (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000) */ @@ -128,3 +128,4 @@ double __ieee754_remainder(double x, double y) } } } +strong_alias (__ieee754_remainder, __remainder_finite) diff --git a/sysdeps/ieee754/dbl-64/e_sinh.c b/sysdeps/ieee754/dbl-64/e_sinh.c index 1701b9bb67..50463c3048 100644 --- a/sysdeps/ieee754/dbl-64/e_sinh.c +++ b/sysdeps/ieee754/dbl-64/e_sinh.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -15,15 +15,15 @@ static char rcsid[] = "$NetBSD: e_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $"; #endif /* __ieee754_sinh(x) - * Method : + * Method : * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2 - * 1. Replace x by |x| (sinh(-x) = -sinh(x)). - * 2. - * E + E/(E+1) + * 1. Replace x by |x| (sinh(-x) = -sinh(x)). + * 2. + * E + E/(E+1) * 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x) - * 2 + * 2 * - * 22 <= x <= lnovft : sinh(x) := exp(x)/2 + * 22 <= x <= lnovft : sinh(x) := exp(x)/2 * lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2) * ln2ovft < x : sinh(x) := x*shuge (overflow) * @@ -35,19 +35,11 @@ static char rcsid[] = "$NetBSD: e_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $"; #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const double one = 1.0, shuge = 1.0e307; -#else -static double one = 1.0, shuge = 1.0e307; -#endif -#ifdef __STDC__ - double __ieee754_sinh(double x) -#else - double __ieee754_sinh(x) - double x; -#endif -{ +double +__ieee754_sinh(double x) +{ double t,w,h; int32_t ix,jx; u_int32_t lx; @@ -57,14 +49,15 @@ static double one = 1.0, shuge = 1.0e307; ix = jx&0x7fffffff; /* x is INF or NaN */ - if(ix>=0x7ff00000) return x+x; + if(__builtin_expect(ix>=0x7ff00000, 0)) return x+x; h = 0.5; if (jx<0) h = -h; /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */ if (ix < 0x40360000) { /* |x|<22 */ - if (ix<0x3e300000) /* |x|<2**-28 */ - if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */ + if (__builtin_expect(ix<0x3e300000, 0)) /* |x|<2**-28 */ + if(shuge+x>one) + return x;/* sinh(tiny) = tiny with inexact */ t = __expm1(fabs(x)); if(ix<0x3ff00000) return h*(2.0*t-t*t/(t+one)); return h*(t+t/(t+one)); @@ -84,3 +77,4 @@ static double one = 1.0, shuge = 1.0e307; /* |x| > overflowthresold, sinh(x) overflow */ return x*shuge; } +strong_alias (__ieee754_sinh, __sinh_finite) diff --git a/sysdeps/ieee754/dbl-64/e_sqrt.c b/sysdeps/ieee754/dbl-64/e_sqrt.c index f7e8055491..05d1e71a0c 100644 --- a/sysdeps/ieee754/dbl-64/e_sqrt.c +++ b/sysdeps/ieee754/dbl-64/e_sqrt.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -86,3 +86,4 @@ double __ieee754_sqrt(double x) { return tm256.x*__ieee754_sqrt(x*t512.x); } } +strong_alias (__ieee754_sqrt, __sqrt_finite) diff --git a/sysdeps/ieee754/dbl-64/halfulp.c b/sysdeps/ieee754/dbl-64/halfulp.c index 478a4bacf6..42b21fb61d 100644 --- a/sysdeps/ieee754/dbl-64/halfulp.c +++ b/sysdeps/ieee754/dbl-64/halfulp.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001, 2005 Free Software Foundation + * Copyright (C) 2001, 2005, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -40,13 +40,11 @@ #include "dla.h" #include "math_private.h" -double __ieee754_sqrt(double x); - static const int4 tab54[32] = { 262143, 11585, 1782, 511, 210, 107, 63, 42, 30, 22, 17, 14, 12, 10, 9, 7, - 7, 6, 5, 5, 5, 4, 4, 4, - 3, 3, 3, 3, 3, 3, 3, 3 }; + 7, 6, 5, 5, 5, 4, 4, 4, + 3, 3, 3, 3, 3, 3, 3, 3 }; double __halfulp(double x, double y) @@ -64,12 +62,12 @@ double __halfulp(double x, double y) z = (double) k; return (z*y == -1075.0)?0: -10.0; } - /* if y > 0 */ + /* if y > 0 */ v.x = y; if (v.i[LOW_HALF] != 0) return -10.0; v.x=x; - /* case where x = 2**n for some integer n */ + /* case where x = 2**n for some integer n */ if (((v.i[HIGH_HALF]&0x000fffff)|v.i[LOW_HALF]) == 0) { k=(v.i[HIGH_HALF]>>20)-1023; return (((double) k)*y == -1075.0)?0:-10.0; @@ -90,7 +88,7 @@ double __halfulp(double x, double y) k = -k; if (k>5) return -10.0; - /* now treat x */ + /* now treat x */ while (k>0) { z = __ieee754_sqrt(x); EMULV(z,z,u,uu,j1,j2,j3,j4,j5); @@ -111,11 +109,11 @@ double __halfulp(double x, double y) m = (k&0x000fffff)|0x00100000; m = m>>(20-l); /* m is the odd integer of x */ - /* now check whether the length of m**n is at most 54 bits */ + /* now check whether the length of m**n is at most 54 bits */ if (m > tab54[n-3]) return -10.0; - /* yes, it is - now compute x**n by simple multiplications */ + /* yes, it is - now compute x**n by simple multiplications */ u = x; for (k=1;k<n;k++) u = u*x; diff --git a/sysdeps/ieee754/dbl-64/s_asinh.c b/sysdeps/ieee754/dbl-64/s_asinh.c index 985cfe32e1..93789fb41e 100644 --- a/sysdeps/ieee754/dbl-64/s_asinh.c +++ b/sysdeps/ieee754/dbl-64/s_asinh.c @@ -10,10 +10,6 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_asinh.c,v 1.9 1995/05/12 04:57:37 jtc Exp $"; -#endif - /* asinh(x) * Method : * Based on @@ -28,40 +24,34 @@ static char rcsid[] = "$NetBSD: s_asinh.c,v 1.9 1995/05/12 04:57:37 jtc Exp $"; #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const double -#else -static double -#endif one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ huge= 1.00000000000000000000e+300; -#ifdef __STDC__ - double __asinh(double x) -#else - double __asinh(x) - double x; -#endif +double +__asinh(double x) { - double t,w; + double w; int32_t hx,ix; GET_HIGH_WORD(hx,x); ix = hx&0x7fffffff; - if(ix>=0x7ff00000) return x+x; /* x is inf or NaN */ - if(ix< 0x3e300000) { /* |x|<2**-28 */ + if(__builtin_expect(ix< 0x3e300000, 0)) { /* |x|<2**-28 */ if(huge+x>one) return x; /* return x inexact except 0 */ } - if(ix>0x41b00000) { /* |x| > 2**28 */ + if(__builtin_expect(ix>0x41b00000, 0)) { /* |x| > 2**28 */ + if(ix>=0x7ff00000) return x+x; /* x is inf or NaN */ w = __ieee754_log(fabs(x))+ln2; - } else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */ - t = fabs(x); - w = __ieee754_log(2.0*t+one/(__ieee754_sqrt(x*x+one)+t)); - } else { /* 2.0 > |x| > 2**-28 */ - t = x*x; - w =__log1p(fabs(x)+t/(one+__ieee754_sqrt(one+t))); + } else { + double xa = fabs(x); + if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */ + w = __ieee754_log(2.0*xa+one/(__ieee754_sqrt(xa*xa+one)+xa)); + } else { /* 2.0 > |x| > 2**-28 */ + double t = xa*xa; + w =__log1p(xa+t/(one+__ieee754_sqrt(one+t))); + } } - if(hx>0) return w; else return -w; + return __copysign(w, x); } weak_alias (__asinh, asinh) #ifdef NO_LONG_DOUBLE |