svn commit: r336362 - in head/lib/msun: ld128 ld80 src
Bruce Evans
bde at FreeBSD.org
Tue Jul 17 07:42:22 UTC 2018
Author: bde
Date: Tue Jul 17 07:42:14 2018
New Revision: 336362
URL: https://svnweb.freebsd.org/changeset/base/336362
Log:
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.
Modified:
head/lib/msun/ld128/e_powl.c
head/lib/msun/ld80/e_powl.c
head/lib/msun/src/catrig.c
head/lib/msun/src/catrigf.c
head/lib/msun/src/catrigl.c
head/lib/msun/src/e_atan2.c
head/lib/msun/src/e_atan2f.c
head/lib/msun/src/e_atan2l.c
head/lib/msun/src/e_fmod.c
head/lib/msun/src/e_fmodf.c
head/lib/msun/src/e_fmodl.c
head/lib/msun/src/e_hypot.c
head/lib/msun/src/e_hypotf.c
head/lib/msun/src/e_hypotl.c
head/lib/msun/src/e_pow.c
head/lib/msun/src/e_powf.c
head/lib/msun/src/e_remainder.c
head/lib/msun/src/e_remainderf.c
head/lib/msun/src/math_private.h
head/lib/msun/src/s_ccosh.c
head/lib/msun/src/s_ccoshf.c
head/lib/msun/src/s_csinh.c
head/lib/msun/src/s_csinhf.c
head/lib/msun/src/s_csqrt.c
head/lib/msun/src/s_csqrtf.c
head/lib/msun/src/s_csqrtl.c
head/lib/msun/src/s_ctanh.c
head/lib/msun/src/s_ctanhf.c
head/lib/msun/src/s_remquo.c
head/lib/msun/src/s_remquof.c
head/lib/msun/src/s_remquol.c
Modified: head/lib/msun/ld128/e_powl.c
==============================================================================
--- head/lib/msun/ld128/e_powl.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/ld128/e_powl.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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
Modified: head/lib/msun/ld80/e_powl.c
==============================================================================
--- head/lib/msun/ld80/e_powl.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/ld80/e_powl.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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 );
Modified: head/lib/msun/src/catrig.c
==============================================================================
--- head/lib/msun/src/catrig.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/catrig.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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)
Modified: head/lib/msun/src/catrigf.c
==============================================================================
--- head/lib/msun/src/catrigf.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/catrigf.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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)
Modified: head/lib/msun/src/catrigl.c
==============================================================================
--- head/lib/msun/src/catrigl.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/catrigl.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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)
Modified: head/lib/msun/src/e_atan2.c
==============================================================================
--- head/lib/msun/src/e_atan2.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/e_atan2.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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) */
Modified: head/lib/msun/src/e_atan2f.c
==============================================================================
--- head/lib/msun/src/e_atan2f.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/e_atan2f.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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) */
Modified: head/lib/msun/src/e_atan2l.c
==============================================================================
--- head/lib/msun/src/e_atan2l.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/e_atan2l.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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) */
Modified: head/lib/msun/src/e_fmod.c
==============================================================================
--- head/lib/msun/src/e_fmod.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/e_fmod.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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)
Modified: head/lib/msun/src/e_fmodf.c
==============================================================================
--- head/lib/msun/src/e_fmodf.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/e_fmodf.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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*/
Modified: head/lib/msun/src/e_fmodl.c
==============================================================================
--- head/lib/msun/src/e_fmodl.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/e_fmodl.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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 &&
Modified: head/lib/msun/src/e_hypot.c
==============================================================================
--- head/lib/msun/src/e_hypot.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/e_hypot.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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);
Modified: head/lib/msun/src/e_hypotf.c
==============================================================================
--- head/lib/msun/src/e_hypotf.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/e_hypotf.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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;
Modified: head/lib/msun/src/e_hypotl.c
==============================================================================
--- head/lib/msun/src/e_hypotl.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/e_hypotl.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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);
Modified: head/lib/msun/src/e_pow.c
==============================================================================
--- head/lib/msun/src/e_pow.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/e_pow.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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
Modified: head/lib/msun/src/e_powf.c
==============================================================================
--- head/lib/msun/src/e_powf.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/e_powf.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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
Modified: head/lib/msun/src/e_remainder.c
==============================================================================
--- head/lib/msun/src/e_remainder.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/e_remainder.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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 */
Modified: head/lib/msun/src/e_remainderf.c
==============================================================================
--- head/lib/msun/src/e_remainderf.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/e_remainderf.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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 */
Modified: head/lib/msun/src/math_private.h
==============================================================================
--- head/lib/msun/src/math_private.h Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/math_private.h Tue Jul 17 07:42:14 2018 (r336362)
@@ -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
/*
Modified: head/lib/msun/src/s_ccosh.c
==============================================================================
--- head/lib/msun/src/s_ccosh.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/s_ccosh.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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
Modified: head/lib/msun/src/s_ccoshf.c
==============================================================================
--- head/lib/msun/src/s_ccoshf.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/s_ccoshf.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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
Modified: head/lib/msun/src/s_csinh.c
==============================================================================
--- head/lib/msun/src/s_csinh.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/s_csinh.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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
Modified: head/lib/msun/src/s_csinhf.c
==============================================================================
--- head/lib/msun/src/s_csinhf.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/s_csinhf.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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
Modified: head/lib/msun/src/s_csqrt.c
==============================================================================
--- head/lib/msun/src/s_csqrt.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/s_csqrt.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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) {
Modified: head/lib/msun/src/s_csqrtf.c
==============================================================================
--- head/lib/msun/src/s_csqrtf.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/s_csqrtf.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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
Modified: head/lib/msun/src/s_csqrtl.c
==============================================================================
--- head/lib/msun/src/s_csqrtl.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/s_csqrtl.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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) {
Modified: head/lib/msun/src/s_ctanh.c
==============================================================================
--- head/lib/msun/src/s_ctanh.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/s_ctanh.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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))));
}
Modified: head/lib/msun/src/s_ctanhf.c
==============================================================================
--- head/lib/msun/src/s_ctanhf.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/s_ctanhf.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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))));
Modified: head/lib/msun/src/s_remquo.c
==============================================================================
--- head/lib/msun/src/s_remquo.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/s_remquo.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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;
Modified: head/lib/msun/src/s_remquof.c
==============================================================================
--- head/lib/msun/src/s_remquof.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/s_remquof.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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 */
Modified: head/lib/msun/src/s_remquol.c
==============================================================================
--- head/lib/msun/src/s_remquol.c Tue Jul 17 04:43:58 2018 (r336361)
+++ head/lib/msun/src/s_remquol.c Tue Jul 17 07:42:14 2018 (r336362)
@@ -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 &&
More information about the svn-src-head
mailing list