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<hy)||(lx<ly)) return x;	/* |x|<|y| return x */
 	    if(lx==ly) 
diff --git a/lib/msun/src/e_fmodf.c b/lib/msun/src/e_fmodf.c
index f337bc854c62..1b6bf36f9cba 100644
--- a/lib/msun/src/e_fmodf.c
+++ b/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*/
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<uy.bits.exp) ||
 	       (ux.bits.manh<=uy.bits.manh &&
diff --git a/lib/msun/src/e_remainder.c b/lib/msun/src/e_remainder.c
index c2ba63cc997f..a4ae0b780fae 100644
--- a/lib/msun/src/e_remainder.c
+++ b/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 */
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<hy)||(lx<ly)) {
 		q = 0;
diff --git a/lib/msun/src/s_remquof.c b/lib/msun/src/s_remquof.c
index cbe618f4647a..f4c101620682 100644
--- a/lib/msun/src/s_remquof.c
+++ b/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 */
diff --git a/lib/msun/src/s_remquol.c b/lib/msun/src/s_remquol.c
index 2c7737949edc..2b03106d6062 100644
--- a/lib/msun/src/s_remquol.c
+++ b/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 &&