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/cosf.c
5//! https://git.musl-libc.org/cgit/musl/tree/src/math/cos.c
6//! https://git.musl-libc.org/cgit/musl/tree/src/math/cosl.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(&__cosh, "__cosh");
24 symbol(&cosf, "cosf");
25 symbol(&cos, "cos");
26 symbol(&__cosx, "__cosx");
27 symbol(&cosq, "cosf128");
28 symbol(&cosl, "cosl");
29 symbol(&cosl, "__cosl"); // required by musl
30}
31
32fn __cosh(x: compiler_rt.f16.Abi) callconv(.c) compiler_rt.f16.Abi {
33 return compiler_rt.f16.toAbi(cos_f16(compiler_rt.f16.fromAbi(x)));
34}
35pub fn cos_f16(x: f16) f16 {
36 // TODO: more efficient implementation
37 return @floatCast(cos_f32(x));
38}
39
40fn cosf(x: compiler_rt.f32.Abi) callconv(.c) compiler_rt.f32.Abi {
41 return compiler_rt.f32.toAbi(cos_f32(compiler_rt.f32.fromAbi(x)));
42}
43pub fn cos_f32(x: f32) f32 {
44 // Small multiples of pi/2 rounded to double precision.
45 const c1pio2: f64 = 1.0 * math.pi / 2.0; // 0x3FF921FB, 0x54442D18
46 const c2pio2: f64 = 2.0 * math.pi / 2.0; // 0x400921FB, 0x54442D18
47 const c3pio2: f64 = 3.0 * math.pi / 2.0; // 0x4012D97C, 0x7F3321D2
48 const c4pio2: 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
57 if (compiler_rt.want_float_exceptions) mem.doNotOptimizeAway(x + 0x1p120);
58 return 1.0;
59 }
60 return trig.cosdf(x);
61 }
62 if (ix <= 0x407b53d1) { // |x| ~<= 5*pi/4
63 if (ix > 0x4016cbe3) { // |x| ~> 3*pi/4
64 return -trig.cosdf(if (sign) x + c2pio2 else x - c2pio2);
65 } else {
66 if (sign) {
67 return trig.sindf(x + c1pio2);
68 } else {
69 return trig.sindf(c1pio2 - x);
70 }
71 }
72 }
73 if (ix <= 0x40e231d5) { // |x| ~<= 9*pi/4
74 if (ix > 0x40afeddf) { // |x| ~> 7*pi/4
75 return trig.cosdf(if (sign) x + c4pio2 else x - c4pio2);
76 } else {
77 if (sign) {
78 return trig.sindf(-x - c3pio2);
79 } else {
80 return trig.sindf(x - c3pio2);
81 }
82 }
83 }
84
85 // cos(Inf or NaN) is NaN
86 if (ix >= 0x7f800000) {
87 return x - x;
88 }
89
90 var y: f64 = undefined;
91 const n = rem_pio2f(x, &y);
92 return switch (n & 3) {
93 0 => trig.cosdf(y),
94 1 => trig.sindf(-y),
95 2 => -trig.cosdf(y),
96 else => trig.sindf(y),
97 };
98}
99
100fn cos(x: compiler_rt.f64.Abi) callconv(.c) compiler_rt.f64.Abi {
101 return compiler_rt.f64.toAbi(cos_f64(compiler_rt.f64.fromAbi(x)));
102}
103pub fn cos_f64(x: f64) f64 {
104 var ix = @as(u64, @bitCast(x)) >> 32;
105 ix &= 0x7fffffff;
106
107 // |x| ~< pi/4
108 if (ix <= 0x3fe921fb) {
109 if (ix < 0x3e46a09e) { // |x| < 2**-27 * sqrt(2)
110 // raise inexact if x!=0
111 if (compiler_rt.want_float_exceptions) mem.doNotOptimizeAway(x + 0x1p120);
112 return 1.0;
113 }
114 return trig.cos(x, 0);
115 }
116
117 // cos(Inf or NaN) is NaN
118 if (ix >= 0x7ff00000) {
119 return x - x;
120 }
121
122 var y: [2]f64 = undefined;
123 const n = rem_pio2(x, &y);
124 return switch (n & 3) {
125 0 => trig.cos(y[0], y[1]),
126 1 => -trig.sin(y[0], y[1], 1),
127 2 => -trig.cos(y[0], y[1]),
128 else => trig.sin(y[0], y[1], 1),
129 };
130}
131
132fn __cosx(x: compiler_rt.f80.Abi) callconv(.c) compiler_rt.f80.Abi {
133 return compiler_rt.f80.toAbi(cos_f80(compiler_rt.f80.fromAbi(x)));
134}
135pub fn cos_f80(x: f80) f80 {
136 const se = ld.signExponent(x) & 0x7fff;
137 if (se == 0x7fff) {
138 return x - x;
139 }
140
141 if (@abs(x) < trig.pi_4) {
142 if (se < 0x3fff - math.floatMantissaBits(f80)) {
143 // raise inexact if x!=0
144 return 1.0 + x;
145 }
146 return trig.cosx(x, 0.0);
147 }
148
149 var y: [2]f80 = undefined;
150 const n = rem_pio2l(f80, x, &y);
151 return switch (n & 3) {
152 0 => trig.cosx(y[0], y[1]),
153 1 => -trig.sinx(y[0], y[1], 1),
154 2 => -trig.cosx(y[0], y[1]),
155 else => trig.sinx(y[0], y[1], 1),
156 };
157}
158
159fn cosq(x: compiler_rt.f128.Abi) callconv(.c) compiler_rt.f128.Abi {
160 return compiler_rt.f128.toAbi(cos_f128(compiler_rt.f128.fromAbi(x)));
161}
162pub fn cos_f128(x: f128) f128 {
163 const se = ld.signExponent(x) & 0x7fff;
164 if (se == 0x7fff) {
165 return x - x;
166 }
167
168 if (@abs(x) < trig.pi_4) {
169 if (se < 0x3fff - math.floatMantissaBits(f128)) {
170 // raise inexact if x!=0
171 return 1.0 + x;
172 }
173 return trig.cosq(x, 0.0);
174 }
175
176 var y: [2]f128 = undefined;
177 const n = rem_pio2l(f128, x, &y);
178 return switch (n & 3) {
179 0 => trig.cosq(y[0], y[1]),
180 1 => -trig.sinq(y[0], y[1], 1),
181 2 => -trig.cosq(y[0], y[1]),
182 else => trig.sinq(y[0], y[1], 1),
183 };
184}
185
186pub fn cosl(x: c_longdouble) callconv(.c) c_longdouble {
187 switch (@typeInfo(c_longdouble).float.bits) {
188 64 => return cos_f64(x),
189 80 => return cos_f80(x),
190 128 => return cos_f128(x),
191 else => comptime unreachable,
192 }
193}
194
195fn testCosSpecial(comptime T: type) !void {
196 const f = switch (T) {
197 f16 => cos_f16,
198 f32 => cos_f32,
199 f64 => cos_f64,
200 f80 => cos_f80,
201 f128 => cos_f128,
202 else => comptime unreachable,
203 };
204
205 try expect(f(0.0) == 1.0);
206 try expect(f(-0.0) == 1.0);
207 try expect(math.isNan(f(math.inf(T))));
208 try expect(math.isNan(f(-math.inf(T))));
209 try expect(math.isNan(f(math.nan(T))));
210}
211
212test "cos32.normal" {
213 const epsilon = math.floatEps(f32);
214 try expectApproxEqAbs(@as(f32, 1.0), cos_f32(0.0), epsilon);
215 try expectApproxEqAbs(@as(f32, 0.9800666), cos_f32(0.2), epsilon);
216 try expectApproxEqAbs(@as(f32, 0.6276231), cos_f32(0.8923), epsilon);
217 try expectApproxEqAbs(@as(f32, 0.0707372), cos_f32(1.5), epsilon);
218 try expectApproxEqAbs(@as(f32, 0.0707372), cos_f32(-1.5), epsilon);
219 try expectApproxEqAbs(@as(f32, 0.96913195), cos_f32(37.45), epsilon);
220 try expectApproxEqAbs(@as(f32, 0.40079966), cos_f32(89.123), epsilon);
221}
222
223test "cos32.special" {
224 try testCosSpecial(f32);
225}
226
227test "cos64.normal" {
228 const epsilon = math.floatEps(f64);
229 try expectApproxEqAbs(@as(f64, 1.0), cos_f64(0.0), epsilon);
230 try expectApproxEqAbs(@as(f64, 0.9800665778412416), cos_f64(0.2), epsilon);
231 try expectApproxEqAbs(@as(f64, 0.6276230983360804), cos_f64(0.8923), epsilon);
232 try expectApproxEqAbs(@as(f64, 0.0707372016677029), cos_f64(1.5), epsilon);
233 try expectApproxEqAbs(@as(f64, 0.0707372016677029), cos_f64(-1.5), epsilon);
234 try expectApproxEqAbs(@as(f64, 0.9691317730707778), cos_f64(37.45), epsilon);
235 try expectApproxEqAbs(@as(f64, 0.4008006809354791), cos_f64(89.123), epsilon);
236}
237
238test "cos64.special" {
239 try testCosSpecial(f64);
240}
241
242test "cos80.normal" {
243 const epsilon = math.floatEps(f80);
244 try expectApproxEqAbs(@as(f80, 1.0), cos_f80(0.0), epsilon);
245 try expectApproxEqAbs(@as(f80, 0.98006657784124163112419651674816888), cos_f80(0.2), epsilon);
246 try expectApproxEqAbs(@as(f80, 0.62762309833608037003563995939286067), cos_f80(0.8923), epsilon);
247 try expectApproxEqAbs(@as(f80, 0.070737201667702910088189851434268747), cos_f80(1.5), epsilon);
248 try expectApproxEqAbs(@as(f80, 0.070737201667702910088189851434268747), cos_f80(-1.5), epsilon);
249 try expectApproxEqAbs(@as(f80, 0.9691317730707771246), cos_f80(37.45), epsilon);
250 try expectApproxEqAbs(@as(f80, 0.4008006809354834001), cos_f80(89.123), epsilon);
251}
252
253test "cos80.special" {
254 try testCosSpecial(f80);
255}
256
257test "cos128.normal" {
258 const epsilon = math.floatEps(f128);
259 try expectApproxEqAbs(@as(f128, 1.0), cos_f128(0.0), epsilon);
260 try expectApproxEqAbs(@as(f128, 0.98006657784124163112419651674816888), cos_f128(0.2), epsilon);
261 try expectApproxEqAbs(@as(f128, 0.62762309833608037003563995939286067), cos_f128(0.8923), epsilon);
262 try expectApproxEqAbs(@as(f128, 0.070737201667702910088189851434268747), cos_f128(1.5), epsilon);
263 try expectApproxEqAbs(@as(f128, 0.070737201667702910088189851434268747), cos_f128(-1.5), epsilon);
264 try expectApproxEqAbs(@as(f128, 0.96913177307077712443149563847233230), cos_f128(37.45), epsilon);
265 try expectApproxEqAbs(@as(f128, 0.40080068093548339848199454493704702), cos_f128(89.123), epsilon);
266}
267
268test "cos128.special" {
269 try testCosSpecial(f128);
270}