diff options
author | das <das@FreeBSD.org> | 2011-10-21 06:29:32 +0000 |
---|---|---|
committer | das <das@FreeBSD.org> | 2011-10-21 06:29:32 +0000 |
commit | 3578083bd5e5c8cc7b2d82f501d5901051537ae4 (patch) | |
tree | e4935310dc992380b0af45c9fd0a57d6e8febb84 /lib | |
parent | f66ae96060daab6b90cb9d24fa049e4c17f3b80f (diff) | |
download | FreeBSD-src-3578083bd5e5c8cc7b2d82f501d5901051537ae4.zip FreeBSD-src-3578083bd5e5c8cc7b2d82f501d5901051537ae4.tar.gz |
Improved handling of large x in ccosh{,f}():
- Handle cases where exp(x) would overflow, but ccosh(x) ~= exp(x) / 2
shouldn't.
- Use the ccosh(x) ~= exp(x) / 2 approximation to simplify the calculation
when x is large.
Similarly for csinh(). Also fixed the return value of csinh(-Inf +- 0i).
Diffstat (limited to 'lib')
-rw-r--r-- | lib/msun/src/s_ccosh.c | 23 | ||||
-rw-r--r-- | lib/msun/src/s_ccoshf.c | 23 | ||||
-rw-r--r-- | lib/msun/src/s_csinh.c | 27 | ||||
-rw-r--r-- | lib/msun/src/s_csinhf.c | 25 |
4 files changed, 83 insertions, 15 deletions
diff --git a/lib/msun/src/s_ccosh.c b/lib/msun/src/s_ccosh.c index 4086394..9ea962b 100644 --- a/lib/msun/src/s_ccosh.c +++ b/lib/msun/src/s_ccosh.c @@ -42,10 +42,12 @@ __FBSDID("$FreeBSD$"); #include "math_private.h" +static const double huge = 0x1p1023; + double complex ccosh(double complex z) { - double x, y; + double x, y, h; int32_t hx, hy, ix, iy, lx, ly; x = creal(z); @@ -61,8 +63,23 @@ ccosh(double complex z) if (ix < 0x7ff00000 && iy < 0x7ff00000) { if ((iy | ly) == 0) return (cpack(cosh(x), x * y)); - /* XXX We don't handle |x| > DBL_MAX ln(2) yet. */ - return (cpack(cosh(x) * cos(y), sinh(x) * sin(y))); + if (ix < 0x40360000) /* small x: normal case */ + return (cpack(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))); + } 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))); + } else { + /* x >= 1455: the result always overflows */ + h = huge * x; + return (cpack(h * h * cos(y), h * sin(y))); + } } /* diff --git a/lib/msun/src/s_ccoshf.c b/lib/msun/src/s_ccoshf.c index bf9b78f..1de9ad4 100644 --- a/lib/msun/src/s_ccoshf.c +++ b/lib/msun/src/s_ccoshf.c @@ -36,10 +36,12 @@ __FBSDID("$FreeBSD$"); #include "math_private.h" +static const float huge = 0x1p127; + float complex ccoshf(float complex z) { - float x, y; + float x, y, h; int32_t hx, hy, ix, iy; x = crealf(z); @@ -54,8 +56,23 @@ ccoshf(float complex z) if (ix < 0x7f800000 && iy < 0x7f800000) { if (iy == 0) return (cpackf(coshf(x), x * y)); - /* XXX We don't handle |x| > FLT_MAX ln(2) yet. */ - return (cpackf(coshf(x) * cosf(y), sinhf(x) * sinf(y))); + if (ix < 0x41100000) /* small x: normal case */ + return (cpackf(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))); + } 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))); + } else { + /* x >= 192.7: the result always overflows */ + h = huge * x; + return (cpackf(h * h * cosf(y), h * sinf(y))); + } } if (ix == 0 && iy >= 0x7f800000) diff --git a/lib/msun/src/s_csinh.c b/lib/msun/src/s_csinh.c index 74a0aff..c192f30 100644 --- a/lib/msun/src/s_csinh.c +++ b/lib/msun/src/s_csinh.c @@ -42,10 +42,12 @@ __FBSDID("$FreeBSD$"); #include "math_private.h" +static const double huge = 0x1p1023; + double complex csinh(double complex z) { - double x, y; + double x, y, h; int32_t hx, hy, ix, iy, lx, ly; x = creal(z); @@ -61,8 +63,23 @@ csinh(double complex z) if (ix < 0x7ff00000 && iy < 0x7ff00000) { if ((iy | ly) == 0) return (cpack(sinh(x), y)); - /* XXX We don't handle |x| > DBL_MAX ln(2) yet. */ - return (cpack(sinh(x) * cos(y), cosh(x) * sin(y))); + if (ix < 0x40360000) /* small x: normal case */ + return (cpack(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))); + } 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))); + } else { + /* x >= 1455: the result always overflows */ + h = huge * x; + return (cpack(h * cos(y), h * h * sin(y))); + } } /* @@ -78,13 +95,13 @@ csinh(double complex z) return (cpack(copysign(0, x * (y - y)), y - y)); /* - * sinh(+-Inf +- I 0) = +-Inf + I (+-)(+-)0. + * sinh(+-Inf +- I 0) = +-Inf + I +-0. * * sinh(NaN +- I 0) = d(NaN) + I +-0. */ if ((iy | ly) == 0 && ix >= 0x7ff00000) { if (((hx & 0xfffff) | lx) == 0) - return (cpack(x, copysign(0, x) * y)); + return (cpack(x, y)); return (cpack(x, copysign(0, y))); } diff --git a/lib/msun/src/s_csinhf.c b/lib/msun/src/s_csinhf.c index e1478ac..c523125 100644 --- a/lib/msun/src/s_csinhf.c +++ b/lib/msun/src/s_csinhf.c @@ -36,10 +36,12 @@ __FBSDID("$FreeBSD$"); #include "math_private.h" +static const float huge = 0x1p127; + float complex csinhf(float complex z) { - float x, y; + float x, y, h; int32_t hx, hy, ix, iy; x = crealf(z); @@ -54,8 +56,23 @@ csinhf(float complex z) if (ix < 0x7f800000 && iy < 0x7f800000) { if (iy == 0) return (cpackf(sinhf(x), y)); - /* XXX We don't handle |x| > FLT_MAX ln(2) yet. */ - return (cpackf(sinhf(x) * cosf(y), coshf(x) * sinf(y))); + if (ix < 0x41100000) /* small x: normal case */ + return (cpackf(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))); + } 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))); + } else { + /* x >= 192.7: the result always overflows */ + h = huge * x; + return (cpackf(h * cosf(y), h * h * sinf(y))); + } } if (ix == 0 && iy >= 0x7f800000) @@ -63,7 +80,7 @@ csinhf(float complex z) if (iy == 0 && ix >= 0x7f800000) { if ((hx & 0x7fffff) == 0) - return (cpackf(x, copysignf(0, x) * y)); + return (cpackf(x, y)); return (cpackf(x, copysignf(0, y))); } |