diff --git a/lib/msun/src/s_ccoshf.c b/lib/msun/src/s_ccoshf.c index e72395c277d5..5d7a09ba5f8d 100644 --- a/lib/msun/src/s_ccoshf.c +++ b/lib/msun/src/s_ccoshf.c @@ -1,104 +1,104 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice unmodified, this list of conditions, and the following * disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ /* * Float version of ccosh(). See s_ccosh.c for details. */ #include __FBSDID("$FreeBSD$"); #include #include #include "math_private.h" static const float huge = 0x1p127; float complex ccoshf(float complex z) { - volatile float x, y, h; + float x, y, h; int32_t hx, hy, ix, iy; x = crealf(z); y = cimagf(z); GET_FLOAT_WORD(hx, x); GET_FLOAT_WORD(hy, y); ix = 0x7fffffff & hx; iy = 0x7fffffff & hy; if (ix < 0x7f800000 && iy < 0x7f800000) { if (iy == 0) return (CMPLXF(coshf(x), x * y)); 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; return (CMPLXF(h * cosf(y), copysignf(h, x) * sinf(y))); } else if (ix < 0x4340b1e7) { /* x < 192.7: scale to avoid overflow */ 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 (CMPLXF(h * h * cosf(y), h * sinf(y))); } } if (ix == 0) /* && iy >= 0x7f800000 */ return (CMPLXF(y - y, x * copysignf(0, y))); if (iy == 0) /* && ix >= 0x7f800000 */ return (CMPLXF(x * x, copysignf(0, x) * y)); if (ix < 0x7f800000) /* && iy >= 0x7f800000 */ return (CMPLXF(y - y, x * (y - y))); if (ix == 0x7f800000) { if (iy >= 0x7f800000) return (CMPLXF(INFINITY, x * (y - y))); return (CMPLXF(INFINITY * cosf(y), x * sinf(y))); } return (CMPLXF(((long double)x * x) * (y - y), ((long double)x + x) * (y - y))); } float complex ccosf(float complex z) { return (ccoshf(CMPLXF(-cimagf(z), crealf(z)))); } diff --git a/lib/msun/src/s_ctanh.c b/lib/msun/src/s_ctanh.c index 93e5ad444501..e5840a1bf67b 100644 --- a/lib/msun/src/s_ctanh.c +++ b/lib/msun/src/s_ctanh.c @@ -1,149 +1,149 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2011 David Schultz * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice unmodified, this list of conditions, and the following * disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ /* * Hyperbolic tangent of a complex argument z = x + I y. * * The algorithm is from: * * W. Kahan. Branch Cuts for Complex Elementary Functions or Much * Ado About Nothing's Sign Bit. In The State of the Art in * Numerical Analysis, pp. 165 ff. Iserles and Powell, eds., 1987. * * Method: * * Let t = tan(x) * beta = 1/cos^2(y) * s = sinh(x) * rho = cosh(x) * * We have: * * tanh(z) = sinh(z) / cosh(z) * * sinh(x) cos(y) + I cosh(x) sin(y) * = --------------------------------- * cosh(x) cos(y) + I sinh(x) sin(y) * * cosh(x) sinh(x) / cos^2(y) + I tan(y) * = ------------------------------------- * 1 + sinh^2(x) / cos^2(y) * * beta rho s + I t * = ---------------- * 1 + beta s^2 * * Modifications: * * I omitted the original algorithm's handling of overflow in tan(x) after * verifying with nearpi.c that this can't happen in IEEE single or double * precision. I also handle large x differently. */ #include __FBSDID("$FreeBSD$"); #include #include #include "math_private.h" double complex ctanh(double complex z) { - volatile double x, y; + double x, y; double t, beta, s, rho, denom; uint32_t hx, ix, lx; x = creal(z); y = cimag(z); EXTRACT_WORDS(hx, lx, x); ix = hx & 0x7fffffff; /* * ctanh(NaN +- I 0) = d(NaN) +- I 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 +- I 0 * * 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 * y is infinite. */ if (ix >= 0x7ff00000) { if ((ix & 0xfffff) | lx) /* x is NaN */ 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)))); } /* * ctanh(+-0 + i NAN) = +-0 + i NaN * ctanh(+-0 +- i Inf) = +-0 + i NaN * ctanh(x + i NAN) = NaN + i NaN * ctanh(x +- i Inf) = NaN + i NaN */ if (!isfinite(y)) return (CMPLX(x ? y - y : x, y - y)); /* * 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 */ double exp_mx = exp(-fabs(x)); return (CMPLX(copysign(1, x), 4 * sin(y) * cos(y) * exp_mx * exp_mx)); } /* Kahan's algorithm */ t = tan(y); beta = 1.0 + t * t; /* = 1 / cos^2(y) */ s = sinh(x); rho = sqrt(1 + s * s); /* = cosh(x) */ denom = 1 + beta * s * s; return (CMPLX((beta * rho * s) / denom, t / denom)); } double complex ctan(double complex 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 164a2c23df9e..c46f86d2e116 100644 --- a/lib/msun/src/s_ctanhf.c +++ b/lib/msun/src/s_ctanhf.c @@ -1,87 +1,87 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2011 David Schultz * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice unmodified, this list of conditions, and the following * disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ /* * Hyperbolic tangent of a complex argument z. See s_ctanh.c for details. */ #include __FBSDID("$FreeBSD$"); #include #include #include "math_private.h" float complex ctanhf(float complex z) { - volatile float x, y; + float x, y; float t, beta, s, rho, denom; uint32_t hx, ix; x = crealf(z); y = cimagf(z); GET_FLOAT_WORD(hx, x); ix = hx & 0x7fffffff; if (ix >= 0x7f800000) { if (ix & 0x7fffff) 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)))); } if (!isfinite(y)) return (CMPLXF(ix ? y - y : x, y - y)); 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)); } t = tanf(y); beta = 1.0 + t * t; s = sinhf(x); rho = sqrtf(1 + s * s); denom = 1 + beta * s * s; return (CMPLXF((beta * rho * s) / denom, t / denom)); } float complex ctanf(float complex z) { z = ctanhf(CMPLXF(cimagf(z), crealf(z))); return (CMPLXF(cimagf(z), crealf(z))); }