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

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;