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/csinhf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/complex/csinh.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 sine of z.
16pub fn sinh(z: anytype) Complex(@TypeOf(z.re, z.im)) {
17 const T = @TypeOf(z.re, z.im);
18 return switch (T) {
19 f32 => sinh32(z),
20 f64 => sinh64(z),
21 else => @compileError("tan not implemented for " ++ @typeName(z)),
22 };
23}
24
25fn sinh32(z: Complex(f32)) Complex(f32) {
26 const x = z.re;
27 const y = z.im;
28
29 const hx = @as(u32, @bitCast(x));
30 const ix = hx & 0x7fffffff;
31
32 const hy = @as(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.sinh(x), y);
38 }
39 // small x: normal case
40 if (ix < 0x41100000) {
41 return Complex(f32).init(math.sinh(x) * @cos(y), math.cosh(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(math.copysign(h, x) * @cos(y), h * @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 * math.copysign(@as(f32, 1.0), x), r.im);
55 }
56 // x >= 192.7: result always overflows
57 else {
58 const h = 0x1p127 * x;
59 return Complex(f32).init(h * @cos(y), h * h * @sin(y));
60 }
61 }
62
63 if (ix == 0 and iy >= 0x7f800000) {
64 return Complex(f32).init(math.copysign(@as(f32, 0.0), x * (y - y)), y - y);
65 }
66
67 if (iy == 0 and ix >= 0x7f800000) {
68 if (hx & 0x7fffff == 0) {
69 return Complex(f32).init(x, y);
70 }
71 return Complex(f32).init(x, math.copysign(@as(f32, 0.0), 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 * @cos(y), math.inf(f32) * @sin(y));
83 }
84
85 return Complex(f32).init((x * x) * (y - y), (x + x) * (y - y));
86}
87
88fn sinh64(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 if (ix < 0x7ff00000 and iy < 0x7ff00000) {
103 if (iy | ly == 0) {
104 return Complex(f64).init(math.sinh(x), y);
105 }
106 // small x: normal case
107 if (ix < 0x40360000) {
108 return Complex(f64).init(math.sinh(x) * @cos(y), math.cosh(x) * @sin(y));
109 }
110
111 // |x|>= 22, so cosh(x) ~= exp(|x|)
112 if (ix < 0x40862e42) {
113 // x < 710: exp(|x|) won't overflow
114 const h = @exp(@abs(x)) * 0.5;
115 return Complex(f64).init(math.copysign(h, x) * @cos(y), h * @sin(y));
116 }
117 // x < 1455: scale to avoid overflow
118 else if (ix < 0x4096bbaa) {
119 const v = Complex(f64).init(@abs(x), y);
120 const r = ldexp_cexp(v, -1);
121 return Complex(f64).init(r.re * math.copysign(@as(f64, 1.0), x), r.im);
122 }
123 // x >= 1455: result always overflows
124 else {
125 const h = 0x1p1023 * x;
126 return Complex(f64).init(h * @cos(y), h * h * @sin(y));
127 }
128 }
129
130 if (ix | lx == 0 and iy >= 0x7ff00000) {
131 return Complex(f64).init(math.copysign(@as(f64, 0.0), x * (y - y)), y - y);
132 }
133
134 if (iy | ly == 0 and ix >= 0x7ff00000) {
135 if ((hx & 0xfffff) | lx == 0) {
136 return Complex(f64).init(x, y);
137 }
138 return Complex(f64).init(x, math.copysign(@as(f64, 0.0), y));
139 }
140
141 if (ix < 0x7ff00000 and iy >= 0x7ff00000) {
142 return Complex(f64).init(y - y, x * (y - y));
143 }
144
145 if (ix >= 0x7ff00000 and (hx & 0xfffff) | lx == 0) {
146 if (iy >= 0x7ff00000) {
147 return Complex(f64).init(x * x, x * (y - y));
148 }
149 return Complex(f64).init(x * @cos(y), math.inf(f64) * @sin(y));
150 }
151
152 return Complex(f64).init((x * x) * (y - y), (x + x) * (y - y));
153}
154
155test sinh32 {
156 const epsilon = math.floatEps(f32);
157 const a = Complex(f32).init(5, 3);
158 const c = sinh(a);
159
160 try testing.expectApproxEqAbs(-73.460617, c.re, epsilon);
161 try testing.expectApproxEqAbs(10.472508, c.im, epsilon);
162}
163
164test sinh64 {
165 const epsilon = math.floatEps(f64);
166 const a = Complex(f64).init(5, 3);
167 const c = sinh(a);
168
169 try testing.expectApproxEqAbs(-73.46062169567367, c.re, epsilon);
170 try testing.expectApproxEqAbs(10.472508533940392, c.im, epsilon);
171}