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/tanf.c
5//! https://git.musl-libc.org/cgit/musl/tree/src/math/tan.c
6//! https://git.musl-libc.org/cgit/musl/tree/src/math/tanl.c
7//! https://golang.org/src/math/tan.go
8
9const std = @import("std");
10const builtin = @import("builtin");
11const math = std.math;
12const ld = math.long_double;
13const mem = std.mem;
14const expect = std.testing.expect;
15const expectApproxEqAbs = std.testing.expectApproxEqAbs;
16
17const kernel = @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
22const arch = builtin.cpu.arch;
23const compiler_rt = @import("../compiler_rt.zig");
24const symbol = compiler_rt.symbol;
25
26comptime {
27 symbol(&__tanh, "__tanh");
28 symbol(&tanf, "tanf");
29 symbol(&tan, "tan");
30 symbol(&__tanx, "__tanx");
31 symbol(&tanq, "tanf128");
32 symbol(&tanl, "tanl");
33}
34
35fn __tanh(x: compiler_rt.f16.Abi) callconv(.c) compiler_rt.f16.Abi {
36 return compiler_rt.f16.toAbi(tan_f16(compiler_rt.f16.fromAbi(x)));
37}
38pub fn tan_f16(x: f16) f16 {
39 // TODO: more efficient implementation
40 return @floatCast(tan_f32(x));
41}
42
43fn tanf(x: compiler_rt.f32.Abi) callconv(.c) compiler_rt.f32.Abi {
44 return compiler_rt.f32.toAbi(tan_f32(compiler_rt.f32.fromAbi(x)));
45}
46pub fn tan_f32(x: f32) f32 {
47 // Small multiples of pi/2 rounded to double precision.
48 const t1pio2: f64 = 1.0 * math.pi / 2.0; // 0x3FF921FB, 0x54442D18
49 const t2pio2: f64 = 2.0 * math.pi / 2.0; // 0x400921FB, 0x54442D18
50 const t3pio2: f64 = 3.0 * math.pi / 2.0; // 0x4012D97C, 0x7F3321D2
51 const t4pio2: f64 = 4.0 * math.pi / 2.0; // 0x401921FB, 0x54442D18
52
53 var ix: u32 = @bitCast(x);
54 const sign = ix >> 31 != 0;
55 ix &= 0x7fffffff;
56
57 if (ix <= 0x3f490fda) { // |x| ~<= pi/4
58 if (ix < 0x39800000) { // |x| < 2**-12
59 // raise inexact if x!=0 and underflow if subnormal
60 if (compiler_rt.want_float_exceptions) {
61 if (ix < 0x00800000) {
62 mem.doNotOptimizeAway(x / 0x1p120);
63 } else {
64 mem.doNotOptimizeAway(x + 0x1p120);
65 }
66 }
67 return x;
68 }
69 return kernel.tandf(x, false);
70 }
71 if (ix <= 0x407b53d1) { // |x| ~<= 5*pi/4
72 if (ix <= 0x4016cbe3) { // |x| ~<= 3pi/4
73 return kernel.tandf((if (sign) x + t1pio2 else x - t1pio2), true);
74 } else {
75 return kernel.tandf((if (sign) x + t2pio2 else x - t2pio2), false);
76 }
77 }
78 if (ix <= 0x40e231d5) { // |x| ~<= 9*pi/4
79 if (ix <= 0x40afeddf) { // |x| ~<= 7*pi/4
80 return kernel.tandf((if (sign) x + t3pio2 else x - t3pio2), true);
81 } else {
82 return kernel.tandf((if (sign) x + t4pio2 else x - t4pio2), false);
83 }
84 }
85
86 // tan(Inf or NaN) is NaN
87 if (ix >= 0x7f800000) {
88 return x - x;
89 }
90
91 var y: f64 = undefined;
92 const n = rem_pio2f(x, &y);
93 return kernel.tandf(y, n & 1 != 0);
94}
95
96fn tan(x: compiler_rt.f64.Abi) callconv(.c) compiler_rt.f64.Abi {
97 return compiler_rt.f64.toAbi(tan_f64(compiler_rt.f64.fromAbi(x)));
98}
99pub fn tan_f64(x: f64) f64 {
100 var ix = @as(u64, @bitCast(x)) >> 32;
101 ix &= 0x7fffffff;
102
103 // |x| ~< pi/4
104 if (ix <= 0x3fe921fb) {
105 if (ix < 0x3e400000) { // |x| < 2**-27
106 // raise inexact if x!=0 and underflow if subnormal
107 if (compiler_rt.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 kernel.tan(x, 0.0, false);
117 }
118
119 // tan(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 kernel.tan(y[0], y[1], n & 1 != 0);
127}
128
129fn __tanx(x: compiler_rt.f80.Abi) callconv(.c) compiler_rt.f80.Abi {
130 return compiler_rt.f80.toAbi(tan_f80(compiler_rt.f80.fromAbi(x)));
131}
132pub fn tan_f80(x: f80) f80 {
133 const se = ld.signExponent(x) & 0x7fff;
134 if (se == 0x7fff) {
135 return x - x;
136 }
137
138 if (@abs(x) < kernel.pi_4) {
139 if (se < 0x3fff - math.floatMantissaBits(f80) / 2) {
140 if (compiler_rt.want_float_exceptions) {
141 mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120);
142 }
143 return x;
144 }
145 return kernel.tanx(x, 0.0, 0);
146 }
147
148 var y: [2]f80 = undefined;
149 const n = rem_pio2l(f80, x, &y);
150 return kernel.tanx(y[0], y[1], n & 1);
151}
152
153fn tanq(x: compiler_rt.f128.Abi) callconv(.c) compiler_rt.f128.Abi {
154 return compiler_rt.f128.toAbi(tan_f128(compiler_rt.f128.fromAbi(x)));
155}
156pub fn tan_f128(x: f128) f128 {
157 const se = ld.signExponent(x) & 0x7fff;
158 if (se == 0x7fff) {
159 return x - x;
160 }
161
162 if (@abs(x) < kernel.pi_4) {
163 if (se < 0x3fff - math.floatMantissaBits(f128) / 2) {
164 if (compiler_rt.want_float_exceptions) {
165 mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120);
166 }
167 return x;
168 }
169 return kernel.tanq(x, 0.0, 0);
170 }
171
172 var y: [2]f128 = undefined;
173 const n = rem_pio2l(f128, x, &y);
174 return kernel.tanq(y[0], y[1], n & 1);
175}
176
177pub fn tanl(x: c_longdouble) callconv(.c) c_longdouble {
178 switch (@typeInfo(c_longdouble).float.bits) {
179 64 => return tan_f64(x),
180 80 => return tan_f80(x),
181 128 => return tan_f128(x),
182 else => comptime unreachable,
183 }
184}
185
186fn testTanNormal(comptime T: type) !void {
187 const f = switch (T) {
188 f16 => tan_f16,
189 f32 => tan_f32,
190 f64 => tan_f64,
191 f80 => tan_f80,
192 f128 => tan_f128,
193 else => comptime unreachable,
194 };
195 const epsilon = 0.00001;
196
197 try expectApproxEqAbs(@as(T, 0.0), f(0.0), epsilon);
198 try expectApproxEqAbs(@as(T, 0.202710), f(0.2), epsilon);
199 try expectApproxEqAbs(@as(T, 1.240422), f(0.8923), epsilon);
200 try expectApproxEqAbs(@as(T, 14.101420), f(1.5), epsilon);
201 try expectApproxEqAbs(@as(T, -0.254397), f(37.45), epsilon);
202 try expectApproxEqAbs(@as(T, 2.285837), f(89.123), epsilon);
203}
204
205fn testTanSpecial(comptime T: type) !void {
206 const f = switch (T) {
207 f16 => tan_f16,
208 f32 => tan_f32,
209 f64 => tan_f64,
210 f80 => tan_f80,
211 f128 => tan_f128,
212 else => comptime unreachable,
213 };
214
215 try expect(math.isPositiveZero(f(0.0)));
216 try expect(math.isNegativeZero(f(-0.0)));
217 try expect(math.isNan(f(math.inf(f32))));
218 try expect(math.isNan(f(-math.inf(f32))));
219 try expect(math.isNan(f(math.nan(f32))));
220}
221
222test "tan32.normal" {
223 try testTanNormal(f32);
224}
225
226test "tan64.normal" {
227 try testTanNormal(f64);
228}
229
230test "tan80.normal" {
231 const epsilon = math.floatEps(f80);
232
233 try expectApproxEqAbs(@as(f80, 0.0), tan_f80(0.0), epsilon);
234 try expectApproxEqAbs(@as(f80, 0.2027100355086724833213582716475345), tan_f80(0.2), epsilon);
235 try expectApproxEqAbs(@as(f80, 1.2404217445497097995561220131857544), tan_f80(0.8923), epsilon);
236 try expectApproxEqAbs(@as(f80, 14.10141994717171938764), tan_f80(1.5), epsilon);
237 try expectApproxEqAbs(@as(f80, -0.25439607116885656232), tan_f80(37.45), epsilon);
238 try expectApproxEqAbs(@as(f80, 2.2858376251355320963), tan_f80(89.123), epsilon);
239}
240
241test "tan128.normal" {
242 const epsilon = math.floatEps(f128);
243
244 try expectApproxEqAbs(@as(f128, 0.0), tan_f128(0.0), epsilon);
245 try expectApproxEqAbs(@as(f128, 0.2027100355086724833213582716475345), tan_f128(0.2), epsilon);
246 try expectApproxEqAbs(@as(f128, 1.2404217445497097995561220131857544), tan_f128(0.8923), epsilon);
247 try expectApproxEqAbs(@as(f128, 14.101419947171719387646083651987755), tan_f128(1.5), epsilon);
248 try expectApproxEqAbs(@as(f128, -0.2543960711688565630469573224504774), tan_f128(37.45), epsilon);
249 try expectApproxEqAbs(@as(f128, 2.2858376251355321074066028114094292), tan_f128(89.123), epsilon);
250}
251
252test "tan32.special" {
253 try testTanSpecial(f32);
254}
255
256test "tan64.special" {
257 try testTanSpecial(f64);
258}
259
260test "tan80.special" {
261 try testTanSpecial(f80);
262}
263
264test "tan128.special" {
265 try testTanSpecial(f128);
266}