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 &&
_______________________________________________
svn-src-head@freebsd.org mailing list
https://lists.freebsd.org/mailman/listinfo/svn-src-head
To unsubscribe, send any mail to "svn-src-head-unsubscr...@freebsd.org"
