d/freebsd-dev
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add a macro nan_mix() and use it to get NaN results that are (bitwise)

Browse Source 
independent of the precision in most cases.  This is mainly to simplify
checking for errors.  r176266 did this for e_pow[f].c using a less
refined expression that often didn't work.  r176276 fixes an error in
the log message for r176266.  The main refinement is to always expand
to long double precision.  See old log messages (especially these 2)
and the comment on the macro for more general details.

Specific details:
- using nan_mix() consistently for the new and old pow*() functions was
  the only thing needed to make my consistency test for powl() vs pow()
  pass on amd64.

- catrig[fl].c already had all the refinements, but open-coded.

- e_atan2[fl].c, e_fmod[fl].c and s_remquo[fl] only had primitive NaN
  mixing.

- e_hypot[fl].c already had a different refined version of r176266.  Refine
  this further.  nan_mix() is not directly usable here since we want to
  clear the sign bit.

- e_remainder[f].c already had an earlier version of r176266.

- s_ccosh[f].c,/s_csinh[f].c already had a version equivalent to r176266.
  Refine this further.  nan_mix() is not directly usable here since the
  expression has to handle some non-NaN cases.

- s_csqrt.[fl]: the mixing was special and mostly wrong.  Partially fix the
  special version.

- s_ctanh[f].c already had a version of r176266.
This commit is contained in:
Bruce Evans 2018-07-17 07:42:14 +00:00
parent 15430057b3
commit 6f1b8a0792
Notes: svn2git 2020-12-20 02:59:44 +00:00 
svn path=/head/; revision=336362
31 changed files with 77 additions and 55 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
lib/msun/ld128/e_powl.c
View File

@ -182,7 +182,7 @@ powl(long double x, long double y)
|| (iy > 0x7fff0000)
|| ((iy == 0x7fff0000)
&& ((q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) != 0)))
return x + y;
return nan_mix(x, y);
/* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer

4
lib/msun/ld80/e_powl.c
View File

@ -216,9 +216,9 @@ if( x == 1.0L )
return( 1.0L );
if( isnan(x) )
return( x );
return ( nan_mix(x, y) );
if( isnan(y) )
return( y );
return ( nan_mix(x, y) );
if( y == 1.0L )
return( x );

6
lib/msun/src/catrig.c
View File

@ -300,7 +300,7 @@ casinh(double complex z)
* C99 leaves it optional whether to raise invalid if one of
* the arguments is not NaN, so we opt not to raise it.
*/
return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLX(nan_mix(x, y), nan_mix(x, y)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@ -384,7 +384,7 @@ cacos(double complex z)
* C99 leaves it optional whether to raise invalid if one of
* the arguments is not NaN, so we opt not to raise it.
*/
return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLX(nan_mix(x, y), nan_mix(x, y)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@ -601,7 +601,7 @@ catanh(double complex z)
* C99 leaves it optional whether to raise invalid if one of
* the arguments is not NaN, so we opt not to raise it.
*/
return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLX(nan_mix(x, y), nan_mix(x, y)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)

6
lib/msun/src/catrigf.c
View File

@ -163,7 +163,7 @@ casinhf(float complex z)
return (CMPLXF(y, x + x));
if (y == 0)
return (CMPLXF(x + x, y));
return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLXF(nan_mix(x, y), nan_mix(x, y)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@ -221,7 +221,7 @@ cacosf(float complex z)
return (CMPLXF(x + x, -y));
if (x == 0)
return (CMPLXF(pio2_hi + pio2_lo, y + y));
return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLXF(nan_mix(x, y), nan_mix(x, y)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@ -359,7 +359,7 @@ catanhf(float complex z)
if (isinf(y))
return (CMPLXF(copysignf(0, x),
copysignf(pio2_hi + pio2_lo, y)));
return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLXF(nan_mix(x, y), nan_mix(x, y)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)

6
lib/msun/src/catrigl.c
View File

@ -182,7 +182,7 @@ casinhl(long double complex z)
return (CMPLXL(y, x + x));
if (y == 0)
return (CMPLXL(x + x, y));
return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLXL(nan_mix(x, y), nan_mix(x, y)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@ -241,7 +241,7 @@ cacosl(long double complex z)
return (CMPLXL(x + x, -y));
if (x == 0)
return (CMPLXL(pio2_hi + pio2_lo, y + y));
return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLXL(nan_mix(x, y), nan_mix(x, y)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@ -380,7 +380,7 @@ catanhl(long double complex z)
if (isinf(y))
return (CMPLXL(copysignl(0, x),
copysignl(pio2_hi + pio2_lo, y)));
return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLXL(nan_mix(x, y), nan_mix(x, y)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)

2
lib/msun/src/e_atan2.c
View File

@ -70,7 +70,7 @@ __ieee754_atan2(double y, double x)
iy = hy&0x7fffffff;
if(((ix|((lx|-lx)>>31))>0x7ff00000)||
((iy|((ly|-ly)>>31))>0x7ff00000)) /* x or y is NaN */
return x+y;
return nan_mix(x, y);
if(hx==0x3ff00000&&lx==0) return atan(y); /* x=1.0 */
m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */

2
lib/msun/src/e_atan2f.c
View File

@ -41,7 +41,7 @@ __ieee754_atan2f(float y, float x)
iy = hy&0x7fffffff;
if((ix>0x7f800000)||
(iy>0x7f800000)) /* x or y is NaN */
return x+y;
return nan_mix(x, y);
if(hx==0x3f800000) return atanf(y); /* x=1.0 */
m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */

2
lib/msun/src/e_atan2l.c
View File

@ -62,7 +62,7 @@ atan2l(long double y, long double x)
((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)!=0) || /* x is NaN */
(expty==BIAS+LDBL_MAX_EXP &&
((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* y is NaN */
return x+y;
return nan_mix(x, y);
if (expsignx==BIAS && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)==0)
return atanl(y); /* x=1.0 */
m = ((expsigny>>15)&1)|((expsignx>>14)&2); /* 2*sign(x)+sign(y) */

2
lib/msun/src/e_fmod.c
View File

@ -42,7 +42,7 @@ __ieee754_fmod(double x, double y)
/* purge off exception values */
if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */
((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */
return (x*y)/(x*y);
return nan_mix(x, y)/nan_mix(x, y);
if(hx<=hy) {
if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
if(lx==ly)

2
lib/msun/src/e_fmodf.c
View File

@ -41,7 +41,7 @@ __ieee754_fmodf(float x, float y)
/* purge off exception values */
if(hy==0||(hx>=0x7f800000)|| /* y=0,or x not finite */
(hy>0x7f800000)) /* or y is NaN */
return (x*y)/(x*y);
return nan_mix(x, y)/nan_mix(x, y);
if(hx<hy) return x; /* |x|<|y| return x */
if(hx==hy)
return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/

2
lib/msun/src/e_fmodl.c
View File

@ -79,7 +79,7 @@ fmodl(long double x, long double y)
(ux.bits.exp == BIAS + LDBL_MAX_EXP) || /* or x not finite */
(uy.bits.exp == BIAS + LDBL_MAX_EXP &&
((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* or y is NaN */
return (x*y)/(x*y);
return nan_mix(x, y)/nan_mix(x, y);
if(ux.bits.exp<=uy.bits.exp) {
if((ux.bits.exp<uy.bits.exp) ||
(ux.bits.manh<=uy.bits.manh &&

2
lib/msun/src/e_hypot.c
View File

@ -70,7 +70,7 @@ __ieee754_hypot(double x, double y)
if(ha >= 0x7ff00000) { /* Inf or NaN */
u_int32_t low;
/* Use original arg order iff result is NaN; quieten sNaNs. */
w = fabs(x+0.0)-fabs(y+0.0);
w = fabsl(x+0.0L)-fabs(y+0);
GET_LOW_WORD(low,a);
if(((ha&0xfffff)|low)==0) w = a;
GET_LOW_WORD(low,b);

2
lib/msun/src/e_hypotf.c
View File

@ -37,7 +37,7 @@ __ieee754_hypotf(float x, float y)
if(ha > 0x58800000) { /* a>2**50 */
if(ha >= 0x7f800000) { /* Inf or NaN */
/* Use original arg order iff result is NaN; quieten sNaNs. */
w = fabsf(x+0.0F)-fabsf(y+0.0F);
w = fabsl(x+0.0L)-fabsf(y+0);
if(ha == 0x7f800000) w = a;
if(hb == 0x7f800000) w = b;
return w;

2
lib/msun/src/e_hypotl.c
View File

@ -64,7 +64,7 @@ hypotl(long double x, long double y)
if(ha >= ESW(MAX_EXP)) { /* Inf or NaN */
man_t manh, manl;
/* Use original arg order iff result is NaN; quieten sNaNs. */
w = fabsl(x+0.0)-fabsl(y+0.0);
w = fabsl(x+0.0L)-fabsl(y+0);
GET_LDBL_MAN(manh,manl,a);
if (manh == LDBL_NBIT && manl == 0) w = a;
GET_LDBL_MAN(manh,manl,b);

2
lib/msun/src/e_pow.c
View File

@ -119,7 +119,7 @@ __ieee754_pow(double x, double y)
/* y!=zero: result is NaN if either arg is NaN */
if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
return (x+0.0)+(y+0.0);
return nan_mix(x, y);
/* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer

2
lib/msun/src/e_powf.c
View File

@ -76,7 +76,7 @@ __ieee754_powf(float x, float y)
/* y!=zero: result is NaN if either arg is NaN */
if(ix > 0x7f800000 ||
iy > 0x7f800000)
return (x+0.0F)+(y+0.0F);
return nan_mix(x, y);
/* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer

6
lib/msun/src/e_remainder.c
View File

@ -45,11 +45,11 @@ __ieee754_remainder(double x, double p)
hx &= 0x7fffffff;
/* purge off exception values */
if((hp|lp)==0) return (x*p)/(x*p); /* p = 0 */
if((hx>=0x7ff00000)|| /* x not finite */
if(((hp|lp)==0)|| /* p = 0 */
(hx>=0x7ff00000)|| /* x not finite */
((hp>=0x7ff00000)&& /* p is NaN */
(((hp-0x7ff00000)|lp)!=0)))
return ((long double)x*p)/((long double)x*p);
return nan_mix(x, p)/nan_mix(x, p);
if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p); /* now x < 2p */

6
lib/msun/src/e_remainderf.c
View File

@ -36,10 +36,10 @@ __ieee754_remainderf(float x, float p)
hx &= 0x7fffffff;
/* purge off exception values */
if(hp==0) return (x*p)/(x*p); /* p = 0 */
if((hx>=0x7f800000)|| /* x not finite */
if((hp==0)|| /* p = 0 */
(hx>=0x7f800000)|| /* x not finite */
((hp>0x7f800000))) /* p is NaN */
return ((long double)x*p)/((long double)x*p);
return nan_mix(x, p)/nan_mix(x, p);
if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p); /* now x < 2p */

18
lib/msun/src/math_private.h
View File

@ -478,6 +478,24 @@ do { \
*/
void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
/*
* Mix 1 or 2 NaNs. First add 0 to each arg. This normally just turns
* signaling NaNs into quiet NaNs by setting a quiet bit. We do this
* because we want to never return a signaling NaN, and also because we
* don't want the quiet bit to affect the result. Then mix the converted
* args using addition. The result is typically the arg whose mantissa
* bits (considered as in integer) are largest.
*
* Technical complications: the result in bits might depend on the precision
* and/or on compiler optimizations, especially when different register sets
* are used for different precisions. Try to make the result not depend on
* at least the precision by always doing the main mixing step in long double
* precision. Try to reduce dependencies on optimizations by adding the
* the 0's in different precisions (unless everything is in long double
* precision).
*/
#define nan_mix(x, y) (((x) + 0.0L) + ((y) + 0))
#ifdef _COMPLEX_H
/*

3
lib/msun/src/s_ccosh.c
View File

@ -146,7 +146,8 @@ ccosh(double complex z)
* Optionally raises the invalid floating-point exception for finite
* nonzero y. Choice = don't raise (except for signaling NaNs).
*/
return (CMPLX((x * x) * (y - y), (x + x) * (y - y)));
return (CMPLX(((long double)x * x) * (y - y),
((long double)x + x) * (y - y)));
}
double complex

3
lib/msun/src/s_ccoshf.c
View File

@ -92,7 +92,8 @@ ccoshf(float complex z)
return (CMPLXF(INFINITY * cosf(y), x * sinf(y)));
}
return (CMPLXF((x * x) * (y - y), (x + x) * (y - y)));
return (CMPLXF(((long double)x * x) * (y - y),
((long double)x + x) * (y - y)));
}
float complex

3
lib/msun/src/s_csinh.c
View File

@ -145,7 +145,8 @@ csinh(double complex z)
* Optionally raises the invalid floating-point exception for finite
* nonzero y. Choice = don't raise (except for signaling NaNs).
*/
return (CMPLX((x + x) * (y - y), (x * x) * (y - y)));
return (CMPLX(((long double)x + x) * (y - y),
((long double)x * x) * (y - y)));
}
double complex

3
lib/msun/src/s_csinhf.c
View File

@ -92,7 +92,8 @@ csinhf(float complex z)
return (CMPLXF(x * cosf(y), INFINITY * sinf(y)));
}
return (CMPLXF((x + x) * (y - y), (x * x) * (y - y)));
return (CMPLXF(((long double)x + x) * (y - y),
((long double)x * x) * (y - y)));
}
float complex

10
lib/msun/src/s_csqrt.c
View File

@ -65,7 +65,7 @@ csqrt(double complex z)
return (CMPLX(INFINITY, b));
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
return (CMPLX(a, t)); /* return NaN + NaN i */
return (CMPLX(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */
}
if (isinf(a)) {
/*
@ -79,10 +79,10 @@ csqrt(double complex z)
else
return (CMPLX(a, copysign(b - b, b)));
}
/*
* The remaining special case (b is NaN) is handled just fine by
* the normal code path below.
*/
if (isnan(b)) {
t = (a - a) / (a - a); /* raise invalid */
return (CMPLX(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */
}
/* Scale to avoid overflow. */
if (fabs(a) >= THRESH || fabs(b) >= THRESH) {

10
lib/msun/src/s_csqrtf.c
View File

@ -56,7 +56,7 @@ csqrtf(float complex z)
return (CMPLXF(INFINITY, b));
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
return (CMPLXF(a, t)); /* return NaN + NaN i */
return (CMPLXF(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */
}
if (isinf(a)) {
/*
@ -70,10 +70,10 @@ csqrtf(float complex z)
else
return (CMPLXF(a, copysignf(b - b, b)));
}
/*
* The remaining special case (b is NaN) is handled just fine by
* the normal code path below.
*/
if (isnan(b)) {
t = (a - a) / (a - a); /* raise invalid */
return (CMPLXF(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */
}
/*
* We compute t in double precision to avoid overflow and to

10
lib/msun/src/s_csqrtl.c
View File

@ -73,7 +73,7 @@ csqrtl(long double complex z)
return (CMPLXL(INFINITY, b));
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
return (CMPLXL(a, t)); /* return NaN + NaN i */
return (CMPLXL(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */
}
if (isinf(a)) {
/*
@ -87,10 +87,10 @@ csqrtl(long double complex z)
else
return (CMPLXL(a, copysignl(b - b, b)));
}
/*
* The remaining special case (b is NaN) is handled just fine by
* the normal code path below.
*/
if (isnan(b)) {
t = (a - a) / (a - a); /* raise invalid */
return (CMPLXL(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */
}
/* Scale to avoid overflow. */
if (fabsl(a) >= THRESH || fabsl(b) >= THRESH) {

4
lib/msun/src/s_ctanh.c
View File

@ -104,8 +104,8 @@ ctanh(double complex z)
*/
if (ix >= 0x7ff00000) {
if ((ix & 0xfffff) | lx) /* x is NaN */
return (CMPLX((x + 0) * (y + 0),
y == 0 ? y : (x + 0) * (y + 0)));
return (CMPLX(nan_mix(x, y),
y == 0 ? y : nan_mix(x, y)));
SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */
return (CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y))));
}

4
lib/msun/src/s_ctanhf.c
View File

@ -53,8 +53,8 @@ ctanhf(float complex z)
if (ix >= 0x7f800000) {
if (ix & 0x7fffff)
return (CMPLXF((x + 0) * (y + 0),
y == 0 ? y : (x + 0) * (y + 0)));
return (CMPLXF(nan_mix(x, y),
y == 0 ? y : nan_mix(x, y)));
SET_FLOAT_WORD(x, hx - 0x40000000);
return (CMPLXF(x,
copysignf(0, isinf(y) ? y : sinf(y) * cosf(y))));

2
lib/msun/src/s_remquo.c
View File

@ -44,7 +44,7 @@ remquo(double x, double y, int *quo)
/* purge off exception values */
if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */
((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */
return (x*y)/(x*y);
return nan_mix(x, y)/nan_mix(x, y);
if(hx<=hy) {
if((hx<hy)||(lx<ly)) {
q = 0;

2
lib/msun/src/s_remquof.c
View File

@ -41,7 +41,7 @@ remquof(float x, float y, int *quo)
/* purge off exception values */
if(hy==0||hx>=0x7f800000||hy>0x7f800000) /* y=0,NaN;or x not finite */
return (x*y)/(x*y);
return nan_mix(x, y)/nan_mix(x, y);
if(hx<hy) {
q = 0;
goto fixup; /* |x|<|y| return x or x-y */

2
lib/msun/src/s_remquol.c
View File

@ -86,7 +86,7 @@ remquol(long double x, long double y, int *quo)
(ux.bits.exp == BIAS + LDBL_MAX_EXP) || /* or x not finite */
(uy.bits.exp == BIAS + LDBL_MAX_EXP &&
((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* or y is NaN */
return (x*y)/(x*y);
return nan_mix(x, y)/nan_mix(x, y);
if(ux.bits.exp<=uy.bits.exp) {
if((ux.bits.exp<uy.bits.exp) ||
(ux.bits.manh<=uy.bits.manh &&