master
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
7const std = @import("std");
8const builtin = @import("builtin");
9const arch = builtin.cpu.arch;
10const math = std.math;
11const mem = std.mem;
12const expect = std.testing.expect;
13const common = @import("common.zig");
14
15const trig = @import("trig.zig");
16const rem_pio2 = @import("rem_pio2.zig").rem_pio2;
17const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f;
18
19pub const panic = common.panic;
20
21comptime {
22 @export(&__sinh, .{ .name = "__sinh", .linkage = common.linkage, .visibility = common.visibility });
23 @export(&sinf, .{ .name = "sinf", .linkage = common.linkage, .visibility = common.visibility });
24 @export(&sin, .{ .name = "sin", .linkage = common.linkage, .visibility = common.visibility });
25 @export(&__sinx, .{ .name = "__sinx", .linkage = common.linkage, .visibility = common.visibility });
26 if (common.want_ppc_abi) {
27 @export(&sinq, .{ .name = "sinf128", .linkage = common.linkage, .visibility = common.visibility });
28 }
29 @export(&sinq, .{ .name = "sinq", .linkage = common.linkage, .visibility = common.visibility });
30 @export(&sinl, .{ .name = "sinl", .linkage = common.linkage, .visibility = common.visibility });
31}
32
33pub fn __sinh(x: f16) callconv(.c) f16 {
34 // TODO: more efficient implementation
35 return @floatCast(sinf(x));
36}
37
38pub fn sinf(x: f32) callconv(.c) f32 {
39 // Small multiples of pi/2 rounded to double precision.
40 const s1pio2: f64 = 1.0 * math.pi / 2.0; // 0x3FF921FB, 0x54442D18
41 const s2pio2: f64 = 2.0 * math.pi / 2.0; // 0x400921FB, 0x54442D18
42 const s3pio2: f64 = 3.0 * math.pi / 2.0; // 0x4012D97C, 0x7F3321D2
43 const s4pio2: f64 = 4.0 * math.pi / 2.0; // 0x401921FB, 0x54442D18
44
45 var ix: u32 = @bitCast(x);
46 const sign = ix >> 31 != 0;
47 ix &= 0x7fffffff;
48
49 if (ix <= 0x3f490fda) { // |x| ~<= pi/4
50 if (ix < 0x39800000) { // |x| < 2**-12
51 // raise inexact if x!=0 and underflow if subnormal
52 if (common.want_float_exceptions) {
53 if (ix < 0x00800000) {
54 mem.doNotOptimizeAway(x / 0x1p120);
55 } else {
56 mem.doNotOptimizeAway(x + 0x1p120);
57 }
58 }
59 return x;
60 }
61 return trig.__sindf(x);
62 }
63 if (ix <= 0x407b53d1) { // |x| ~<= 5*pi/4
64 if (ix <= 0x4016cbe3) { // |x| ~<= 3pi/4
65 if (sign) {
66 return -trig.__cosdf(x + s1pio2);
67 } else {
68 return trig.__cosdf(x - s1pio2);
69 }
70 }
71 return trig.__sindf(if (sign) -(x + s2pio2) else -(x - s2pio2));
72 }
73 if (ix <= 0x40e231d5) { // |x| ~<= 9*pi/4
74 if (ix <= 0x40afeddf) { // |x| ~<= 7*pi/4
75 if (sign) {
76 return trig.__cosdf(x + s3pio2);
77 } else {
78 return -trig.__cosdf(x - s3pio2);
79 }
80 }
81 return trig.__sindf(if (sign) x + s4pio2 else x - s4pio2);
82 }
83
84 // sin(Inf or NaN) is NaN
85 if (ix >= 0x7f800000) {
86 return x - x;
87 }
88
89 var y: f64 = undefined;
90 const n = rem_pio2f(x, &y);
91 return switch (n & 3) {
92 0 => trig.__sindf(y),
93 1 => trig.__cosdf(y),
94 2 => trig.__sindf(-y),
95 else => -trig.__cosdf(y),
96 };
97}
98
99pub fn sin(x: f64) callconv(.c) f64 {
100 var ix = @as(u64, @bitCast(x)) >> 32;
101 ix &= 0x7fffffff;
102
103 // |x| ~< pi/4
104 if (ix <= 0x3fe921fb) {
105 if (ix < 0x3e500000) { // |x| < 2**-26
106 // raise inexact if x != 0 and underflow if subnormal
107 if (common.want_float_exceptions) {
108 if (ix < 0x00100000) {
109 mem.doNotOptimizeAway(x / 0x1p120);
110 } else {
111 mem.doNotOptimizeAway(x + 0x1p120);
112 }
113 }
114 return x;
115 }
116 return trig.__sin(x, 0.0, 0);
117 }
118
119 // sin(Inf or NaN) is NaN
120 if (ix >= 0x7ff00000) {
121 return x - x;
122 }
123
124 var y: [2]f64 = undefined;
125 const n = rem_pio2(x, &y);
126 return switch (n & 3) {
127 0 => trig.__sin(y[0], y[1], 1),
128 1 => trig.__cos(y[0], y[1]),
129 2 => -trig.__sin(y[0], y[1], 1),
130 else => -trig.__cos(y[0], y[1]),
131 };
132}
133
134pub fn __sinx(x: f80) callconv(.c) f80 {
135 // TODO: more efficient implementation
136 return @floatCast(sinq(x));
137}
138
139pub fn sinq(x: f128) callconv(.c) f128 {
140 // TODO: more correct implementation
141 return sin(@floatCast(x));
142}
143
144pub fn sinl(x: c_longdouble) callconv(.c) c_longdouble {
145 switch (@typeInfo(c_longdouble).float.bits) {
146 16 => return __sinh(x),
147 32 => return sinf(x),
148 64 => return sin(x),
149 80 => return __sinx(x),
150 128 => return sinq(x),
151 else => @compileError("unreachable"),
152 }
153}
154
155test "sin32" {
156 const epsilon = 0.00001;
157
158 try expect(math.approxEqAbs(f32, sinf(0.0), 0.0, epsilon));
159 try expect(math.approxEqAbs(f32, sinf(0.2), 0.198669, epsilon));
160 try expect(math.approxEqAbs(f32, sinf(0.8923), 0.778517, epsilon));
161 try expect(math.approxEqAbs(f32, sinf(1.5), 0.997495, epsilon));
162 try expect(math.approxEqAbs(f32, sinf(-1.5), -0.997495, epsilon));
163 try expect(math.approxEqAbs(f32, sinf(37.45), -0.246544, epsilon));
164 try expect(math.approxEqAbs(f32, sinf(89.123), 0.916166, epsilon));
165}
166
167test "sin64" {
168 const epsilon = 0.000001;
169
170 try expect(math.approxEqAbs(f64, sin(0.0), 0.0, epsilon));
171 try expect(math.approxEqAbs(f64, sin(0.2), 0.198669, epsilon));
172 try expect(math.approxEqAbs(f64, sin(0.8923), 0.778517, epsilon));
173 try expect(math.approxEqAbs(f64, sin(1.5), 0.997495, epsilon));
174 try expect(math.approxEqAbs(f64, sin(-1.5), -0.997495, epsilon));
175 try expect(math.approxEqAbs(f64, sin(37.45), -0.246543, epsilon));
176 try expect(math.approxEqAbs(f64, sin(89.123), 0.916166, epsilon));
177}
178
179test "sin32.special" {
180 try expect(sinf(0.0) == 0.0);
181 try expect(sinf(-0.0) == -0.0);
182 try expect(math.isNan(sinf(math.inf(f32))));
183 try expect(math.isNan(sinf(-math.inf(f32))));
184 try expect(math.isNan(sinf(math.nan(f32))));
185}
186
187test "sin64.special" {
188 try expect(sin(0.0) == 0.0);
189 try expect(sin(-0.0) == -0.0);
190 try expect(math.isNan(sin(math.inf(f64))));
191 try expect(math.isNan(sin(-math.inf(f64))));
192 try expect(math.isNan(sin(math.nan(f64))));
193}
194
195test "sin32 #9901" {
196 const float: f32 = @bitCast(@as(u32, 0b11100011111111110000000000000000));
197 _ = sinf(float);
198}
199
200test "sin64 #9901" {
201 const float: f64 = @bitCast(@as(u64, 0b1111111101000001000000001111110111111111100000000000000000000001));
202 _ = sin(float);
203}