diff --git a/lib/msun/ld128/e_powl.c b/lib/msun/ld128/e_powl.c
index e91d6d9335d9..d6e72333e4b9 100644
--- a/lib/msun/ld128/e_powl.c
+++ b/lib/msun/ld128/e_powl.c
@@ -182,7 +182,7 @@ powl(long double x, long double y)
       || (iy > 0x7fff0000)
       || ((iy == 0x7fff0000)
 	  && ((q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) != 0)))
-    return x + y;
+    return nan_mix(x, y);
 
   /* determine if y is an odd int when x < 0
    * yisint = 0       ... y is not an integer
diff --git a/lib/msun/ld80/e_powl.c b/lib/msun/ld80/e_powl.c
index 2bae94b58cde..87c75688d61d 100644
--- a/lib/msun/ld80/e_powl.c
+++ b/lib/msun/ld80/e_powl.c
@@ -216,9 +216,9 @@ if( x == 1.0L )
 	return( 1.0L );
 
 if( isnan(x) )
-	return( x );
+	return ( nan_mix(x, y) );
 if( isnan(y) )
-	return( y );
+	return ( nan_mix(x, y) );
 
 if( y == 1.0L )
 	return( x );
diff --git a/lib/msun/src/catrig.c b/lib/msun/src/catrig.c
index 6ca35af15ff9..431d8e6742fe 100644
--- a/lib/msun/src/catrig.c
+++ b/lib/msun/src/catrig.c
@@ -300,7 +300,7 @@ casinh(double complex z)
 		 * C99 leaves it optional whether to raise invalid if one of
 		 * the arguments is not NaN, so we opt not to raise it.
 		 */
-		return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+		return (CMPLX(nan_mix(x, y), nan_mix(x, y)));
 	}
 
 	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@@ -384,7 +384,7 @@ cacos(double complex z)
 		 * C99 leaves it optional whether to raise invalid if one of
 		 * the arguments is not NaN, so we opt not to raise it.
 		 */
-		return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+		return (CMPLX(nan_mix(x, y), nan_mix(x, y)));
 	}
 
 	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@@ -601,7 +601,7 @@ catanh(double complex z)
 		 * C99 leaves it optional whether to raise invalid if one of
 		 * the arguments is not NaN, so we opt not to raise it.
 		 */
-		return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+		return (CMPLX(nan_mix(x, y), nan_mix(x, y)));
 	}
 
 	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
diff --git a/lib/msun/src/catrigf.c b/lib/msun/src/catrigf.c
index 9533aab71751..62bcd3995603 100644
--- a/lib/msun/src/catrigf.c
+++ b/lib/msun/src/catrigf.c
@@ -163,7 +163,7 @@ casinhf(float complex z)
 			return (CMPLXF(y, x + x));
 		if (y == 0)
 			return (CMPLXF(x + x, y));
-		return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+		return (CMPLXF(nan_mix(x, y), nan_mix(x, y)));
 	}
 
 	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@@ -221,7 +221,7 @@ cacosf(float complex z)
 			return (CMPLXF(x + x, -y));
 		if (x == 0)
 			return (CMPLXF(pio2_hi + pio2_lo, y + y));
-		return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+		return (CMPLXF(nan_mix(x, y), nan_mix(x, y)));
 	}
 
 	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@@ -359,7 +359,7 @@ catanhf(float complex z)
 		if (isinf(y))
 			return (CMPLXF(copysignf(0, x),
 			    copysignf(pio2_hi + pio2_lo, y)));
-		return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+		return (CMPLXF(nan_mix(x, y), nan_mix(x, y)));
 	}
 
 	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
diff --git a/lib/msun/src/catrigl.c b/lib/msun/src/catrigl.c
index b79857082700..e66f87a24a4b 100644
--- a/lib/msun/src/catrigl.c
+++ b/lib/msun/src/catrigl.c
@@ -182,7 +182,7 @@ casinhl(long double complex z)
 			return (CMPLXL(y, x + x));
 		if (y == 0)
 			return (CMPLXL(x + x, y));
-		return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+		return (CMPLXL(nan_mix(x, y), nan_mix(x, y)));
 	}
 
 	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@@ -241,7 +241,7 @@ cacosl(long double complex z)
 			return (CMPLXL(x + x, -y));
 		if (x == 0)
 			return (CMPLXL(pio2_hi + pio2_lo, y + y));
-		return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+		return (CMPLXL(nan_mix(x, y), nan_mix(x, y)));
 	}
 
 	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@@ -380,7 +380,7 @@ catanhl(long double complex z)
 		if (isinf(y))
 			return (CMPLXL(copysignl(0, x),
 			    copysignl(pio2_hi + pio2_lo, y)));
-		return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+		return (CMPLXL(nan_mix(x, y), nan_mix(x, y)));
 	}
 
 	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
diff --git a/lib/msun/src/e_atan2.c b/lib/msun/src/e_atan2.c
index 352db171d219..231a1611ee12 100644
--- a/lib/msun/src/e_atan2.c
+++ b/lib/msun/src/e_atan2.c
@@ -70,7 +70,7 @@ __ieee754_atan2(double y, double x)
 	iy = hy&0x7fffffff;
 	if(((ix|((lx|-lx)>>31))>0x7ff00000)||
 	   ((iy|((ly|-ly)>>31))>0x7ff00000))	/* x or y is NaN */
-	   return x+y;
+	    return nan_mix(x, y);
 	if(hx==0x3ff00000&&lx==0) return atan(y);   /* x=1.0 */
 	m = ((hy>>31)&1)|((hx>>30)&2);	/* 2*sign(x)+sign(y) */
 
diff --git a/lib/msun/src/e_atan2f.c b/lib/msun/src/e_atan2f.c
index fc77bffac273..346d76746c05 100644
--- a/lib/msun/src/e_atan2f.c
+++ b/lib/msun/src/e_atan2f.c
@@ -41,7 +41,7 @@ __ieee754_atan2f(float y, float x)
 	iy = hy&0x7fffffff;
 	if((ix>0x7f800000)||
 	   (iy>0x7f800000))	/* x or y is NaN */
-	   return x+y;
+	    return nan_mix(x, y);
 	if(hx==0x3f800000) return atanf(y);   /* x=1.0 */
 	m = ((hy>>31)&1)|((hx>>30)&2);	/* 2*sign(x)+sign(y) */
 
diff --git a/lib/msun/src/e_atan2l.c b/lib/msun/src/e_atan2l.c
index 032648248999..94ebdec54322 100644
--- a/lib/msun/src/e_atan2l.c
+++ b/lib/msun/src/e_atan2l.c
@@ -62,7 +62,7 @@ atan2l(long double y, long double x)
 	     ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)!=0) ||	/* x is NaN */
 	    (expty==BIAS+LDBL_MAX_EXP &&
 	     ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0))	/* y is NaN */
-	    return x+y;
+	    return nan_mix(x, y);
 	if (expsignx==BIAS && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)==0)
 	    return atanl(y);					/* x=1.0 */
 	m = ((expsigny>>15)&1)|((expsignx>>14)&2);	/* 2*sign(x)+sign(y) */
diff --git a/lib/msun/src/e_fmod.c b/lib/msun/src/e_fmod.c
index 970278ee850f..6c5c865b29ee 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 (x*y)/(x*y);
+	    return nan_mix(x, y)/nan_mix(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 52ce373ab8ee..f337bc854c62 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 (x*y)/(x*y);
+	    return nan_mix(x, y)/nan_mix(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 e315f761d051..6ec899727af3 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 (x*y)/(x*y);
+	    return nan_mix(x, y)/nan_mix(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_hypot.c b/lib/msun/src/e_hypot.c
index 2398e9837fec..e67e6935575d 100644
--- a/lib/msun/src/e_hypot.c
+++ b/lib/msun/src/e_hypot.c
@@ -70,7 +70,7 @@ __ieee754_hypot(double x, double y)
 	   if(ha >= 0x7ff00000) {	/* Inf or NaN */
 	       u_int32_t low;
 	       /* Use original arg order iff result is NaN; quieten sNaNs. */
-	       w = fabs(x+0.0)-fabs(y+0.0);
+	       w = fabsl(x+0.0L)-fabs(y+0);
 	       GET_LOW_WORD(low,a);
 	       if(((ha&0xfffff)|low)==0) w = a;
 	       GET_LOW_WORD(low,b);
diff --git a/lib/msun/src/e_hypotf.c b/lib/msun/src/e_hypotf.c
index 6d083e4d2c3f..bd45c6195128 100644
--- a/lib/msun/src/e_hypotf.c
+++ b/lib/msun/src/e_hypotf.c
@@ -37,7 +37,7 @@ __ieee754_hypotf(float x, float y)
 	if(ha > 0x58800000) {	/* a>2**50 */
 	   if(ha >= 0x7f800000) {	/* Inf or NaN */
 	       /* Use original arg order iff result is NaN; quieten sNaNs. */
-	       w = fabsf(x+0.0F)-fabsf(y+0.0F);
+	       w = fabsl(x+0.0L)-fabsf(y+0);
 	       if(ha == 0x7f800000) w = a;
 	       if(hb == 0x7f800000) w = b;
 	       return w;
diff --git a/lib/msun/src/e_hypotl.c b/lib/msun/src/e_hypotl.c
index 7b5ab89562fe..9189b1fab54d 100644
--- a/lib/msun/src/e_hypotl.c
+++ b/lib/msun/src/e_hypotl.c
@@ -64,7 +64,7 @@ hypotl(long double x, long double y)
 	   if(ha >= ESW(MAX_EXP)) {	/* Inf or NaN */
 	       man_t manh, manl;
 	       /* Use original arg order iff result is NaN; quieten sNaNs. */
-	       w = fabsl(x+0.0)-fabsl(y+0.0);
+	       w = fabsl(x+0.0L)-fabsl(y+0);
 	       GET_LDBL_MAN(manh,manl,a);
 	       if (manh == LDBL_NBIT && manl == 0) w = a;
 	       GET_LDBL_MAN(manh,manl,b);
diff --git a/lib/msun/src/e_pow.c b/lib/msun/src/e_pow.c
index 51e913a6684a..69ddb7ffcb48 100644
--- a/lib/msun/src/e_pow.c
+++ b/lib/msun/src/e_pow.c
@@ -119,7 +119,7 @@ __ieee754_pow(double x, double y)
     /* y!=zero: result is NaN if either arg is NaN */
 	if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
 	   iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
-		return (x+0.0)+(y+0.0);
+	    return nan_mix(x, y);
 
     /* determine if y is an odd int when x < 0
      * yisint = 0	... y is not an integer
diff --git a/lib/msun/src/e_powf.c b/lib/msun/src/e_powf.c
index d7611e571731..53f1d37d6bec 100644
--- a/lib/msun/src/e_powf.c
+++ b/lib/msun/src/e_powf.c
@@ -76,7 +76,7 @@ __ieee754_powf(float x, float y)
     /* y!=zero: result is NaN if either arg is NaN */
 	if(ix > 0x7f800000 ||
 	   iy > 0x7f800000)
-		return (x+0.0F)+(y+0.0F);
+	    return nan_mix(x, y);
 
     /* determine if y is an odd int when x < 0
      * yisint = 0	... y is not an integer
diff --git a/lib/msun/src/e_remainder.c b/lib/msun/src/e_remainder.c
index 9be513b19022..c2ba63cc997f 100644
--- a/lib/msun/src/e_remainder.c
+++ b/lib/msun/src/e_remainder.c
@@ -45,11 +45,11 @@ __ieee754_remainder(double x, double p)
 	hx &= 0x7fffffff;
 
     /* purge off exception values */
-	if((hp|lp)==0) return (x*p)/(x*p); 	/* p = 0 */
-	if((hx>=0x7ff00000)||			/* x not finite */
+	if(((hp|lp)==0)||		 	/* p = 0 */
+	  (hx>=0x7ff00000)||			/* x not finite */
 	  ((hp>=0x7ff00000)&&			/* p is NaN */
 	  (((hp-0x7ff00000)|lp)!=0)))
-	    return ((long double)x*p)/((long double)x*p);
+	    return nan_mix(x, p)/nan_mix(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 b0014ae214f7..2b83fe0ec967 100644
--- a/lib/msun/src/e_remainderf.c
+++ b/lib/msun/src/e_remainderf.c
@@ -36,10 +36,10 @@ __ieee754_remainderf(float x, float p)
 	hx &= 0x7fffffff;
 
     /* purge off exception values */
-	if(hp==0) return (x*p)/(x*p);	 	/* p = 0 */
-	if((hx>=0x7f800000)||			/* x not finite */
+	if((hp==0)||			 	/* p = 0 */
+	  (hx>=0x7f800000)||			/* x not finite */
 	  ((hp>0x7f800000)))			/* p is NaN */
-	    return ((long double)x*p)/((long double)x*p);
+	    return nan_mix(x, p)/nan_mix(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 c87db67bf359..1b1105e1c3ac 100644
--- a/lib/msun/src/math_private.h
+++ b/lib/msun/src/math_private.h
@@ -478,6 +478,24 @@ 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
+ * 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.
+ *
+ * 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
+ * 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))
+
 #ifdef _COMPLEX_H
 
 /*
diff --git a/lib/msun/src/s_ccosh.c b/lib/msun/src/s_ccosh.c
index db6f108a429e..7652b3c49886 100644
--- a/lib/msun/src/s_ccosh.c
+++ b/lib/msun/src/s_ccosh.c
@@ -146,7 +146,8 @@ ccosh(double complex z)
 	 * Optionally raises the invalid floating-point exception for finite
 	 * nonzero y.  Choice = don't raise (except for signaling NaNs).
 	 */
-	return (CMPLX((x * x) * (y - y), (x + x) * (y - y)));
+	return (CMPLX(((long double)x * x) * (y - y),
+	    ((long double)x + x) * (y - y)));
 }
 
 double complex
diff --git a/lib/msun/src/s_ccoshf.c b/lib/msun/src/s_ccoshf.c
index 3e3d242e346a..5d7a09ba5f8d 100644
--- a/lib/msun/src/s_ccoshf.c
+++ b/lib/msun/src/s_ccoshf.c
@@ -92,7 +92,8 @@ ccoshf(float complex z)
 		return (CMPLXF(INFINITY * cosf(y), x * sinf(y)));
 	}
 
-	return (CMPLXF((x * x) * (y - y), (x + x) * (y - y)));
+	return (CMPLXF(((long double)x * x) * (y - y),
+	    ((long double)x + x) * (y - y)));
 }
 
 float complex
diff --git a/lib/msun/src/s_csinh.c b/lib/msun/src/s_csinh.c
index 1a32bf99b812..09b2dcbecb7b 100644
--- a/lib/msun/src/s_csinh.c
+++ b/lib/msun/src/s_csinh.c
@@ -145,7 +145,8 @@ csinh(double complex z)
 	 * Optionally raises the invalid floating-point exception for finite
 	 * nonzero y.  Choice = don't raise (except for signaling NaNs).
 	 */
-	return (CMPLX((x + x) * (y - y), (x * x) * (y - y)));
+	return (CMPLX(((long double)x + x) * (y - y),
+	    ((long double)x * x) * (y - y)));
 }
 
 double complex
diff --git a/lib/msun/src/s_csinhf.c b/lib/msun/src/s_csinhf.c
index 03c2082df207..ba1a79497fc6 100644
--- a/lib/msun/src/s_csinhf.c
+++ b/lib/msun/src/s_csinhf.c
@@ -92,7 +92,8 @@ csinhf(float complex z)
 		return (CMPLXF(x * cosf(y), INFINITY * sinf(y)));
 	}
 
-	return (CMPLXF((x + x) * (y - y), (x * x) * (y - y)));
+	return (CMPLXF(((long double)x + x) * (y - y),
+	    ((long double)x * x) * (y - y)));
 }
 
 float complex
diff --git a/lib/msun/src/s_csqrt.c b/lib/msun/src/s_csqrt.c
index 807532d2ae09..a2f3441ab7de 100644
--- a/lib/msun/src/s_csqrt.c
+++ b/lib/msun/src/s_csqrt.c
@@ -65,7 +65,7 @@ csqrt(double complex z)
 		return (CMPLX(INFINITY, b));
 	if (isnan(a)) {
 		t = (b - b) / (b - b);	/* raise invalid if b is not a NaN */
-		return (CMPLX(a, t));	/* return NaN + NaN i */
+		return (CMPLX(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */
 	}
 	if (isinf(a)) {
 		/*
@@ -79,10 +79,10 @@ csqrt(double complex z)
 		else
 			return (CMPLX(a, copysign(b - b, b)));
 	}
-	/*
-	 * The remaining special case (b is NaN) is handled just fine by
-	 * the normal code path below.
-	 */
+	if (isnan(b)) {
+		t = (a - a) / (a - a);	/* raise invalid */
+		return (CMPLX(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */
+	}
 
 	/* Scale to avoid overflow. */
 	if (fabs(a) >= THRESH || fabs(b) >= THRESH) {
diff --git a/lib/msun/src/s_csqrtf.c b/lib/msun/src/s_csqrtf.c
index e7a2f68b4283..0f36aa0ffb51 100644
--- a/lib/msun/src/s_csqrtf.c
+++ b/lib/msun/src/s_csqrtf.c
@@ -56,7 +56,7 @@ csqrtf(float complex z)
 		return (CMPLXF(INFINITY, b));
 	if (isnan(a)) {
 		t = (b - b) / (b - b);	/* raise invalid if b is not a NaN */
-		return (CMPLXF(a, t));	/* return NaN + NaN i */
+		return (CMPLXF(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */
 	}
 	if (isinf(a)) {
 		/*
@@ -70,10 +70,10 @@ csqrtf(float complex z)
 		else
 			return (CMPLXF(a, copysignf(b - b, b)));
 	}
-	/*
-	 * The remaining special case (b is NaN) is handled just fine by
-	 * the normal code path below.
-	 */
+	if (isnan(b)) {
+		t = (a - a) / (a - a);	/* raise invalid */
+		return (CMPLXF(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */
+	}
 
 	/*
 	 * We compute t in double precision to avoid overflow and to
diff --git a/lib/msun/src/s_csqrtl.c b/lib/msun/src/s_csqrtl.c
index d46fca7f4857..a34ad12235ca 100644
--- a/lib/msun/src/s_csqrtl.c
+++ b/lib/msun/src/s_csqrtl.c
@@ -73,7 +73,7 @@ csqrtl(long double complex z)
 		return (CMPLXL(INFINITY, b));
 	if (isnan(a)) {
 		t = (b - b) / (b - b);	/* raise invalid if b is not a NaN */
-		return (CMPLXL(a, t));	/* return NaN + NaN i */
+		return (CMPLXL(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */
 	}
 	if (isinf(a)) {
 		/*
@@ -87,10 +87,10 @@ csqrtl(long double complex z)
 		else
 			return (CMPLXL(a, copysignl(b - b, b)));
 	}
-	/*
-	 * The remaining special case (b is NaN) is handled just fine by
-	 * the normal code path below.
-	 */
+	if (isnan(b)) {
+		t = (a - a) / (a - a);	/* raise invalid */
+		return (CMPLXL(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */
+	}
 
 	/* Scale to avoid overflow. */
 	if (fabsl(a) >= THRESH || fabsl(b) >= THRESH) {
diff --git a/lib/msun/src/s_ctanh.c b/lib/msun/src/s_ctanh.c
index c8ecd5d713f7..88afeb50e26e 100644
--- a/lib/msun/src/s_ctanh.c
+++ b/lib/msun/src/s_ctanh.c
@@ -104,8 +104,8 @@ ctanh(double complex z)
 	 */
 	if (ix >= 0x7ff00000) {
 		if ((ix & 0xfffff) | lx)	/* x is NaN */
-			return (CMPLX((x + 0) * (y + 0),
-			    y == 0 ? y : (x + 0) * (y + 0)));
+			return (CMPLX(nan_mix(x, y),
+			    y == 0 ? y : nan_mix(x, y)));
 		SET_HIGH_WORD(x, hx - 0x40000000);	/* x = copysign(1, x) */
 		return (CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y))));
 	}
diff --git a/lib/msun/src/s_ctanhf.c b/lib/msun/src/s_ctanhf.c
index cc4ab41c861f..d2bd0b6786f0 100644
--- a/lib/msun/src/s_ctanhf.c
+++ b/lib/msun/src/s_ctanhf.c
@@ -53,8 +53,8 @@ ctanhf(float complex z)
 
 	if (ix >= 0x7f800000) {
 		if (ix & 0x7fffff)
-			return (CMPLXF((x + 0) * (y + 0),
-			    y == 0 ? y : (x + 0) * (y + 0)));
+			return (CMPLXF(nan_mix(x, y),
+			    y == 0 ? y : nan_mix(x, y)));
 		SET_FLOAT_WORD(x, hx - 0x40000000);
 		return (CMPLXF(x,
 		    copysignf(0, isinf(y) ? y : sinf(y) * cosf(y))));
diff --git a/lib/msun/src/s_remquo.c b/lib/msun/src/s_remquo.c
index d811c69da90c..9c61fa7487e5 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 (x*y)/(x*y);
+	    return nan_mix(x, y)/nan_mix(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 f7b4c0065ce8..cbe618f4647a 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 (x*y)/(x*y);
+	    return nan_mix(x, y)/nan_mix(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 712651c94ee6..2c7737949edc 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 (x*y)/(x*y);
+	    return nan_mix(x, y)/nan_mix(x, y);
 	if(ux.bits.exp<=uy.bits.exp) {
 	    if((ux.bits.exp<uy.bits.exp) ||
 	       (ux.bits.manh<=uy.bits.manh &&