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/asinhf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/asinh.c
6
7const std = @import("../std.zig");
8const math = std.math;
9const mem = std.mem;
10const expect = std.testing.expect;
11const maxInt = std.math.maxInt;
12
13/// Returns the hyperbolic arc-sin of x.
14///
15/// Special Cases:
16/// - asinh(+-0) = +-0
17/// - asinh(+-inf) = +-inf
18/// - asinh(nan) = nan
19pub fn asinh(x: anytype) @TypeOf(x) {
20 const T = @TypeOf(x);
21 return switch (T) {
22 f32 => asinh32(x),
23 f64 => asinh64(x),
24 else => @compileError("asinh not implemented for " ++ @typeName(T)),
25 };
26}
27
28// asinh(x) = sign(x) * log(|x| + sqrt(x * x + 1)) ~= x - x^3/6 + o(x^5)
29fn asinh32(x: f32) f32 {
30 const u = @as(u32, @bitCast(x));
31 const i = u & 0x7FFFFFFF;
32 const s = u >> 31;
33
34 var rx = @as(f32, @bitCast(i)); // |x|
35
36 // |x| >= 0x1p12 or inf or nan
37 if (i >= 0x3F800000 + (12 << 23)) {
38 rx = @log(rx) + 0.69314718055994530941723212145817656;
39 }
40 // |x| >= 2
41 else if (i >= 0x3F800000 + (1 << 23)) {
42 rx = @log(2 * rx + 1 / (@sqrt(rx * rx + 1) + rx));
43 }
44 // |x| >= 0x1p-12, up to 1.6ulp error
45 else if (i >= 0x3F800000 - (12 << 23)) {
46 rx = math.log1p(rx + rx * rx / (@sqrt(rx * rx + 1) + 1));
47 }
48 // |x| < 0x1p-12, inexact if x != 0
49 else {
50 mem.doNotOptimizeAway(rx + 0x1.0p120);
51 }
52
53 return if (s != 0) -rx else rx;
54}
55
56fn asinh64(x: f64) f64 {
57 const u = @as(u64, @bitCast(x));
58 const e = (u >> 52) & 0x7FF;
59 const s = u >> 63;
60
61 var rx = @as(f64, @bitCast(u & (maxInt(u64) >> 1))); // |x|
62
63 // |x| >= 0x1p26 or inf or nan
64 if (e >= 0x3FF + 26) {
65 rx = @log(rx) + 0.693147180559945309417232121458176568;
66 }
67 // |x| >= 2
68 else if (e >= 0x3FF + 1) {
69 rx = @log(2 * rx + 1 / (@sqrt(rx * rx + 1) + rx));
70 }
71 // |x| >= 0x1p-12, up to 1.6ulp error
72 else if (e >= 0x3FF - 26) {
73 rx = math.log1p(rx + rx * rx / (@sqrt(rx * rx + 1) + 1));
74 }
75 // |x| < 0x1p-12, inexact if x != 0
76 else {
77 mem.doNotOptimizeAway(rx + 0x1.0p120);
78 }
79
80 return if (s != 0) -rx else rx;
81}
82
83test asinh {
84 try expect(asinh(@as(f32, 0.0)) == asinh32(0.0));
85 try expect(asinh(@as(f64, 0.0)) == asinh64(0.0));
86}
87
88test asinh32 {
89 const epsilon = 0.000001;
90
91 try expect(math.approxEqAbs(f32, asinh32(0.0), 0.0, epsilon));
92 try expect(math.approxEqAbs(f32, asinh32(-0.2), -0.198690, epsilon));
93 try expect(math.approxEqAbs(f32, asinh32(0.2), 0.198690, epsilon));
94 try expect(math.approxEqAbs(f32, asinh32(0.8923), 0.803133, epsilon));
95 try expect(math.approxEqAbs(f32, asinh32(1.5), 1.194763, epsilon));
96 try expect(math.approxEqAbs(f32, asinh32(37.45), 4.316332, epsilon));
97 try expect(math.approxEqAbs(f32, asinh32(89.123), 5.183196, epsilon));
98 try expect(math.approxEqAbs(f32, asinh32(123123.234375), 12.414088, epsilon));
99}
100
101test asinh64 {
102 const epsilon = 0.000001;
103
104 try expect(math.approxEqAbs(f64, asinh64(0.0), 0.0, epsilon));
105 try expect(math.approxEqAbs(f64, asinh64(-0.2), -0.198690, epsilon));
106 try expect(math.approxEqAbs(f64, asinh64(0.2), 0.198690, epsilon));
107 try expect(math.approxEqAbs(f64, asinh64(0.8923), 0.803133, epsilon));
108 try expect(math.approxEqAbs(f64, asinh64(1.5), 1.194763, epsilon));
109 try expect(math.approxEqAbs(f64, asinh64(37.45), 4.316332, epsilon));
110 try expect(math.approxEqAbs(f64, asinh64(89.123), 5.183196, epsilon));
111 try expect(math.approxEqAbs(f64, asinh64(123123.234375), 12.414088, epsilon));
112}
113
114test "asinh32.special" {
115 try expect(math.isPositiveZero(asinh32(0.0)));
116 try expect(math.isNegativeZero(asinh32(-0.0)));
117 try expect(math.isPositiveInf(asinh32(math.inf(f32))));
118 try expect(math.isNegativeInf(asinh32(-math.inf(f32))));
119 try expect(math.isNan(asinh32(math.nan(f32))));
120}
121
122test "asinh64.special" {
123 try expect(math.isPositiveZero(asinh64(0.0)));
124 try expect(math.isNegativeZero(asinh64(-0.0)));
125 try expect(math.isPositiveInf(asinh64(math.inf(f64))));
126 try expect(math.isNegativeInf(asinh64(-math.inf(f64))));
127 try expect(math.isNan(asinh64(math.nan(f64))));
128}