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/sinf.c
5//! https://git.musl-libc.org/cgit/musl/tree/src/math/sin.c
6//! https://git.musl-libc.org/cgit/musl/tree/src/math/sinl.c
7
8const std = @import("std");
9const math = std.math;
10const ld = math.long_double;
11const mem = std.mem;
12const expect = std.testing.expect;
13const expectApproxEqAbs = std.testing.expectApproxEqAbs;
14
15const compiler_rt = @import("../compiler_rt.zig");
16const symbol = compiler_rt.symbol;
17const trig = @import("trig.zig");
18const rem_pio2 = @import("rem_pio2.zig").rem_pio2;
19const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f;
20const rem_pio2l = @import("rem_pio2l.zig").rem_pio2l;
21
22comptime {
23 symbol(&__sinh, "__sinh");
24 symbol(&sinf, "sinf");
25 symbol(&sin, "sin");
26 symbol(&__sinx, "__sinx");
27 symbol(&sinq, "sinf128");
28 symbol(&sinl, "sinl");
29 symbol(&sinl, "__sinl"); // required by musl
30}
31
32fn __sinh(x: compiler_rt.f16.Abi) callconv(.c) compiler_rt.f16.Abi {
33 return compiler_rt.f16.toAbi(sin_f16(compiler_rt.f16.fromAbi(x)));
34}
35pub fn sin_f16(x: f16) f16 {
36 // TODO: more efficient implementation
37 return @floatCast(sin_f32(x));
38}
39
40fn sinf(x: compiler_rt.f32.Abi) callconv(.c) compiler_rt.f32.Abi {
41 return compiler_rt.f32.toAbi(sin_f32(compiler_rt.f32.fromAbi(x)));
42}
43pub fn sin_f32(x: f32) f32 {
44 // Small multiples of pi/2 rounded to double precision.
45 const s1pio2: f64 = 1.0 * math.pi / 2.0; // 0x3FF921FB, 0x54442D18
46 const s2pio2: f64 = 2.0 * math.pi / 2.0; // 0x400921FB, 0x54442D18
47 const s3pio2: f64 = 3.0 * math.pi / 2.0; // 0x4012D97C, 0x7F3321D2
48 const s4pio2: f64 = 4.0 * math.pi / 2.0; // 0x401921FB, 0x54442D18
49
50 var ix: u32 = @bitCast(x);
51 const sign = ix >> 31 != 0;
52 ix &= 0x7fffffff;
53
54 if (ix <= 0x3f490fda) { // |x| ~<= pi/4
55 if (ix < 0x39800000) { // |x| < 2**-12
56 // raise inexact if x!=0 and underflow if subnormal
57 if (compiler_rt.want_float_exceptions) {
58 if (ix < 0x00800000) {
59 mem.doNotOptimizeAway(x / 0x1p120);
60 } else {
61 mem.doNotOptimizeAway(x + 0x1p120);
62 }
63 }
64 return x;
65 }
66 return trig.sindf(x);
67 }
68 if (ix <= 0x407b53d1) { // |x| ~<= 5*pi/4
69 if (ix <= 0x4016cbe3) { // |x| ~<= 3pi/4
70 if (sign) {
71 return -trig.cosdf(x + s1pio2);
72 } else {
73 return trig.cosdf(x - s1pio2);
74 }
75 }
76 return trig.sindf(if (sign) -(x + s2pio2) else -(x - s2pio2));
77 }
78 if (ix <= 0x40e231d5) { // |x| ~<= 9*pi/4
79 if (ix <= 0x40afeddf) { // |x| ~<= 7*pi/4
80 if (sign) {
81 return trig.cosdf(x + s3pio2);
82 } else {
83 return -trig.cosdf(x - s3pio2);
84 }
85 }
86 return trig.sindf(if (sign) x + s4pio2 else x - s4pio2);
87 }
88
89 // sin(Inf or NaN) is NaN
90 if (ix >= 0x7f800000) {
91 return x - x;
92 }
93
94 var y: f64 = undefined;
95 const n = rem_pio2f(x, &y);
96 return switch (n & 3) {
97 0 => trig.sindf(y),
98 1 => trig.cosdf(y),
99 2 => trig.sindf(-y),
100 else => -trig.cosdf(y),
101 };
102}
103
104fn sin(x: compiler_rt.f64.Abi) callconv(.c) compiler_rt.f64.Abi {
105 return compiler_rt.f64.toAbi(sin_f64(compiler_rt.f64.fromAbi(x)));
106}
107pub fn sin_f64(x: f64) f64 {
108 var ix = @as(u64, @bitCast(x)) >> 32;
109 ix &= 0x7fffffff;
110
111 // |x| ~< pi/4
112 if (ix <= 0x3fe921fb) {
113 if (ix < 0x3e500000) { // |x| < 2**-26
114 // raise inexact if x != 0 and underflow if subnormal
115 if (compiler_rt.want_float_exceptions) {
116 if (ix < 0x00100000) {
117 mem.doNotOptimizeAway(x / 0x1p120);
118 } else {
119 mem.doNotOptimizeAway(x + 0x1p120);
120 }
121 }
122 return x;
123 }
124 return trig.sin(x, 0.0, 0);
125 }
126
127 // sin(Inf or NaN) is NaN
128 if (ix >= 0x7ff00000) {
129 return x - x;
130 }
131
132 var y: [2]f64 = undefined;
133 const n = rem_pio2(x, &y);
134 return switch (n & 3) {
135 0 => trig.sin(y[0], y[1], 1),
136 1 => trig.cos(y[0], y[1]),
137 2 => -trig.sin(y[0], y[1], 1),
138 else => -trig.cos(y[0], y[1]),
139 };
140}
141
142fn __sinx(x: compiler_rt.f80.Abi) callconv(.c) compiler_rt.f80.Abi {
143 return compiler_rt.f80.toAbi(sin_f80(compiler_rt.f80.fromAbi(x)));
144}
145pub fn sin_f80(x: f80) f80 {
146 const se = ld.signExponent(x) & 0x7fff;
147 if (se == 0x7fff) {
148 return x - x;
149 }
150
151 if (@abs(x) < trig.pi_4) {
152 if (se < 0x3fff - (math.floatMantissaBits(f80) / 2)) {
153 // raise inexact if x!=0 and underflow if subnormal
154 if (compiler_rt.want_float_exceptions) {
155 mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120);
156 }
157 return x;
158 }
159 return trig.sinx(x, 0.0, 0);
160 }
161
162 var y: [2]f80 = undefined;
163 const n = rem_pio2l(f80, x, &y);
164 return switch (n & 3) {
165 0 => trig.sinx(y[0], y[1], 1),
166 1 => trig.cosx(y[0], y[1]),
167 2 => -trig.sinx(y[0], y[1], 1),
168 else => -trig.cosx(y[0], y[1]),
169 };
170}
171
172fn sinq(x: compiler_rt.f128.Abi) callconv(.c) compiler_rt.f128.Abi {
173 return compiler_rt.f128.toAbi(sin_f128(compiler_rt.f128.fromAbi(x)));
174}
175pub fn sin_f128(x: f128) f128 {
176 const se = ld.signExponent(x) & 0x7fff;
177 if (se == 0x7fff) {
178 return x - x;
179 }
180
181 if (@abs(x) < trig.pi_4) {
182 if (se < 0x3fff - (math.floatMantissaBits(f128) / 2)) {
183 // raise inexact if x!=0 and underflow if subnormal
184 if (compiler_rt.want_float_exceptions) {
185 mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120);
186 }
187 return x;
188 }
189 return trig.sinq(x, 0.0, 0);
190 }
191
192 var y: [2]f128 = undefined;
193 const n = rem_pio2l(f128, x, &y);
194 return switch (n & 3) {
195 0 => trig.sinq(y[0], y[1], 1),
196 1 => trig.cosq(y[0], y[1]),
197 2 => -trig.sinq(y[0], y[1], 1),
198 else => -trig.cosq(y[0], y[1]),
199 };
200}
201
202pub fn sinl(x: c_longdouble) callconv(.c) c_longdouble {
203 switch (@typeInfo(c_longdouble).float.bits) {
204 64 => return sin_f64(x),
205 80 => return sin_f80(x),
206 128 => return sin_f128(x),
207 else => comptime unreachable,
208 }
209}
210
211fn testSinSpecial(comptime T: type) !void {
212 const f = switch (T) {
213 f16 => sin_f16,
214 f32 => sin_f32,
215 f64 => sin_f64,
216 f80 => sin_f80,
217 f128 => sin_f128,
218 else => comptime unreachable,
219 };
220
221 try expect(math.isPositiveZero(f(0.0)));
222 try expect(math.isNegativeZero(f(-0.0)));
223 try expect(math.isNan(f(math.inf(T))));
224 try expect(math.isNan(f(-math.inf(T))));
225 try expect(math.isNan(f(math.nan(T))));
226}
227
228test "sin32.normal" {
229 const epsilon = math.floatEps(f32);
230 try expectApproxEqAbs(@as(f32, 0.0), sin_f32(0.0), epsilon);
231 try expectApproxEqAbs(@as(f32, 0.19866933), sin_f32(0.2), epsilon);
232 try expectApproxEqAbs(@as(f32, 0.77851737), sin_f32(0.8923), epsilon);
233 try expectApproxEqAbs(@as(f32, 0.997495), sin_f32(1.5), epsilon);
234 try expectApproxEqAbs(@as(f32, -0.997495), sin_f32(-1.5), epsilon);
235 try expectApproxEqAbs(@as(f32, -0.24654257), sin_f32(37.45), epsilon);
236 try expectApproxEqAbs(@as(f32, 0.9161657), sin_f32(89.123), epsilon);
237}
238
239test "sin32.special" {
240 try testSinSpecial(f32);
241}
242
243test "sin64.normal" {
244 const epsilon = math.floatEps(f64);
245 try expectApproxEqAbs(@as(f64, 0.0), sin_f64(0.0), epsilon);
246 try expectApproxEqAbs(@as(f64, 0.19866933079506122), sin_f64(0.2), epsilon);
247 try expectApproxEqAbs(@as(f64, 0.7785173385577349), sin_f64(0.8923), epsilon);
248 try expectApproxEqAbs(@as(f64, 0.9974949866040544), sin_f64(1.5), epsilon);
249 try expectApproxEqAbs(@as(f64, -0.9974949866040544), sin_f64(-1.5), epsilon);
250 try expectApproxEqAbs(@as(f64, -0.24654331551411082), sin_f64(37.45), epsilon);
251 try expectApproxEqAbs(@as(f64, 0.9161652766622714), sin_f64(89.123), epsilon);
252}
253
254test "sin64.special" {
255 try testSinSpecial(f64);
256}
257
258test "sin80.normal" {
259 const epsilon = math.floatEps(f80);
260 try expectApproxEqAbs(@as(f80, 0.0), sin_f80(0.0), epsilon);
261 try expectApproxEqAbs(@as(f80, 0.19866933079506121545941262711838975), sin_f80(0.2), epsilon);
262 try expectApproxEqAbs(@as(f80, 0.77851733855773487830689285621486050), sin_f80(0.8923), epsilon);
263 try expectApproxEqAbs(@as(f80, 0.99749498660405443094172337114148732), sin_f80(1.5), epsilon);
264 try expectApproxEqAbs(@as(f80, -0.99749498660405443094172337114148732), sin_f80(-1.5), epsilon);
265 try expectApproxEqAbs(@as(f80, -0.24654331551411356504), sin_f80(37.45), epsilon);
266 try expectApproxEqAbs(@as(f80, 0.91616527666226951006), sin_f80(89.123), epsilon);
267}
268
269test "sin80.special" {
270 try testSinSpecial(f80);
271}
272
273test "sin128.normal" {
274 const epsilon = math.floatEps(f128);
275 try expectApproxEqAbs(@as(f128, 0.0), sin_f128(0.0), epsilon);
276 try expectApproxEqAbs(@as(f128, 0.19866933079506121545941262711838975), sin_f128(0.2), epsilon);
277 try expectApproxEqAbs(@as(f128, 0.77851733855773487830689285621486050), sin_f128(0.8923), epsilon);
278 try expectApproxEqAbs(@as(f128, 0.99749498660405443094172337114148732), sin_f128(1.5), epsilon);
279 try expectApproxEqAbs(@as(f128, -0.99749498660405443094172337114148732), sin_f128(-1.5), epsilon);
280 try expectApproxEqAbs(@as(f128, -0.24654331551411356571238581321661085), sin_f128(37.45), epsilon);
281 try expectApproxEqAbs(@as(f128, 0.91616527666226951075019849560482170), sin_f128(89.123), epsilon);
282}
283
284test "sin128.special" {
285 try testSinSpecial(f128);
286}
287
288test "sin32 #9901" {
289 const float: f32 = @bitCast(@as(u32, 0b11100011111111110000000000000000));
290 _ = sin_f32(float);
291}
292
293test "sin64 #9901" {
294 const float: f64 = @bitCast(@as(u64, 0b1111111101000001000000001111110111111111100000000000000000000001));
295 _ = sin_f64(float);
296}