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
This commit is contained in:
parent
752ba93078
commit
733331540a
Notes:
svn2git
2020-12-20 02:59:44 +00:00
svn path=/head/; revision=270947
@ -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;
|
||||
|
@ -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;
|
||||
|
Loading…
Reference in New Issue
Block a user