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/sinhf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/sinh.c
6
7const std = @import("../std.zig");
8const math = std.math;
9const expect = std.testing.expect;
10const expo2 = @import("expo2.zig").expo2;
11const maxInt = std.math.maxInt;
12
13/// Returns the hyperbolic sine of x.
14///
15/// Special Cases:
16/// - sinh(+-0) = +-0
17/// - sinh(+-inf) = +-inf
18/// - sinh(nan) = nan
19pub fn sinh(x: anytype) @TypeOf(x) {
20 const T = @TypeOf(x);
21 return switch (T) {
22 f32 => sinh32(x),
23 f64 => sinh64(x),
24 else => @compileError("sinh not implemented for " ++ @typeName(T)),
25 };
26}
27
28// sinh(x) = (exp(x) - 1 / exp(x)) / 2
29// = (exp(x) - 1 + (exp(x) - 1) / exp(x)) / 2
30// = x + x^3 / 6 + o(x^5)
31fn sinh32(x: f32) f32 {
32 const u = @as(u32, @bitCast(x));
33 const ux = u & 0x7FFFFFFF;
34 const ax = @as(f32, @bitCast(ux));
35
36 if (x == 0.0 or math.isNan(x)) {
37 return x;
38 }
39
40 var h: f32 = 0.5;
41 if (u >> 31 != 0) {
42 h = -h;
43 }
44
45 // |x| < log(FLT_MAX)
46 if (ux < 0x42B17217) {
47 const t = math.expm1(ax);
48 if (ux < 0x3F800000) {
49 if (ux < 0x3F800000 - (12 << 23)) {
50 return x;
51 } else {
52 return h * (2 * t - t * t / (t + 1));
53 }
54 }
55 return h * (t + t / (t + 1));
56 }
57
58 // |x| > log(FLT_MAX) or nan
59 return 2 * h * expo2(ax);
60}
61
62fn sinh64(x: f64) f64 {
63 const u = @as(u64, @bitCast(x));
64 const w = @as(u32, @intCast(u >> 32)) & (maxInt(u32) >> 1);
65 const ax = @as(f64, @bitCast(u & (maxInt(u64) >> 1)));
66
67 if (x == 0.0 or math.isNan(x)) {
68 return x;
69 }
70
71 var h: f32 = 0.5;
72 if (u >> 63 != 0) {
73 h = -h;
74 }
75
76 // |x| < log(FLT_MAX)
77 if (w < 0x40862E42) {
78 const t = math.expm1(ax);
79 if (w < 0x3FF00000) {
80 if (w < 0x3FF00000 - (26 << 20)) {
81 return x;
82 } else {
83 return h * (2 * t - t * t / (t + 1));
84 }
85 }
86 // NOTE: |x| > log(0x1p26) + eps could be h * exp(x)
87 return h * (t + t / (t + 1));
88 }
89
90 // |x| > log(DBL_MAX) or nan
91 return 2 * h * expo2(ax);
92}
93
94test sinh {
95 try expect(sinh(@as(f32, 1.5)) == sinh32(1.5));
96 try expect(sinh(@as(f64, 1.5)) == sinh64(1.5));
97}
98
99test sinh32 {
100 const epsilon = 0.000001;
101
102 try expect(math.approxEqAbs(f32, sinh32(0.0), 0.0, epsilon));
103 try expect(math.approxEqAbs(f32, sinh32(0.2), 0.201336, epsilon));
104 try expect(math.approxEqAbs(f32, sinh32(0.8923), 1.015512, epsilon));
105 try expect(math.approxEqAbs(f32, sinh32(1.5), 2.129279, epsilon));
106 try expect(math.approxEqAbs(f32, sinh32(-0.0), -0.0, epsilon));
107 try expect(math.approxEqAbs(f32, sinh32(-0.2), -0.201336, epsilon));
108 try expect(math.approxEqAbs(f32, sinh32(-0.8923), -1.015512, epsilon));
109 try expect(math.approxEqAbs(f32, sinh32(-1.5), -2.129279, epsilon));
110}
111
112test sinh64 {
113 const epsilon = 0.000001;
114
115 try expect(math.approxEqAbs(f64, sinh64(0.0), 0.0, epsilon));
116 try expect(math.approxEqAbs(f64, sinh64(0.2), 0.201336, epsilon));
117 try expect(math.approxEqAbs(f64, sinh64(0.8923), 1.015512, epsilon));
118 try expect(math.approxEqAbs(f64, sinh64(1.5), 2.129279, epsilon));
119 try expect(math.approxEqAbs(f64, sinh64(-0.0), -0.0, epsilon));
120 try expect(math.approxEqAbs(f64, sinh64(-0.2), -0.201336, epsilon));
121 try expect(math.approxEqAbs(f64, sinh64(-0.8923), -1.015512, epsilon));
122 try expect(math.approxEqAbs(f64, sinh64(-1.5), -2.129279, epsilon));
123}
124
125test "sinh32.special" {
126 try expect(math.isPositiveZero(sinh32(0.0)));
127 try expect(math.isNegativeZero(sinh32(-0.0)));
128 try expect(math.isPositiveInf(sinh32(math.inf(f32))));
129 try expect(math.isNegativeInf(sinh32(-math.inf(f32))));
130 try expect(math.isNan(sinh32(math.nan(f32))));
131}
132
133test "sinh64.special" {
134 try expect(math.isPositiveZero(sinh64(0.0)));
135 try expect(math.isNegativeZero(sinh64(-0.0)));
136 try expect(math.isPositiveInf(sinh64(math.inf(f64))));
137 try expect(math.isNegativeInf(sinh64(-math.inf(f64))));
138 try expect(math.isNan(sinh64(math.nan(f64))));
139}