From 6f1b8a07926a4a00060628d9ed70412b4842ed69 Mon Sep 17 00:00:00 2001 From: Bruce Evans Date: Tue, 17 Jul 2018 07:42:14 +0000 Subject: [PATCH] Add a macro nan_mix() and use it to get NaN results that are (bitwise) independent of the precision in most cases. This is mainly to simplify checking for errors. r176266 did this for e_pow[f].c using a less refined expression that often didn't work. r176276 fixes an error in the log message for r176266. The main refinement is to always expand to long double precision. See old log messages (especially these 2) and the comment on the macro for more general details. Specific details: - using nan_mix() consistently for the new and old pow*() functions was the only thing needed to make my consistency test for powl() vs pow() pass on amd64. - catrig[fl].c already had all the refinements, but open-coded. - e_atan2[fl].c, e_fmod[fl].c and s_remquo[fl] only had primitive NaN mixing. - e_hypot[fl].c already had a different refined version of r176266. Refine this further. nan_mix() is not directly usable here since we want to clear the sign bit. - e_remainder[f].c already had an earlier version of r176266. - s_ccosh[f].c,/s_csinh[f].c already had a version equivalent to r176266. Refine this further. nan_mix() is not directly usable here since the expression has to handle some non-NaN cases. - s_csqrt.[fl]: the mixing was special and mostly wrong. Partially fix the special version. - s_ctanh[f].c already had a version of r176266. --- lib/msun/ld128/e_powl.c | 2 +- lib/msun/ld80/e_powl.c | 4 ++-- lib/msun/src/catrig.c | 6 +++--- lib/msun/src/catrigf.c | 6 +++--- lib/msun/src/catrigl.c | 6 +++--- lib/msun/src/e_atan2.c | 2 +- lib/msun/src/e_atan2f.c | 2 +- lib/msun/src/e_atan2l.c | 2 +- lib/msun/src/e_fmod.c | 2 +- lib/msun/src/e_fmodf.c | 2 +- lib/msun/src/e_fmodl.c | 2 +- lib/msun/src/e_hypot.c | 2 +- lib/msun/src/e_hypotf.c | 2 +- lib/msun/src/e_hypotl.c | 2 +- lib/msun/src/e_pow.c | 2 +- lib/msun/src/e_powf.c | 2 +- lib/msun/src/e_remainder.c | 6 +++--- lib/msun/src/e_remainderf.c | 6 +++--- lib/msun/src/math_private.h | 18 ++++++++++++++++++ lib/msun/src/s_ccosh.c | 3 ++- lib/msun/src/s_ccoshf.c | 3 ++- lib/msun/src/s_csinh.c | 3 ++- lib/msun/src/s_csinhf.c | 3 ++- lib/msun/src/s_csqrt.c | 10 +++++----- lib/msun/src/s_csqrtf.c | 10 +++++----- lib/msun/src/s_csqrtl.c | 10 +++++----- lib/msun/src/s_ctanh.c | 4 ++-- lib/msun/src/s_ctanhf.c | 4 ++-- lib/msun/src/s_remquo.c | 2 +- lib/msun/src/s_remquof.c | 2 +- lib/msun/src/s_remquol.c | 2 +- 31 files changed, 77 insertions(+), 55 deletions(-) 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=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>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= 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=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