Date: Tue, 24 Jul 2018 10:10:17 +0000 (UTC) From: Bruce Evans <bde@FreeBSD.org> To: src-committers@freebsd.org, svn-src-all@freebsd.org, svn-src-head@freebsd.org Subject: svn commit: r336663 - head/lib/msun/src Message-ID: <201807241010.w6OAAHpe031564@repo.freebsd.org>
next in thread | raw e-mail | index | archive | help
Author: bde Date: Tue Jul 24 10:10:16 2018 New Revision: 336663 URL: https://svnweb.freebsd.org/changeset/base/336663 Log: 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 Modified: head/lib/msun/src/e_fmod.c head/lib/msun/src/e_fmodf.c head/lib/msun/src/e_fmodl.c head/lib/msun/src/e_remainder.c head/lib/msun/src/e_remainderf.c head/lib/msun/src/math_private.h head/lib/msun/src/s_remquo.c head/lib/msun/src/s_remquof.c head/lib/msun/src/s_remquol.c Modified: head/lib/msun/src/e_fmod.c ============================================================================== --- head/lib/msun/src/e_fmod.c Tue Jul 24 08:15:02 2018 (r336662) +++ head/lib/msun/src/e_fmod.c Tue Jul 24 10:10:16 2018 (r336663) @@ -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<hy)||(lx<ly)) return x; /* |x|<|y| return x */ if(lx==ly) Modified: head/lib/msun/src/e_fmodf.c ============================================================================== --- head/lib/msun/src/e_fmodf.c Tue Jul 24 08:15:02 2018 (r336662) +++ head/lib/msun/src/e_fmodf.c Tue Jul 24 10:10:16 2018 (r336663) @@ -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 nan_mix(x, y)/nan_mix(x, y); + return nan_mix_op(x, y, *)/nan_mix_op(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*/ Modified: head/lib/msun/src/e_fmodl.c ============================================================================== --- head/lib/msun/src/e_fmodl.c Tue Jul 24 08:15:02 2018 (r336662) +++ head/lib/msun/src/e_fmodl.c Tue Jul 24 10:10:16 2018 (r336663) @@ -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<uy.bits.exp) || (ux.bits.manh<=uy.bits.manh && Modified: head/lib/msun/src/e_remainder.c ============================================================================== --- head/lib/msun/src/e_remainder.c Tue Jul 24 08:15:02 2018 (r336662) +++ head/lib/msun/src/e_remainder.c Tue Jul 24 10:10:16 2018 (r336663) @@ -49,7 +49,7 @@ __ieee754_remainder(double x, double p) (hx>=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 */ Modified: head/lib/msun/src/e_remainderf.c ============================================================================== --- head/lib/msun/src/e_remainderf.c Tue Jul 24 08:15:02 2018 (r336662) +++ head/lib/msun/src/e_remainderf.c Tue Jul 24 10:10:16 2018 (r336663) @@ -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 */ Modified: head/lib/msun/src/math_private.h ============================================================================== --- head/lib/msun/src/math_private.h Tue Jul 24 08:15:02 2018 (r336662) +++ head/lib/msun/src/math_private.h Tue Jul 24 10:10:16 2018 (r336663) @@ -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 Modified: head/lib/msun/src/s_remquo.c ============================================================================== --- head/lib/msun/src/s_remquo.c Tue Jul 24 08:15:02 2018 (r336662) +++ head/lib/msun/src/s_remquo.c Tue Jul 24 10:10:16 2018 (r336663) @@ -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<hy)||(lx<ly)) { q = 0; Modified: head/lib/msun/src/s_remquof.c ============================================================================== --- head/lib/msun/src/s_remquof.c Tue Jul 24 08:15:02 2018 (r336662) +++ head/lib/msun/src/s_remquof.c Tue Jul 24 10:10:16 2018 (r336663) @@ -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 nan_mix(x, y)/nan_mix(x, y); + return nan_mix_op(x, y, *)/nan_mix_op(x, y, *); if(hx<hy) { q = 0; goto fixup; /* |x|<|y| return x or x-y */ Modified: head/lib/msun/src/s_remquol.c ============================================================================== --- head/lib/msun/src/s_remquol.c Tue Jul 24 08:15:02 2018 (r336662) +++ head/lib/msun/src/s_remquol.c Tue Jul 24 10:10:16 2018 (r336663) @@ -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 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<uy.bits.exp) || (ux.bits.manh<=uy.bits.manh &&
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?201807241010.w6OAAHpe031564>