13daee1d6SDavid Schultz /*-
2*4d846d26SWarner Losh * SPDX-License-Identifier: BSD-2-Clause
35e53a4f9SPedro F. Giffuni *
43daee1d6SDavid Schultz * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl
53daee1d6SDavid Schultz * All rights reserved.
63daee1d6SDavid Schultz *
73daee1d6SDavid Schultz * Redistribution and use in source and binary forms, with or without
83daee1d6SDavid Schultz * modification, are permitted provided that the following conditions
93daee1d6SDavid Schultz * are met:
103daee1d6SDavid Schultz * 1. Redistributions of source code must retain the above copyright
113daee1d6SDavid Schultz * notice unmodified, this list of conditions, and the following
123daee1d6SDavid Schultz * disclaimer.
133daee1d6SDavid Schultz * 2. Redistributions in binary form must reproduce the above copyright
143daee1d6SDavid Schultz * notice, this list of conditions and the following disclaimer in the
153daee1d6SDavid Schultz * documentation and/or other materials provided with the distribution.
163daee1d6SDavid Schultz *
173daee1d6SDavid Schultz * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
183daee1d6SDavid Schultz * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
193daee1d6SDavid Schultz * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
203daee1d6SDavid Schultz * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
213daee1d6SDavid Schultz * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
223daee1d6SDavid Schultz * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
233daee1d6SDavid Schultz * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
243daee1d6SDavid Schultz * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
253daee1d6SDavid Schultz * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
263daee1d6SDavid Schultz * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
273daee1d6SDavid Schultz */
283daee1d6SDavid Schultz
293daee1d6SDavid Schultz /*
30db7548d0STijl Coosemans * Float version of ccosh(). See s_ccosh.c for details.
313daee1d6SDavid Schultz */
323daee1d6SDavid Schultz
333daee1d6SDavid Schultz #include <complex.h>
343daee1d6SDavid Schultz #include <math.h>
353daee1d6SDavid Schultz
363daee1d6SDavid Schultz #include "math_private.h"
373daee1d6SDavid Schultz
38c6df46baSDavid Schultz static const float huge = 0x1p127;
39c6df46baSDavid Schultz
403daee1d6SDavid Schultz float complex
ccoshf(float complex z)413daee1d6SDavid Schultz ccoshf(float complex z)
423daee1d6SDavid Schultz {
43e2157cd0SDimitry Andric float x, y, h;
443daee1d6SDavid Schultz int32_t hx, hy, ix, iy;
453daee1d6SDavid Schultz
463daee1d6SDavid Schultz x = crealf(z);
473daee1d6SDavid Schultz y = cimagf(z);
483daee1d6SDavid Schultz
493daee1d6SDavid Schultz GET_FLOAT_WORD(hx, x);
503daee1d6SDavid Schultz GET_FLOAT_WORD(hy, y);
513daee1d6SDavid Schultz
523daee1d6SDavid Schultz ix = 0x7fffffff & hx;
533daee1d6SDavid Schultz iy = 0x7fffffff & hy;
543daee1d6SDavid Schultz
553daee1d6SDavid Schultz if (ix < 0x7f800000 && iy < 0x7f800000) {
563daee1d6SDavid Schultz if (iy == 0)
572cec876aSEd Schouten return (CMPLXF(coshf(x), x * y));
58db7548d0STijl Coosemans if (ix < 0x41100000) /* |x| < 9: normal case */
592cec876aSEd Schouten return (CMPLXF(coshf(x) * cosf(y), sinhf(x) * sinf(y)));
60c6df46baSDavid Schultz
61c6df46baSDavid Schultz /* |x| >= 9, so cosh(x) ~= exp(|x|) */
62c6df46baSDavid Schultz if (ix < 0x42b17218) {
63c6df46baSDavid Schultz /* x < 88.7: expf(|x|) won't overflow */
64db7548d0STijl Coosemans h = expf(fabsf(x)) * 0.5F;
652cec876aSEd Schouten return (CMPLXF(h * cosf(y), copysignf(h, x) * sinf(y)));
66c6df46baSDavid Schultz } else if (ix < 0x4340b1e7) {
67c6df46baSDavid Schultz /* x < 192.7: scale to avoid overflow */
682cec876aSEd Schouten z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1);
692cec876aSEd Schouten return (CMPLXF(crealf(z), cimagf(z) * copysignf(1, x)));
70c6df46baSDavid Schultz } else {
71c6df46baSDavid Schultz /* x >= 192.7: the result always overflows */
72c6df46baSDavid Schultz h = huge * x;
732cec876aSEd Schouten return (CMPLXF(h * h * cosf(y), h * sinf(y)));
74c6df46baSDavid Schultz }
753daee1d6SDavid Schultz }
763daee1d6SDavid Schultz
77db7548d0STijl Coosemans if (ix == 0) /* && iy >= 0x7f800000 */
78db7548d0STijl Coosemans return (CMPLXF(y - y, x * copysignf(0, y)));
793daee1d6SDavid Schultz
80db7548d0STijl Coosemans if (iy == 0) /* && ix >= 0x7f800000 */
812cec876aSEd Schouten return (CMPLXF(x * x, copysignf(0, x) * y));
823daee1d6SDavid Schultz
83db7548d0STijl Coosemans if (ix < 0x7f800000) /* && iy >= 0x7f800000 */
842cec876aSEd Schouten return (CMPLXF(y - y, x * (y - y)));
853daee1d6SDavid Schultz
86db7548d0STijl Coosemans if (ix == 0x7f800000) {
873daee1d6SDavid Schultz if (iy >= 0x7f800000)
88db7548d0STijl Coosemans return (CMPLXF(INFINITY, x * (y - y)));
89db7548d0STijl Coosemans return (CMPLXF(INFINITY * cosf(y), x * sinf(y)));
903daee1d6SDavid Schultz }
913daee1d6SDavid Schultz
926f1b8a07SBruce Evans return (CMPLXF(((long double)x * x) * (y - y),
936f1b8a07SBruce Evans ((long double)x + x) * (y - y)));
943daee1d6SDavid Schultz }
953daee1d6SDavid Schultz
963daee1d6SDavid Schultz float complex
ccosf(float complex z)973daee1d6SDavid Schultz ccosf(float complex z)
983daee1d6SDavid Schultz {
993daee1d6SDavid Schultz
1002cec876aSEd Schouten return (ccoshf(CMPLXF(-cimagf(z), crealf(z))));
1013daee1d6SDavid Schultz }
102