commit 9a53be5fb444873ee561eb54a0ea1aaeaaad844d (tree)
parent a6d6af947d09fbd52662c96c6f0a12510ef15f70
Author: mihael <hi@mihaelm.com>
Date: Thu, 26 Mar 2026 21:43:55 +0100
`libzigc/math`: Implement more precise `sincosl` in `compiler_rt`
The changes were tested by running:
```bash
$ ./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=sincosl -fqemu -fwasmtime --summary line
Build Summary: 369/369 steps succeeded
```
Additionally, I've added unit tests for `sincos`.
Diffstat:
4 files changed, 223 insertions(+), 119 deletions(-)
diff --git a/lib/compiler_rt/sincos.zig b/lib/compiler_rt/sincos.zig
@@ -3,9 +3,13 @@ const builtin = @import("builtin");
const arch = builtin.cpu.arch;
const math = std.math;
const mem = std.mem;
+const expect = std.testing.expect;
+const expectApproxEqAbs = std.testing.expectApproxEqAbs;
const trig = @import("trig.zig");
const rem_pio2 = @import("rem_pio2.zig").rem_pio2;
const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f;
+const rem_pio2l = @import("rem_pio2l.zig").rem_pio2l;
+const utils = @import("math_utils.zig");
const compiler_rt = @import("../compiler_rt.zig");
const symbol = compiler_rt.symbol;
@@ -188,50 +192,12 @@ pub fn sincos(x: f64, r_sin: *f64, r_cos: *f64) callconv(.c) void {
}
}
-pub fn __sincosx(x: f80, r_sin: *f80, r_cos: *f80) callconv(.c) void {
- // TODO: more efficient implementation
- //return sincos_generic(f80, x, r_sin, r_cos);
- var big_sin: f128 = undefined;
- var big_cos: f128 = undefined;
- sincosq(x, &big_sin, &big_cos);
- r_sin.* = @as(f80, @floatCast(big_sin));
- r_cos.* = @as(f80, @floatCast(big_cos));
-}
-
-pub fn sincosq(x: f128, r_sin: *f128, r_cos: *f128) callconv(.c) void {
- // TODO: more correct implementation
- //return sincos_generic(f128, x, r_sin, r_cos);
- var small_sin: f64 = undefined;
- var small_cos: f64 = undefined;
- sincos(@as(f64, @floatCast(x)), &small_sin, &small_cos);
- r_sin.* = small_sin;
- r_cos.* = small_cos;
-}
-
-pub fn sincosl(x: c_longdouble, r_sin: *c_longdouble, r_cos: *c_longdouble) callconv(.c) void {
- switch (@typeInfo(c_longdouble).float.bits) {
- 16 => return __sincosh(x, r_sin, r_cos),
- 32 => return sincosf(x, r_sin, r_cos),
- 64 => return sincos(x, r_sin, r_cos),
- 80 => return __sincosx(x, r_sin, r_cos),
- 128 => return sincosq(x, r_sin, r_cos),
- else => @compileError("unreachable"),
+fn sincoslGeneric(comptime T: type, x: T, r_sin: *T, r_cos: *T) void {
+ if (T != f80 and T != f128) {
+ @compileError("`sincoslGeneric` implemented only for `f80` and `f128`, got: " ++ @typeName(T));
}
-}
-
-pub const rem_pio2_generic = @compileError("TODO");
-
-/// Ported from musl sincosl.c. Needs the following dependencies to be complete:
-/// * rem_pio2_generic ported from __rem_pio2l.c
-/// * trig.sin_generic ported from __sinl.c
-/// * trig.cos_generic ported from __cosl.c
-inline fn sincos_generic(comptime F: type, x: F, r_sin: *F, r_cos: *F) void {
- const sc1pio4: F = 1.0 * math.pi / 4.0;
- const bits = @typeInfo(F).float.bits;
- const I = std.meta.Int(.unsigned, bits);
- const ix = @as(I, @bitCast(x)) & (math.maxInt(I) >> 1);
- const se: u16 = @truncate(ix >> (bits - 16));
+ const se = utils.ldSignExponent(x) & 0x7fff;
if (se == 0x7fff) {
const result = x - x;
r_sin.* = result;
@@ -239,26 +205,26 @@ inline fn sincos_generic(comptime F: type, x: F, r_sin: *F, r_cos: *F) void {
return;
}
- if (@as(F, @bitCast(ix)) < sc1pio4) {
- if (se < 0x3fff - math.floatFractionalBits(F) - 1) {
+ if (@abs(x) < utils.pi_4) {
+ if (se < 0x3fff - math.floatMantissaBits(T)) {
// raise underflow if subnormal
- if (se == 0) {
- if (compiler_rt.want_float_exceptions) mem.doNotOptimizeAway(x * 0x1p-120);
+ if (compiler_rt.want_float_exceptions and se == 0) {
+ mem.doNotOptimizeAway(x * 0x1p-120);
}
r_sin.* = x;
// raise inexact if x!=0
r_cos.* = 1.0 + x;
return;
}
- r_sin.* = trig.sin_generic(F, x, 0, 0);
- r_cos.* = trig.cos_generic(F, x, 0);
+ r_sin.* = trig.__sinl(T, x, 0.0, 0);
+ r_cos.* = trig.__cosl(T, x, 0.0);
return;
}
- var y: [2]F = undefined;
- const n = rem_pio2_generic(F, x, &y);
- const s = trig.sin_generic(F, y[0], y[1], 1);
- const c = trig.cos_generic(F, y[0], y[1]);
+ var y: [2]T = undefined;
+ const n = rem_pio2l(T, x, &y);
+ const s = trig.__sinl(T, y[0], y[1], 1);
+ const c = trig.__cosl(T, y[0], y[1]);
switch (n & 3) {
0 => {
r_sin.* = s;
@@ -278,3 +244,207 @@ inline fn sincos_generic(comptime F: type, x: F, r_sin: *F, r_cos: *F) void {
},
}
}
+
+pub fn __sincosx(x: f80, r_sin: *f80, r_cos: *f80) callconv(.c) void {
+ return sincoslGeneric(f80, x, r_sin, r_cos);
+}
+
+pub fn sincosq(x: f128, r_sin: *f128, r_cos: *f128) callconv(.c) void {
+ return sincoslGeneric(f128, x, r_sin, r_cos);
+}
+
+pub fn sincosl(x: c_longdouble, r_sin: *c_longdouble, r_cos: *c_longdouble) callconv(.c) void {
+ switch (@typeInfo(c_longdouble).float.bits) {
+ 16 => return __sincosh(x, r_sin, r_cos),
+ 32 => return sincosf(x, r_sin, r_cos),
+ 64 => return sincos(x, r_sin, r_cos),
+ 80 => return __sincosx(x, r_sin, r_cos),
+ 128 => return sincosq(x, r_sin, r_cos),
+ else => @compileError("unreachable"),
+ }
+}
+
+fn testSincosSpecial(comptime T: type) !void {
+ const f = switch (T) {
+ f32 => sincosf,
+ f64 => sincos,
+ f80 => __sincosx,
+ f128 => sincosq,
+ else => @compileError("unimplemented"),
+ };
+
+ var s: T = undefined;
+ var c: T = undefined;
+
+ f(0.0, &s, &c);
+ try expect(math.isPositiveZero(s));
+ try expect(c == 1.0);
+
+ f(-0.0, &s, &c);
+ try expect(math.isNegativeZero(s));
+ try expect(c == 1.0);
+
+ f(math.inf(T), &s, &c);
+ try expect(math.isNan(s));
+ try expect(math.isNan(c));
+
+ f(-math.inf(T), &s, &c);
+ try expect(math.isNan(s));
+ try expect(math.isNan(c));
+
+ f(math.nan(T), &s, &c);
+ try expect(math.isNan(s));
+ try expect(math.isNan(c));
+}
+
+test "sincos32.normal" {
+ const epsilon = math.floatEps(f32);
+ var s: f32 = undefined;
+ var c: f32 = undefined;
+
+ sincosf(0.0, &s, &c);
+ try expectApproxEqAbs(@as(f32, 0.0), s, epsilon);
+ try expectApproxEqAbs(@as(f32, 1.0), c, epsilon);
+
+ sincosf(0.2, &s, &c);
+ try expectApproxEqAbs(@as(f32, 0.19866933), s, epsilon);
+ try expectApproxEqAbs(@as(f32, 0.9800666), c, epsilon);
+
+ sincosf(0.8923, &s, &c);
+ try expectApproxEqAbs(@as(f32, 0.77851737), s, epsilon);
+ try expectApproxEqAbs(@as(f32, 0.6276231), c, epsilon);
+
+ sincosf(1.5, &s, &c);
+ try expectApproxEqAbs(@as(f32, 0.997495), s, epsilon);
+ try expectApproxEqAbs(@as(f32, 0.0707372), c, epsilon);
+
+ sincosf(-1.5, &s, &c);
+ try expectApproxEqAbs(@as(f32, -0.997495), s, epsilon);
+ try expectApproxEqAbs(@as(f32, 0.0707372), c, epsilon);
+
+ sincosf(37.45, &s, &c);
+ try expectApproxEqAbs(@as(f32, -0.24654257), s, epsilon);
+ try expectApproxEqAbs(@as(f32, 0.96913195), c, epsilon);
+
+ sincosf(89.123, &s, &c);
+ try expectApproxEqAbs(@as(f32, 0.9161657), s, epsilon);
+ try expectApproxEqAbs(@as(f32, 0.40079966), c, epsilon);
+}
+
+test "sincos32.special" {
+ try testSincosSpecial(f32);
+}
+
+test "sincos64.normal" {
+ const epsilon = math.floatEps(f64);
+ var s: f64 = undefined;
+ var c: f64 = undefined;
+
+ sincos(0.0, &s, &c);
+ try expectApproxEqAbs(@as(f64, 0.0), s, epsilon);
+ try expectApproxEqAbs(@as(f64, 1.0), c, epsilon);
+
+ sincos(0.2, &s, &c);
+ try expectApproxEqAbs(@as(f64, 0.19866933079506122), s, epsilon);
+ try expectApproxEqAbs(@as(f64, 0.9800665778412416), c, epsilon);
+
+ sincos(0.8923, &s, &c);
+ try expectApproxEqAbs(@as(f64, 0.7785173385577349), s, epsilon);
+ try expectApproxEqAbs(@as(f64, 0.6276230983360804), c, epsilon);
+
+ sincos(1.5, &s, &c);
+ try expectApproxEqAbs(@as(f64, 0.9974949866040544), s, epsilon);
+ try expectApproxEqAbs(@as(f64, 0.0707372016677029), c, epsilon);
+
+ sincos(-1.5, &s, &c);
+ try expectApproxEqAbs(@as(f64, -0.9974949866040544), s, epsilon);
+ try expectApproxEqAbs(@as(f64, 0.0707372016677029), c, epsilon);
+
+ sincos(37.45, &s, &c);
+ try expectApproxEqAbs(@as(f64, -0.24654331551411082), s, epsilon);
+ try expectApproxEqAbs(@as(f64, 0.9691317730707778), c, epsilon);
+
+ sincos(89.123, &s, &c);
+ try expectApproxEqAbs(@as(f64, 0.9161652766622714), s, epsilon);
+ try expectApproxEqAbs(@as(f64, 0.4008006809354791), c, epsilon);
+}
+
+test "sincos64.special" {
+ try testSincosSpecial(f64);
+}
+
+test "sincos80.normal" {
+ const epsilon = math.floatEps(f80);
+ var s: f80 = undefined;
+ var c: f80 = undefined;
+
+ __sincosx(0.0, &s, &c);
+ try expectApproxEqAbs(@as(f80, 0.0), s, epsilon);
+ try expectApproxEqAbs(@as(f80, 1.0), c, epsilon);
+
+ __sincosx(0.2, &s, &c);
+ try expectApproxEqAbs(@as(f80, 0.19866933079506121545941262711838975), s, epsilon);
+ try expectApproxEqAbs(@as(f80, 0.98006657784124163112419651674816888), c, epsilon);
+
+ __sincosx(0.8923, &s, &c);
+ try expectApproxEqAbs(@as(f80, 0.77851733855773487830689285621486050), s, epsilon);
+ try expectApproxEqAbs(@as(f80, 0.62762309833608037003563995939286067), c, epsilon);
+
+ __sincosx(1.5, &s, &c);
+ try expectApproxEqAbs(@as(f80, 0.99749498660405443094172337114148732), s, epsilon);
+ try expectApproxEqAbs(@as(f80, 0.070737201667702910088189851434268747), c, epsilon);
+
+ __sincosx(-1.5, &s, &c);
+ try expectApproxEqAbs(@as(f80, -0.99749498660405443094172337114148732), s, epsilon);
+ try expectApproxEqAbs(@as(f80, 0.070737201667702910088189851434268747), c, epsilon);
+
+ __sincosx(37.45, &s, &c);
+ try expectApproxEqAbs(@as(f80, -0.24654331551411356504), s, epsilon);
+ try expectApproxEqAbs(@as(f80, 0.9691317730707771246), c, epsilon);
+
+ __sincosx(89.123, &s, &c);
+ try expectApproxEqAbs(@as(f80, 0.91616527666226951006), s, epsilon);
+ try expectApproxEqAbs(@as(f80, 0.4008006809354834001), c, epsilon);
+}
+
+test "sincos80.special" {
+ try testSincosSpecial(f80);
+}
+
+test "sincos128.normal" {
+ const epsilon = math.floatEps(f128);
+ var s: f128 = undefined;
+ var c: f128 = undefined;
+
+ sincosq(0.0, &s, &c);
+ try expectApproxEqAbs(@as(f128, 0.0), s, epsilon);
+ try expectApproxEqAbs(@as(f128, 1.0), c, epsilon);
+
+ sincosq(0.2, &s, &c);
+ try expectApproxEqAbs(@as(f128, 0.19866933079506121545941262711838975), s, epsilon);
+ try expectApproxEqAbs(@as(f128, 0.98006657784124163112419651674816888), c, epsilon);
+
+ sincosq(0.8923, &s, &c);
+ try expectApproxEqAbs(@as(f128, 0.77851733855773487830689285621486050), s, epsilon);
+ try expectApproxEqAbs(@as(f128, 0.62762309833608037003563995939286067), c, epsilon);
+
+ sincosq(1.5, &s, &c);
+ try expectApproxEqAbs(@as(f128, 0.99749498660405443094172337114148732), s, epsilon);
+ try expectApproxEqAbs(@as(f128, 0.070737201667702910088189851434268747), c, epsilon);
+
+ sincosq(-1.5, &s, &c);
+ try expectApproxEqAbs(@as(f128, -0.99749498660405443094172337114148732), s, epsilon);
+ try expectApproxEqAbs(@as(f128, 0.070737201667702910088189851434268747), c, epsilon);
+
+ sincosq(37.45, &s, &c);
+ try expectApproxEqAbs(@as(f128, -0.24654331551411356571238581321661085), s, epsilon);
+ try expectApproxEqAbs(@as(f128, 0.96913177307077712443149563847233230), c, epsilon);
+
+ sincosq(89.123, &s, &c);
+ try expectApproxEqAbs(@as(f128, 0.91616527666226951075019849560482170), s, epsilon);
+ try expectApproxEqAbs(@as(f128, 0.40080068093548339848199454493704702), c, epsilon);
+}
+
+test "sincos128.special" {
+ try testSincosSpecial(f128);
+}
diff --git a/lib/libc/mingw/math/arm/sincos.S b/lib/libc/mingw/math/arm/sincos.S
@@ -1,30 +0,0 @@
-/**
- * This file has no copyright assigned and is placed in the Public Domain.
- * This file is part of the mingw-w64 runtime package.
- * No warranty is given; refer to the file DISCLAIMER.PD within this package.
- */
-#include <_mingw_mac.h>
-
- .file "sincos.S"
- .text
- .align 2
- /* zig patch: remove sincos symbol because sincos in compiler_rt is used instead */
- .globl __MINGW_USYMBOL(sincosl)
- .def __MINGW_USYMBOL(sincosl); .scl 2; .type 32; .endef
-__MINGW_USYMBOL(sincosl):
- push {r4, r5, r11, lr}
- add r11, sp, #8
- vpush {d8}
-
- mov r4, r0
- mov r5, r1
- vmov.f64 d8, d0
- bl sin
- vstr d0, [r4]
-
- vmov.f64 d0, d8
- bl cos
- vstr d0, [r5]
-
- vpop {d8}
- pop {r4, r5, r11, pc}
diff --git a/lib/libc/mingw/math/arm64/sincos.S b/lib/libc/mingw/math/arm64/sincos.S
@@ -1,32 +0,0 @@
-/**
- * This file has no copyright assigned and is placed in the Public Domain.
- * This file is part of the mingw-w64 runtime package.
- * No warranty is given; refer to the file DISCLAIMER.PD within this package.
- */
-#include <_mingw_mac.h>
-
- .file "sincos.S"
- .text
- .align 2
- /* zig patch: remove sincos symbol because sincos in compiler_rt is used instead */
- .globl __MINGW_USYMBOL(sincosl)
- .def __MINGW_USYMBOL(sincosl); .scl 2; .type 32; .endef
-__MINGW_USYMBOL(sincosl):
- str d8, [sp, #-32]!
- str x30, [sp, #8]
- stp x19, x20, [sp, #16]
-
- mov x19, x0
- mov x20, x1
- fmov d8, d0
- bl sin
- str d0, [x19]
-
- fmov d0, d8
- bl cos
- str d0, [x20]
-
- ldp x19, x20, [sp, #16]
- ldr x30, [sp, #8]
- ldr d8, [sp], #32
- ret
diff --git a/src/libs/mingw.zig b/src/libs/mingw.zig
@@ -986,13 +986,9 @@ const mingw32_arm32_src = [_][]const u8{
// mingwex
"math" ++ path.sep_str ++ "arm" ++ path.sep_str ++ "s_rint.c",
"math" ++ path.sep_str ++ "arm" ++ path.sep_str ++ "s_rintf.c",
- "math" ++ path.sep_str ++ "arm" ++ path.sep_str ++ "sincos.S",
};
-const mingw32_arm64_src = [_][]const u8{
- // mingwex
- "math" ++ path.sep_str ++ "arm64" ++ path.sep_str ++ "sincos.S",
-};
+const mingw32_arm64_src = [_][]const u8{};
const mingw32_winpthreads_src = [_][]const u8{
// winpthreads
