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}