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

@@ -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) 

lib/msun/src/e_fmodf.c

@@ -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*/

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<uy.bits.exp) ||
 	       (ux.bits.manh<=uy.bits.manh &&

lib/msun/src/e_remainder.c

@@ -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 */

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 */

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
 

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<hy)||(lx<ly)) {
 		q = 0;

lib/msun/src/s_remquof.c

@@ -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 */

lib/msun/src/s_remquol.c

@@ -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 &&