1const std = @import("std");
2const builtin = @import("builtin");
3const arch = builtin.cpu.arch;
4const math = std.math;
5const ld = math.long_double;
6const mem = std.mem;
7const expect = std.testing.expect;
8const expectApproxEqAbs = std.testing.expectApproxEqAbs;
9const trig = @import("trig.zig");
10const rem_pio2 = @import("rem_pio2.zig").rem_pio2;
11const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f;
12const rem_pio2l = @import("rem_pio2l.zig").rem_pio2l;
13const compiler_rt = @import("../compiler_rt.zig");
14const symbol = compiler_rt.symbol;
15
16comptime {
17 symbol(&sincosh, "__sincosh");
18 symbol(&sincosf, "sincosf");
19 symbol(&sincos, "sincos");
20 symbol(&sincosx, "__sincosx");
21 symbol(&sincosq, "sincosf128");
22 symbol(&sincosl, "sincosl");
23}
24
25fn sincosh(x: compiler_rt.f16.Abi, r_sin: *compiler_rt.f16.Abi, r_cos: *compiler_rt.f16.Abi) callconv(.c) void {
26 const s, const c = sincos_f16(compiler_rt.f16.fromAbi(x));
27 r_sin.* = compiler_rt.f16.toAbi(s);
28 r_cos.* = compiler_rt.f16.toAbi(c);
29}
30pub fn sincos_f16(x: f16) struct { f16, f16 } {
31 // TODO: more efficient implementation
32 const s, const c = sincos_f32(x);
33 return .{ @floatCast(s), @floatCast(c) };
34}
35
36fn sincosf(x: compiler_rt.f32.Abi, r_sin: *compiler_rt.f32.Abi, r_cos: *compiler_rt.f32.Abi) callconv(.c) void {
37 const s, const c = sincos_f32(compiler_rt.f32.fromAbi(x));
38 r_sin.* = compiler_rt.f32.toAbi(s);
39 r_cos.* = compiler_rt.f32.toAbi(c);
40}
41pub fn sincos_f32(x: f32) struct { f32, f32 } {
42 const sc1pio2: f64 = 1.0 * math.pi / 2.0; // 0x3FF921FB, 0x54442D18
43 const sc2pio2: f64 = 2.0 * math.pi / 2.0; // 0x400921FB, 0x54442D18
44 const sc3pio2: f64 = 3.0 * math.pi / 2.0; // 0x4012D97C, 0x7F3321D2
45 const sc4pio2: f64 = 4.0 * math.pi / 2.0; // 0x401921FB, 0x54442D18
46
47 const pre_ix = @as(u32, @bitCast(x));
48 const sign = pre_ix >> 31 != 0;
49 const ix = pre_ix & 0x7fffffff;
50
51 // |x| ~<= pi/4
52 if (ix <= 0x3f490fda) {
53 // |x| < 2**-12
54 if (ix < 0x39800000) {
55 // raise inexact if x!=0 and underflow if subnormal
56 if (compiler_rt.want_float_exceptions) {
57 if (ix < 0x00100000) {
58 mem.doNotOptimizeAway(x / 0x1p120);
59 } else {
60 mem.doNotOptimizeAway(x + 0x1p120);
61 }
62 }
63 return .{ x, 1.0 };
64 }
65 return .{ trig.sindf(x), trig.cosdf(x) };
66 }
67
68 // |x| ~<= 5*pi/4
69 if (ix <= 0x407b53d1) {
70 // |x| ~<= 3pi/4
71 if (ix <= 0x4016cbe3) {
72 if (sign) {
73 return .{ -trig.cosdf(x + sc1pio2), trig.sindf(x + sc1pio2) };
74 } else {
75 return .{ trig.cosdf(sc1pio2 - x), trig.sindf(sc1pio2 - x) };
76 }
77 }
78 // -sin(x+c) is not correct if x+c could be 0: -0 vs +0
79 return .{
80 -trig.sindf(if (sign) x + sc2pio2 else x - sc2pio2),
81 -trig.cosdf(if (sign) x + sc2pio2 else x - sc2pio2),
82 };
83 }
84
85 // |x| ~<= 9*pi/4
86 if (ix <= 0x40e231d5) {
87 // |x| ~<= 7*pi/4
88 if (ix <= 0x40afeddf) {
89 if (sign) {
90 return .{ trig.cosdf(x + sc3pio2), -trig.sindf(x + sc3pio2) };
91 } else {
92 return .{ -trig.cosdf(x - sc3pio2), trig.sindf(x - sc3pio2) };
93 }
94 }
95 return .{
96 trig.sindf(if (sign) x + sc4pio2 else x - sc4pio2),
97 trig.cosdf(if (sign) x + sc4pio2 else x - sc4pio2),
98 };
99 }
100
101 // sin(Inf or NaN) is NaN
102 if (ix >= 0x7f800000) {
103 const result = x - x;
104 return .{ result, result };
105 }
106
107 // general argument reduction needed
108 var y: f64 = undefined;
109 const n = rem_pio2f(x, &y);
110 const s = trig.sindf(y);
111 const c = trig.cosdf(y);
112 return switch (@as(u2, @truncate(@as(u32, @bitCast(n))))) {
113 0 => .{ s, c },
114 1 => .{ c, -s },
115 2 => .{ -s, -c },
116 3 => .{ -c, s },
117 };
118}
119
120fn sincos(x: compiler_rt.f64.Abi, r_sin: *compiler_rt.f64.Abi, r_cos: *compiler_rt.f64.Abi) callconv(.c) void {
121 const s, const c = sincos_f64(compiler_rt.f64.fromAbi(x));
122 r_sin.* = compiler_rt.f64.toAbi(s);
123 r_cos.* = compiler_rt.f64.toAbi(c);
124}
125pub fn sincos_f64(x: f64) struct { f64, f64 } {
126 const ix = @as(u32, @truncate(@as(u64, @bitCast(x)) >> 32)) & 0x7fffffff;
127
128 // |x| ~< pi/4
129 if (ix <= 0x3fe921fb) {
130 // if |x| < 2**-27 * sqrt(2)
131 if (ix < 0x3e46a09e) {
132 // raise inexact if x != 0 and underflow if subnormal
133 if (compiler_rt.want_float_exceptions) {
134 if (ix < 0x00100000) {
135 mem.doNotOptimizeAway(x / 0x1p120);
136 } else {
137 mem.doNotOptimizeAway(x + 0x1p120);
138 }
139 }
140 return .{ x, 1.0 };
141 }
142 return .{ trig.sin(x, 0.0, 0), trig.cos(x, 0.0) };
143 }
144
145 // sincos(Inf or NaN) is NaN
146 if (ix >= 0x7ff00000) {
147 const result = x - x;
148 return .{ result, result };
149 }
150
151 // argument reduction needed
152 var y: [2]f64 = undefined;
153 const n = rem_pio2(x, &y);
154 const s = trig.sin(y[0], y[1], 1);
155 const c = trig.cos(y[0], y[1]);
156 return switch (@as(u2, @truncate(@as(u32, @bitCast(n))))) {
157 0 => .{ s, c },
158 1 => .{ c, -s },
159 2 => .{ -s, -c },
160 3 => .{ -c, s },
161 };
162}
163
164fn sincosx(x: compiler_rt.f80.Abi, r_sin: *compiler_rt.f80.Abi, r_cos: *compiler_rt.f80.Abi) callconv(.c) void {
165 const s, const c = sincos_f80(compiler_rt.f80.fromAbi(x));
166 r_sin.* = compiler_rt.f80.toAbi(s);
167 r_cos.* = compiler_rt.f80.toAbi(c);
168}
169pub fn sincos_f80(x: f80) struct { f80, f80 } {
170 const se = ld.signExponent(x) & 0x7fff;
171 if (se == 0x7fff) {
172 const result = x - x;
173 return .{ result, result };
174 }
175
176 if (@abs(x) < trig.pi_4) {
177 if (se < 0x3fff - math.floatMantissaBits(f80)) {
178 // raise underflow if subnormal
179 if (compiler_rt.want_float_exceptions and se == 0) {
180 mem.doNotOptimizeAway(x * 0x1p-120);
181 }
182 // raise inexact if x!=0
183 return .{ x, 1.0 + x };
184 }
185 return .{ trig.sinx(x, 0.0, 0), trig.cosx(x, 0.0) };
186 }
187
188 var y: [2]f80 = undefined;
189 const n = rem_pio2l(f80, x, &y);
190 const s = trig.sinx(y[0], y[1], 1);
191 const c = trig.cosx(y[0], y[1]);
192 return switch (@as(u2, @truncate(@as(u32, @bitCast(n))))) {
193 0 => .{ s, c },
194 1 => .{ c, -s },
195 2 => .{ -s, -c },
196 3 => .{ -c, s },
197 };
198}
199
200fn sincosq(x: compiler_rt.f128.Abi, r_sin: *compiler_rt.f128.Abi, r_cos: *compiler_rt.f128.Abi) callconv(.c) void {
201 const s, const c = sincos_f128(compiler_rt.f128.fromAbi(x));
202 r_sin.* = compiler_rt.f128.toAbi(s);
203 r_cos.* = compiler_rt.f128.toAbi(c);
204}
205pub fn sincos_f128(x: f128) struct { f128, f128 } {
206 const se = ld.signExponent(x) & 0x7fff;
207 if (se == 0x7fff) {
208 const result = x - x;
209 return .{ result, result };
210 }
211
212 if (@abs(x) < trig.pi_4) {
213 if (se < 0x3fff - math.floatMantissaBits(f128)) {
214 // raise underflow if subnormal
215 if (compiler_rt.want_float_exceptions and se == 0) {
216 mem.doNotOptimizeAway(x * 0x1p-120);
217 }
218 // raise inexact if x!=0
219 return .{ x, 1.0 + x };
220 }
221 return .{ trig.sinq(x, 0.0, 0), trig.cosq(x, 0.0) };
222 }
223
224 var y: [2]f128 = undefined;
225 const n = rem_pio2l(f128, x, &y);
226 const s = trig.sinq(y[0], y[1], 1);
227 const c = trig.cosq(y[0], y[1]);
228 return switch (@as(u2, @truncate(@as(u32, @bitCast(n))))) {
229 0 => .{ s, c },
230 1 => .{ c, -s },
231 2 => .{ -s, -c },
232 3 => .{ -c, s },
233 };
234}
235
236pub fn sincosl(x: c_longdouble, r_sin: *c_longdouble, r_cos: *c_longdouble) callconv(.c) void {
237 r_sin.*, r_cos.* = switch (@typeInfo(c_longdouble).float.bits) {
238 64 => sincos_f64(x),
239 80 => sincos_f80(x),
240 128 => sincos_f128(x),
241 else => comptime unreachable,
242 };
243}
244
245fn testSincosSpecial(comptime T: type) !void {
246 const f = switch (T) {
247 f16 => sincos_f16,
248 f32 => sincos_f32,
249 f64 => sincos_f64,
250 f80 => sincos_f80,
251 f128 => sincos_f128,
252 else => @compileError("unimplemented"),
253 };
254
255 var s: T = undefined;
256 var c: T = undefined;
257
258 s, c = f(0.0);
259 try expect(math.isPositiveZero(s));
260 try expect(c == 1.0);
261
262 s, c = f(-0.0);
263 try expect(math.isNegativeZero(s));
264 try expect(c == 1.0);
265
266 s, c = f(math.inf(T));
267 try expect(math.isNan(s));
268 try expect(math.isNan(c));
269
270 s, c = f(-math.inf(T));
271 try expect(math.isNan(s));
272 try expect(math.isNan(c));
273
274 s, c = f(math.nan(T));
275 try expect(math.isNan(s));
276 try expect(math.isNan(c));
277}
278
279test "sincos32.normal" {
280 const epsilon = math.floatEps(f32);
281 var s: f32 = undefined;
282 var c: f32 = undefined;
283
284 s, c = sincos_f32(0.0);
285 try expectApproxEqAbs(@as(f32, 0.0), s, epsilon);
286 try expectApproxEqAbs(@as(f32, 1.0), c, epsilon);
287
288 s, c = sincos_f32(0.2);
289 try expectApproxEqAbs(@as(f32, 0.19866933), s, epsilon);
290 try expectApproxEqAbs(@as(f32, 0.9800666), c, epsilon);
291
292 s, c = sincos_f32(0.8923);
293 try expectApproxEqAbs(@as(f32, 0.77851737), s, epsilon);
294 try expectApproxEqAbs(@as(f32, 0.6276231), c, epsilon);
295
296 s, c = sincos_f32(1.5);
297 try expectApproxEqAbs(@as(f32, 0.997495), s, epsilon);
298 try expectApproxEqAbs(@as(f32, 0.0707372), c, epsilon);
299
300 s, c = sincos_f32(-1.5);
301 try expectApproxEqAbs(@as(f32, -0.997495), s, epsilon);
302 try expectApproxEqAbs(@as(f32, 0.0707372), c, epsilon);
303
304 s, c = sincos_f32(37.45);
305 try expectApproxEqAbs(@as(f32, -0.24654257), s, epsilon);
306 try expectApproxEqAbs(@as(f32, 0.96913195), c, epsilon);
307
308 s, c = sincos_f32(89.123);
309 try expectApproxEqAbs(@as(f32, 0.9161657), s, epsilon);
310 try expectApproxEqAbs(@as(f32, 0.40079966), c, epsilon);
311}
312
313test "sincos32.special" {
314 try testSincosSpecial(f32);
315}
316
317test "sincos64.normal" {
318 const epsilon = math.floatEps(f64);
319 var s: f64 = undefined;
320 var c: f64 = undefined;
321
322 s, c = sincos_f64(0.0);
323 try expectApproxEqAbs(@as(f64, 0.0), s, epsilon);
324 try expectApproxEqAbs(@as(f64, 1.0), c, epsilon);
325
326 s, c = sincos_f64(0.2);
327 try expectApproxEqAbs(@as(f64, 0.19866933079506122), s, epsilon);
328 try expectApproxEqAbs(@as(f64, 0.9800665778412416), c, epsilon);
329
330 s, c = sincos_f64(0.8923);
331 try expectApproxEqAbs(@as(f64, 0.7785173385577349), s, epsilon);
332 try expectApproxEqAbs(@as(f64, 0.6276230983360804), c, epsilon);
333
334 s, c = sincos_f64(1.5);
335 try expectApproxEqAbs(@as(f64, 0.9974949866040544), s, epsilon);
336 try expectApproxEqAbs(@as(f64, 0.0707372016677029), c, epsilon);
337
338 s, c = sincos_f64(-1.5);
339 try expectApproxEqAbs(@as(f64, -0.9974949866040544), s, epsilon);
340 try expectApproxEqAbs(@as(f64, 0.0707372016677029), c, epsilon);
341
342 s, c = sincos_f64(37.45);
343 try expectApproxEqAbs(@as(f64, -0.24654331551411082), s, epsilon);
344 try expectApproxEqAbs(@as(f64, 0.9691317730707778), c, epsilon);
345
346 s, c = sincos_f64(89.123);
347 try expectApproxEqAbs(@as(f64, 0.9161652766622714), s, epsilon);
348 try expectApproxEqAbs(@as(f64, 0.4008006809354791), c, epsilon);
349}
350
351test "sincos64.special" {
352 try testSincosSpecial(f64);
353}
354
355test "sincos80.normal" {
356 const epsilon = math.floatEps(f80);
357 var s: f80 = undefined;
358 var c: f80 = undefined;
359
360 s, c = sincos_f80(0.0);
361 try expectApproxEqAbs(@as(f80, 0.0), s, epsilon);
362 try expectApproxEqAbs(@as(f80, 1.0), c, epsilon);
363
364 s, c = sincos_f80(0.2);
365 try expectApproxEqAbs(@as(f80, 0.19866933079506121545941262711838975), s, epsilon);
366 try expectApproxEqAbs(@as(f80, 0.98006657784124163112419651674816888), c, epsilon);
367
368 s, c = sincos_f80(0.8923);
369 try expectApproxEqAbs(@as(f80, 0.77851733855773487830689285621486050), s, epsilon);
370 try expectApproxEqAbs(@as(f80, 0.62762309833608037003563995939286067), c, epsilon);
371
372 s, c = sincos_f80(1.5);
373 try expectApproxEqAbs(@as(f80, 0.99749498660405443094172337114148732), s, epsilon);
374 try expectApproxEqAbs(@as(f80, 0.070737201667702910088189851434268747), c, epsilon);
375
376 s, c = sincos_f80(-1.5);
377 try expectApproxEqAbs(@as(f80, -0.99749498660405443094172337114148732), s, epsilon);
378 try expectApproxEqAbs(@as(f80, 0.070737201667702910088189851434268747), c, epsilon);
379
380 s, c = sincos_f80(37.45);
381 try expectApproxEqAbs(@as(f80, -0.24654331551411356504), s, epsilon);
382 try expectApproxEqAbs(@as(f80, 0.9691317730707771246), c, epsilon);
383
384 s, c = sincos_f80(89.123);
385 try expectApproxEqAbs(@as(f80, 0.91616527666226951006), s, epsilon);
386 try expectApproxEqAbs(@as(f80, 0.4008006809354834001), c, epsilon);
387}
388
389test "sincos80.special" {
390 try testSincosSpecial(f80);
391}
392
393test "sincos128.normal" {
394 const epsilon = math.floatEps(f128);
395 var s: f128 = undefined;
396 var c: f128 = undefined;
397
398 s, c = sincos_f128(0.0);
399 try expectApproxEqAbs(@as(f128, 0.0), s, epsilon);
400 try expectApproxEqAbs(@as(f128, 1.0), c, epsilon);
401
402 s, c = sincos_f128(0.2);
403 try expectApproxEqAbs(@as(f128, 0.19866933079506121545941262711838975), s, epsilon);
404 try expectApproxEqAbs(@as(f128, 0.98006657784124163112419651674816888), c, epsilon);
405
406 s, c = sincos_f128(0.8923);
407 try expectApproxEqAbs(@as(f128, 0.77851733855773487830689285621486050), s, epsilon);
408 try expectApproxEqAbs(@as(f128, 0.62762309833608037003563995939286067), c, epsilon);
409
410 s, c = sincos_f128(1.5);
411 try expectApproxEqAbs(@as(f128, 0.99749498660405443094172337114148732), s, epsilon);
412 try expectApproxEqAbs(@as(f128, 0.070737201667702910088189851434268747), c, epsilon);
413
414 s, c = sincos_f128(-1.5);
415 try expectApproxEqAbs(@as(f128, -0.99749498660405443094172337114148732), s, epsilon);
416 try expectApproxEqAbs(@as(f128, 0.070737201667702910088189851434268747), c, epsilon);
417
418 s, c = sincos_f128(37.45);
419 try expectApproxEqAbs(@as(f128, -0.24654331551411356571238581321661085), s, epsilon);
420 try expectApproxEqAbs(@as(f128, 0.96913177307077712443149563847233230), c, epsilon);
421
422 s, c = sincos_f128(89.123);
423 try expectApproxEqAbs(@as(f128, 0.91616527666226951075019849560482170), s, epsilon);
424 try expectApproxEqAbs(@as(f128, 0.40080068093548339848199454493704702), c, epsilon);
425}
426
427test "sincos128.special" {
428 try testSincosSpecial(f128);
429}