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