13daee1d6SDavid Schultz /*- 23daee1d6SDavid Schultz * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl 33daee1d6SDavid Schultz * All rights reserved. 43daee1d6SDavid Schultz * 53daee1d6SDavid Schultz * Redistribution and use in source and binary forms, with or without 63daee1d6SDavid Schultz * modification, are permitted provided that the following conditions 73daee1d6SDavid Schultz * are met: 83daee1d6SDavid Schultz * 1. Redistributions of source code must retain the above copyright 93daee1d6SDavid Schultz * notice unmodified, this list of conditions, and the following 103daee1d6SDavid Schultz * disclaimer. 113daee1d6SDavid Schultz * 2. Redistributions in binary form must reproduce the above copyright 123daee1d6SDavid Schultz * notice, this list of conditions and the following disclaimer in the 133daee1d6SDavid Schultz * documentation and/or other materials provided with the distribution. 143daee1d6SDavid Schultz * 153daee1d6SDavid Schultz * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR 163daee1d6SDavid Schultz * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 173daee1d6SDavid Schultz * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 183daee1d6SDavid Schultz * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, 193daee1d6SDavid Schultz * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 203daee1d6SDavid Schultz * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 213daee1d6SDavid Schultz * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 223daee1d6SDavid Schultz * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 233daee1d6SDavid Schultz * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 243daee1d6SDavid Schultz * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 253daee1d6SDavid Schultz */ 263daee1d6SDavid Schultz 273daee1d6SDavid Schultz /* 28*db7548d0STijl Coosemans * Float version of ccosh(). See s_ccosh.c for details. 293daee1d6SDavid Schultz */ 303daee1d6SDavid Schultz 313daee1d6SDavid Schultz #include <sys/cdefs.h> 323daee1d6SDavid Schultz __FBSDID("$FreeBSD$"); 333daee1d6SDavid Schultz 343daee1d6SDavid Schultz #include <complex.h> 353daee1d6SDavid Schultz #include <math.h> 363daee1d6SDavid Schultz 373daee1d6SDavid Schultz #include "math_private.h" 383daee1d6SDavid Schultz 39c6df46baSDavid Schultz static const float huge = 0x1p127; 40c6df46baSDavid Schultz 413daee1d6SDavid Schultz float complex 423daee1d6SDavid Schultz ccoshf(float complex z) 433daee1d6SDavid Schultz { 44c6df46baSDavid Schultz float x, y, h; 453daee1d6SDavid Schultz int32_t hx, hy, ix, iy; 463daee1d6SDavid Schultz 473daee1d6SDavid Schultz x = crealf(z); 483daee1d6SDavid Schultz y = cimagf(z); 493daee1d6SDavid Schultz 503daee1d6SDavid Schultz GET_FLOAT_WORD(hx, x); 513daee1d6SDavid Schultz GET_FLOAT_WORD(hy, y); 523daee1d6SDavid Schultz 533daee1d6SDavid Schultz ix = 0x7fffffff & hx; 543daee1d6SDavid Schultz iy = 0x7fffffff & hy; 553daee1d6SDavid Schultz 563daee1d6SDavid Schultz if (ix < 0x7f800000 && iy < 0x7f800000) { 573daee1d6SDavid Schultz if (iy == 0) 582cec876aSEd Schouten return (CMPLXF(coshf(x), x * y)); 59*db7548d0STijl Coosemans if (ix < 0x41100000) /* |x| < 9: normal case */ 602cec876aSEd Schouten return (CMPLXF(coshf(x) * cosf(y), sinhf(x) * sinf(y))); 61c6df46baSDavid Schultz 62c6df46baSDavid Schultz /* |x| >= 9, so cosh(x) ~= exp(|x|) */ 63c6df46baSDavid Schultz if (ix < 0x42b17218) { 64c6df46baSDavid Schultz /* x < 88.7: expf(|x|) won't overflow */ 65*db7548d0STijl Coosemans h = expf(fabsf(x)) * 0.5F; 662cec876aSEd Schouten return (CMPLXF(h * cosf(y), copysignf(h, x) * sinf(y))); 67c6df46baSDavid Schultz } else if (ix < 0x4340b1e7) { 68c6df46baSDavid Schultz /* x < 192.7: scale to avoid overflow */ 692cec876aSEd Schouten z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1); 702cec876aSEd Schouten return (CMPLXF(crealf(z), cimagf(z) * copysignf(1, x))); 71c6df46baSDavid Schultz } else { 72c6df46baSDavid Schultz /* x >= 192.7: the result always overflows */ 73c6df46baSDavid Schultz h = huge * x; 742cec876aSEd Schouten return (CMPLXF(h * h * cosf(y), h * sinf(y))); 75c6df46baSDavid Schultz } 763daee1d6SDavid Schultz } 773daee1d6SDavid Schultz 78*db7548d0STijl Coosemans if (ix == 0) /* && iy >= 0x7f800000 */ 79*db7548d0STijl Coosemans return (CMPLXF(y - y, x * copysignf(0, y))); 803daee1d6SDavid Schultz 81*db7548d0STijl Coosemans if (iy == 0) /* && ix >= 0x7f800000 */ 822cec876aSEd Schouten return (CMPLXF(x * x, copysignf(0, x) * y)); 833daee1d6SDavid Schultz 84*db7548d0STijl Coosemans if (ix < 0x7f800000) /* && iy >= 0x7f800000 */ 852cec876aSEd Schouten return (CMPLXF(y - y, x * (y - y))); 863daee1d6SDavid Schultz 87*db7548d0STijl Coosemans if (ix == 0x7f800000) { 883daee1d6SDavid Schultz if (iy >= 0x7f800000) 89*db7548d0STijl Coosemans return (CMPLXF(INFINITY, x * (y - y))); 90*db7548d0STijl Coosemans return (CMPLXF(INFINITY * cosf(y), x * sinf(y))); 913daee1d6SDavid Schultz } 923daee1d6SDavid Schultz 932cec876aSEd Schouten return (CMPLXF((x * x) * (y - y), (x + x) * (y - y))); 943daee1d6SDavid Schultz } 953daee1d6SDavid Schultz 963daee1d6SDavid Schultz float complex 973daee1d6SDavid Schultz ccosf(float complex z) 983daee1d6SDavid Schultz { 993daee1d6SDavid Schultz 1002cec876aSEd Schouten return (ccoshf(CMPLXF(-cimagf(z), crealf(z)))); 1013daee1d6SDavid Schultz } 102