From daa1e391104cea6c63abbb3f1b2b6785f4c986fa Mon Sep 17 00:00:00 2001 From: Bruce Evans Date: Tue, 24 Jul 2018 10:10:16 +0000 Subject: [PATCH] Fix the conversion to use nan_mix() in r336362. fmod*(x, y), remainder*(x, y) and remquo*(x, y, quo) were broken for y = 0 by changing multiplication by y to addition of y. (When y is 0, the result should be NaN but became 1 for finite x.) Use a new macro nan_mix_op() to give more control over the mixing, and expand comments. Recent re-testing missed finding this bug since I only tested the macro version on amd64 and i386 and these arches don't use the C versions (they use either asm versions or builtins). Reported by: enh via freebsd-numerics --- lib/msun/src/e_fmod.c | 2 +- lib/msun/src/e_fmodf.c | 2 +- lib/msun/src/e_fmodl.c | 2 +- lib/msun/src/e_remainder.c | 2 +- lib/msun/src/e_remainderf.c | 2 +- lib/msun/src/math_private.h | 23 +++++++++++++++-------- lib/msun/src/s_remquo.c | 2 +- lib/msun/src/s_remquof.c | 2 +- lib/msun/src/s_remquol.c | 2 +- 9 files changed, 23 insertions(+), 16 deletions(-) diff --git a/lib/msun/src/e_fmod.c b/lib/msun/src/e_fmod.c index 6c5c865b29ee..3a28dc4ff1f3 100644 --- a/lib/msun/src/e_fmod.c +++ b/lib/msun/src/e_fmod.c @@ -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 nan_mix(x, y)/nan_mix(x, y); + return nan_mix_op(x, y, *)/nan_mix_op(x, y, *); if(hx<=hy) { if((hx=0x7f800000)|| /* y=0,or x not finite */ (hy>0x7f800000)) /* or y is NaN */ - return nan_mix(x, y)/nan_mix(x, y); + return nan_mix_op(x, y, *)/nan_mix_op(x, y, *); if(hx>31]; /* |x|=|y| return x*0*/ diff --git a/lib/msun/src/e_fmodl.c b/lib/msun/src/e_fmodl.c index 6ec899727af3..ad3bcc34c9a2 100644 --- a/lib/msun/src/e_fmodl.c +++ b/lib/msun/src/e_fmodl.c @@ -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 nan_mix(x, y)/nan_mix(x, y); + return nan_mix_op(x, y, *)/nan_mix_op(x, y, *); if(ux.bits.exp<=uy.bits.exp) { if((ux.bits.exp=0x7ff00000)|| /* x not finite */ ((hp>=0x7ff00000)&& /* p is NaN */ (((hp-0x7ff00000)|lp)!=0))) - return nan_mix(x, p)/nan_mix(x, p); + return nan_mix_op(x, p, *)/nan_mix_op(x, p, *); if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p); /* now x < 2p */ diff --git a/lib/msun/src/e_remainderf.c b/lib/msun/src/e_remainderf.c index 2b83fe0ec967..8004493de77b 100644 --- a/lib/msun/src/e_remainderf.c +++ b/lib/msun/src/e_remainderf.c @@ -39,7 +39,7 @@ __ieee754_remainderf(float x, float p) if((hp==0)|| /* p = 0 */ (hx>=0x7f800000)|| /* x not finite */ ((hp>0x7f800000))) /* p is NaN */ - return nan_mix(x, p)/nan_mix(x, p); + return nan_mix_op(x, p, *)/nan_mix_op(x, p, *); if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p); /* now x < 2p */ diff --git a/lib/msun/src/math_private.h b/lib/msun/src/math_private.h index 4006654cf3fe..b91b54cea689 100644 --- a/lib/msun/src/math_private.h +++ b/lib/msun/src/math_private.h @@ -479,22 +479,29 @@ 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 + * Mix 0, 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. + * args using the specified operation. * - * 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 + * When one arg is NaN, the result is typically that arg quieted. When both + * args are NaNs, the result is typically the quietening of the arg whose + * mantissa is largest after quietening. When neither arg is NaN, the + * result may be NaN because it is indeterminate, or finite for subsequent + * construction of a NaN as the indeterminate 0.0L/0.0L. + * + * Technical complications: the result in bits after rounding to the final + * precision might depend on the runtime 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 + * runtime 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)) +#define nan_mix(x, y) (nan_mix_op((x), (y), +)) +#define nan_mix_op(x, y, op) (((x) + 0.0L) op ((y) + 0)) #ifdef _COMPLEX_H diff --git a/lib/msun/src/s_remquo.c b/lib/msun/src/s_remquo.c index 9c61fa7487e5..6a111dfc6e10 100644 --- a/lib/msun/src/s_remquo.c +++ b/lib/msun/src/s_remquo.c @@ -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 nan_mix(x, y)/nan_mix(x, y); + return nan_mix_op(x, y, *)/nan_mix_op(x, y, *); if(hx<=hy) { if((hx=0x7f800000||hy>0x7f800000) /* y=0,NaN;or x not finite */ - return nan_mix(x, y)/nan_mix(x, y); + return nan_mix_op(x, y, *)/nan_mix_op(x, y, *); if(hx