Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 1 | /* origin: FreeBSD /usr/src/lib/msun/src/s_ccosh.c */ |
| 2 | /*- |
| 3 | * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl |
| 4 | * All rights reserved. |
| 5 | * |
| 6 | * Redistribution and use in source and binary forms, with or without |
| 7 | * modification, are permitted provided that the following conditions |
| 8 | * are met: |
| 9 | * 1. Redistributions of source code must retain the above copyright |
| 10 | * notice unmodified, this list of conditions, and the following |
| 11 | * disclaimer. |
| 12 | * 2. Redistributions in binary form must reproduce the above copyright |
| 13 | * notice, this list of conditions and the following disclaimer in the |
| 14 | * documentation and/or other materials provided with the distribution. |
| 15 | * |
| 16 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR |
| 17 | * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES |
| 18 | * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. |
| 19 | * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, |
| 20 | * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT |
| 21 | * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, |
| 22 | * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY |
| 23 | * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT |
| 24 | * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF |
| 25 | * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
| 26 | */ |
| 27 | /* |
| 28 | * Hyperbolic cosine of a complex argument z = x + i y. |
| 29 | * |
| 30 | * cosh(z) = cosh(x+iy) |
| 31 | * = cosh(x) cos(y) + i sinh(x) sin(y). |
| 32 | * |
| 33 | * Exceptional values are noted in the comments within the source code. |
| 34 | * These values and the return value were taken from n1124.pdf. |
| 35 | */ |
| 36 | |
Szabolcs Nagy | 2d72b58 | 2018-11-29 22:09:53 +0000 | [diff] [blame^] | 37 | #include "complex_impl.h" |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 38 | |
| 39 | static const double huge = 0x1p1023; |
| 40 | |
| 41 | double complex ccosh(double complex z) |
| 42 | { |
| 43 | double x, y, h; |
| 44 | int32_t hx, hy, ix, iy, lx, ly; |
| 45 | |
| 46 | x = creal(z); |
| 47 | y = cimag(z); |
| 48 | |
| 49 | EXTRACT_WORDS(hx, lx, x); |
| 50 | EXTRACT_WORDS(hy, ly, y); |
| 51 | |
| 52 | ix = 0x7fffffff & hx; |
| 53 | iy = 0x7fffffff & hy; |
| 54 | |
| 55 | /* Handle the nearly-non-exceptional cases where x and y are finite. */ |
| 56 | if (ix < 0x7ff00000 && iy < 0x7ff00000) { |
| 57 | if ((iy | ly) == 0) |
Szabolcs Nagy | cfbaba7 | 2012-11-13 01:31:49 +0100 | [diff] [blame] | 58 | return CMPLX(cosh(x), x * y); |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 59 | if (ix < 0x40360000) /* small x: normal case */ |
Szabolcs Nagy | cfbaba7 | 2012-11-13 01:31:49 +0100 | [diff] [blame] | 60 | return CMPLX(cosh(x) * cos(y), sinh(x) * sin(y)); |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 61 | |
| 62 | /* |x| >= 22, so cosh(x) ~= exp(|x|) */ |
| 63 | if (ix < 0x40862e42) { |
| 64 | /* x < 710: exp(|x|) won't overflow */ |
| 65 | h = exp(fabs(x)) * 0.5; |
Szabolcs Nagy | cfbaba7 | 2012-11-13 01:31:49 +0100 | [diff] [blame] | 66 | return CMPLX(h * cos(y), copysign(h, x) * sin(y)); |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 67 | } else if (ix < 0x4096bbaa) { |
| 68 | /* x < 1455: scale to avoid overflow */ |
Szabolcs Nagy | cfbaba7 | 2012-11-13 01:31:49 +0100 | [diff] [blame] | 69 | z = __ldexp_cexp(CMPLX(fabs(x), y), -1); |
| 70 | return CMPLX(creal(z), cimag(z) * copysign(1, x)); |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 71 | } else { |
| 72 | /* x >= 1455: the result always overflows */ |
| 73 | h = huge * x; |
Szabolcs Nagy | cfbaba7 | 2012-11-13 01:31:49 +0100 | [diff] [blame] | 74 | return CMPLX(h * h * cos(y), h * sin(y)); |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 75 | } |
| 76 | } |
| 77 | |
| 78 | /* |
| 79 | * cosh(+-0 +- I Inf) = dNaN + I sign(d(+-0, dNaN))0. |
| 80 | * The sign of 0 in the result is unspecified. Choice = normally |
| 81 | * the same as dNaN. Raise the invalid floating-point exception. |
| 82 | * |
| 83 | * cosh(+-0 +- I NaN) = d(NaN) + I sign(d(+-0, NaN))0. |
| 84 | * The sign of 0 in the result is unspecified. Choice = normally |
| 85 | * the same as d(NaN). |
| 86 | */ |
| 87 | if ((ix | lx) == 0 && iy >= 0x7ff00000) |
Szabolcs Nagy | cfbaba7 | 2012-11-13 01:31:49 +0100 | [diff] [blame] | 88 | return CMPLX(y - y, copysign(0, x * (y - y))); |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 89 | |
| 90 | /* |
| 91 | * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0. |
| 92 | * |
| 93 | * cosh(NaN +- I 0) = d(NaN) + I sign(d(NaN, +-0))0. |
| 94 | * The sign of 0 in the result is unspecified. |
| 95 | */ |
| 96 | if ((iy | ly) == 0 && ix >= 0x7ff00000) { |
| 97 | if (((hx & 0xfffff) | lx) == 0) |
Szabolcs Nagy | cfbaba7 | 2012-11-13 01:31:49 +0100 | [diff] [blame] | 98 | return CMPLX(x * x, copysign(0, x) * y); |
| 99 | return CMPLX(x * x, copysign(0, (x + x) * y)); |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 100 | } |
| 101 | |
| 102 | /* |
| 103 | * cosh(x +- I Inf) = dNaN + I dNaN. |
| 104 | * Raise the invalid floating-point exception for finite nonzero x. |
| 105 | * |
| 106 | * cosh(x + I NaN) = d(NaN) + I d(NaN). |
| 107 | * Optionally raises the invalid floating-point exception for finite |
| 108 | * nonzero x. Choice = don't raise (except for signaling NaNs). |
| 109 | */ |
| 110 | if (ix < 0x7ff00000 && iy >= 0x7ff00000) |
Szabolcs Nagy | cfbaba7 | 2012-11-13 01:31:49 +0100 | [diff] [blame] | 111 | return CMPLX(y - y, x * (y - y)); |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 112 | |
| 113 | /* |
| 114 | * cosh(+-Inf + I NaN) = +Inf + I d(NaN). |
| 115 | * |
| 116 | * cosh(+-Inf +- I Inf) = +Inf + I dNaN. |
| 117 | * The sign of Inf in the result is unspecified. Choice = always +. |
| 118 | * Raise the invalid floating-point exception. |
| 119 | * |
| 120 | * cosh(+-Inf + I y) = +Inf cos(y) +- I Inf sin(y) |
| 121 | */ |
| 122 | if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) { |
| 123 | if (iy >= 0x7ff00000) |
Szabolcs Nagy | cfbaba7 | 2012-11-13 01:31:49 +0100 | [diff] [blame] | 124 | return CMPLX(x * x, x * (y - y)); |
| 125 | return CMPLX((x * x) * cos(y), x * sin(y)); |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 126 | } |
| 127 | |
| 128 | /* |
| 129 | * cosh(NaN + I NaN) = d(NaN) + I d(NaN). |
| 130 | * |
| 131 | * cosh(NaN +- I Inf) = d(NaN) + I d(NaN). |
| 132 | * Optionally raises the invalid floating-point exception. |
| 133 | * Choice = raise. |
| 134 | * |
| 135 | * cosh(NaN + I y) = d(NaN) + I d(NaN). |
| 136 | * Optionally raises the invalid floating-point exception for finite |
| 137 | * nonzero y. Choice = don't raise (except for signaling NaNs). |
| 138 | */ |
Szabolcs Nagy | cfbaba7 | 2012-11-13 01:31:49 +0100 | [diff] [blame] | 139 | return CMPLX((x * x) * (y - y), (x + x) * (y - y)); |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 140 | } |