1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/complex/ccoshf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/complex/ccosh.c
6
7const std = @import("../../std.zig");
8const testing = std.testing;
9const math = std.math;
10const cmath = math.complex;
11const Complex = cmath.Complex;
12
13const ldexp_cexp = @import("ldexp.zig").ldexp_cexp;
14
15/// Returns the hyperbolic arc-cosine of z.
16pub fn cosh(z: anytype) Complex(@TypeOf(z.re, z.im)) {
17 const T = @TypeOf(z.re, z.im);
18 return switch (T) {
19 f32 => cosh32(z),
20 f64 => cosh64(z),
21 else => @compileError("cosh not implemented for " ++ @typeName(z)),
22 };
23}
24
25fn cosh32(z: Complex(f32)) Complex(f32) {
26 const x = z.re;
27 const y = z.im;
28
29 const hx: u32 = @bitCast(x);
30 const ix = hx & 0x7fffffff;
31
32 const hy: u32 = @bitCast(y);
33 const iy = hy & 0x7fffffff;
34
35 if (ix < 0x7f800000 and iy < 0x7f800000) {
36 if (iy == 0) {
37 return Complex(f32).init(math.cosh(x), x * y);
38 }
39 // small x: normal case
40 if (ix < 0x41100000) {
41 return Complex(f32).init(math.cosh(x) * @cos(y), math.sinh(x) * @sin(y));
42 }
43
44 // |x|>= 9, so cosh(x) ~= exp(|x|)
45 if (ix < 0x42b17218) {
46 // x < 88.7: exp(|x|) won't overflow
47 const h = @exp(@abs(x)) * 0.5;
48 return Complex(f32).init(h * @cos(y), math.copysign(h, x) * @sin(y));
49 }
50 // x < 192.7: scale to avoid overflow
51 else if (ix < 0x4340b1e7) {
52 const v = Complex(f32).init(@abs(x), y);
53 const r = ldexp_cexp(v, -1);
54 return Complex(f32).init(r.re, r.im * math.copysign(@as(f32, 1.0), x));
55 }
56 // x >= 192.7: result always overflows
57 else {
58 const h = 0x1p127 * x;
59 return Complex(f32).init(h * h * @cos(y), h * @sin(y));
60 }
61 }
62
63 if (ix == 0 and iy >= 0x7f800000) {
64 return Complex(f32).init(y - y, math.copysign(@as(f32, 0.0), x * (y - y)));
65 }
66
67 if (iy == 0 and ix >= 0x7f800000) {
68 if (hx & 0x7fffff == 0) {
69 return Complex(f32).init(x * x, math.copysign(@as(f32, 0.0), x) * y);
70 }
71 return Complex(f32).init(x * x, math.copysign(@as(f32, 0.0), (x + x) * y));
72 }
73
74 if (ix < 0x7f800000 and iy >= 0x7f800000) {
75 return Complex(f32).init(y - y, x * (y - y));
76 }
77
78 if (ix >= 0x7f800000 and (hx & 0x7fffff) == 0) {
79 if (iy >= 0x7f800000) {
80 return Complex(f32).init(x * x, x * (y - y));
81 }
82 return Complex(f32).init((x * x) * @cos(y), x * @sin(y));
83 }
84
85 return Complex(f32).init((x * x) * (y - y), (x + x) * (y - y));
86}
87
88fn cosh64(z: Complex(f64)) Complex(f64) {
89 const x = z.re;
90 const y = z.im;
91
92 const fx: u64 = @bitCast(x);
93 const hx: u32 = @intCast(fx >> 32);
94 const lx: u32 = @truncate(fx);
95 const ix = hx & 0x7fffffff;
96
97 const fy: u64 = @bitCast(y);
98 const hy: u32 = @intCast(fy >> 32);
99 const ly: u32 = @truncate(fy);
100 const iy = hy & 0x7fffffff;
101
102 // nearly non-exceptional case where x, y are finite
103 if (ix < 0x7ff00000 and iy < 0x7ff00000) {
104 if (iy | ly == 0) {
105 return Complex(f64).init(math.cosh(x), x * y);
106 }
107 // small x: normal case
108 if (ix < 0x40360000) {
109 return Complex(f64).init(math.cosh(x) * @cos(y), math.sinh(x) * @sin(y));
110 }
111
112 // |x|>= 22, so cosh(x) ~= exp(|x|)
113 if (ix < 0x40862e42) {
114 // x < 710: exp(|x|) won't overflow
115 const h = @exp(@abs(x)) * 0.5;
116 return Complex(f64).init(h * @cos(y), math.copysign(h, x) * @sin(y));
117 }
118 // x < 1455: scale to avoid overflow
119 else if (ix < 0x4096bbaa) {
120 const v = Complex(f64).init(@abs(x), y);
121 const r = ldexp_cexp(v, -1);
122 return Complex(f64).init(r.re, r.im * math.copysign(@as(f64, 1.0), x));
123 }
124 // x >= 1455: result always overflows
125 else {
126 const h = 0x1p1023 * x;
127 return Complex(f64).init(h * h * @cos(y), h * @sin(y));
128 }
129 }
130
131 if (ix | lx == 0 and iy >= 0x7ff00000) {
132 return Complex(f64).init(y - y, math.copysign(@as(f64, 0.0), x * (y - y)));
133 }
134
135 if (iy | ly == 0 and ix >= 0x7ff00000) {
136 if ((hx & 0xfffff) | lx == 0) {
137 return Complex(f64).init(x * x, math.copysign(@as(f64, 0.0), x) * y);
138 }
139 return Complex(f64).init(x * x, math.copysign(@as(f64, 0.0), (x + x) * y));
140 }
141
142 if (ix < 0x7ff00000 and iy >= 0x7ff00000) {
143 return Complex(f64).init(y - y, x * (y - y));
144 }
145
146 if (ix >= 0x7ff00000 and (hx & 0xfffff) | lx == 0) {
147 if (iy >= 0x7ff00000) {
148 return Complex(f64).init(x * x, x * (y - y));
149 }
150 return Complex(f64).init(x * x * @cos(y), x * @sin(y));
151 }
152
153 return Complex(f64).init((x * x) * (y - y), (x + x) * (y - y));
154}
155
156test cosh32 {
157 const epsilon = math.floatEps(f32);
158 const a = Complex(f32).init(5, 3);
159 const c = cosh(a);
160
161 try testing.expectApproxEqAbs(-73.467300, c.re, epsilon);
162 try testing.expectApproxEqAbs(10.471557, c.im, epsilon);
163}
164
165test cosh64 {
166 const epsilon = math.floatEps(f64);
167 const a = Complex(f64).init(5, 3);
168 const c = cosh(a);
169
170 try testing.expectApproxEqAbs(-73.46729221264526, c.re, epsilon);
171 try testing.expectApproxEqAbs(10.471557674805572, c.im, epsilon);
172}
173
174test "cosh64 musl" {
175 const epsilon = math.floatEps(f64);
176 const a = Complex(f64).init(7.44648873421389e17, 1.6008058402057622e19);
177 const c = cosh(a);
178
179 try testing.expectApproxEqAbs(std.math.inf(f64), c.re, epsilon);
180 try testing.expectApproxEqAbs(std.math.inf(f64), c.im, epsilon);
181}