From c76644caf3051dd2f837a6c40987a64e01b09c54 Mon Sep 17 00:00:00 2001 From: mihael Date: Mon, 23 Mar 2026 22:50:52 +0100 Subject: [PATCH] `libzigc/math`: Port `rem_pio2l` from musl to `compiler_rt` Additionally, add helper functions for fetching the sign+exponent and top 16 bits of a f80/f128's mantissa. We'll need these functions to implement `cosl`, `sinl`, `tanl`, and, transitively, `sincosl`. --- lib/compiler_rt/math_utils.zig | 35 +++++++ lib/compiler_rt/rem_pio2l.zig | 173 +++++++++++++++++++++++++++++++++ 2 files changed, 208 insertions(+) create mode 100644 lib/compiler_rt/math_utils.zig create mode 100644 lib/compiler_rt/rem_pio2l.zig diff --git a/lib/compiler_rt/math_utils.zig b/lib/compiler_rt/math_utils.zig new file mode 100644 index 0000000000000000000000000000000000000000..76cbcec468a74602c2a2281d1a07c9a7e1c83039 --- /dev/null +++ b/lib/compiler_rt/math_utils.zig @@ -0,0 +1,35 @@ +const std = @import("std"); + +pub const U80 = std.meta.Int(.unsigned, 80); + +/// Returns the sign + exponent bits of a `long double` +pub fn ldSignExponent(x: anytype) u16 { + const T = @TypeOf(x); + switch (T) { + f80 => { + const bits: U80 = @bitCast(x); + return @intCast(bits >> 64); + }, + f128 => { + const bits: u128 = @bitCast(x); + return @intCast(bits >> 112); + }, + else => @compileError("`ldSignExponent` supports only `f80` and `f128`, got: " ++ @typeName(T)), + } +} + +/// Takes the top 16 bits of a `long double`'s mantissa +pub fn ldMantissaTop(x: anytype) u16 { + const T = @TypeOf(x); + switch (T) { + f80 => { + const bits: U80 = @bitCast(x); + return @intCast((bits >> 48) & 0xFFFF); + }, + f128 => { + const bits: u128 = @bitCast(x); + return @intCast((bits >> 96) & 0xFFFF); + }, + else => @compileError("`ldMantissaTop` supports only `f80` and `f128`, got: " ++ @typeName(T)), + } +} diff --git a/lib/compiler_rt/rem_pio2l.zig b/lib/compiler_rt/rem_pio2l.zig new file mode 100644 index 0000000000000000000000000000000000000000..1ca9d170bcdcbbfbc065f23367581d9afc63e980 --- /dev/null +++ b/lib/compiler_rt/rem_pio2l.zig @@ -0,0 +1,173 @@ +// Ported from musl, which is licensed under the MIT license: +// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT +// +// https://git.musl-libc.org/cgit/musl/tree/src/math/__rem_pio2l.c + +const std = @import("std"); +const math = std.math; + +const utils = @import("math_utils.zig"); +const rem_pio2_large = @import("rem_pio2_large.zig").rem_pio2_large; + +pub fn rem_pio2l(comptime T: type, x: T, y: *[2]T) i32 { + const impl = switch (T) { + f80 => struct { + const round1: i8 = 22; + const round2: i8 = 61; + const nx: i8 = 3; + const ny: i8 = 2; + + const pio4: T = 0x1.921fb54442d1846ap-1; + // 64 bits of 2/pi + const invpio2: T = 6.36619772367581343076e-01; // 0xa2f9836e4e44152a.0p-64 + // first 39 bits of pi/2 + const pio2_1: f64 = 1.57079632679597125389e+00; // 0x3FF921FB, 0x54444000 + // pi/2 - pio2_1 + const pio2_1t: T = -1.07463465549719416346e-12; // -0x973dcb3b399d747f.0p-103 + // second 39 bits of pi/2 + const pio2_2: f64 = -1.07463465549783099519e-12; // -0x12e7b967674000.0p-92 + // pi/2 - (pio2_1+pio2_2) + const pio2_2t: T = 6.36831716351095013979e-25; // 0xc51701b839a25205.0p-144 + // pi/2 - (pio2_1+pio2_2+pio2_3) + const pio2_3t: T = -2.75299651904407171810e-37; // -0xbb5bf6c7ddd660ce.0p-185 + // third 39 bits of pi/2 + const pio2_3: f64 = 6.36831716351370313614e-25; // 0x18a2e037074000.0p-133 + + fn small(x_val: T) bool { + const se = utils.ldSignExponent(x_val); + const top = utils.ldMantissaTop(x_val); + const lhs = (@as(u32, se & 0x7fff) << 16) | top; + const rhs: u32 = ((0x3fff + 25) << 16) | 0x921f >> 1 | 0x8000; + return lhs < rhs; + } + + fn quobits(v: T) i32 { + const q: i32 = @intFromFloat(v); + return @intCast(@as(u32, @bitCast(q)) & 0x7fffffff); + } + }, + f128 => struct { + const round1: i8 = 51; + const round2: i8 = 119; + const nx: i8 = 5; + const ny: i8 = 3; + + const pio4: T = 0x1.921fb54442d18469898cc51701b8p-1; + const invpio2: T = 6.3661977236758134307553505349005747e-01; + const pio2_1: T = 1.5707963267948966192292994253909555e+00; + const pio2_1t: T = 2.0222662487959507323996846200947577e-21; + const pio2_2: T = 2.0222662487959507323994779168837751e-21; + const pio2_2t: T = 2.0670321098263988236496903051604844e-43; + const pio2_3: T = 2.0670321098263988236499468110329591e-43; + const pio2_3t: T = -2.5650587247459238361625433492959285e-65; + + fn small(x_val: T) bool { + const se = utils.ldSignExponent(x_val); + const top = utils.ldMantissaTop(x_val); + const lhs = (@as(u32, se & 0x7fff) << 16) | top; + const rhs: u32 = ((0x3fff + 45) << 16) | 0x921f; + return lhs < rhs; + } + + fn quobits(fn_val: T) i32 { + const q: i64 = @intFromFloat(fn_val); + return @intCast(@as(u64, @bitCast(q)) & 0x7fffffff); + } + }, + else => @compileError("rem_pio2l supports only f80 and f128, got: " ++ @typeName(T)), + }; + + const x_se = utils.ldSignExponent(x); + const ex: i32 = @intCast(x_se & 0x7fff); + + if (impl.small(x)) { + // rint(x/(pi/2)) + const toint: T = 1.5 / math.floatEps(T); + var fn_ = x * impl.invpio2 + toint - toint; + var n = impl.quobits(fn_); + var r = x - fn_ * @as(T, impl.pio2_1); + var w = fn_ * impl.pio2_1t; // 1st round good to 102/180 bits + + // Matters with directed rounding. + if (r - w < -impl.pio4) { + @branchHint(.unlikely); + n -= 1; + fn_ -= 1; + r = x - fn_ * @as(T, impl.pio2_1); + w = fn_ * impl.pio2_1t; + } else if (r - w > impl.pio4) { + @branchHint(.unlikely); + n += 1; + fn_ += 1; + r = x - fn_ * @as(T, impl.pio2_1); + w = fn_ * impl.pio2_1t; + } + + y[0] = r - w; + + const ey: i32 = @intCast(utils.ldSignExponent(y[0]) & 0x7fff); + if (ex - ey > impl.round1) { + var t = r; + w = fn_ * impl.pio2_2; + r = t - w; + w = fn_ * impl.pio2_2t - ((t - r) - w); + y[0] = r - w; + const ey2: i32 = @intCast(utils.ldSignExponent(y[0]) & 0x7fff); + if (ex - ey2 > impl.round2) { + t = r; + w = fn_ * impl.pio2_3; + r = t - w; + w = fn_ * impl.pio2_3t - ((t - r) - w); + y[0] = r - w; + } + } + y[1] = (r - y[0]) - w; + return n; + } + + // all other (large) arguments + if (ex == 0x7fff) { // x is inf or NaN + y[0] = x - x; + y[1] = y[0]; + return 0; + } + + var z: T = math.scalbn(@abs(x), -math.ilogb(x) + 23); + var tx: [impl.nx]f64 = undefined; + var ty: [impl.ny]f64 = undefined; + var i: usize = 0; + + while (i < impl.nx - 1) : (i += 1) { + tx[i] = @floatFromInt(@as(i32, @intFromFloat(z))); + z = (z - @as(T, tx[i])) * 0x1p24; + } + + tx[i] = @floatCast(z); + while (tx[i] == 0.0) { + i -= 1; + } + + const n = rem_pio2_large( + tx[0..(i + 1)], + ty[0..impl.ny], + ex - 0x3fff - 23, + @intCast(i + 1), + impl.ny, + ); + var w: f64 = ty[1]; + if (impl.ny == 3) { + w += ty[2]; + } + const r = ty[0] + w; + w -= r - ty[0]; + + if (x_se >> 15 != 0) { + y[0] = -@as(T, r); + y[1] = -@as(T, w); + return -n; + } + + y[0] = @as(T, r); + y[1] = @as(T, w); + return n; +} -- 2.54.0