authorgravatar for hi@mihaelm.commihael <hi@mihaelm.com> 2026-03-24 01:11:56+01:00
committergravatar for hi@mihaelm.commihael <hi@mihaelm.com> 2026-04-02 23:54:20+02:00
logffd6f6cc6e05e84a4e029b02ff13a3e84cbf1aa4
tree532cfa00e01532ba84d41e5375682c645be06e90
parentc76644caf3051dd2f837a6c40987a64e01b09c54
signaturebadge-check Signed by SSH key SHA256:aoFoShdYLdrqMichqKXSSieTKUfACUIDJHsKc4V2tQg

`libzigc/math`: Implement more precise `tanl` in `compiler_rt`

The logic was more or less ported from `musl`, with small adjustments where it was convenient. The 'internal' `__tanl` function was implemented in the `trig.zig` module along with other 'internal' trigonometric functions. Now, the `tanl` implementation is precise enough to pass all the relevant `libc-test` suite tests: ``` $ ./build/stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib $ stage4/bin/zig build test-libc -Dlibc-test-path=<LIBC-TEST-PATH> -Dtest-filter=tanl -fqemu -fwasmtime --summary line Build Summary: 553/553 steps succeeded ``` The unit tests were also extended to include cases for `f80` and `f128`, and they're passing.

2 files changed, 223 insertions(+), 33 deletions(-)

lib/compiler_rt/tan.zig+96-33
......@@ -3,6 +3,7 @@
33//!
44//! https://git.musl-libc.org/cgit/musl/tree/src/math/tanf.c
55//! 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
67//! https://golang.org/src/math/tan.go
78
89const std = @import("std");
......@@ -10,10 +11,13 @@ const builtin = @import("builtin");
1011const math = std.math;
1112const mem = std.mem;
1213const expect = std.testing.expect;
14const expectApproxEqAbs = std.testing.expectApproxEqAbs;
1315
1416const kernel = @import("trig.zig");
1517const rem_pio2 = @import("rem_pio2.zig").rem_pio2;
1618const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f;
19const rem_pio2l = @import("rem_pio2l.zig").rem_pio2l;
20const utils = @import("math_utils.zig");
1721
1822const arch = builtin.cpu.arch;
1923const compiler_rt = @import("../compiler_rt.zig");
......@@ -116,14 +120,38 @@ pub fn tan(x: f64) callconv(.c) f64 {
116120 return kernel.__tan(y[0], y[1], n & 1 != 0);
117121}
118122
123fn tanlGeneric(comptime T: type, x: T) T {
124 if (!(T == f80 or T == f128)) {
125 @compileError("`tanlGeneric` implemented only for `f80` and `f128`, got: " ++ T);
126 }
127
128 const se = utils.ldSignExponent(x) & 0x7fff;
129 if (se == 0x7fff) {
130 return x - x;
131 }
132
133 const pi_4 = 0.78539816339744830962;
134 if (@abs(x) < pi_4) {
135 if (se < 0x3fff - math.floatMantissaBits(T) / 2) {
136 if (compiler_rt.want_float_exceptions) {
137 mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120);
138 }
139 return x;
140 }
141 return kernel.__tanl(T, x, 0.0, 0);
142 }
143
144 var y: [2]T = undefined;
145 const n = rem_pio2l(T, x, &y);
146 return kernel.__tanl(T, y[0], y[1], n & 1);
147}
148
119149pub fn __tanx(x: f80) callconv(.c) f80 {
120 // TODO: more efficient implementation
121 return @floatCast(tanq(x));
150 return tanlGeneric(f80, x);
122151}
123152
124153pub fn tanq(x: f128) callconv(.c) f128 {
125 // TODO: more correct implementation
126 return tan(@floatCast(x));
154 return tanlGeneric(f128, x);
127155}
128156
129157pub fn tanl(x: c_longdouble) callconv(.c) c_longdouble {
......@@ -137,45 +165,80 @@ pub fn tanl(x: c_longdouble) callconv(.c) c_longdouble {
137165 }
138166}
139167
140test "tan" {
141 try expect(tan(@as(f32, 0.0)) == tanf(0.0));
142 try expect(tan(@as(f64, 0.0)) == tan(0.0));
168fn testTanNormal(comptime T: type) !void {
169 const f = switch (T) {
170 f32 => tanf,
171 f64 => tan,
172 else => @compileError("unimplemented"),
173 };
174 const epsilon = 0.00001;
175
176 try expectApproxEqAbs(@as(T, 0.0), f(0.0), epsilon);
177 try expectApproxEqAbs(@as(T, 0.202710), f(0.2), epsilon);
178 try expectApproxEqAbs(@as(T, 1.240422), f(0.8923), epsilon);
179 try expectApproxEqAbs(@as(T, 14.101420), f(1.5), epsilon);
180 try expectApproxEqAbs(@as(T, -0.254397), f(37.45), epsilon);
181 try expectApproxEqAbs(@as(T, 2.285837), f(89.123), epsilon);
143182}
144183
145test "tan32" {
146 const epsilon = 0.00001;
184fn testTanSpecial(comptime T: type) !void {
185 const f = switch (T) {
186 f32 => tanf,
187 f64 => tan,
188 f80 => __tanx,
189 f128 => tanq,
190 else => @compileError("unimplemented"),
191 };
192
193 try expect(math.isPositiveZero(f(0.0)));
194 try expect(math.isNegativeZero(f(-0.0)));
195 try expect(math.isNan(f(math.inf(f32))));
196 try expect(math.isNan(f(-math.inf(f32))));
197 try expect(math.isNan(f(math.nan(f32))));
198}
147199
148 try expect(math.approxEqAbs(f32, tanf(0.0), 0.0, epsilon));
149 try expect(math.approxEqAbs(f32, tanf(0.2), 0.202710, epsilon));
150 try expect(math.approxEqAbs(f32, tanf(0.8923), 1.240422, epsilon));
151 try expect(math.approxEqAbs(f32, tanf(1.5), 14.101420, epsilon));
152 try expect(math.approxEqAbs(f32, tanf(37.45), -0.254397, epsilon));
153 try expect(math.approxEqAbs(f32, tanf(89.123), 2.285852, epsilon));
200test "tan32.normal" {
201 try testTanNormal(f32);
154202}
155203
156test "tan64" {
157 const epsilon = 0.000001;
204test "tan64.normal" {
205 try testTanNormal(f64);
206}
207
208test "tan80.normal" {
209 const epsilon = math.floatEps(f80);
158210
159 try expect(math.approxEqAbs(f64, tan(0.0), 0.0, epsilon));
160 try expect(math.approxEqAbs(f64, tan(0.2), 0.202710, epsilon));
161 try expect(math.approxEqAbs(f64, tan(0.8923), 1.240422, epsilon));
162 try expect(math.approxEqAbs(f64, tan(1.5), 14.101420, epsilon));
163 try expect(math.approxEqAbs(f64, tan(37.45), -0.254397, epsilon));
164 try expect(math.approxEqAbs(f64, tan(89.123), 2.2858376, epsilon));
211 try expectApproxEqAbs(@as(f80, 0.0), __tanx(0.0), epsilon);
212 try expectApproxEqAbs(@as(f80, 0.2027100355086724833213582716475345), __tanx(0.2), epsilon);
213 try expectApproxEqAbs(@as(f80, 1.2404217445497097995561220131857544), __tanx(0.8923), epsilon);
214 try expectApproxEqAbs(@as(f80, 14.10141994717171938764), __tanx(1.5), epsilon);
215 try expectApproxEqAbs(@as(f80, -0.25439607116885656232), __tanx(37.45), epsilon);
216 try expectApproxEqAbs(@as(f80, 2.2858376251355320963), __tanx(89.123), epsilon);
217}
218
219test "tan128.normal" {
220 const epsilon = math.floatEps(f128);
221
222 try expectApproxEqAbs(@as(f128, 0.0), tanq(0.0), epsilon);
223 try expectApproxEqAbs(@as(f128, 0.2027100355086724833213582716475345), tanq(0.2), epsilon);
224 try expectApproxEqAbs(@as(f128, 1.2404217445497097995561220131857544), tanq(0.8923), epsilon);
225 try expectApproxEqAbs(@as(f128, 14.101419947171719387646083651987755), tanq(1.5), epsilon);
226 try expectApproxEqAbs(@as(f128, -0.2543960711688565630469573224504774), tanq(37.45), epsilon);
227 try expectApproxEqAbs(@as(f128, 2.2858376251355321074066028114094292), tanq(89.123), epsilon);
165228}
166229
167230test "tan32.special" {
168 try expect(tanf(0.0) == 0.0);
169 try expect(tanf(-0.0) == -0.0);
170 try expect(math.isNan(tanf(math.inf(f32))));
171 try expect(math.isNan(tanf(-math.inf(f32))));
172 try expect(math.isNan(tanf(math.nan(f32))));
231 try testTanSpecial(f32);
173232}
174233
175234test "tan64.special" {
176 try expect(tan(0.0) == 0.0);
177 try expect(tan(-0.0) == -0.0);
178 try expect(math.isNan(tan(math.inf(f64))));
179 try expect(math.isNan(tan(-math.inf(f64))));
180 try expect(math.isNan(tan(math.nan(f64))));
235 try testTanSpecial(f64);
236}
237
238test "tan80.special" {
239 try testTanSpecial(f80);
240}
241
242test "tan128.special" {
243 try testTanSpecial(f128);
181244}
lib/compiler_rt/trig.zig+127
......@@ -7,6 +7,7 @@
77// https://git.musl-libc.org/cgit/musl/tree/src/math/__sindf.c
88// https://git.musl-libc.org/cgit/musl/tree/src/math/__tand.c
99// https://git.musl-libc.org/cgit/musl/tree/src/math/__tandf.c
10// https://git.musl-libc.org/cgit/musl/tree/src/math/__tanl.c
1011
1112/// kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
1213/// Input x is assumed to be bounded by ~pi/4 in magnitude.
......@@ -271,3 +272,129 @@ pub fn __tandf(x: f64, odd: bool) f32 {
271272 const r0 = (x + s * u) + (s * w) * (t + w * r);
272273 return @floatCast(if (odd) -1.0 / r0 else r0);
273274}
275
276pub fn __tanl(comptime T: type, x_: T, y_: T, odd: i32) T {
277 var x = x_;
278 var y = y_;
279 const impl = switch (T) {
280 f80 => struct {
281 const pio4: T = 0.785398163397448309628;
282 const pio4lo: T = -1.25413940316708300586e-20;
283
284 const T3: T = 0.333333333333333333180;
285 const T5: T = 0.133333333333333372290;
286 const T7: T = 0.0539682539682504975744;
287 const T9: f64 = 0.021869488536312216;
288 const T11: f64 = 0.0088632355256619590;
289 const T13: f64 = 0.0035921281113786528;
290 const T15: f64 = 0.0014558334756312418;
291 const T17: f64 = 0.00059003538700862256;
292 const T19: f64 = 0.00023907843576635544;
293 const T21: f64 = 0.000097154625656538905;
294 const T23: f64 = 0.000038440165747303162;
295 const T25: f64 = 0.000018082171885432524;
296 const T27: f64 = 0.0000024196006108814377;
297 const T29: f64 = 0.0000078293456938132840;
298 const T31: f64 = -0.0000032609076735050182;
299 const T33: f64 = 0.0000023261313142559411;
300
301 inline fn rpoly(w: T) T {
302 return T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 +
303 w * (T25 + w * (T29 + w * T33))))));
304 }
305
306 inline fn vpoly(w: T) T {
307 return T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 +
308 w * (T27 + w * T31)))));
309 }
310 },
311 f128 => struct {
312 const pio4: T = 0x1.921fb54442d18469898cc51701b8p-1;
313 const pio4lo: T = 0x1.cd129024e088a67cc74020bbea60p-116;
314
315 const T3: T = 0x1.5555555555555555555555555553p-2;
316 const T5: T = 0x1.1111111111111111111111111eb5p-3;
317 const T7: T = 0x1.ba1ba1ba1ba1ba1ba1ba1b694cd6p-5;
318 const T9: T = 0x1.664f4882c10f9f32d6bbe09d8bcdp-6;
319 const T11: T = 0x1.226e355e6c23c8f5b4f5762322eep-7;
320 const T13: T = 0x1.d6d3d0e157ddfb5fed8e84e27b37p-9;
321 const T15: T = 0x1.7da36452b75e2b5fce9ee7c2c92ep-10;
322 const T17: T = 0x1.355824803674477dfcf726649efep-11;
323 const T19: T = 0x1.f57d7734d1656e0aceb716f614c2p-13;
324 const T21: T = 0x1.967e18afcb180ed942dfdc518d6cp-14;
325 const T23: T = 0x1.497d8eea21e95bc7e2aa79b9f2cdp-15;
326 const T25: T = 0x1.0b132d39f055c81be49eff7afd50p-16;
327 const T27: T = 0x1.b0f72d33eff7bfa2fbc1059d90b6p-18;
328 const T29: T = 0x1.5ef2daf21d1113df38d0fbc00267p-19;
329 const T31: T = 0x1.1c77d6eac0234988cdaa04c96626p-20;
330 const T33: T = 0x1.cd2a5a292b180e0bdd701057dfe3p-22;
331 const T35: T = 0x1.75c7357d0298c01a31d0a6f7d518p-23;
332 const T37: T = 0x1.2f3190f4718a9a520f98f50081fcp-24;
333 const T39: f64 = 0.000000028443389121318352;
334 const T41: f64 = 0.000000011981013102001973;
335 const T43: f64 = 0.0000000038303578044958070;
336 const T45: f64 = 0.0000000034664378216909893;
337 const T47: f64 = -0.0000000015090641701997785;
338 const T49: f64 = 0.0000000029449552300483952;
339 const T51: f64 = -0.0000000022006995706097711;
340 const T53: f64 = 0.0000000015468200913196612;
341 const T55: f64 = -0.00000000061311613386849674;
342 const T57: f64 = 1.4912469681508012e-10;
343
344 inline fn rpoly(w: T) T {
345 return T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 +
346 w * (T25 + w * (T29 + w * (T33 + w * (T37 + w * (T41 +
347 w * (T45 + w * (T49 + w * (T53 + w * T57))))))))))));
348 }
349
350 inline fn vpoly(w: T) T {
351 return T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 +
352 w * (T27 + w * (T31 + w * (T35 + w * (T39 + w * (T43 +
353 w * (T47 + w * (T51 + w * T55)))))))))));
354 }
355 },
356 else => @compileError("__tanl supports only f80 and f128, got: " ++ @typeName(T)),
357 };
358
359 const big = @abs(x) >= 0.67434;
360 var sign: i8 = 0;
361
362 if (big) {
363 if (x < 0) {
364 sign = -1;
365 x = -x;
366 y = -y;
367 }
368 x = (impl.pio4 - x) + (impl.pio4lo - y);
369 y = 0.0;
370 }
371
372 var z = x * x;
373 var w = z * z;
374 var r = impl.rpoly(w);
375 var v = z * impl.vpoly(w);
376 var s = z * x;
377 r = y + z * (s * (r + v) + y) + impl.T3 * s;
378 w = x + r;
379
380 if (big) {
381 s = @as(T, @floatFromInt(1 - 2 * odd));
382 v = s - 2.0 * (x + (r - w * w / (w + s)));
383 return if (sign == -1) -v else v;
384 }
385
386 if (odd == 0) {
387 return w;
388 }
389
390 // if allow error up to 2 ulp, simply return
391 // -1.0 / (x+r) here
392 //
393 // compute -1.0 / (x+r) accurately
394 z = w + 0x1p32 - 0x1p32;
395 v = r - (z - x);
396 const a = -1.0 / w;
397 const t = a + 0x1p32 - 0x1p32;
398 s = 1.0 + t * z;
399 return t + a * (s + t * v);
400}