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/ctanhf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/complex/ctanh.c
6
7const std = @import("../../std.zig");
8const testing = std.testing;
9const math = std.math;
10const cmath = math.complex;
11const Complex = cmath.Complex;
12
13/// Returns the hyperbolic tangent of z.
14pub fn tanh(z: anytype) Complex(@TypeOf(z.re, z.im)) {
15 const T = @TypeOf(z.re, z.im);
16 return switch (T) {
17 f32 => tanh32(z),
18 f64 => tanh64(z),
19 else => @compileError("tan not implemented for " ++ @typeName(z)),
20 };
21}
22
23fn tanh32(z: Complex(f32)) Complex(f32) {
24 const x = z.re;
25 const y = z.im;
26
27 const hx = @as(u32, @bitCast(x));
28 const ix = hx & 0x7fffffff;
29
30 if (ix >= 0x7f800000) {
31 if (ix & 0x7fffff != 0) {
32 const r = if (y == 0) y else x * y;
33 return Complex(f32).init(x, r);
34 }
35 const xx = @as(f32, @bitCast(hx - 0x40000000));
36 const r = if (math.isInf(y)) y else @sin(y) * @cos(y);
37 return Complex(f32).init(xx, math.copysign(@as(f32, 0.0), r));
38 }
39
40 if (!math.isFinite(y)) {
41 const r = if (ix != 0) y - y else x;
42 return Complex(f32).init(r, y - y);
43 }
44
45 // x >= 11
46 if (ix >= 0x41300000) {
47 const exp_mx = @exp(-@abs(x));
48 return Complex(f32).init(math.copysign(@as(f32, 1.0), x), 4 * @sin(y) * @cos(y) * exp_mx * exp_mx);
49 }
50
51 // Kahan's algorithm
52 const t = @tan(y);
53 const beta = 1.0 + t * t;
54 const s = math.sinh(x);
55 const rho = @sqrt(1 + s * s);
56 const den = 1 + beta * s * s;
57
58 return Complex(f32).init((beta * rho * s) / den, t / den);
59}
60
61fn tanh64(z: Complex(f64)) Complex(f64) {
62 const x = z.re;
63 const y = z.im;
64
65 const fx: u64 = @bitCast(x);
66 // TODO: zig should allow this conversion implicitly because it can notice that the value necessarily
67 // fits in range.
68 const hx: u32 = @intCast(fx >> 32);
69 const lx: u32 = @truncate(fx);
70 const ix = hx & 0x7fffffff;
71
72 if (ix >= 0x7ff00000) {
73 if ((ix & 0xfffff) | lx != 0) {
74 const r = if (y == 0) y else x * y;
75 return Complex(f64).init(x, r);
76 }
77
78 const xx: f64 = @bitCast((@as(u64, hx - 0x40000000) << 32) | lx);
79 const r = if (math.isInf(y)) y else @sin(y) * @cos(y);
80 return Complex(f64).init(xx, math.copysign(@as(f64, 0.0), r));
81 }
82
83 if (!math.isFinite(y)) {
84 const r = if (ix != 0) y - y else x;
85 return Complex(f64).init(r, y - y);
86 }
87
88 // x >= 22
89 if (ix >= 0x40360000) {
90 const exp_mx = @exp(-@abs(x));
91 return Complex(f64).init(math.copysign(@as(f64, 1.0), x), 4 * @sin(y) * @cos(y) * exp_mx * exp_mx);
92 }
93
94 // Kahan's algorithm
95 const t = @tan(y);
96 const beta = 1.0 + t * t;
97 const s = math.sinh(x);
98 const rho = @sqrt(1 + s * s);
99 const den = 1 + beta * s * s;
100
101 return Complex(f64).init((beta * rho * s) / den, t / den);
102}
103
104test tanh32 {
105 const epsilon = math.floatEps(f32);
106 const a = Complex(f32).init(5, 3);
107 const c = tanh(a);
108
109 try testing.expectApproxEqAbs(0.99991274, c.re, epsilon);
110 try testing.expectApproxEqAbs(-0.00002536878, c.im, epsilon);
111}
112
113test tanh64 {
114 const epsilon = math.floatEps(f64);
115 const a = Complex(f64).init(5, 3);
116 const c = tanh(a);
117
118 try testing.expectApproxEqAbs(0.9999128201513536, c.re, epsilon);
119 try testing.expectApproxEqAbs(-0.00002536867620767604, c.im, epsilon);
120}
121
122test "tanh64 musl" {
123 const epsilon = math.floatEps(f64);
124 const a = Complex(f64).init(std.math.inf(f64), std.math.inf(f64));
125 const c = tanh(a);
126
127 try testing.expectApproxEqAbs(1, c.re, epsilon);
128 try testing.expectApproxEqAbs(0, c.im, epsilon);
129}