diff options
87 files changed, 339 insertions, 372 deletions
diff --git a/include/complex.h b/include/complex.h index b4bcb22e..8206e026 100644 --- a/include/complex.h +++ b/include/complex.h @@ -112,6 +112,15 @@ long double creall(long double complex); #define cimagf(x) __CIMAG(x, float) #define cimagl(x) __CIMAG(x, long double) +#define __CMPLX(x, y, t) \ + ((union { _Complex t __z; t __xy[2]; }){.__xy = {(x),(y)}}.__z) + +#if __STDC_VERSION__ >= 201112L +#define CMPLX(x, y) __CMPLX(x, y, double) +#define CMPLXF(x, y) __CMPLX(x, y, float) +#define CMPLXL(x, y) __CMPLX(x, y, long double) +#endif + #ifdef __cplusplus } #endif diff --git a/src/complex/__cexp.c b/src/complex/__cexp.c index f603e2be..05ac28c7 100644 --- a/src/complex/__cexp.c +++ b/src/complex/__cexp.c @@ -83,5 +83,5 @@ double complex __ldexp_cexp(double complex z, int expt) half_expt = expt - half_expt; INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0); - return cpack(cos(y) * exp_x * scale1 * scale2, sin(y) * exp_x * scale1 * scale2); + return CMPLX(cos(y) * exp_x * scale1 * scale2, sin(y) * exp_x * scale1 * scale2); } diff --git a/src/complex/__cexpf.c b/src/complex/__cexpf.c index 47168e8f..69b54045 100644 --- a/src/complex/__cexpf.c +++ b/src/complex/__cexpf.c @@ -63,6 +63,6 @@ float complex __ldexp_cexpf(float complex z, int expt) half_expt = expt - half_expt; SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23); - return cpackf(cosf(y) * exp_x * scale1 * scale2, + return CMPLXF(cosf(y) * exp_x * scale1 * scale2, sinf(y) * exp_x * scale1 * scale2); } diff --git a/src/complex/cacos.c b/src/complex/cacos.c index 3aca0519..27c35636 100644 --- a/src/complex/cacos.c +++ b/src/complex/cacos.c @@ -7,5 +7,5 @@ double complex cacos(double complex z) { z = casin(z); - return cpack(M_PI_2 - creal(z), -cimag(z)); + return CMPLX(M_PI_2 - creal(z), -cimag(z)); } diff --git a/src/complex/cacosf.c b/src/complex/cacosf.c index 563766e7..11852659 100644 --- a/src/complex/cacosf.c +++ b/src/complex/cacosf.c @@ -5,5 +5,5 @@ float complex cacosf(float complex z) { z = casinf(z); - return cpackf((float)M_PI_2 - crealf(z), -cimagf(z)); + return CMPLXF((float)M_PI_2 - crealf(z), -cimagf(z)); } diff --git a/src/complex/cacosh.c b/src/complex/cacosh.c index c2dfc1ba..8c68cb01 100644 --- a/src/complex/cacosh.c +++ b/src/complex/cacosh.c @@ -5,5 +5,5 @@ double complex cacosh(double complex z) { z = cacos(z); - return cpack(-cimag(z), creal(z)); + return CMPLX(-cimag(z), creal(z)); } diff --git a/src/complex/cacoshf.c b/src/complex/cacoshf.c index 37ff8800..ade01c09 100644 --- a/src/complex/cacoshf.c +++ b/src/complex/cacoshf.c @@ -3,5 +3,5 @@ float complex cacoshf(float complex z) { z = cacosf(z); - return cpackf(-cimagf(z), crealf(z)); + return CMPLXF(-cimagf(z), crealf(z)); } diff --git a/src/complex/cacoshl.c b/src/complex/cacoshl.c index 2a04e27b..65342557 100644 --- a/src/complex/cacoshl.c +++ b/src/complex/cacoshl.c @@ -9,6 +9,6 @@ long double complex cacoshl(long double complex z) long double complex cacoshl(long double complex z) { z = cacosl(z); - return cpackl(-cimagl(z), creall(z)); + return CMPLXL(-cimagl(z), creall(z)); } #endif diff --git a/src/complex/cacosl.c b/src/complex/cacosl.c index 5992e056..7fd4a2f6 100644 --- a/src/complex/cacosl.c +++ b/src/complex/cacosl.c @@ -11,6 +11,6 @@ long double complex cacosl(long double complex z) long double complex cacosl(long double complex z) { z = casinl(z); - return cpackl(PI_2 - creall(z), -cimagl(z)); + return CMPLXL(PI_2 - creall(z), -cimagl(z)); } #endif diff --git a/src/complex/casin.c b/src/complex/casin.c index 79aff278..dfdda988 100644 --- a/src/complex/casin.c +++ b/src/complex/casin.c @@ -11,6 +11,6 @@ double complex casin(double complex z) x = creal(z); y = cimag(z); - w = cpack(1.0 - (x - y)*(x + y), -2.0*x*y); - return clog(cpack(-y, x) + csqrt(w)); + w = CMPLX(1.0 - (x - y)*(x + y), -2.0*x*y); + return clog(CMPLX(-y, x) + csqrt(w)); } diff --git a/src/complex/casinf.c b/src/complex/casinf.c index cb9863f6..93f0e335 100644 --- a/src/complex/casinf.c +++ b/src/complex/casinf.c @@ -9,6 +9,6 @@ float complex casinf(float complex z) x = crealf(z); y = cimagf(z); - w = cpackf(1.0 - (x - y)*(x + y), -2.0*x*y); - return clogf(cpackf(-y, x) + csqrtf(w)); + w = CMPLXF(1.0 - (x - y)*(x + y), -2.0*x*y); + return clogf(CMPLXF(-y, x) + csqrtf(w)); } diff --git a/src/complex/casinh.c b/src/complex/casinh.c index f2b3fef8..b57fe8c4 100644 --- a/src/complex/casinh.c +++ b/src/complex/casinh.c @@ -4,6 +4,6 @@ double complex casinh(double complex z) { - z = casin(cpack(-cimag(z), creal(z))); - return cpack(cimag(z), -creal(z)); + z = casin(CMPLX(-cimag(z), creal(z))); + return CMPLX(cimag(z), -creal(z)); } diff --git a/src/complex/casinhf.c b/src/complex/casinhf.c index ed4af643..a11bf902 100644 --- a/src/complex/casinhf.c +++ b/src/complex/casinhf.c @@ -2,6 +2,6 @@ float complex casinhf(float complex z) { - z = casinf(cpackf(-cimagf(z), crealf(z))); - return cpackf(cimagf(z), -crealf(z)); + z = casinf(CMPLXF(-cimagf(z), crealf(z))); + return CMPLXF(cimagf(z), -crealf(z)); } diff --git a/src/complex/casinhl.c b/src/complex/casinhl.c index e5d80cef..714f1893 100644 --- a/src/complex/casinhl.c +++ b/src/complex/casinhl.c @@ -8,7 +8,7 @@ long double complex casinhl(long double complex z) #else long double complex casinhl(long double complex z) { - z = casinl(cpackl(-cimagl(z), creall(z))); - return cpackl(cimagl(z), -creall(z)); + z = casinl(CMPLXL(-cimagl(z), creall(z))); + return CMPLXL(cimagl(z), -creall(z)); } #endif diff --git a/src/complex/casinl.c b/src/complex/casinl.c index f9aa8ded..0916c60f 100644 --- a/src/complex/casinl.c +++ b/src/complex/casinl.c @@ -14,7 +14,7 @@ long double complex casinl(long double complex z) x = creall(z); y = cimagl(z); - w = cpackl(1.0 - (x - y)*(x + y), -2.0*x*y); - return clogl(cpackl(-y, x) + csqrtl(w)); + w = CMPLXL(1.0 - (x - y)*(x + y), -2.0*x*y); + return clogl(CMPLXL(-y, x) + csqrtl(w)); } #endif diff --git a/src/complex/catanh.c b/src/complex/catanh.c index b1628022..e248d9b9 100644 --- a/src/complex/catanh.c +++ b/src/complex/catanh.c @@ -4,6 +4,6 @@ double complex catanh(double complex z) { - z = catan(cpack(-cimag(z), creal(z))); - return cpack(cimag(z), -creal(z)); + z = catan(CMPLX(-cimag(z), creal(z))); + return CMPLX(cimag(z), -creal(z)); } diff --git a/src/complex/catanhf.c b/src/complex/catanhf.c index e1d1e648..4a5eb040 100644 --- a/src/complex/catanhf.c +++ b/src/complex/catanhf.c @@ -2,6 +2,6 @@ float complex catanhf(float complex z) { - z = catanf(cpackf(-cimagf(z), crealf(z))); - return cpackf(cimagf(z), -crealf(z)); + z = catanf(CMPLXF(-cimagf(z), crealf(z))); + return CMPLXF(cimagf(z), -crealf(z)); } diff --git a/src/complex/catanhl.c b/src/complex/catanhl.c index 0a9374a3..a5dd538e 100644 --- a/src/complex/catanhl.c +++ b/src/complex/catanhl.c @@ -8,7 +8,7 @@ long double complex catanhl(long double complex z) #else long double complex catanhl(long double complex z) { - z = catanl(cpackl(-cimagl(z), creall(z))); - return cpackl(cimagl(z), -creall(z)); + z = catanl(CMPLXL(-cimagl(z), creall(z))); + return CMPLXL(cimagl(z), -creall(z)); } #endif diff --git a/src/complex/ccos.c b/src/complex/ccos.c index 5754c238..645aec29 100644 --- a/src/complex/ccos.c +++ b/src/complex/ccos.c @@ -4,5 +4,5 @@ double complex ccos(double complex z) { - return ccosh(cpack(-cimag(z), creal(z))); + return ccosh(CMPLX(-cimag(z), creal(z))); } diff --git a/src/complex/ccosf.c b/src/complex/ccosf.c index 9b72c4f4..9a67241f 100644 --- a/src/complex/ccosf.c +++ b/src/complex/ccosf.c @@ -2,5 +2,5 @@ float complex ccosf(float complex z) { - return ccoshf(cpackf(-cimagf(z), crealf(z))); + return ccoshf(CMPLXF(-cimagf(z), crealf(z))); } diff --git a/src/complex/ccosh.c b/src/complex/ccosh.c index 81f29432..401f3c60 100644 --- a/src/complex/ccosh.c +++ b/src/complex/ccosh.c @@ -55,23 +55,23 @@ double complex ccosh(double complex z) /* Handle the nearly-non-exceptional cases where x and y are finite. */ if (ix < 0x7ff00000 && iy < 0x7ff00000) { if ((iy | ly) == 0) - return cpack(cosh(x), x * y); + return CMPLX(cosh(x), x * y); if (ix < 0x40360000) /* small x: normal case */ - return cpack(cosh(x) * cos(y), sinh(x) * sin(y)); + return CMPLX(cosh(x) * cos(y), sinh(x) * sin(y)); /* |x| >= 22, so cosh(x) ~= exp(|x|) */ if (ix < 0x40862e42) { /* x < 710: exp(|x|) won't overflow */ h = exp(fabs(x)) * 0.5; - return cpack(h * cos(y), copysign(h, x) * sin(y)); + return CMPLX(h * cos(y), copysign(h, x) * sin(y)); } else if (ix < 0x4096bbaa) { /* x < 1455: scale to avoid overflow */ - z = __ldexp_cexp(cpack(fabs(x), y), -1); - return cpack(creal(z), cimag(z) * copysign(1, x)); + z = __ldexp_cexp(CMPLX(fabs(x), y), -1); + return CMPLX(creal(z), cimag(z) * copysign(1, x)); } else { /* x >= 1455: the result always overflows */ h = huge * x; - return cpack(h * h * cos(y), h * sin(y)); + return CMPLX(h * h * cos(y), h * sin(y)); } } @@ -85,7 +85,7 @@ double complex ccosh(double complex z) * the same as d(NaN). */ if ((ix | lx) == 0 && iy >= 0x7ff00000) - return cpack(y - y, copysign(0, x * (y - y))); + return CMPLX(y - y, copysign(0, x * (y - y))); /* * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0. @@ -95,8 +95,8 @@ double complex ccosh(double complex z) */ if ((iy | ly) == 0 && ix >= 0x7ff00000) { if (((hx & 0xfffff) | lx) == 0) - return cpack(x * x, copysign(0, x) * y); - return cpack(x * x, copysign(0, (x + x) * y)); + return CMPLX(x * x, copysign(0, x) * y); + return CMPLX(x * x, copysign(0, (x + x) * y)); } /* @@ -108,7 +108,7 @@ double complex ccosh(double complex z) * nonzero x. Choice = don't raise (except for signaling NaNs). */ if (ix < 0x7ff00000 && iy >= 0x7ff00000) - return cpack(y - y, x * (y - y)); + return CMPLX(y - y, x * (y - y)); /* * cosh(+-Inf + I NaN) = +Inf + I d(NaN). @@ -121,8 +121,8 @@ double complex ccosh(double complex z) */ if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) { if (iy >= 0x7ff00000) - return cpack(x * x, x * (y - y)); - return cpack((x * x) * cos(y), x * sin(y)); + return CMPLX(x * x, x * (y - y)); + return CMPLX((x * x) * cos(y), x * sin(y)); } /* @@ -136,5 +136,5 @@ double complex ccosh(double complex z) * Optionally raises the invalid floating-point exception for finite * nonzero y. Choice = don't raise (except for signaling NaNs). */ - return cpack((x * x) * (y - y), (x + x) * (y - y)); + return CMPLX((x * x) * (y - y), (x + x) * (y - y)); } diff --git a/src/complex/ccoshf.c b/src/complex/ccoshf.c index 683e77fa..90acfe05 100644 --- a/src/complex/ccoshf.c +++ b/src/complex/ccoshf.c @@ -48,43 +48,43 @@ float complex ccoshf(float complex z) if (ix < 0x7f800000 && iy < 0x7f800000) { if (iy == 0) - return cpackf(coshf(x), x * y); + return CMPLXF(coshf(x), x * y); if (ix < 0x41100000) /* small x: normal case */ - return cpackf(coshf(x) * cosf(y), sinhf(x) * sinf(y)); + return CMPLXF(coshf(x) * cosf(y), sinhf(x) * sinf(y)); /* |x| >= 9, so cosh(x) ~= exp(|x|) */ if (ix < 0x42b17218) { /* x < 88.7: expf(|x|) won't overflow */ h = expf(fabsf(x)) * 0.5f; - return cpackf(h * cosf(y), copysignf(h, x) * sinf(y)); + return CMPLXF(h * cosf(y), copysignf(h, x) * sinf(y)); } else if (ix < 0x4340b1e7) { /* x < 192.7: scale to avoid overflow */ - z = __ldexp_cexpf(cpackf(fabsf(x), y), -1); - return cpackf(crealf(z), cimagf(z) * copysignf(1, x)); + z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1); + return CMPLXF(crealf(z), cimagf(z) * copysignf(1, x)); } else { /* x >= 192.7: the result always overflows */ h = huge * x; - return cpackf(h * h * cosf(y), h * sinf(y)); + return CMPLXF(h * h * cosf(y), h * sinf(y)); } } if (ix == 0 && iy >= 0x7f800000) - return cpackf(y - y, copysignf(0, x * (y - y))); + return CMPLXF(y - y, copysignf(0, x * (y - y))); if (iy == 0 && ix >= 0x7f800000) { if ((hx & 0x7fffff) == 0) - return cpackf(x * x, copysignf(0, x) * y); - return cpackf(x * x, copysignf(0, (x + x) * y)); + return CMPLXF(x * x, copysignf(0, x) * y); + return CMPLXF(x * x, copysignf(0, (x + x) * y)); } if (ix < 0x7f800000 && iy >= 0x7f800000) - return cpackf(y - y, x * (y - y)); + return CMPLXF(y - y, x * (y - y)); if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) { if (iy >= 0x7f800000) - return cpackf(x * x, x * (y - y)); - return cpackf((x * x) * cosf(y), x * sinf(y)); + return CMPLXF(x * x, x * (y - y)); + return CMPLXF((x * x) * cosf(y), x * sinf(y)); } - return cpackf((x * x) * (y - y), (x + x) * (y - y)); + return CMPLXF((x * x) * (y - y), (x + x) * (y - y)); } diff --git a/src/complex/ccosl.c b/src/complex/ccosl.c index e37825a9..d787047f 100644 --- a/src/complex/ccosl.c +++ b/src/complex/ccosl.c @@ -8,6 +8,6 @@ long double complex ccosl(long double complex z) #else long double complex ccosl(long double complex z) { - return ccoshl(cpackl(-cimagl(z), creall(z))); + return ccoshl(CMPLXL(-cimagl(z), creall(z))); } #endif diff --git a/src/complex/cexp.c b/src/complex/cexp.c index 3b8bb752..5118e00e 100644 --- a/src/complex/cexp.c +++ b/src/complex/cexp.c @@ -44,22 +44,22 @@ double complex cexp(double complex z) /* cexp(x + I 0) = exp(x) + I 0 */ if ((hy | ly) == 0) - return cpack(exp(x), y); + return CMPLX(exp(x), y); EXTRACT_WORDS(hx, lx, x); /* cexp(0 + I y) = cos(y) + I sin(y) */ if (((hx & 0x7fffffff) | lx) == 0) - return cpack(cos(y), sin(y)); + return CMPLX(cos(y), sin(y)); if (hy >= 0x7ff00000) { if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) { /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ - return cpack(y - y, y - y); + return CMPLX(y - y, y - y); } else if (hx & 0x80000000) { /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ - return cpack(0.0, 0.0); + return CMPLX(0.0, 0.0); } else { /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ - return cpack(x, y - y); + return CMPLX(x, y - y); } } @@ -78,6 +78,6 @@ double complex cexp(double complex z) * - x = NaN (spurious inexact exception from y) */ exp_x = exp(x); - return cpack(exp_x * cos(y), exp_x * sin(y)); + return CMPLX(exp_x * cos(y), exp_x * sin(y)); } } diff --git a/src/complex/cexpf.c b/src/complex/cexpf.c index 0cf13a3d..1a09964c 100644 --- a/src/complex/cexpf.c +++ b/src/complex/cexpf.c @@ -44,22 +44,22 @@ float complex cexpf(float complex z) /* cexp(x + I 0) = exp(x) + I 0 */ if (hy == 0) - return cpackf(expf(x), y); + return CMPLXF(expf(x), y); GET_FLOAT_WORD(hx, x); /* cexp(0 + I y) = cos(y) + I sin(y) */ if ((hx & 0x7fffffff) == 0) - return cpackf(cosf(y), sinf(y)); + return CMPLXF(cosf(y), sinf(y)); if (hy >= 0x7f800000) { if ((hx & 0x7fffffff) != 0x7f800000) { /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ - return cpackf(y - y, y - y); + return CMPLXF(y - y, y - y); } else if (hx & 0x80000000) { /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ - return cpackf(0.0, 0.0); + return CMPLXF(0.0, 0.0); } else { /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ - return cpackf(x, y - y); + return CMPLXF(x, y - y); } } @@ -78,6 +78,6 @@ float complex cexpf(float complex z) * - x = NaN (spurious inexact exception from y) */ exp_x = expf(x); - return cpackf(exp_x * cosf(y), exp_x * sinf(y)); + return CMPLXF(exp_x * cosf(y), exp_x * sinf(y)); } } diff --git a/src/complex/clog.c b/src/complex/clog.c index 6f10a115..12aae9c7 100644 --- a/src/complex/clog.c +++ b/src/complex/clog.c @@ -10,5 +10,5 @@ double complex clog(double complex z) r = cabs(z); phi = carg(z); - return cpack(log(r), phi); + return CMPLX(log(r), phi); } diff --git a/src/complex/clogf.c b/src/complex/clogf.c index f3aec54d..e9b32e60 100644 --- a/src/complex/clogf.c +++ b/src/complex/clogf.c @@ -8,5 +8,5 @@ float complex clogf(float complex z) r = cabsf(z); phi = cargf(z); - return cpackf(logf(r), phi); + return CMPLXF(logf(r), phi); } diff --git a/src/complex/clogl.c b/src/complex/clogl.c index 5b84ba59..18f16088 100644 --- a/src/complex/clogl.c +++ b/src/complex/clogl.c @@ -13,6 +13,6 @@ long double complex clogl(long double complex z) r = cabsl(z); phi = cargl(z); - return cpackl(logl(r), phi); + return CMPLXL(logl(r), phi); } #endif diff --git a/src/complex/conj.c b/src/complex/conj.c index 4aceea7b..0b3f5f46 100644 --- a/src/complex/conj.c +++ b/src/complex/conj.c @@ -2,5 +2,5 @@ double complex conj(double complex z) { - return cpack(creal(z), -cimag(z)); + return CMPLX(creal(z), -cimag(z)); } diff --git a/src/complex/conjf.c b/src/complex/conjf.c index 31556800..9af6b2c3 100644 --- a/src/complex/conjf.c +++ b/src/complex/conjf.c @@ -2,5 +2,5 @@ float complex conjf(float complex z) { - return cpackf(crealf(z), -cimagf(z)); + return CMPLXF(crealf(z), -cimagf(z)); } diff --git a/src/complex/conjl.c b/src/complex/conjl.c index 01332262..67f11b9d 100644 --- a/src/complex/conjl.c +++ b/src/complex/conjl.c @@ -2,5 +2,5 @@ long double complex conjl(long double complex z) { - return cpackl(creall(z), -cimagl(z)); + return CMPLXL(creall(z), -cimagl(z)); } diff --git a/src/complex/cproj.c b/src/complex/cproj.c index 1cf9bb94..15f358a1 100644 --- a/src/complex/cproj.c +++ b/src/complex/cproj.c @@ -3,6 +3,6 @@ double complex cproj(double complex z) { if (isinf(creal(z)) || isinf(cimag(z))) - return cpack(INFINITY, copysign(0.0, creal(z))); + return CMPLX(INFINITY, copysign(0.0, creal(z))); return z; } diff --git a/src/complex/cprojf.c b/src/complex/cprojf.c index 71129743..653be5e8 100644 --- a/src/complex/cprojf.c +++ b/src/complex/cprojf.c @@ -3,6 +3,6 @@ float complex cprojf(float complex z) { if (isinf(crealf(z)) || isinf(cimagf(z))) - return cpackf(INFINITY, copysignf(0.0, crealf(z))); + return CMPLXF(INFINITY, copysignf(0.0, crealf(z))); return z; } diff --git a/src/complex/cprojl.c b/src/complex/cprojl.c index 72e50cf5..6731aaa2 100644 --- a/src/complex/cprojl.c +++ b/src/complex/cprojl.c @@ -9,7 +9,7 @@ long double complex cprojl(long double complex z) long double complex cprojl(long double complex z) { if (isinf(creall(z)) || isinf(cimagl(z))) - return cpackl(INFINITY, copysignl(0.0, creall(z))); + return CMPLXL(INFINITY, copysignl(0.0, creall(z))); return z; } #endif diff --git a/src/complex/csin.c b/src/complex/csin.c index 246a4595..ad8ae67a 100644 --- a/src/complex/csin.c +++ b/src/complex/csin.c @@ -4,6 +4,6 @@ double complex csin(double complex z) { - z = csinh(cpack(-cimag(z), creal(z))); - return cpack(cimag(z), -creal(z)); + z = csinh(CMPLX(-cimag(z), creal(z))); + return CMPLX(cimag(z), -creal(z)); } diff --git a/src/complex/csinf.c b/src/complex/csinf.c index 3aabe8f8..60b3cbaa 100644 --- a/src/complex/csinf.c +++ b/src/complex/csinf.c @@ -2,6 +2,6 @@ float complex csinf(float complex z) { - z = csinhf(cpackf(-cimagf(z), crealf(z))); - return cpackf(cimagf(z), -crealf(z)); + z = csinhf(CMPLXF(-cimagf(z), crealf(z))); + return CMPLXF(cimagf(z), -crealf(z)); } diff --git a/src/complex/csinh.c b/src/complex/csinh.c index fe16f06b..0f8035d1 100644 --- a/src/complex/csinh.c +++ b/src/complex/csinh.c @@ -55,23 +55,23 @@ double complex csinh(double complex z) /* Handle the nearly-non-exceptional cases where x and y are finite. */ if (ix < 0x7ff00000 && iy < 0x7ff00000) { if ((iy | ly) == 0) - return cpack(sinh(x), y); + return CMPLX(sinh(x), y); if (ix < 0x40360000) /* small x: normal case */ - return cpack(sinh(x) * cos(y), cosh(x) * sin(y)); + return CMPLX(sinh(x) * cos(y), cosh(x) * sin(y)); /* |x| >= 22, so cosh(x) ~= exp(|x|) */ if (ix < 0x40862e42) { /* x < 710: exp(|x|) won't overflow */ h = exp(fabs(x)) * 0.5; - return cpack(copysign(h, x) * cos(y), h * sin(y)); + return CMPLX(copysign(h, x) * cos(y), h * sin(y)); } else if (ix < 0x4096bbaa) { /* x < 1455: scale to avoid overflow */ - z = __ldexp_cexp(cpack(fabs(x), y), -1); - return cpack(creal(z) * copysign(1, x), cimag(z)); + z = __ldexp_cexp(CMPLX(fabs(x), y), -1); + return CMPLX(creal(z) * copysign(1, x), cimag(z)); } else { /* x >= 1455: the result always overflows */ h = huge * x; - return cpack(h * cos(y), h * h * sin(y)); + return CMPLX(h * cos(y), h * h * sin(y)); } } @@ -85,7 +85,7 @@ double complex csinh(double complex z) * the same as d(NaN). */ if ((ix | lx) == 0 && iy >= 0x7ff00000) - return cpack(copysign(0, x * (y - y)), y - y); + return CMPLX(copysign(0, x * (y - y)), y - y); /* * sinh(+-Inf +- I 0) = +-Inf + I +-0. @@ -94,8 +94,8 @@ double complex csinh(double complex z) */ if ((iy | ly) == 0 && ix >= 0x7ff00000) { if (((hx & 0xfffff) | lx) == 0) - return cpack(x, y); - return cpack(x, copysign(0, y)); + return CMPLX(x, y); + return CMPLX(x, copysign(0, y)); } /* @@ -107,7 +107,7 @@ double complex csinh(double complex z) * nonzero x. Choice = don't raise (except for signaling NaNs). */ if (ix < 0x7ff00000 && iy >= 0x7ff00000) - return cpack(y - y, x * (y - y)); + return CMPLX(y - y, x * (y - y)); /* * sinh(+-Inf + I NaN) = +-Inf + I d(NaN). @@ -122,8 +122,8 @@ double complex csinh(double complex z) */ if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) { if (iy >= 0x7ff00000) - return cpack(x * x, x * (y - y)); - return cpack(x * cos(y), INFINITY * sin(y)); + return CMPLX(x * x, x * (y - y)); + return CMPLX(x * cos(y), INFINITY * sin(y)); } /* @@ -137,5 +137,5 @@ double complex csinh(double complex z) * Optionally raises the invalid floating-point exception for finite * nonzero y. Choice = don't raise (except for signaling NaNs). */ - return cpack((x * x) * (y - y), (x + x) * (y - y)); + return CMPLX((x * x) * (y - y), (x + x) * (y - y)); } diff --git a/src/complex/csinhf.c b/src/complex/csinhf.c index bbb116c2..49697f02 100644 --- a/src/complex/csinhf.c +++ b/src/complex/csinhf.c @@ -48,43 +48,43 @@ float complex csinhf(float complex z) if (ix < 0x7f800000 && iy < 0x7f800000) { if (iy == 0) - return cpackf(sinhf(x), y); + return CMPLXF(sinhf(x), y); if (ix < 0x41100000) /* small x: normal case */ - return cpackf(sinhf(x) * cosf(y), coshf(x) * sinf(y)); + return CMPLXF(sinhf(x) * cosf(y), coshf(x) * sinf(y)); /* |x| >= 9, so cosh(x) ~= exp(|x|) */ if (ix < 0x42b17218) { /* x < 88.7: expf(|x|) won't overflow */ h = expf(fabsf(x)) * 0.5f; - return cpackf(copysignf(h, x) * cosf(y), h * sinf(y)); + return CMPLXF(copysignf(h, x) * cosf(y), h * sinf(y)); } else if (ix < 0x4340b1e7) { /* x < 192.7: scale to avoid overflow */ - z = __ldexp_cexpf(cpackf(fabsf(x), y), -1); - return cpackf(crealf(z) * copysignf(1, x), cimagf(z)); + z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1); + return CMPLXF(crealf(z) * copysignf(1, x), cimagf(z)); } else { /* x >= 192.7: the result always overflows */ h = huge * x; - return cpackf(h * cosf(y), h * h * sinf(y)); + return CMPLXF(h * cosf(y), h * h * sinf(y)); } } if (ix == 0 && iy >= 0x7f800000) - return cpackf(copysignf(0, x * (y - y)), y - y); + return CMPLXF(copysignf(0, x * (y - y)), y - y); if (iy == 0 && ix >= 0x7f800000) { if ((hx & 0x7fffff) == 0) - return cpackf(x, y); - return cpackf(x, copysignf(0, y)); + return CMPLXF(x, y); + return CMPLXF(x, copysignf(0, y)); } if (ix < 0x7f800000 && iy >= 0x7f800000) - return cpackf(y - y, x * (y - y)); + return CMPLXF(y - y, x * (y - y)); if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) { if (iy >= 0x7f800000) - return cpackf(x * x, x * (y - y)); - return cpackf(x * cosf(y), INFINITY * sinf(y)); + return CMPLXF(x * x, x * (y - y)); + return CMPLXF(x * cosf(y), INFINITY * sinf(y)); } - return cpackf((x * x) * (y - y), (x + x) * (y - y)); + return CMPLXF((x * x) * (y - y), (x + x) * (y - y)); } diff --git a/src/complex/csinl.c b/src/complex/csinl.c index 4ad86745..4e9f86c3 100644 --- a/src/complex/csinl.c +++ b/src/complex/csinl.c @@ -8,7 +8,7 @@ long double complex csinl(long double complex z) #else long double complex csinl(long double complex z) { - z = csinhl(cpackl(-cimagl(z), creall(z))); - return cpackl(cimagl(z), -creall(z)); + z = csinhl(CMPLXL(-cimagl(z), creall(z))); + return CMPLXL(cimagl(z), -creall(z)); } #endif diff --git a/src/complex/csqrt.c b/src/complex/csqrt.c index 21fb879d..8a2ba608 100644 --- a/src/complex/csqrt.c +++ b/src/complex/csqrt.c @@ -51,12 +51,12 @@ double complex csqrt(double complex z) /* Handle special cases. */ if (z == 0) - return cpack(0, b); + return CMPLX(0, b); if (isinf(b)) - return cpack(INFINITY, b); + return CMPLX(INFINITY, b); if (isnan(a)) { t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ - return cpack(a, t); /* return NaN + NaN i */ + return CMPLX(a, t); /* return NaN + NaN i */ } if (isinf(a)) { /* @@ -66,9 +66,9 @@ double complex csqrt(double complex z) * csqrt(-inf + y i) = 0 + inf i */ if (signbit(a)) - return cpack(fabs(b - b), copysign(a, b)); + return CMPLX(fabs(b - b), copysign(a, b)); else - return cpack(a, copysign(b - b, b)); + return CMPLX(a, copysign(b - b, b)); } /* * The remaining special case (b is NaN) is handled just fine by @@ -87,10 +87,10 @@ double complex csqrt(double complex z) /* Algorithm 312, CACM vol 10, Oct 1967. */ if (a >= 0) { t = sqrt((a + hypot(a, b)) * 0.5); - result = cpack(t, b / (2 * t)); + result = CMPLX(t, b / (2 * t)); } else { t = sqrt((-a + hypot(a, b)) * 0.5); - result = cpack(fabs(b) / (2 * t), copysign(t, b)); + result = CMPLX(fabs(b) / (2 * t), copysign(t, b)); } /* Rescale. */ diff --git a/src/complex/csqrtf.c b/src/complex/csqrtf.c index 16487c23..ab5102f0 100644 --- a/src/complex/csqrtf.c +++ b/src/complex/csqrtf.c @@ -43,12 +43,12 @@ float complex csqrtf(float complex z) /* Handle special cases. */ if (z == 0) - return cpackf(0, b); + return CMPLXF(0, b); if (isinf(b)) - return cpackf(INFINITY, b); + return CMPLXF(INFINITY, b); if (isnan(a)) { t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ - return cpackf(a, t); /* return NaN + NaN i */ + return CMPLXF(a, t); /* return NaN + NaN i */ } if (isinf(a)) { /* @@ -58,9 +58,9 @@ float complex csqrtf(float complex z) * csqrtf(-inf + y i) = 0 + inf i */ if (signbit(a)) - return cpackf(fabsf(b - b), copysignf(a, b)); + return CMPLXF(fabsf(b - b), copysignf(a, b)); else - return cpackf(a, copysignf(b - b, b)); + return CMPLXF(a, copysignf(b - b, b)); } /* * The remaining special case (b is NaN) is handled just fine by @@ -74,9 +74,9 @@ float complex csqrtf(float complex z) */ if (a >= 0) { t = sqrt((a + hypot(a, b)) * 0.5); - return cpackf(t, b / (2.0 * t)); + return CMPLXF(t, b / (2.0 * t)); } else { t = sqrt((-a + hypot(a, b)) * 0.5); - return cpackf(fabsf(b) / (2.0 * t), copysignf(t, b)); + return CMPLXF(fabsf(b) / (2.0 * t), copysignf(t, b)); } } diff --git a/src/complex/ctan.c b/src/complex/ctan.c index 4741a4df..c0926374 100644 --- a/src/complex/ctan.c +++ b/src/complex/ctan.c @@ -4,6 +4,6 @@ double complex ctan(double complex z) { - z = ctanh(cpack(-cimag(z), creal(z))); - return cpack(cimag(z), -creal(z)); + z = ctanh(CMPLX(-cimag(z), creal(z))); + return CMPLX(cimag(z), -creal(z)); } diff --git a/src/complex/ctanf.c b/src/complex/ctanf.c index 9bbeb051..009b1921 100644 --- a/src/complex/ctanf.c +++ b/src/complex/ctanf.c @@ -2,6 +2,6 @@ float complex ctanf(float complex z) { - z = ctanhf(cpackf(-cimagf(z), crealf(z))); - return cpackf(cimagf(z), -crealf(z)); + z = ctanhf(CMPLXF(-cimagf(z), crealf(z))); + return CMPLXF(cimagf(z), -crealf(z)); } diff --git a/src/complex/ctanh.c b/src/complex/ctanh.c index dd569fc3..0461050d 100644 --- a/src/complex/ctanh.c +++ b/src/complex/ctanh.c @@ -95,9 +95,9 @@ double complex ctanh(double complex z) */ if (ix >= 0x7ff00000) { if ((ix & 0xfffff) | lx) /* x is NaN */ - return cpack(x, (y == 0 ? y : x * y)); + return CMPLX(x, (y == 0 ? y : x * y)); SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */ - return cpack(x, copysign(0, isinf(y) ? y : sin(y) * cos(y))); + return CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y))); } /* @@ -105,7 +105,7 @@ double complex ctanh(double complex z) * ctanh(x +- i Inf) = NaN + i NaN */ if (!isfinite(y)) - return cpack(y - y, y - y); + return CMPLX(y - y, y - y); /* * ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the @@ -114,7 +114,7 @@ double complex ctanh(double complex z) */ if (ix >= 0x40360000) { /* x >= 22 */ double exp_mx = exp(-fabs(x)); - return cpack(copysign(1, x), 4 * sin(y) * cos(y) * exp_mx * exp_mx); + return CMPLX(copysign(1, x), 4 * sin(y) * cos(y) * exp_mx * exp_mx); } /* Kahan's algorithm */ @@ -123,5 +123,5 @@ double complex ctanh(double complex z) s = sinh(x); rho = sqrt(1 + s * s); /* = cosh(x) */ denom = 1 + beta * s * s; - return cpack((beta * rho * s) / denom, t / denom); + return CMPLX((beta * rho * s) / denom, t / denom); } diff --git a/src/complex/ctanhf.c b/src/complex/ctanhf.c index 7d746134..a7e1a5fc 100644 --- a/src/complex/ctanhf.c +++ b/src/complex/ctanhf.c @@ -44,17 +44,17 @@ float complex ctanhf(float complex z) if (ix >= 0x7f800000) { if (ix & 0x7fffff) - return cpackf(x, (y == 0 ? y : x * y)); + return CMPLXF(x, (y == 0 ? y : x * y)); SET_FLOAT_WORD(x, hx - 0x40000000); - return cpackf(x, copysignf(0, isinf(y) ? y : sinf(y) * cosf(y))); + return CMPLXF(x, copysignf(0, isinf(y) ? y : sinf(y) * cosf(y))); } if (!isfinite(y)) - return cpackf(y - y, y - y); + return CMPLXF(y - y, y - y); if (ix >= 0x41300000) { /* x >= 11 */ float exp_mx = expf(-fabsf(x)); - return cpackf(copysignf(1, x), 4 * sinf(y) * cosf(y) * exp_mx * exp_mx); + return CMPLXF(copysignf(1, x), 4 * sinf(y) * cosf(y) * exp_mx * exp_mx); } t = tanf(y); @@ -62,5 +62,5 @@ float complex ctanhf(float complex z) s = sinhf(x); rho = sqrtf(1 + s * s); denom = 1 + beta * s * s; - return cpackf((beta * rho * s) / denom, t / denom); + return CMPLXF((beta * rho * s) / denom, t / denom); } diff --git a/src/complex/ctanl.c b/src/complex/ctanl.c index 4b4c99b6..ac1c3e0a 100644 --- a/src/complex/ctanl.c +++ b/src/complex/ctanl.c @@ -8,7 +8,7 @@ long double complex ctanl(long double complex z) #else long double complex ctanl(long double complex z) { - z = ctanhl(cpackl(-cimagl(z), creall(z))); - return cpackl(cimagl(z), -creall(z)); + z = ctanhl(CMPLXL(-cimagl(z), creall(z))); + return CMPLXL(cimagl(z), -creall(z)); } #endif diff --git a/src/fenv/fenv.c b/src/fenv/fenv.c index f77599bc..b2dfd254 100644 --- a/src/fenv/fenv.c +++ b/src/fenv/fenv.c @@ -19,7 +19,7 @@ int fetestexcept(int mask) int fegetround(void) { - return 0; + return FE_TONEAREST; } int fesetround(int r) diff --git a/src/fenv/feupdateenv.c b/src/fenv/feupdateenv.c index 50cef8e5..f45ed7c0 100644 --- a/src/fenv/feupdateenv.c +++ b/src/fenv/feupdateenv.c @@ -2,6 +2,7 @@ int feupdateenv(const fenv_t *envp) { + #pragma STDC FENV_ACCESS ON int ex = fetestexcept(FE_ALL_EXCEPT); fesetenv(envp); feraiseexcept(ex); diff --git a/src/internal/libm.h b/src/internal/libm.h index a71c4c05..46c4b564 100644 --- a/src/internal/libm.h +++ b/src/internal/libm.h @@ -157,38 +157,23 @@ long double __tanl(long double, long double, int); long double __polevll(long double, const long double *, int); long double __p1evll(long double, const long double *, int); -// FIXME: not needed when -fexcess-precision=standard is supported (>=gcc4.5) -/* - * Attempt to get strict C99 semantics for assignment with non-C99 compilers. - */ -#if 1 +#if 0 +/* Attempt to get strict C99 semantics for assignment with non-C99 compilers. */ #define STRICT_ASSIGN(type, lval, rval) do { \ volatile type __v = (rval); \ (lval) = __v; \ } while (0) #else +/* Should work with -fexcess-precision=standard (>=gcc-4.5) or -ffloat-store */ #define STRICT_ASSIGN(type, lval, rval) ((lval) = (type)(rval)) #endif - /* complex */ -union dcomplex { - double complex z; - double a[2]; -}; -union fcomplex { - float complex z; - float a[2]; -}; -union lcomplex { - long double complex z; - long double a[2]; -}; - -/* x + y*I is not supported properly by gcc */ -#define cpack(x,y) ((union dcomplex){.a={(x),(y)}}.z) -#define cpackf(x,y) ((union fcomplex){.a={(x),(y)}}.z) -#define cpackl(x,y) ((union lcomplex){.a={(x),(y)}}.z) +#ifndef CMPLX +#define CMPLX(x, y) __CMPLX(x, y, double) +#define CMPLXF(x, y) __CMPLX(x, y, float) +#define CMPLXL(x, y) __CMPLX(x, y, long double) +#endif #endif diff --git a/src/math/__invtrigl.c b/src/math/__invtrigl.c index a821842c..a51330e2 100644 --- a/src/math/__invtrigl.c +++ b/src/math/__invtrigl.c @@ -28,9 +28,8 @@ #include "__invtrigl.h" #if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 -/* - * asinl() and acosl() - */ + +/* coefficients used in asinl() and acosl() */ const long double pS0 = 1.66666666666666666631e-01L, pS1 = -4.16313987993683104320e-01L, @@ -45,38 +44,9 @@ qS3 = -1.68285799854822427013e+00L, qS4 = 3.90699412641738801874e-01L, qS5 = -3.14365703596053263322e-02L; -/* - * atanl() - */ -const long double atanhi[] = { - 4.63647609000806116202e-01L, - 7.85398163397448309628e-01L, - 9.82793723247329067960e-01L, - 1.57079632679489661926e+00L, -}; - -const long double atanlo[] = { - 1.18469937025062860669e-20L, - -1.25413940316708300586e-20L, - 2.55232234165405176172e-20L, - -2.50827880633416601173e-20L, -}; - -const long double aT[] = { - 3.33333333333333333017e-01L, - -1.99999999999999632011e-01L, - 1.42857142857046531280e-01L, - -1.11111111100562372733e-01L, - 9.09090902935647302252e-02L, - -7.69230552476207730353e-02L, - 6.66661718042406260546e-02L, - -5.88158892835030888692e-02L, - 5.25499891539726639379e-02L, - -4.70119845393155721494e-02L, - 4.03539201366454414072e-02L, - -2.91303858419364158725e-02L, - 1.24822046299269234080e-02L, -}; - +const long double pi_hi = 3.1415926535897932384626433832795L; const long double pi_lo = -5.01655761266833202345e-20L; +const long double pio2_hi = 1.57079632679489661926L; +const long double pio2_lo = -2.50827880633416601173e-20L; + #endif diff --git a/src/math/__invtrigl.h b/src/math/__invtrigl.h index c3ad3c49..22748b68 100644 --- a/src/math/__invtrigl.h +++ b/src/math/__invtrigl.h @@ -32,15 +32,6 @@ #define BIAS (LDBL_MAX_EXP - 1) #define MANH_SIZE LDBL_MANH_SIZE -/* Approximation thresholds. */ -#define ASIN_LINEAR (BIAS - 32) /* 2**-32 */ -#define ACOS_CONST (BIAS - 65) /* 2**-65 */ -#define ATAN_CONST (BIAS + 65) /* 2**65 */ -#define ATAN_LINEAR (BIAS - 32) /* 2**-32 */ - -/* 0.95 */ -#define THRESH ((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|LDBL_NBIT) - /* Constants shared by the long double inverse trig functions. */ #define pS0 __pS0 #define pS1 __pS1 @@ -54,56 +45,24 @@ #define qS3 __qS3 #define qS4 __qS4 #define qS5 __qS5 -#define atanhi __atanhi -#define atanlo __atanlo -#define aT __aT +#define pi_hi __pi_hi #define pi_lo __pi_lo +#define pio2_hi __pio2_hi +#define pio2_lo __pio2_lo -#define pio2_hi atanhi[3] -#define pio2_lo atanlo[3] -#define pio4_hi atanhi[1] - -#ifdef STRUCT_DECLS -typedef struct longdouble { - uint64_t mant; - uint16_t expsign; -} LONGDOUBLE; -#else -typedef long double LONGDOUBLE; -#endif - -extern const LONGDOUBLE pS0, pS1, pS2, pS3, pS4, pS5, pS6; -extern const LONGDOUBLE qS1, qS2, qS3, qS4, qS5; -extern const LONGDOUBLE atanhi[], atanlo[], aT[]; -extern const LONGDOUBLE pi_lo; - -#ifndef STRUCT_DECLS -static inline long double -P(long double x) -{ - return (x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 + x * \ - (pS4 + x * (pS5 + x * pS6))))))); -} - -static inline long double -Q(long double x) -{ - return (1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * qS5))))); -} +extern const long double pS0, pS1, pS2, pS3, pS4, pS5, pS6; +extern const long double qS1, qS2, qS3, qS4, qS5; +extern const long double pi_hi, pi_lo, pio2_hi, pio2_lo; -static inline long double -T_even(long double x) +static long double P(long double x) { - return (aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + x * \ - (aT[8] + x * (aT[10] + x * aT[12])))))); + return x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 + + x * (pS4 + x * (pS5 + x * pS6)))))); } -static inline long double -T_odd(long double x) +static long double Q(long double x) { - return (aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + x * \ - (aT[9] + x * aT[11]))))); + return 1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * qS5)))); } -#endif #endif diff --git a/src/math/acos.c b/src/math/acos.c index 54d266ee..0eb15bed 100644 --- a/src/math/acos.c +++ b/src/math/acos.c @@ -38,6 +38,7 @@ static const double pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ pio2_hi = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */ +// FIXME static const volatile double pio2_lo = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */ static const double diff --git a/src/math/acosl.c b/src/math/acosl.c index cc565336..83857d49 100644 --- a/src/math/acosl.c +++ b/src/math/acosl.c @@ -23,9 +23,7 @@ long double acosl(long double x) } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 #include "__invtrigl.h" - -static const long double -pi = 3.14159265358979323846264338327950280e+00L; +#define ACOS_CONST (BIAS - 65) /* 2**-65 */ long double acosl(long double x) { @@ -41,7 +39,8 @@ long double acosl(long double x) if (expsign > 0) return 0.0; /* acos(1) = 0 */ else - return pi + 2.0 * pio2_lo; /* acos(-1)= pi */ + // FIXME + return pi_hi + 2.0 * pio2_lo; /* acos(-1)= pi */ } return (x - x) / (x - x); /* acos(|x|>1) is NaN */ } @@ -60,7 +59,7 @@ long double acosl(long double x) s = sqrtl(z); r = p / q; w = r * s - pio2_lo; - return pi - 2.0 * (s + w); + return pi_hi - 2.0 * (s + w); } else { /* x > 0.5 */ z = (1.0 - x) * 0.5; s = sqrtl(z); diff --git a/src/math/asinl.c b/src/math/asinl.c index ddd807e2..7572767c 100644 --- a/src/math/asinl.c +++ b/src/math/asinl.c @@ -24,6 +24,10 @@ long double asinl(long double x) #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 #include "__invtrigl.h" static const long double huge = 1.000e+300; +static const long double pio4_hi = 7.85398163397448309628e-01L; +#define ASIN_LINEAR (BIAS - 32) /* 2**-32 */ +/* 0.95 */ +#define THRESH ((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|LDBL_NBIT) long double asinl(long double x) { diff --git a/src/math/atan2.c b/src/math/atan2.c index 143c3834..825a38fe 100644 --- a/src/math/atan2.c +++ b/src/math/atan2.c @@ -39,6 +39,7 @@ #include "libm.h" +// FIXME static const volatile double tiny = 1.0e-300; static const double diff --git a/src/math/atan2l.c b/src/math/atan2l.c index 0fc901c8..45cbfcc1 100644 --- a/src/math/atan2l.c +++ b/src/math/atan2l.c @@ -24,10 +24,8 @@ long double atan2l(long double y, long double x) } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 #include "__invtrigl.h" -static const volatile long double -tiny = 1.0e-300; -static const long double -pi = 3.14159265358979323846264338327950280e+00L; +// FIXME: +static const volatile long double tiny = 1.0e-300; long double atan2l(long double y, long double x) { @@ -55,9 +53,9 @@ long double atan2l(long double y, long double x) if (expty==0 && ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)==0) { switch(m) { case 0: - case 1: return y; /* atan(+-0,+anything)=+-0 */ - case 2: return pi+tiny; /* atan(+0,-anything) = pi */ - case 3: return -pi-tiny; /* atan(-0,-anything) =-pi */ + case 1: return y; /* atan(+-0,+anything)=+-0 */ + case 2: return pi_hi+tiny; /* atan(+0,-anything) = pi */ + case 3: return -pi_hi-tiny; /* atan(-0,-anything) =-pi */ } } /* when x = 0 */ @@ -69,15 +67,15 @@ long double atan2l(long double y, long double x) switch(m) { case 0: return pio2_hi*0.5+tiny; /* atan(+INF,+INF) */ case 1: return -pio2_hi*0.5-tiny; /* atan(-INF,+INF) */ - case 2: return 1.5*pio2_hi+tiny; /*atan(+INF,-INF)*/ - case 3: return -1.5*pio2_hi-tiny; /*atan(-INF,-INF)*/ + case 2: return 1.5*pio2_hi+tiny; /* atan(+INF,-INF) */ + case 3: return -1.5*pio2_hi-tiny; /* atan(-INF,-INF) */ } } else { switch(m) { - case 0: return 0.0; /* atan(+...,+INF) */ - case 1: return -0.0; /* atan(-...,+INF) */ - case 2: return pi+tiny; /* atan(+...,-INF) */ - case 3: return -pi-tiny; /* atan(-...,-INF) */ + case 0: return 0.0; /* atan(+...,+INF) */ + case 1: return -0.0; /* atan(-...,+INF) */ + case 2: return pi_hi+tiny; /* atan(+...,-INF) */ + case 3: return -pi_hi-tiny; /* atan(-...,-INF) */ } } } @@ -95,11 +93,11 @@ long double atan2l(long double y, long double x) else /* safe to do y/x */ z = atanl(fabsl(y/x)); switch (m) { - case 0: return z; /* atan(+,+) */ - case 1: return -z; /* atan(-,+) */ - case 2: return pi - (z-pi_lo); /* atan(+,-) */ + case 0: return z; /* atan(+,+) */ + case 1: return -z; /* atan(-,+) */ + case 2: return pi_hi-(z-pi_lo); /* atan(+,-) */ default: /* case 3 */ - return (z-pi_lo) - pi; /* atan(-,-) */ + return (z-pi_lo)-pi_hi; /* atan(-,-) */ } } #endif diff --git a/src/math/atanl.c b/src/math/atanl.c index 36072c17..33ecff0d 100644 --- a/src/math/atanl.c +++ b/src/math/atanl.c @@ -23,8 +23,53 @@ long double atanl(long double x) } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 #include "__invtrigl.h" + +#define ATAN_CONST (BIAS + 65) /* 2**65 */ +#define ATAN_LINEAR (BIAS - 32) /* 2**-32 */ static const long double huge = 1.0e300; +static const long double atanhi[] = { + 4.63647609000806116202e-01L, + 7.85398163397448309628e-01L, + 9.82793723247329067960e-01L, + 1.57079632679489661926e+00L, +}; + +static const long double atanlo[] = { + 1.18469937025062860669e-20L, + -1.25413940316708300586e-20L, + 2.55232234165405176172e-20L, + -2.50827880633416601173e-20L, +}; + +static const long double aT[] = { + 3.33333333333333333017e-01L, + -1.99999999999999632011e-01L, + 1.42857142857046531280e-01L, + -1.11111111100562372733e-01L, + 9.09090902935647302252e-02L, + -7.69230552476207730353e-02L, + 6.66661718042406260546e-02L, + -5.88158892835030888692e-02L, + 5.25499891539726639379e-02L, + -4.70119845393155721494e-02L, + 4.03539201366454414072e-02L, + -2.91303858419364158725e-02L, + 1.24822046299269234080e-02L, +}; + +static long double T_even(long double x) +{ + return aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + + x * (aT[8] + x * (aT[10] + x * aT[12]))))); +} + +static long double T_odd(long double x) +{ + return aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + + x * (aT[9] + x * aT[11])))); +} + long double atanl(long double x) { union IEEEl2bits u; diff --git a/src/math/exp10l.c b/src/math/exp10l.c index 4b2b7ef6..22a4636a 100644 --- a/src/math/exp10l.c +++ b/src/math/exp10l.c @@ -5,8 +5,8 @@ long double exp10l(long double x) { static const long double p10[] = { - 1e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-10, - 1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, + 1e-15L, 1e-14L, 1e-13L, 1e-12L, 1e-11L, 1e-10L, + 1e-9L, 1e-8L, 1e-7L, 1e-6L, 1e-5L, 1e-4L, 1e-3L, 1e-2L, 1e-1L, 1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, 1e12, 1e13, 1e14, 1e15 }; diff --git a/src/math/expm1l.c b/src/math/expm1l.c index ca0d141f..21a86c00 100644 --- a/src/math/expm1l.c +++ b/src/math/expm1l.c @@ -44,13 +44,7 @@ * * Relative error: * arithmetic domain # trials peak rms - * IEEE -45,+MAXLOG 200,000 1.2e-19 2.5e-20 - * - * ERROR MESSAGES: - * - * message condition value returned - * expm1l overflow x > MAXLOG MAXNUM - * + * IEEE -45,+maxarg 200,000 1.2e-19 2.5e-20 */ #include "libm.h" @@ -61,7 +55,6 @@ long double expm1l(long double x) return expm1(x); } #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 -static const long double MAXLOGL = 1.1356523406294143949492E4L; /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x) -.5 ln 2 < x < .5 ln 2 @@ -83,19 +76,20 @@ C1 = 6.93145751953125E-1L, C2 = 1.428606820309417232121458176568075500134E-6L, /* ln 2^-65 */ minarg = -4.5054566736396445112120088E1L, -huge = 0x1p10000L; +/* ln 2^16384 */ +maxarg = 1.1356523406294143949492E4L; long double expm1l(long double x) { long double px, qx, xx; int k; - /* Overflow. */ - if (x > MAXLOGL) - return huge*huge; /* overflow */ + if (isnan(x)) + return x; + if (x > maxarg) + return x*0x1p16383L; /* overflow, unless x==inf */ if (x == 0.0) return x; - /* Minimum value.*/ if (x < minarg) return -1.0; diff --git a/src/math/fma.c b/src/math/fma.c index 7076d4b9..17f1cdcc 100644 --- a/src/math/fma.c +++ b/src/math/fma.c @@ -89,6 +89,7 @@ static int getexp(long double x) double fma(double x, double y, double z) { + #pragma STDC FENV_ACCESS ON long double hi, lo1, lo2, xy; int round, ez, exy; @@ -306,6 +307,7 @@ static inline struct dd dd_mul(double a, double b) */ double fma(double x, double y, double z) { + #pragma STDC FENV_ACCESS ON double xs, ys, zs, adj; struct dd xy, r; int oround; diff --git a/src/math/fmaf.c b/src/math/fmaf.c index cc427fb5..a1c7f4f8 100644 --- a/src/math/fmaf.c +++ b/src/math/fmaf.c @@ -37,6 +37,7 @@ */ float fmaf(float x, float y, float z) { + #pragma STDC FENV_ACCESS ON double xy, result; uint32_t hr, lr; diff --git a/src/math/fmal.c b/src/math/fmal.c index be64f145..ccbe434d 100644 --- a/src/math/fmal.c +++ b/src/math/fmal.c @@ -162,6 +162,7 @@ static inline struct dd dd_mul(long double a, long double b) */ long double fmal(long double x, long double y, long double z) { + #pragma STDC FENV_ACCESS ON long double xs, ys, zs, adj; struct dd xy, r; int oround; diff --git a/src/math/hypot.c b/src/math/hypot.c index ba4c7575..9a4cbdb3 100644 --- a/src/math/hypot.c +++ b/src/math/hypot.c @@ -117,12 +117,7 @@ double hypot(double x, double y) t2 = a - t1; w = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b))); } - if (k != 0) { - uint32_t high; - t1 = 1.0; - GET_HIGH_WORD(high, t1); - SET_HIGH_WORD(t1, high+(k<<20)); - return t1*w; - } + if (k) + w = scalbn(w, k); return w; } diff --git a/src/math/hypotf.c b/src/math/hypotf.c index 9fd77e6a..4d80178d 100644 --- a/src/math/hypotf.c +++ b/src/math/hypotf.c @@ -80,9 +80,7 @@ float hypotf(float x, float y) t2 = a - t1; w = sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b))); } - if (k != 0) { - SET_FLOAT_WORD(t1, 0x3f800000+(k<<23)); - return t1*w; - } + if (k) + w = scalbnf(w, k); return w; } diff --git a/src/math/ilogb.c b/src/math/ilogb.c index 0a3a6a46..5a1819d8 100644 --- a/src/math/ilogb.c +++ b/src/math/ilogb.c @@ -8,13 +8,17 @@ int ilogb(double x) if (!e) { u.bits <<= 12; - if (u.bits == 0) + if (u.bits == 0) { + FORCE_EVAL(0/0.0f); return FP_ILOGB0; + } /* subnormal x */ for (e = -0x3ff; u.bits < (uint64_t)1<<63; e--, u.bits<<=1); return e; } - if (e == 0x7ff) + if (e == 0x7ff) { + FORCE_EVAL(0/0.0f); return u.bits<<12 ? FP_ILOGBNAN : INT_MAX; + } return e - 0x3ff; } diff --git a/src/math/ilogbf.c b/src/math/ilogbf.c index 272cbdac..42cd62e2 100644 --- a/src/math/ilogbf.c +++ b/src/math/ilogbf.c @@ -8,13 +8,17 @@ int ilogbf(float x) if (!e) { u.bits <<= 9; - if (u.bits == 0) + if (u.bits == 0) { + FORCE_EVAL(0/0.0f); return FP_ILOGB0; + } /* subnormal x */ for (e = -0x7f; u.bits < (uint32_t)1<<31; e--, u.bits<<=1); return e; } - if (e == 0xff) + if (e == 0xff) { + FORCE_EVAL(0/0.0f); return u.bits<<9 ? FP_ILOGBNAN : INT_MAX; + } return e - 0x7f; } diff --git a/src/math/ilogbl.c b/src/math/ilogbl.c index ed9ddcbc..1512934f 100644 --- a/src/math/ilogbl.c +++ b/src/math/ilogbl.c @@ -14,15 +14,19 @@ int ilogbl(long double x) int e = u.bits.exp; if (!e) { - if (m == 0) + if (m == 0) { + FORCE_EVAL(0/0.0f); return FP_ILOGB0; + } /* subnormal x */ for (e = -0x3fff+1; m < (uint64_t)1<<63; e--, m<<=1); return e; } - if (e == 0x7fff) + if (e == 0x7fff) { + FORCE_EVAL(0/0.0f); /* in ld80 msb is set in inf */ return m & (uint64_t)-1>>1 ? FP_ILOGBNAN : INT_MAX; + } return e - 0x3fff; } #endif diff --git a/src/math/llrintl.c b/src/math/llrintl.c index 32bb8b03..3449f6f2 100644 --- a/src/math/llrintl.c +++ b/src/math/llrintl.c @@ -18,6 +18,7 @@ raises inexact (with tonearest or upward rounding mode) */ long long llrintl(long double x) { + #pragma STDC FENV_ACCESS ON int e; e = fetestexcept(FE_INEXACT); diff --git a/src/math/log1pl.c b/src/math/log1pl.c index 1400d365..edb48df1 100644 --- a/src/math/log1pl.c +++ b/src/math/log1pl.c @@ -46,11 +46,6 @@ * Relative error: * arithmetic domain # trials peak rms * IEEE -1.0, 9.0 100000 8.2e-20 2.5e-20 - * - * ERROR MESSAGES: - * - * log singularity: x-1 = 0; returns -INFINITY - * log domain: x-1 < 0; returns NAN */ #include "libm.h" @@ -123,8 +118,8 @@ long double log1pl(long double xm1) /* Test for domain errors. */ if (x <= 0.0) { if (x == 0.0) - return -INFINITY; - return NAN; + return -1/x; /* -inf with divbyzero */ + return 0/0.0f; /* nan with invalid */ } /* Separate mantissa from exponent. diff --git a/src/math/log2l.c b/src/math/log2l.c index 8ebce9c4..345b395d 100644 --- a/src/math/log2l.c +++ b/src/math/log2l.c @@ -50,11 +50,6 @@ * In the tests over the interval exp(+-10000), the logarithms * of the random arguments were uniformly distributed over * [-10000, +10000]. - * - * ERROR MESSAGES: - * - * log singularity: x = 0; returns -INFINITY - * log domain: x < 0; returns NAN */ #include "libm.h" @@ -113,8 +108,7 @@ static const long double S[4] = { long double log2l(long double x) { - volatile long double z; - long double y; + long double y, z; int e; if (isnan(x)) @@ -123,8 +117,8 @@ long double log2l(long double x) return x; if (x <= 0.0) { if (x == 0.0) - return -INFINITY; - return NAN; + return -1/(x+0); /* -inf with divbyzero */ + return 0/0.0f; /* nan with invalid */ } /* separate mantissa from exponent */ diff --git a/src/math/logb.c b/src/math/logb.c index f7cd7613..624425a8 100644 --- a/src/math/logb.c +++ b/src/math/logb.c @@ -1,20 +1,17 @@ -#include <limits.h> #include "libm.h" /* special cases: - logb(+-0) = -inf + logb(+-0) = -inf, and raise divbyzero logb(+-inf) = +inf logb(nan) = nan -these are calculated at runtime to raise fp exceptions */ -double logb(double x) { - int i = ilogb(x); - - if (i == FP_ILOGB0) - return -1.0/fabs(x); - if (i == FP_ILOGBNAN || i == INT_MAX) +double logb(double x) +{ + if (!isfinite(x)) return x * x; - return i; + if (x == 0) + return -1/(x+0); + return ilogb(x); } diff --git a/src/math/logbf.c b/src/math/logbf.c index 934827f8..950d3569 100644 --- a/src/math/logbf.c +++ b/src/math/logbf.c @@ -1,12 +1,10 @@ -#include <limits.h> #include "libm.h" -float logbf(float x) { - int i = ilogbf(x); - - if (i == FP_ILOGB0) - return -1.0f/fabsf(x); - if (i == FP_ILOGBNAN || i == INT_MAX) +float logbf(float x) +{ + if (!isfinite(x)) return x * x; - return i; + if (x == 0) + return -1/(x+0); + return ilogbf(x); } diff --git a/src/math/logbl.c b/src/math/logbl.c index 5d04abd7..f3850745 100644 --- a/src/math/logbl.c +++ b/src/math/logbl.c @@ -1,4 +1,3 @@ -#include <limits.h> #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 long double logbl(long double x) @@ -8,12 +7,10 @@ long double logbl(long double x) #else long double logbl(long double x) { - int i = ilogbl(x); - - if (i == FP_ILOGB0) - return -1.0/fabsl(x); - if (i == FP_ILOGBNAN || i == INT_MAX) + if (!isfinite(x)) return x * x; - return i; + if (x == 0) + return -1/(x+0); + return ilogbl(x); } #endif diff --git a/src/math/logl.c b/src/math/logl.c index ffd81038..ef2b5515 100644 --- a/src/math/logl.c +++ b/src/math/logl.c @@ -50,11 +50,6 @@ * In the tests over the interval exp(+-10000), the logarithms * of the random arguments were uniformly distributed over * [-10000, +10000]. - * - * ERROR MESSAGES: - * - * log singularity: x = 0; returns -INFINITY - * log domain: x < 0; returns NAN */ #include "libm.h" @@ -121,8 +116,8 @@ long double logl(long double x) return x; if (x <= 0.0) { if (x == 0.0) - return -INFINITY; - return NAN; + return -1/(x+0); /* -inf with divbyzero */ + return 0/0.0f; /* nan with invalid */ } /* separate mantissa from exponent */ diff --git a/src/math/lrint.c b/src/math/lrint.c index 67091835..bdca8b7c 100644 --- a/src/math/lrint.c +++ b/src/math/lrint.c @@ -28,6 +28,7 @@ as a double. #if LONG_MAX < 1U<<53 && defined(FE_INEXACT) long lrint(double x) { + #pragma STDC FENV_ACCESS ON int e; e = fetestexcept(FE_INEXACT); diff --git a/src/math/lrintl.c b/src/math/lrintl.c index 5eb1ba7e..b2a8106d 100644 --- a/src/math/lrintl.c +++ b/src/math/lrintl.c @@ -18,6 +18,7 @@ raises inexact (with tonearest or upward rounding mode) */ long lrintl(long double x) { + #pragma STDC FENV_ACCESS ON int e; e = fetestexcept(FE_INEXACT); diff --git a/src/math/modf.c b/src/math/modf.c index de45069f..8f337ef0 100644 --- a/src/math/modf.c +++ b/src/math/modf.c @@ -1,5 +1,4 @@ -#include <math.h> -#include <stdint.h> +#include "libm.h" double modf(double x, double *iptr) { @@ -33,5 +32,6 @@ double modf(double x, double *iptr) } u.n &= ~mask; *iptr = u.x; - return x - *iptr; + STRICT_ASSIGN(double, x, x - *iptr); + return x; } diff --git a/src/math/modff.c b/src/math/modff.c index 84d0b82a..90caf316 100644 --- a/src/math/modff.c +++ b/src/math/modff.c @@ -1,5 +1,4 @@ -#include <math.h> -#include <stdint.h> +#include "libm.h" float modff(float x, float *iptr) { @@ -33,5 +32,6 @@ float modff(float x, float *iptr) } u.n &= ~mask; *iptr = u.x; - return x - *iptr; + STRICT_ASSIGN(float, x, x - *iptr); + return x; } diff --git a/src/math/nearbyint.c b/src/math/nearbyint.c index 7a4c58cf..f4e8aac4 100644 --- a/src/math/nearbyint.c +++ b/src/math/nearbyint.c @@ -6,6 +6,7 @@ double nearbyint(double x) { #ifdef FE_INEXACT + #pragma STDC FENV_ACCESS ON int e; e = fetestexcept(FE_INEXACT); diff --git a/src/math/nearbyintf.c b/src/math/nearbyintf.c index 39c3d73b..092e9ffa 100644 --- a/src/math/nearbyintf.c +++ b/src/math/nearbyintf.c @@ -4,6 +4,7 @@ float nearbyintf(float x) { #ifdef FE_INEXACT + #pragma STDC FENV_ACCESS ON int e; e = fetestexcept(FE_INEXACT); diff --git a/src/math/nearbyintl.c b/src/math/nearbyintl.c index 0ff4b1f9..82852492 100644 --- a/src/math/nearbyintl.c +++ b/src/math/nearbyintl.c @@ -11,6 +11,7 @@ long double nearbyintl(long double x) long double nearbyintl(long double x) { #ifdef FE_INEXACT + #pragma STDC FENV_ACCESS ON int e; e = fetestexcept(FE_INEXACT); diff --git a/src/math/nextafter.c b/src/math/nextafter.c index a3b42c99..9ee82518 100644 --- a/src/math/nextafter.c +++ b/src/math/nextafter.c @@ -27,7 +27,7 @@ double nextafter(double x, double y) e = ux.bits >> 52 & 0x7ff; /* raise overflow if ux.value is infinite and x is finite */ if (e == 0x7ff) - return x + x; + FORCE_EVAL(x+x); /* raise underflow if ux.value is subnormal or zero */ if (e == 0) FORCE_EVAL(x*x + ux.value*ux.value); diff --git a/src/math/nextafterf.c b/src/math/nextafterf.c index b703487b..22b61dce 100644 --- a/src/math/nextafterf.c +++ b/src/math/nextafterf.c @@ -26,7 +26,7 @@ float nextafterf(float x, float y) e = ux.bits & 0x7f800000; /* raise overflow if ux.value is infinite and x is finite */ if (e == 0x7f800000) - return x + x; + FORCE_EVAL(x+x); /* raise underflow if ux.value is subnormal or zero */ if (e == 0) FORCE_EVAL(x*x + ux.value*ux.value); diff --git a/src/math/nexttoward.c b/src/math/nexttoward.c index 7355f2f1..6f32eca4 100644 --- a/src/math/nexttoward.c +++ b/src/math/nexttoward.c @@ -36,7 +36,7 @@ double nexttoward(double x, long double y) e = ux.bits>>52 & 0x7ff; /* raise overflow if ux.value is infinite and x is finite */ if (e == 0x7ff) - return x + x; + FORCE_EVAL(x+x); /* raise underflow if ux.value is subnormal or zero */ if (e == 0) FORCE_EVAL(x*x + ux.value*ux.value); diff --git a/src/math/nexttowardf.c b/src/math/nexttowardf.c index 8648be6a..9a693b1a 100644 --- a/src/math/nexttowardf.c +++ b/src/math/nexttowardf.c @@ -28,7 +28,7 @@ float nexttowardf(float x, long double y) e = ux.bits & 0x7f800000; /* raise overflow if ux.value is infinite and x is finite */ if (e == 0x7f800000) - return x + x; + FORCE_EVAL(x+x); /* raise underflow if ux.value is subnormal or zero */ if (e == 0) FORCE_EVAL(x*x + ux.value*ux.value); diff --git a/src/math/scalbn.c b/src/math/scalbn.c index c9c7af80..003141e3 100644 --- a/src/math/scalbn.c +++ b/src/math/scalbn.c @@ -10,8 +10,10 @@ double scalbn(double x, int n) if (n > 1023) { x *= 0x1p1023; n -= 1023; - if (n > 1023) - return x * 0x1p1023; + if (n > 1023) { + STRICT_ASSIGN(double, x, x * 0x1p1023); + return x; + } } } else if (n < -1022) { x *= 0x1p-1022; @@ -19,10 +21,13 @@ double scalbn(double x, int n) if (n < -1022) { x *= 0x1p-1022; n += 1022; - if (n < -1022) - return x * 0x1p-1022; + if (n < -1022) { + STRICT_ASSIGN(double, x, x * 0x1p-1022); + return x; + } } } INSERT_WORDS(scale, (uint32_t)(0x3ff+n)<<20, 0); - return x * scale; + STRICT_ASSIGN(double, x, x * scale); + return x; } diff --git a/src/math/scalbnf.c b/src/math/scalbnf.c index 243dafdf..f94b5d59 100644 --- a/src/math/scalbnf.c +++ b/src/math/scalbnf.c @@ -10,8 +10,10 @@ float scalbnf(float x, int n) if (n > 127) { x *= 0x1p127f; n -= 127; - if (n > 127) - return x * 0x1p127f; + if (n > 127) { + STRICT_ASSIGN(float, x, x * 0x1p127f); + return x; + } } } else if (n < -126) { x *= 0x1p-126f; @@ -19,10 +21,13 @@ float scalbnf(float x, int n) if (n < -126) { x *= 0x1p-126f; n += 126; - if (n < -126) - return x * 0x1p-126f; + if (n < -126) { + STRICT_ASSIGN(float, x, x * 0x1p-126f); + return x; + } } } SET_FLOAT_WORD(scale, (uint32_t)(0x7f+n)<<23); - return x * scale; + STRICT_ASSIGN(float, x, x * scale); + return x; } |