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/math/tanhf.c
  5// https://git.musl-libc.org/cgit/musl/tree/src/math/tanh.c
  6
  7const std = @import("../std.zig");
  8const math = std.math;
  9const mem = std.mem;
 10const expect = std.testing.expect;
 11const expo2 = @import("expo2.zig").expo2;
 12const maxInt = std.math.maxInt;
 13
 14/// Returns the hyperbolic tangent of x.
 15///
 16/// Special Cases:
 17///  - tanh(+-0)   = +-0
 18///  - tanh(+-inf) = +-1
 19///  - tanh(nan)   = nan
 20pub fn tanh(x: anytype) @TypeOf(x) {
 21    const T = @TypeOf(x);
 22    return switch (T) {
 23        f32 => tanh32(x),
 24        f64 => tanh64(x),
 25        else => @compileError("tanh not implemented for " ++ @typeName(T)),
 26    };
 27}
 28
 29// tanh(x) = (exp(x) - exp(-x)) / (exp(x) + exp(-x))
 30//         = (exp(2x) - 1) / (exp(2x) - 1 + 2)
 31//         = (1 - exp(-2x)) / (exp(-2x) - 1 + 2)
 32fn tanh32(x: f32) f32 {
 33    const u = @as(u32, @bitCast(x));
 34    const ux = u & 0x7FFFFFFF;
 35    const ax = @as(f32, @bitCast(ux));
 36    const sign = (u >> 31) != 0;
 37
 38    var t: f32 = undefined;
 39
 40    // |x| < log(3) / 2 ~= 0.5493 or nan
 41    if (ux > 0x3F0C9F54) {
 42        // |x| > 10
 43        if (ux > 0x41200000) {
 44            t = 1.0 + 0 / x;
 45        } else {
 46            t = math.expm1(2 * ax);
 47            t = 1 - 2 / (t + 2);
 48        }
 49    }
 50    // |x| > log(5 / 3) / 2 ~= 0.2554
 51    else if (ux > 0x3E82C578) {
 52        t = math.expm1(2 * ax);
 53        t = t / (t + 2);
 54    }
 55    // |x| >= 0x1.0p-126
 56    else if (ux >= 0x00800000) {
 57        t = math.expm1(-2 * ax);
 58        t = -t / (t + 2);
 59    }
 60    // |x| is subnormal
 61    else {
 62        mem.doNotOptimizeAway(ax * ax);
 63        t = ax;
 64    }
 65
 66    return if (sign) -t else t;
 67}
 68
 69fn tanh64(x: f64) f64 {
 70    const u = @as(u64, @bitCast(x));
 71    const ux = u & 0x7FFFFFFFFFFFFFFF;
 72    const w = @as(u32, @intCast(ux >> 32));
 73    const ax = @as(f64, @bitCast(ux));
 74    const sign = (u >> 63) != 0;
 75
 76    var t: f64 = undefined;
 77
 78    // |x| < log(3) / 2 ~= 0.5493 or nan
 79    if (w > 0x3FE193EA) {
 80        // |x| > 20 or nan
 81        if (w > 0x40340000) {
 82            t = 1.0 - 0 / ax;
 83        } else {
 84            t = math.expm1(2 * ax);
 85            t = 1 - 2 / (t + 2);
 86        }
 87    }
 88    // |x| > log(5 / 3) / 2 ~= 0.2554
 89    else if (w > 0x3FD058AE) {
 90        t = math.expm1(2 * ax);
 91        t = t / (t + 2);
 92    }
 93    // |x| >= 0x1.0p-1022
 94    else if (w >= 0x00100000) {
 95        t = math.expm1(-2 * ax);
 96        t = -t / (t + 2);
 97    }
 98    // |x| is subnormal
 99    else {
100        mem.doNotOptimizeAway(@as(f32, @floatCast(ax)));
101        t = ax;
102    }
103
104    return if (sign) -t else t;
105}
106
107test tanh {
108    try expect(tanh(@as(f32, 1.5)) == tanh32(1.5));
109    try expect(tanh(@as(f64, 1.5)) == tanh64(1.5));
110}
111
112test tanh32 {
113    const epsilon = 0.000001;
114
115    try expect(math.approxEqAbs(f32, tanh32(0.0), 0.0, epsilon));
116    try expect(math.approxEqAbs(f32, tanh32(0.2), 0.197375, epsilon));
117    try expect(math.approxEqAbs(f32, tanh32(0.8923), 0.712528, epsilon));
118    try expect(math.approxEqAbs(f32, tanh32(1.5), 0.905148, epsilon));
119    try expect(math.approxEqAbs(f32, tanh32(37.45), 1.0, epsilon));
120    try expect(math.approxEqAbs(f32, tanh32(-0.8923), -0.712528, epsilon));
121    try expect(math.approxEqAbs(f32, tanh32(-1.5), -0.905148, epsilon));
122    try expect(math.approxEqAbs(f32, tanh32(-37.45), -1.0, epsilon));
123}
124
125test tanh64 {
126    const epsilon = 0.000001;
127
128    try expect(math.approxEqAbs(f64, tanh64(0.0), 0.0, epsilon));
129    try expect(math.approxEqAbs(f64, tanh64(0.2), 0.197375, epsilon));
130    try expect(math.approxEqAbs(f64, tanh64(0.8923), 0.712528, epsilon));
131    try expect(math.approxEqAbs(f64, tanh64(1.5), 0.905148, epsilon));
132    try expect(math.approxEqAbs(f64, tanh64(37.45), 1.0, epsilon));
133    try expect(math.approxEqAbs(f64, tanh64(-0.8923), -0.712528, epsilon));
134    try expect(math.approxEqAbs(f64, tanh64(-1.5), -0.905148, epsilon));
135    try expect(math.approxEqAbs(f64, tanh64(-37.45), -1.0, epsilon));
136}
137
138test "tanh32.special" {
139    try expect(math.isPositiveZero(tanh32(0.0)));
140    try expect(math.isNegativeZero(tanh32(-0.0)));
141    try expect(tanh32(math.inf(f32)) == 1.0);
142    try expect(tanh32(-math.inf(f32)) == -1.0);
143    try expect(math.isNan(tanh32(math.nan(f32))));
144}
145
146test "tanh64.special" {
147    try expect(math.isPositiveZero(tanh64(0.0)));
148    try expect(math.isNegativeZero(tanh64(-0.0)));
149    try expect(tanh64(math.inf(f64)) == 1.0);
150    try expect(tanh64(-math.inf(f64)) == -1.0);
151    try expect(math.isNan(tanh64(math.nan(f64))));
152}