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