From 9a53be5fb444873ee561eb54a0ea1aaeaaad844d Mon Sep 17 00:00:00 2001 From: mihael Date: Thu, 26 Mar 2026 21:43:55 +0100 Subject: [PATCH] `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= -Dtest-filter=sincosl -fqemu -fwasmtime --summary line Build Summary: 369/369 steps succeeded ``` Additionally, I've added unit tests for `sincos`. --- lib/compiler_rt/sincos.zig | 314 ++++++++++++++++++++++------- lib/libc/mingw/math/arm/sincos.S | 30 --- lib/libc/mingw/math/arm64/sincos.S | 32 --- src/libs/mingw.zig | 6 +- 4 files changed, 243 insertions(+), 139 deletions(-) delete mode 100644 lib/libc/mingw/math/arm/sincos.S delete mode 100644 lib/libc/mingw/math/arm64/sincos.S diff --git a/lib/compiler_rt/sincos.zig b/lib/compiler_rt/sincos.zig index 6aeb1105e59f7f8433eb631fab27b5c20d3fccd4..bc43525fbe73143978bbb172aa1b7408d3fa6470 100644 --- 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,24 +192,65 @@ pub fn sincos(x: f64, r_sin: *f64, r_cos: *f64) callconv(.c) void { } } +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)); + } + + const se = utils.ldSignExponent(x) & 0x7fff; + if (se == 0x7fff) { + const result = x - x; + r_sin.* = result; + r_cos.* = result; + return; + } + + if (@abs(x) < utils.pi_4) { + if (se < 0x3fff - math.floatMantissaBits(T)) { + // raise underflow if subnormal + 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.__sinl(T, x, 0.0, 0); + r_cos.* = trig.__cosl(T, x, 0.0); + return; + } + + 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; + r_cos.* = c; + }, + 1 => { + r_sin.* = c; + r_cos.* = -s; + }, + 2 => { + r_sin.* = -s; + r_cos.* = -c; + }, + else => { + r_sin.* = -c; + r_cos.* = s; + }, + } +} + 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)); + return sincoslGeneric(f80, x, r_sin, r_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; + 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 { @@ -219,62 +264,187 @@ pub fn sincosl(x: c_longdouble, r_sin: *c_longdouble, r_cos: *c_longdouble) call } } -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)); - - if (se == 0x7fff) { - const result = x - x; - r_sin.* = result; - r_cos.* = result; - return; - } - - if (@as(F, @bitCast(ix)) < sc1pio4) { - if (se < 0x3fff - math.floatFractionalBits(F) - 1) { - // raise underflow if subnormal - if (se == 0) { - if (compiler_rt.want_float_exceptions) 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); - 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]); - switch (n & 3) { - 0 => { - r_sin.* = s; - r_cos.* = c; - }, - 1 => { - r_sin.* = c; - r_cos.* = -s; - }, - 2 => { - r_sin.* = -s; - r_cos.* = -c; - }, - else => { - r_sin.* = -c; - r_cos.* = s; - }, - } +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 deleted file mode 100644 index ac22d45f25dbc19e9312f7e2a1ffb1cccfe45f6b..0000000000000000000000000000000000000000 --- a/lib/libc/mingw/math/arm/sincos.S +++ /dev/null @@ -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 deleted file mode 100644 index d6e92e3bba5bb385f47ee6f2d382dc101c8ee237..0000000000000000000000000000000000000000 --- a/lib/libc/mingw/math/arm64/sincos.S +++ /dev/null @@ -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 index f0a44ca4be1082811abd85af6026567ed8a80d7c..e229db09d652cf1459498da17eeb527d3fe1ec1a 100644 --- 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 -- 2.54.0