diff options
author | kargl <kargl@FreeBSD.org> | 2014-09-01 18:57:13 +0000 |
---|---|---|
committer | kargl <kargl@FreeBSD.org> | 2014-09-01 18:57:13 +0000 |
commit | 5bc23acdd22e6981a51a17e96d6102b6597fdae5 (patch) | |
tree | 5e64d62d8933e8a88fe7ebde51c2258b4a7138f0 /lib/msun/src | |
parent | b149a9e7aec5e222b7dba23fdebc5717058e1c33 (diff) | |
download | FreeBSD-src-5bc23acdd22e6981a51a17e96d6102b6597fdae5.zip FreeBSD-src-5bc23acdd22e6981a51a17e96d6102b6597fdae5.tar.gz |
According to the ISO C standard, lgamma(-integer) returns
inf and raises the divided-by-zero exception. Compilers
constant fold one/zero to inf but do not raise the exception.
Introduce a volatile vzero to prevent the constant folding.
Move the declaration of zero into the main declaration block.
While here, fix a nearby disordering of 'lx,ix'
Discussed with: bde
Diffstat (limited to 'lib/msun/src')
-rw-r--r-- | lib/msun/src/e_lgamma_r.c | 13 | ||||
-rw-r--r-- | lib/msun/src/e_lgammaf_r.c | 11 |
2 files changed, 13 insertions, 11 deletions
diff --git a/lib/msun/src/e_lgamma_r.c b/lib/msun/src/e_lgamma_r.c index 87b8ea2..6af9d2c 100644 --- a/lib/msun/src/e_lgamma_r.c +++ b/lib/msun/src/e_lgamma_r.c @@ -86,7 +86,10 @@ __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" +static const volatile double vzero = 0; + static const double +zero= 0.00000000000000000000e+00, two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ @@ -154,8 +157,6 @@ w4 = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */ w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */ w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */ -static const double zero= 0.00000000000000000000e+00; - /* * Compute sin(pi*x) without actually doing the pi*x multiplication. * sin_pi(x) is only called for x < 0 and |x| < 2**(p-1) where p is @@ -204,7 +205,7 @@ __ieee754_lgamma_r(double x, int *signgamp) { double t,y,z,nadj,p,p1,p2,p3,q,r,w; int32_t hx; - int i,lx,ix; + int i,ix,lx; EXTRACT_WORDS(hx,lx,x); @@ -212,7 +213,7 @@ __ieee754_lgamma_r(double x, int *signgamp) *signgamp = 1; ix = hx&0x7fffffff; if(ix>=0x7ff00000) return x*x; - if((ix|lx)==0) return one/zero; + if((ix|lx)==0) return one/vzero; if(ix<0x3b900000) { /* |x|<2**-70, return -log(|x|) */ if(hx<0) { *signgamp = -1; @@ -221,9 +222,9 @@ __ieee754_lgamma_r(double x, int *signgamp) } if(hx<0) { if(ix>=0x43300000) /* |x|>=2**52, must be -integer */ - return one/zero; + return one/vzero; t = sin_pi(x); - if(t==zero) return one/zero; /* -integer */ + if(t==zero) return one/vzero; /* -integer */ nadj = __ieee754_log(pi/fabs(t*x)); if(t<zero) *signgamp = -1; x = -x; diff --git a/lib/msun/src/e_lgammaf_r.c b/lib/msun/src/e_lgammaf_r.c index d13484e..2aaa656 100644 --- a/lib/msun/src/e_lgammaf_r.c +++ b/lib/msun/src/e_lgammaf_r.c @@ -19,7 +19,10 @@ __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" +static const volatile float vzero = 0; + static const float +zero= 0.0000000000e+00, two23= 8.3886080000e+06, /* 0x4b000000 */ half= 5.0000000000e-01, /* 0x3f000000 */ one = 1.0000000000e+00, /* 0x3f800000 */ @@ -87,8 +90,6 @@ w4 = -5.9518753551e-04, /* 0xba1c065c */ w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */ w6 = -1.6309292987e-03; /* 0xbad5c4e8 */ -static const float zero= 0.0000000000e+00; - static float sin_pif(float x) { @@ -140,7 +141,7 @@ __ieee754_lgammaf_r(float x, int *signgamp) *signgamp = 1; ix = hx&0x7fffffff; if(ix>=0x7f800000) return x*x; - if(ix==0) return one/zero; + if(ix==0) return one/vzero; if(ix<0x35000000) { /* |x|<2**-21, return -log(|x|) */ if(hx<0) { *signgamp = -1; @@ -149,9 +150,9 @@ __ieee754_lgammaf_r(float x, int *signgamp) } if(hx<0) { if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */ - return one/zero; + return one/vzero; t = sin_pif(x); - if(t==zero) return one/zero; /* -integer */ + if(t==zero) return one/vzero; /* -integer */ nadj = __ieee754_logf(pi/fabsf(t*x)); if(t<zero) *signgamp = -1; x = -x; |