diff options
Diffstat (limited to 'lib/msun/src')
| -rw-r--r-- | lib/msun/src/s_ccosh.c | 41 | ||||
| -rw-r--r-- | lib/msun/src/s_ccoshf.c | 25 | ||||
| -rw-r--r-- | lib/msun/src/s_csinh.c | 59 | ||||
| -rw-r--r-- | lib/msun/src/s_csinhf.c | 31 | ||||
| -rw-r--r-- | lib/msun/src/s_ctanh.c | 35 | ||||
| -rw-r--r-- | lib/msun/src/s_ctanhf.c | 9 |
6 files changed, 98 insertions, 102 deletions
diff --git a/lib/msun/src/s_ccosh.c b/lib/msun/src/s_ccosh.c index fbe15fc..2843198 100644 --- a/lib/msun/src/s_ccosh.c +++ b/lib/msun/src/s_ccosh.c @@ -32,6 +32,8 @@ * * Exceptional values are noted in the comments within the source code. * These values and the return value were taken from n1124.pdf. + * The sign of the result for some exceptional values is unspecified but + * must satisfy both cosh(conj(z)) == conj(cosh(z)) and cosh(-z) == cosh(z). */ #include <sys/cdefs.h> @@ -63,7 +65,7 @@ ccosh(double complex z) if (ix < 0x7ff00000 && iy < 0x7ff00000) { if ((iy | ly) == 0) return (CMPLX(cosh(x), x * y)); - if (ix < 0x40360000) /* small x: normal case */ + if (ix < 0x40360000) /* |x| < 22: normal case */ return (CMPLX(cosh(x) * cos(y), sinh(x) * sin(y))); /* |x| >= 22, so cosh(x) ~= exp(|x|) */ @@ -83,28 +85,27 @@ ccosh(double complex z) } /* - * cosh(+-0 +- I Inf) = dNaN + I sign(d(+-0, dNaN))0. - * The sign of 0 in the result is unspecified. Choice = normally - * the same as dNaN. Raise the invalid floating-point exception. + * cosh(+-0 +- I Inf) = dNaN + I (+-)(+-)0. + * The sign of 0 in the result is unspecified. Choice = product + * of the signs of the argument. Raise the invalid floating-point + * exception. * - * cosh(+-0 +- I NaN) = d(NaN) + I sign(d(+-0, NaN))0. - * The sign of 0 in the result is unspecified. Choice = normally - * the same as d(NaN). + * cosh(+-0 +- I NaN) = d(NaN) + I (+-)(+-)0. + * The sign of 0 in the result is unspecified. Choice = product + * of the signs of the argument. */ - if ((ix | lx) == 0 && iy >= 0x7ff00000) - return (CMPLX(y - y, copysign(0, x * (y - y)))); + if ((ix | lx) == 0) /* && iy >= 0x7ff00000 */ + return (CMPLX(y - y, x * copysign(0, y))); /* * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0. * - * cosh(NaN +- I 0) = d(NaN) + I sign(d(NaN, +-0))0. - * The sign of 0 in the result is unspecified. + * cosh(NaN +- I 0) = d(NaN) + I (+-)(+-)0. + * The sign of 0 in the result is unspecified. Choice = product + * of the signs of the argument. */ - if ((iy | ly) == 0 && ix >= 0x7ff00000) { - if (((hx & 0xfffff) | lx) == 0) - return (CMPLX(x * x, copysign(0, x) * y)); - return (CMPLX(x * x, copysign(0, (x + x) * y))); - } + if ((iy | ly) == 0) /* && ix >= 0x7ff00000 */ + return (CMPLX(x * x, copysign(0, x) * y)); /* * cosh(x +- I Inf) = dNaN + I dNaN. @@ -114,7 +115,7 @@ ccosh(double complex z) * Optionally raises the invalid floating-point exception for finite * nonzero x. Choice = don't raise (except for signaling NaNs). */ - if (ix < 0x7ff00000 && iy >= 0x7ff00000) + if (ix < 0x7ff00000) /* && iy >= 0x7ff00000 */ return (CMPLX(y - y, x * (y - y))); /* @@ -126,10 +127,10 @@ ccosh(double complex z) * * cosh(+-Inf + I y) = +Inf cos(y) +- I Inf sin(y) */ - if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) { + if (ix == 0x7ff00000 && lx == 0) { if (iy >= 0x7ff00000) - return (CMPLX(x * x, x * (y - y))); - return (CMPLX((x * x) * cos(y), x * sin(y))); + return (CMPLX(INFINITY, x * (y - y))); + return (CMPLX(INFINITY * cos(y), x * sin(y))); } /* diff --git a/lib/msun/src/s_ccoshf.c b/lib/msun/src/s_ccoshf.c index fe8cf89..eeed92f 100644 --- a/lib/msun/src/s_ccoshf.c +++ b/lib/msun/src/s_ccoshf.c @@ -25,7 +25,7 @@ */ /* - * Hyperbolic cosine of a complex argument. See s_ccosh.c for details. + * Float version of ccosh(). See s_ccosh.c for details. */ #include <sys/cdefs.h> @@ -56,13 +56,13 @@ ccoshf(float complex z) if (ix < 0x7f800000 && iy < 0x7f800000) { if (iy == 0) return (CMPLXF(coshf(x), x * y)); - if (ix < 0x41100000) /* small x: normal case */ + if (ix < 0x41100000) /* |x| < 9: normal case */ 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; + h = expf(fabsf(x)) * 0.5F; return (CMPLXF(h * cosf(y), copysignf(h, x) * sinf(y))); } else if (ix < 0x4340b1e7) { /* x < 192.7: scale to avoid overflow */ @@ -75,22 +75,19 @@ ccoshf(float complex z) } } - if (ix == 0 && iy >= 0x7f800000) - return (CMPLXF(y - y, copysignf(0, x * (y - y)))); + if (ix == 0) /* && iy >= 0x7f800000 */ + return (CMPLXF(y - y, x * copysignf(0, y))); - if (iy == 0 && ix >= 0x7f800000) { - if ((hx & 0x7fffff) == 0) - return (CMPLXF(x * x, copysignf(0, x) * y)); - return (CMPLXF(x * x, copysignf(0, (x + x) * y))); - } + if (iy == 0) /* && ix >= 0x7f800000 */ + return (CMPLXF(x * x, copysignf(0, x) * y)); - if (ix < 0x7f800000 && iy >= 0x7f800000) + if (ix < 0x7f800000) /* && iy >= 0x7f800000 */ return (CMPLXF(y - y, x * (y - y))); - if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) { + if (ix == 0x7f800000) { if (iy >= 0x7f800000) - return (CMPLXF(x * x, x * (y - y))); - return (CMPLXF((x * x) * cosf(y), x * sinf(y))); + return (CMPLXF(INFINITY, x * (y - y))); + return (CMPLXF(INFINITY * cosf(y), x * sinf(y))); } return (CMPLXF((x * x) * (y - y), (x + x) * (y - y))); diff --git a/lib/msun/src/s_csinh.c b/lib/msun/src/s_csinh.c index 37f46da..7af3634 100644 --- a/lib/msun/src/s_csinh.c +++ b/lib/msun/src/s_csinh.c @@ -32,6 +32,8 @@ * * Exceptional values are noted in the comments within the source code. * These values and the return value were taken from n1124.pdf. + * The sign of the result for some exceptional values is unspecified but + * must satisfy both sinh(conj(z)) == conj(sinh(z)) and sinh(-z) == -sinh(z). */ #include <sys/cdefs.h> @@ -63,7 +65,7 @@ csinh(double complex z) if (ix < 0x7ff00000 && iy < 0x7ff00000) { if ((iy | ly) == 0) return (CMPLX(sinh(x), y)); - if (ix < 0x40360000) /* small x: normal case */ + if (ix < 0x40360000) /* |x| < 22: normal case */ return (CMPLX(sinh(x) * cos(y), cosh(x) * sin(y))); /* |x| >= 22, so cosh(x) ~= exp(|x|) */ @@ -83,27 +85,24 @@ csinh(double complex z) } /* - * sinh(+-0 +- I Inf) = sign(d(+-0, dNaN))0 + I dNaN. - * The sign of 0 in the result is unspecified. Choice = normally - * the same as dNaN. Raise the invalid floating-point exception. + * sinh(+-0 +- I Inf) = +-0 + I dNaN. + * The sign of 0 in the result is unspecified. Choice = same sign + * as the argument. Raise the invalid floating-point exception. * - * sinh(+-0 +- I NaN) = sign(d(+-0, NaN))0 + I d(NaN). - * The sign of 0 in the result is unspecified. Choice = normally - * the same as d(NaN). + * sinh(+-0 +- I NaN) = +-0 + I d(NaN). + * The sign of 0 in the result is unspecified. Choice = same sign + * as the argument. */ - if ((ix | lx) == 0 && iy >= 0x7ff00000) - return (CMPLX(copysign(0, x * (y - y)), y - y)); + if ((ix | lx) == 0) /* && iy >= 0x7ff00000 */ + return (CMPLX(x, y - y)); /* * 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 (CMPLX(x, y)); - return (CMPLX(x, copysign(0, y))); - } + if ((iy | ly) == 0) /* && ix >= 0x7ff00000 */ + return (CMPLX(x + x, y)); /* * sinh(x +- I Inf) = dNaN + I dNaN. @@ -113,45 +112,45 @@ csinh(double complex z) * Optionally raises the invalid floating-point exception for finite * nonzero x. Choice = don't raise (except for signaling NaNs). */ - if (ix < 0x7ff00000 && iy >= 0x7ff00000) - return (CMPLX(y - y, x * (y - y))); + if (ix < 0x7ff00000) /* && iy >= 0x7ff00000 */ + return (CMPLX(y - y, y - y)); /* * sinh(+-Inf + I NaN) = +-Inf + I d(NaN). - * The sign of Inf in the result is unspecified. Choice = normally - * the same as d(NaN). + * The sign of Inf in the result is unspecified. Choice = same sign + * as the argument. * - * sinh(+-Inf +- I Inf) = +Inf + I dNaN. - * The sign of Inf in the result is unspecified. Choice = always +. - * Raise the invalid floating-point exception. + * sinh(+-Inf +- I Inf) = +-Inf + I dNaN. + * The sign of Inf in the result is unspecified. Choice = same sign + * as the argument. Raise the invalid floating-point exception. * * sinh(+-Inf + I y) = +-Inf cos(y) + I Inf sin(y) */ - if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) { + if (ix == 0x7ff00000 && lx == 0) { if (iy >= 0x7ff00000) - return (CMPLX(x * x, x * (y - y))); + return (CMPLX(x, y - y)); return (CMPLX(x * cos(y), INFINITY * sin(y))); } /* - * sinh(NaN + I NaN) = d(NaN) + I d(NaN). + * sinh(NaN1 + I NaN2) = d(NaN1, NaN2) + I d(NaN1, NaN2). * - * sinh(NaN +- I Inf) = d(NaN) + I d(NaN). + * sinh(NaN +- I Inf) = d(NaN, dNaN) + I d(NaN, dNaN). * Optionally raises the invalid floating-point exception. * Choice = raise. * - * sinh(NaN + I y) = d(NaN) + I d(NaN). + * sinh(NaN + I y) = d(NaN) + I d(NaN). * 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((x + x) * (y - y), (x * x) * (y - y))); } double complex csin(double complex z) { - /* csin(z) = -I * csinh(I * z) */ - z = csinh(CMPLX(-cimag(z), creal(z))); - return (CMPLX(cimag(z), -creal(z))); + /* csin(z) = -I * csinh(I * z) = I * conj(csinh(I * conj(z))). */ + z = csinh(CMPLX(cimag(z), creal(z))); + return (CMPLX(cimag(z), creal(z))); } diff --git a/lib/msun/src/s_csinhf.c b/lib/msun/src/s_csinhf.c index 6348272..c9d9d3c 100644 --- a/lib/msun/src/s_csinhf.c +++ b/lib/msun/src/s_csinhf.c @@ -25,7 +25,7 @@ */ /* - * Hyperbolic sine of a complex argument z. See s_csinh.c for details. + * Float version of csinh(). See s_csinh.c for details. */ #include <sys/cdefs.h> @@ -56,13 +56,13 @@ csinhf(float complex z) if (ix < 0x7f800000 && iy < 0x7f800000) { if (iy == 0) return (CMPLXF(sinhf(x), y)); - if (ix < 0x41100000) /* small x: normal case */ + if (ix < 0x41100000) /* |x| < 9: normal case */ 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; + h = expf(fabsf(x)) * 0.5F; return (CMPLXF(copysignf(h, x) * cosf(y), h * sinf(y))); } else if (ix < 0x4340b1e7) { /* x < 192.7: scale to avoid overflow */ @@ -75,31 +75,28 @@ csinhf(float complex z) } } - if (ix == 0 && iy >= 0x7f800000) - return (CMPLXF(copysignf(0, x * (y - y)), y - y)); + if (ix == 0) /* && iy >= 0x7f800000 */ + return (CMPLXF(x, y - y)); - if (iy == 0 && ix >= 0x7f800000) { - if ((hx & 0x7fffff) == 0) - return (CMPLXF(x, y)); - return (CMPLXF(x, copysignf(0, y))); - } + if (iy == 0) /* && ix >= 0x7f800000 */ + return (CMPLXF(x + x, y)); - if (ix < 0x7f800000 && iy >= 0x7f800000) - return (CMPLXF(y - y, x * (y - y))); + if (ix < 0x7f800000) /* && iy >= 0x7f800000 */ + return (CMPLXF(y - y, y - y)); - if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) { + if (ix == 0x7f800000) { if (iy >= 0x7f800000) - return (CMPLXF(x * x, x * (y - y))); + return (CMPLXF(x, y - y)); return (CMPLXF(x * cosf(y), INFINITY * sinf(y))); } - return (CMPLXF((x * x) * (y - y), (x + x) * (y - y))); + return (CMPLXF((x + x) * (y - y), (x * x) * (y - y))); } float complex csinf(float complex z) { - z = csinhf(CMPLXF(-cimagf(z), crealf(z))); - return (CMPLXF(cimagf(z), -crealf(z))); + z = csinhf(CMPLXF(cimagf(z), crealf(z))); + return (CMPLXF(cimagf(z), crealf(z))); } diff --git a/lib/msun/src/s_ctanh.c b/lib/msun/src/s_ctanh.c index b15c814..f5b9bdd 100644 --- a/lib/msun/src/s_ctanh.c +++ b/lib/msun/src/s_ctanh.c @@ -25,7 +25,7 @@ */ /* - * Hyperbolic tangent of a complex argument z = x + i y. + * Hyperbolic tangent of a complex argument z = x + I y. * * The algorithm is from: * @@ -44,15 +44,15 @@ * * tanh(z) = sinh(z) / cosh(z) * - * sinh(x) cos(y) + i cosh(x) sin(y) + * sinh(x) cos(y) + I cosh(x) sin(y) * = --------------------------------- - * cosh(x) cos(y) + i sinh(x) sin(y) + * cosh(x) cos(y) + I sinh(x) sin(y) * - * cosh(x) sinh(x) / cos^2(y) + i tan(y) + * cosh(x) sinh(x) / cos^2(y) + I tan(y) * = ------------------------------------- * 1 + sinh^2(x) / cos^2(y) * - * beta rho s + i t + * beta rho s + I t * = ---------------- * 1 + beta s^2 * @@ -85,16 +85,16 @@ ctanh(double complex z) ix = hx & 0x7fffffff; /* - * ctanh(NaN + i 0) = NaN + i 0 + * ctanh(NaN +- I 0) = d(NaN) +- I 0 * - * ctanh(NaN + i y) = NaN + i NaN for y != 0 + * ctanh(NaN + I y) = d(NaN,y) + I d(NaN,y) for y != 0 * * The imaginary part has the sign of x*sin(2*y), but there's no * special effort to get this right. * - * ctanh(+-Inf +- i Inf) = +-1 +- 0 + * ctanh(+-Inf +- I Inf) = +-1 +- I 0 * - * ctanh(+-Inf + i y) = +-1 + 0 sin(2y) for y finite + * ctanh(+-Inf + I y) = +-1 + I 0 sin(2y) for y finite * * The imaginary part of the sign is unspecified. This special * case is only needed to avoid a spurious invalid exception when @@ -102,24 +102,25 @@ ctanh(double complex z) */ if (ix >= 0x7ff00000) { if ((ix & 0xfffff) | lx) /* x is NaN */ - return (CMPLX(x, (y == 0 ? y : x * y))); + return (CMPLX((x + 0) * (y + 0), + y == 0 ? y : (x + 0) * (y + 0))); SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */ return (CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)))); } /* - * ctanh(x + i NAN) = NaN + i NaN - * ctanh(x +- i Inf) = NaN + i NaN + * ctanh(x + I NaN) = d(NaN) + I d(NaN) + * ctanh(x +- I Inf) = dNaN + I dNaN */ if (!isfinite(y)) return (CMPLX(y - y, y - y)); /* - * ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the + * ctanh(+-huge +- I y) ~= +-1 +- I 2sin(2y)/exp(2x), using the * approximation sinh^2(huge) ~= exp(2*huge) / 4. * We use a modified formula to avoid spurious overflow. */ - if (ix >= 0x40360000) { /* x >= 22 */ + if (ix >= 0x40360000) { /* |x| >= 22 */ double exp_mx = exp(-fabs(x)); return (CMPLX(copysign(1, x), 4 * sin(y) * cos(y) * exp_mx * exp_mx)); @@ -138,7 +139,7 @@ double complex ctan(double complex z) { - /* ctan(z) = -I * ctanh(I * z) */ - z = ctanh(CMPLX(-cimag(z), creal(z))); - return (CMPLX(cimag(z), -creal(z))); + /* ctan(z) = -I * ctanh(I * z) = I * conj(ctanh(I * conj(z))) */ + z = ctanh(CMPLX(cimag(z), creal(z))); + return (CMPLX(cimag(z), creal(z))); } diff --git a/lib/msun/src/s_ctanhf.c b/lib/msun/src/s_ctanhf.c index 9b0cda1..520bf77 100644 --- a/lib/msun/src/s_ctanhf.c +++ b/lib/msun/src/s_ctanhf.c @@ -51,7 +51,8 @@ ctanhf(float complex z) if (ix >= 0x7f800000) { if (ix & 0x7fffff) - return (CMPLXF(x, (y == 0 ? y : x * y))); + return (CMPLXF((x + 0) * (y + 0), + y == 0 ? y : (x + 0) * (y + 0))); SET_FLOAT_WORD(x, hx - 0x40000000); return (CMPLXF(x, copysignf(0, isinf(y) ? y : sinf(y) * cosf(y)))); @@ -60,7 +61,7 @@ ctanhf(float complex z) if (!isfinite(y)) return (CMPLXF(y - y, y - y)); - if (ix >= 0x41300000) { /* x >= 11 */ + if (ix >= 0x41300000) { /* |x| >= 11 */ float exp_mx = expf(-fabsf(x)); return (CMPLXF(copysignf(1, x), 4 * sinf(y) * cosf(y) * exp_mx * exp_mx)); @@ -78,7 +79,7 @@ float complex ctanf(float complex z) { - z = ctanhf(CMPLXF(-cimagf(z), crealf(z))); - return (CMPLXF(cimagf(z), -crealf(z))); + z = ctanhf(CMPLXF(cimagf(z), crealf(z))); + return (CMPLXF(cimagf(z), crealf(z))); } |
