master
1const std = @import("std");
2const builtin = @import("builtin");
3const arch = builtin.cpu.arch;
4const math = std.math;
5const mem = std.mem;
6const trig = @import("trig.zig");
7const rem_pio2 = @import("rem_pio2.zig").rem_pio2;
8const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f;
9const common = @import("common.zig");
10
11pub const panic = common.panic;
12
13comptime {
14 @export(&__sincosh, .{ .name = "__sincosh", .linkage = common.linkage, .visibility = common.visibility });
15 @export(&sincosf, .{ .name = "sincosf", .linkage = common.linkage, .visibility = common.visibility });
16 @export(&sincos, .{ .name = "sincos", .linkage = common.linkage, .visibility = common.visibility });
17 @export(&__sincosx, .{ .name = "__sincosx", .linkage = common.linkage, .visibility = common.visibility });
18 if (common.want_ppc_abi) {
19 @export(&sincosq, .{ .name = "sincosf128", .linkage = common.linkage, .visibility = common.visibility });
20 }
21 @export(&sincosq, .{ .name = "sincosq", .linkage = common.linkage, .visibility = common.visibility });
22 @export(&sincosl, .{ .name = "sincosl", .linkage = common.linkage, .visibility = common.visibility });
23}
24
25pub fn __sincosh(x: f16, r_sin: *f16, r_cos: *f16) callconv(.c) void {
26 // TODO: more efficient implementation
27 var big_sin: f32 = undefined;
28 var big_cos: f32 = undefined;
29 sincosf(x, &big_sin, &big_cos);
30 r_sin.* = @as(f16, @floatCast(big_sin));
31 r_cos.* = @as(f16, @floatCast(big_cos));
32}
33
34pub fn sincosf(x: f32, r_sin: *f32, r_cos: *f32) callconv(.c) void {
35 const sc1pio2: f64 = 1.0 * math.pi / 2.0; // 0x3FF921FB, 0x54442D18
36 const sc2pio2: f64 = 2.0 * math.pi / 2.0; // 0x400921FB, 0x54442D18
37 const sc3pio2: f64 = 3.0 * math.pi / 2.0; // 0x4012D97C, 0x7F3321D2
38 const sc4pio2: f64 = 4.0 * math.pi / 2.0; // 0x401921FB, 0x54442D18
39
40 const pre_ix = @as(u32, @bitCast(x));
41 const sign = pre_ix >> 31 != 0;
42 const ix = pre_ix & 0x7fffffff;
43
44 // |x| ~<= pi/4
45 if (ix <= 0x3f490fda) {
46 // |x| < 2**-12
47 if (ix < 0x39800000) {
48 // raise inexact if x!=0 and underflow if subnormal
49 if (common.want_float_exceptions) {
50 if (ix < 0x00100000) {
51 mem.doNotOptimizeAway(x / 0x1p120);
52 } else {
53 mem.doNotOptimizeAway(x + 0x1p120);
54 }
55 }
56 r_sin.* = x;
57 r_cos.* = 1.0;
58 return;
59 }
60 r_sin.* = trig.__sindf(x);
61 r_cos.* = trig.__cosdf(x);
62 return;
63 }
64
65 // |x| ~<= 5*pi/4
66 if (ix <= 0x407b53d1) {
67 // |x| ~<= 3pi/4
68 if (ix <= 0x4016cbe3) {
69 if (sign) {
70 r_sin.* = -trig.__cosdf(x + sc1pio2);
71 r_cos.* = trig.__sindf(x + sc1pio2);
72 } else {
73 r_sin.* = trig.__cosdf(sc1pio2 - x);
74 r_cos.* = trig.__sindf(sc1pio2 - x);
75 }
76 return;
77 }
78 // -sin(x+c) is not correct if x+c could be 0: -0 vs +0
79 r_sin.* = -trig.__sindf(if (sign) x + sc2pio2 else x - sc2pio2);
80 r_cos.* = -trig.__cosdf(if (sign) x + sc2pio2 else x - sc2pio2);
81 return;
82 }
83
84 // |x| ~<= 9*pi/4
85 if (ix <= 0x40e231d5) {
86 // |x| ~<= 7*pi/4
87 if (ix <= 0x40afeddf) {
88 if (sign) {
89 r_sin.* = trig.__cosdf(x + sc3pio2);
90 r_cos.* = -trig.__sindf(x + sc3pio2);
91 } else {
92 r_sin.* = -trig.__cosdf(x - sc3pio2);
93 r_cos.* = trig.__sindf(x - sc3pio2);
94 }
95 return;
96 }
97 r_sin.* = trig.__sindf(if (sign) x + sc4pio2 else x - sc4pio2);
98 r_cos.* = trig.__cosdf(if (sign) x + sc4pio2 else x - sc4pio2);
99 return;
100 }
101
102 // sin(Inf or NaN) is NaN
103 if (ix >= 0x7f800000) {
104 const result = x - x;
105 r_sin.* = result;
106 r_cos.* = result;
107 return;
108 }
109
110 // general argument reduction needed
111 var y: f64 = undefined;
112 const n = rem_pio2f(x, &y);
113 const s = trig.__sindf(y);
114 const c = trig.__cosdf(y);
115 switch (n & 3) {
116 0 => {
117 r_sin.* = s;
118 r_cos.* = c;
119 },
120 1 => {
121 r_sin.* = c;
122 r_cos.* = -s;
123 },
124 2 => {
125 r_sin.* = -s;
126 r_cos.* = -c;
127 },
128 else => {
129 r_sin.* = -c;
130 r_cos.* = s;
131 },
132 }
133}
134
135pub fn sincos(x: f64, r_sin: *f64, r_cos: *f64) callconv(.c) void {
136 const ix = @as(u32, @truncate(@as(u64, @bitCast(x)) >> 32)) & 0x7fffffff;
137
138 // |x| ~< pi/4
139 if (ix <= 0x3fe921fb) {
140 // if |x| < 2**-27 * sqrt(2)
141 if (ix < 0x3e46a09e) {
142 // raise inexact if x != 0 and underflow if subnormal
143 if (common.want_float_exceptions) {
144 if (ix < 0x00100000) {
145 mem.doNotOptimizeAway(x / 0x1p120);
146 } else {
147 mem.doNotOptimizeAway(x + 0x1p120);
148 }
149 }
150 r_sin.* = x;
151 r_cos.* = 1.0;
152 return;
153 }
154 r_sin.* = trig.__sin(x, 0.0, 0);
155 r_cos.* = trig.__cos(x, 0.0);
156 return;
157 }
158
159 // sincos(Inf or NaN) is NaN
160 if (ix >= 0x7ff00000) {
161 const result = x - x;
162 r_sin.* = result;
163 r_cos.* = result;
164 return;
165 }
166
167 // argument reduction needed
168 var y: [2]f64 = undefined;
169 const n = rem_pio2(x, &y);
170 const s = trig.__sin(y[0], y[1], 1);
171 const c = trig.__cos(y[0], y[1]);
172 switch (n & 3) {
173 0 => {
174 r_sin.* = s;
175 r_cos.* = c;
176 },
177 1 => {
178 r_sin.* = c;
179 r_cos.* = -s;
180 },
181 2 => {
182 r_sin.* = -s;
183 r_cos.* = -c;
184 },
185 else => {
186 r_sin.* = -c;
187 r_cos.* = s;
188 },
189 }
190}
191
192pub fn __sincosx(x: f80, r_sin: *f80, r_cos: *f80) callconv(.c) void {
193 // TODO: more efficient implementation
194 //return sincos_generic(f80, x, r_sin, r_cos);
195 var big_sin: f128 = undefined;
196 var big_cos: f128 = undefined;
197 sincosq(x, &big_sin, &big_cos);
198 r_sin.* = @as(f80, @floatCast(big_sin));
199 r_cos.* = @as(f80, @floatCast(big_cos));
200}
201
202pub fn sincosq(x: f128, r_sin: *f128, r_cos: *f128) callconv(.c) void {
203 // TODO: more correct implementation
204 //return sincos_generic(f128, x, r_sin, r_cos);
205 var small_sin: f64 = undefined;
206 var small_cos: f64 = undefined;
207 sincos(@as(f64, @floatCast(x)), &small_sin, &small_cos);
208 r_sin.* = small_sin;
209 r_cos.* = small_cos;
210}
211
212pub fn sincosl(x: c_longdouble, r_sin: *c_longdouble, r_cos: *c_longdouble) callconv(.c) void {
213 switch (@typeInfo(c_longdouble).float.bits) {
214 16 => return __sincosh(x, r_sin, r_cos),
215 32 => return sincosf(x, r_sin, r_cos),
216 64 => return sincos(x, r_sin, r_cos),
217 80 => return __sincosx(x, r_sin, r_cos),
218 128 => return sincosq(x, r_sin, r_cos),
219 else => @compileError("unreachable"),
220 }
221}
222
223pub const rem_pio2_generic = @compileError("TODO");
224
225/// Ported from musl sincosl.c. Needs the following dependencies to be complete:
226/// * rem_pio2_generic ported from __rem_pio2l.c
227/// * trig.sin_generic ported from __sinl.c
228/// * trig.cos_generic ported from __cosl.c
229inline fn sincos_generic(comptime F: type, x: F, r_sin: *F, r_cos: *F) void {
230 const sc1pio4: F = 1.0 * math.pi / 4.0;
231 const bits = @typeInfo(F).float.bits;
232 const I = std.meta.Int(.unsigned, bits);
233 const ix = @as(I, @bitCast(x)) & (math.maxInt(I) >> 1);
234 const se: u16 = @truncate(ix >> (bits - 16));
235
236 if (se == 0x7fff) {
237 const result = x - x;
238 r_sin.* = result;
239 r_cos.* = result;
240 return;
241 }
242
243 if (@as(F, @bitCast(ix)) < sc1pio4) {
244 if (se < 0x3fff - math.floatFractionalBits(F) - 1) {
245 // raise underflow if subnormal
246 if (se == 0) {
247 if (common.want_float_exceptions) mem.doNotOptimizeAway(x * 0x1p-120);
248 }
249 r_sin.* = x;
250 // raise inexact if x!=0
251 r_cos.* = 1.0 + x;
252 return;
253 }
254 r_sin.* = trig.sin_generic(F, x, 0, 0);
255 r_cos.* = trig.cos_generic(F, x, 0);
256 return;
257 }
258
259 var y: [2]F = undefined;
260 const n = rem_pio2_generic(F, x, &y);
261 const s = trig.sin_generic(F, y[0], y[1], 1);
262 const c = trig.cos_generic(F, y[0], y[1]);
263 switch (n & 3) {
264 0 => {
265 r_sin.* = s;
266 r_cos.* = c;
267 },
268 1 => {
269 r_sin.* = c;
270 r_cos.* = -s;
271 },
272 2 => {
273 r_sin.* = -s;
274 r_cos.* = -c;
275 },
276 else => {
277 r_sin.* = -c;
278 r_cos.* = s;
279 },
280 }
281}