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
|
// 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);
}
|
