From a8b710697eeb7d83364cdca303167d14571e75e3 Mon Sep 17 00:00:00 2001 From: mihael Date: Wed, 18 Mar 2026 18:50:16 +0100 Subject: [PATCH 01/20] `libzigc/math`: Switch to `compiler_rt` impl of `sincosl` --- lib/libc/mingw/math/arm-common/sincosl.c | 13 ----- lib/libc/musl/src/math/sincosl.c | 60 ------------------------ src/libs/mingw.zig | 1 - src/libs/musl.zig | 1 - src/libs/wasi_libc.zig | 1 - 5 files changed, 76 deletions(-) delete mode 100644 lib/libc/mingw/math/arm-common/sincosl.c delete mode 100644 lib/libc/musl/src/math/sincosl.c diff --git a/lib/libc/mingw/math/arm-common/sincosl.c b/lib/libc/mingw/math/arm-common/sincosl.c deleted file mode 100644 index 8caeaec5a9804f4fe32c1b24eb70d37570ac21a9..0000000000000000000000000000000000000000 --- a/lib/libc/mingw/math/arm-common/sincosl.c +++ /dev/null @@ -1,13 +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 - -void sincosl(long double x, long double *s, long double *c) -{ - *s = sinl(x); - *c = cosl(x); -} diff --git a/lib/libc/musl/src/math/sincosl.c b/lib/libc/musl/src/math/sincosl.c deleted file mode 100644 index d3ac1c4c8c29e11067d55a3a99fa82eb888c5dd5..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/sincosl.c +++ /dev/null @@ -1,60 +0,0 @@ -#define _GNU_SOURCE -#include "libm.h" - -#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 -void sincosl(long double x, long double *sin, long double *cos) -{ - double sind, cosd; - sincos(x, &sind, &cosd); - *sin = sind; - *cos = cosd; -} -#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 -void sincosl(long double x, long double *sin, long double *cos) -{ - union ldshape u = {x}; - unsigned n; - long double y[2], s, c; - - u.i.se &= 0x7fff; - if (u.i.se == 0x7fff) { - *sin = *cos = x - x; - return; - } - if (u.f < M_PI_4) { - if (u.i.se < 0x3fff - LDBL_MANT_DIG) { - /* raise underflow if subnormal */ - if (u.i.se == 0) FORCE_EVAL(x*0x1p-120f); - *sin = x; - /* raise inexact if x!=0 */ - *cos = 1.0 + x; - return; - } - *sin = __sinl(x, 0, 0); - *cos = __cosl(x, 0); - return; - } - n = __rem_pio2l(x, y); - s = __sinl(y[0], y[1], 1); - c = __cosl(y[0], y[1]); - switch (n & 3) { - case 0: - *sin = s; - *cos = c; - break; - case 1: - *sin = c; - *cos = -s; - break; - case 2: - *sin = -s; - *cos = -c; - break; - case 3: - default: - *sin = -c; - *cos = s; - break; - } -} -#endif diff --git a/src/libs/mingw.zig b/src/libs/mingw.zig index 07e67ea468f859f174c390158f9ee141a4f52062..99fad69a06b0fa04ac83d62f88d26fc2f66028d3 100644 --- a/src/libs/mingw.zig +++ b/src/libs/mingw.zig @@ -987,7 +987,6 @@ const mingw32_x86_32_src = [_][]const u8{ const mingw32_arm_src = [_][]const u8{ // mingwex "math" ++ path.sep_str ++ "arm-common" ++ path.sep_str ++ "ldexpl.c", - "math" ++ path.sep_str ++ "arm-common" ++ path.sep_str ++ "sincosl.c", }; const mingw32_arm32_src = [_][]const u8{ diff --git a/src/libs/musl.zig b/src/libs/musl.zig index e206d8636988ea0553652987a841871cf6793065..649bef26a3b2ffe4c54ca2349963f51bd416f7c5 100644 --- a/src/libs/musl.zig +++ b/src/libs/musl.zig @@ -992,7 +992,6 @@ const src_files = [_][]const u8{ "musl/src/math/significand.c", "musl/src/math/significandf.c", "musl/src/math/__sin.c", - "musl/src/math/sincosl.c", "musl/src/math/__sindf.c", "musl/src/math/sinh.c", "musl/src/math/sinhf.c", diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig index 12edf9b783a9a1954c008f229178756a937bfb21..953c3cbff5f17ac4929e9f7f51009beccd1d74b7 100644 --- a/src/libs/wasi_libc.zig +++ b/src/libs/wasi_libc.zig @@ -785,7 +785,6 @@ const libc_top_half_src_files = [_][]const u8{ "musl/src/math/significand.c", "musl/src/math/significandf.c", "musl/src/math/__sin.c", - "musl/src/math/sincosl.c", "musl/src/math/__sindf.c", "musl/src/math/sinhl.c", "musl/src/math/__sinl.c", -- 2.54.0 From ba4af7aa8a33d19bbf7b0225a2d1407d5a756c06 Mon Sep 17 00:00:00 2001 From: mihael Date: Wed, 18 Mar 2026 19:41:29 +0100 Subject: [PATCH 02/20] `libzigc/math`: Switch to `compiler_rt` impl of `tanl` --- lib/libc/mingw/math/x86/tanl.S | 62 -------------- lib/libc/musl/src/math/__tanl.c | 143 -------------------------------- lib/libc/musl/src/math/tanl.c | 29 ------- src/libs/mingw.zig | 1 - src/libs/musl.zig | 2 - src/libs/wasi_libc.zig | 2 - 6 files changed, 239 deletions(-) delete mode 100644 lib/libc/mingw/math/x86/tanl.S delete mode 100644 lib/libc/musl/src/math/__tanl.c delete mode 100644 lib/libc/musl/src/math/tanl.c diff --git a/lib/libc/mingw/math/x86/tanl.S b/lib/libc/mingw/math/x86/tanl.S deleted file mode 100644 index f11b53920beddd16ab407a9e4cc75dd7d595414b..0000000000000000000000000000000000000000 --- a/lib/libc/mingw/math/x86/tanl.S +++ /dev/null @@ -1,62 +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 "tanl.S" - .text -#ifdef __x86_64__ - .align 8 -#else - .align 4 -#endif -.globl __MINGW_USYMBOL(tanl) - .def __MINGW_USYMBOL(tanl); .scl 2; .type 32; .endef -__MINGW_USYMBOL(tanl): -#ifdef __x86_64__ - fldt (%rdx) - fptan - fnstsw %ax - testl $0x400,%eax - jnz 1f - fstp %st(0) - movq %rcx,%rax - movq $0,8(%rcx) - fstpt (%rcx) - ret -1: fldpi - fadd %st(0) - fxch %st(1) -2: fprem1 - fstsw %ax - testl $0x400,%eax - jnz 2b - fstp %st(1) - fptan - fstp %st(0) - movq %rcx,%rax - movq $0,8(%rcx) - fstpt (%rcx) - ret -#else - fldt 4(%esp) - fptan - fnstsw %ax - testl $0x400,%eax - jnz 1f - fstp %st(0) - ret -1: fldpi - fadd %st(0) - fxch %st(1) -2: fprem1 - fstsw %ax - testl $0x400,%eax - jnz 2b - fstp %st(1) - fptan - fstp %st(0) - ret -#endif diff --git a/lib/libc/musl/src/math/__tanl.c b/lib/libc/musl/src/math/__tanl.c deleted file mode 100644 index 54abc3daf1ad126306b363151aaa8dcca5928423..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/__tanl.c +++ /dev/null @@ -1,143 +0,0 @@ -/* origin: FreeBSD /usr/src/lib/msun/ld80/k_tanl.c */ -/* origin: FreeBSD /usr/src/lib/msun/ld128/k_tanl.c */ -/* - * ==================================================== - * Copyright 2004 Sun Microsystems, Inc. All Rights Reserved. - * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans. - * - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#include "libm.h" - -#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 -#if LDBL_MANT_DIG == 64 -/* - * ld80 version of __tan.c. See __tan.c for most comments. - */ -/* - * Domain [-0.67434, 0.67434], range ~[-2.25e-22, 1.921e-22] - * |tan(x)/x - t(x)| < 2**-71.9 - * - * See __cosl.c for more details about the polynomial. - */ -static const long double -T3 = 0.333333333333333333180L, /* 0xaaaaaaaaaaaaaaa5.0p-65 */ -T5 = 0.133333333333333372290L, /* 0x88888888888893c3.0p-66 */ -T7 = 0.0539682539682504975744L, /* 0xdd0dd0dd0dc13ba2.0p-68 */ -pio4 = 0.785398163397448309628L, /* 0xc90fdaa22168c235.0p-64 */ -pio4lo = -1.25413940316708300586e-20L; /* -0xece675d1fc8f8cbb.0p-130 */ -static const double -T9 = 0.021869488536312216, /* 0x1664f4882cc1c2.0p-58 */ -T11 = 0.0088632355256619590, /* 0x1226e355c17612.0p-59 */ -T13 = 0.0035921281113786528, /* 0x1d6d3d185d7ff8.0p-61 */ -T15 = 0.0014558334756312418, /* 0x17da354aa3f96b.0p-62 */ -T17 = 0.00059003538700862256, /* 0x13559358685b83.0p-63 */ -T19 = 0.00023907843576635544, /* 0x1f56242026b5be.0p-65 */ -T21 = 0.000097154625656538905, /* 0x1977efc26806f4.0p-66 */ -T23 = 0.000038440165747303162, /* 0x14275a09b3ceac.0p-67 */ -T25 = 0.000018082171885432524, /* 0x12f5e563e5487e.0p-68 */ -T27 = 0.0000024196006108814377, /* 0x144c0d80cc6896.0p-71 */ -T29 = 0.0000078293456938132840, /* 0x106b59141a6cb3.0p-69 */ -T31 = -0.0000032609076735050182, /* -0x1b5abef3ba4b59.0p-71 */ -T33 = 0.0000023261313142559411; /* 0x13835436c0c87f.0p-71 */ -#define RPOLY(w) (T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + \ - w * (T25 + w * (T29 + w * T33))))))) -#define VPOLY(w) (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + \ - w * (T27 + w * T31)))))) -#elif LDBL_MANT_DIG == 113 -/* - * ld128 version of __tan.c. See __tan.c for most comments. - */ -/* - * Domain [-0.67434, 0.67434], range ~[-3.37e-36, 1.982e-37] - * |tan(x)/x - t(x)| < 2**-117.8 (XXX should be ~1e-37) - * - * See __cosl.c for more details about the polynomial. - */ -static const long double -T3 = 0x1.5555555555555555555555555553p-2L, -T5 = 0x1.1111111111111111111111111eb5p-3L, -T7 = 0x1.ba1ba1ba1ba1ba1ba1ba1b694cd6p-5L, -T9 = 0x1.664f4882c10f9f32d6bbe09d8bcdp-6L, -T11 = 0x1.226e355e6c23c8f5b4f5762322eep-7L, -T13 = 0x1.d6d3d0e157ddfb5fed8e84e27b37p-9L, -T15 = 0x1.7da36452b75e2b5fce9ee7c2c92ep-10L, -T17 = 0x1.355824803674477dfcf726649efep-11L, -T19 = 0x1.f57d7734d1656e0aceb716f614c2p-13L, -T21 = 0x1.967e18afcb180ed942dfdc518d6cp-14L, -T23 = 0x1.497d8eea21e95bc7e2aa79b9f2cdp-15L, -T25 = 0x1.0b132d39f055c81be49eff7afd50p-16L, -T27 = 0x1.b0f72d33eff7bfa2fbc1059d90b6p-18L, -T29 = 0x1.5ef2daf21d1113df38d0fbc00267p-19L, -T31 = 0x1.1c77d6eac0234988cdaa04c96626p-20L, -T33 = 0x1.cd2a5a292b180e0bdd701057dfe3p-22L, -T35 = 0x1.75c7357d0298c01a31d0a6f7d518p-23L, -T37 = 0x1.2f3190f4718a9a520f98f50081fcp-24L, -pio4 = 0x1.921fb54442d18469898cc51701b8p-1L, -pio4lo = 0x1.cd129024e088a67cc74020bbea60p-116L; -static const double -T39 = 0.000000028443389121318352, /* 0x1e8a7592977938.0p-78 */ -T41 = 0.000000011981013102001973, /* 0x19baa1b1223219.0p-79 */ -T43 = 0.0000000038303578044958070, /* 0x107385dfb24529.0p-80 */ -T45 = 0.0000000034664378216909893, /* 0x1dc6c702a05262.0p-81 */ -T47 = -0.0000000015090641701997785, /* -0x19ecef3569ebb6.0p-82 */ -T49 = 0.0000000029449552300483952, /* 0x194c0668da786a.0p-81 */ -T51 = -0.0000000022006995706097711, /* -0x12e763b8845268.0p-81 */ -T53 = 0.0000000015468200913196612, /* 0x1a92fc98c29554.0p-82 */ -T55 = -0.00000000061311613386849674, /* -0x151106cbc779a9.0p-83 */ -T57 = 1.4912469681508012e-10; /* 0x147edbdba6f43a.0p-85 */ -#define RPOLY(w) (T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + \ - w * (T25 + w * (T29 + w * (T33 + w * (T37 + w * (T41 + \ - w * (T45 + w * (T49 + w * (T53 + w * T57))))))))))))) -#define VPOLY(w) (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + \ - w * (T27 + w * (T31 + w * (T35 + w * (T39 + w * (T43 + \ - w * (T47 + w * (T51 + w * T55)))))))))))) -#endif - -long double __tanl(long double x, long double y, int odd) { - long double z, r, v, w, s, a, t; - int big, sign; - - big = fabsl(x) >= 0.67434; - if (big) { - sign = 0; - if (x < 0) { - sign = 1; - x = -x; - y = -y; - } - x = (pio4 - x) + (pio4lo - y); - y = 0.0; - } - z = x * x; - w = z * z; - r = RPOLY(w); - v = z * VPOLY(w); - s = z * x; - r = y + z * (s * (r + v) + y) + T3 * s; - w = x + r; - if (big) { - s = 1 - 2*odd; - v = s - 2.0 * (x + (r - w * w / (w + s))); - return sign ? -v : v; - } - if (!odd) - return w; - /* - * if allow error up to 2 ulp, simply return - * -1.0 / (x+r) here - */ - /* compute -1.0 / (x+r) accurately */ - z = w; - z = z + 0x1p32 - 0x1p32; - v = r - (z - x); /* z+v = r+x */ - t = a = -1.0 / w; /* a = -1.0/w */ - t = t + 0x1p32 - 0x1p32; - s = 1.0 + t * z; - return t + a * (s + t * v); -} -#endif diff --git a/lib/libc/musl/src/math/tanl.c b/lib/libc/musl/src/math/tanl.c deleted file mode 100644 index 6af067127ade990b98d26696594e5737fbf50ade..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/tanl.c +++ /dev/null @@ -1,29 +0,0 @@ -#include "libm.h" - -#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 -long double tanl(long double x) -{ - return tan(x); -} -#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 -long double tanl(long double x) -{ - union ldshape u = {x}; - long double y[2]; - unsigned n; - - u.i.se &= 0x7fff; - if (u.i.se == 0x7fff) - return x - x; - if (u.f < M_PI_4) { - if (u.i.se < 0x3fff - LDBL_MANT_DIG/2) { - /* raise inexact if x!=0 and underflow if subnormal */ - FORCE_EVAL(u.i.se == 0 ? x*0x1p-120f : x+0x1p120f); - return x; - } - return __tanl(x, 0, 0); - } - n = __rem_pio2l(x, y); - return __tanl(y[0], y[1], n&1); -} -#endif diff --git a/src/libs/mingw.zig b/src/libs/mingw.zig index 99fad69a06b0fa04ac83d62f88d26fc2f66028d3..f7b73cc96e8db83c8012be02c561b98c8b082ede 100644 --- a/src/libs/mingw.zig +++ b/src/libs/mingw.zig @@ -967,7 +967,6 @@ const mingw32_x86_src = [_][]const u8{ "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "scalbnl.S", "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "sinl.c", "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "sinl_internal.S", - "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "tanl.S", // ucrtbase "math" ++ path.sep_str ++ "nextafterl.c", "math" ++ path.sep_str ++ "nexttoward.c", diff --git a/src/libs/musl.zig b/src/libs/musl.zig index 649bef26a3b2ffe4c54ca2349963f51bd416f7c5..a3806437d73d7f4734ddd3fde3d69740e8f72e14 100644 --- a/src/libs/musl.zig +++ b/src/libs/musl.zig @@ -1001,8 +1001,6 @@ const src_files = [_][]const u8{ "musl/src/math/__tan.c", "musl/src/math/__tandf.c", "musl/src/math/tanhl.c", - "musl/src/math/__tanl.c", - "musl/src/math/tanl.c", "musl/src/math/tgamma.c", "musl/src/math/tgammaf.c", "musl/src/math/tgammal.c", diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig index 953c3cbff5f17ac4929e9f7f51009beccd1d74b7..aaf896806014dc96e21c99701af9427d8a2400d0 100644 --- a/src/libs/wasi_libc.zig +++ b/src/libs/wasi_libc.zig @@ -792,8 +792,6 @@ const libc_top_half_src_files = [_][]const u8{ "musl/src/math/__tan.c", "musl/src/math/__tandf.c", "musl/src/math/tanhl.c", - "musl/src/math/__tanl.c", - "musl/src/math/tanl.c", "musl/src/math/tgamma.c", "musl/src/math/tgammaf.c", "musl/src/math/tgammal.c", -- 2.54.0 From dcffee067290c03f89c5afcacab1ba543ee28e7c Mon Sep 17 00:00:00 2001 From: mihael Date: Wed, 18 Mar 2026 20:32:18 +0100 Subject: [PATCH 03/20] `libzigc/math`: Implement `frexp` The changes were tested by running: ``` $ ./build/stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib $ stage4/bin/zig build test-libc -Dlibc-test-path= -Dtest-filter=frexp -fqemu -fwasmtime --summary line Build Summary: 737/737 steps succeeded ``` The tests were passing even when it was a straightforward calling of Zig std library, but I wanted the `x is NaN` special case to match the behaviour described in `libc` manpages, and for it to be consistent with how infinities as arguments are handled in Zig. --- lib/c/math.zig | 20 ++++++++++++++++++++ lib/libc/musl/src/math/frexp.c | 23 ----------------------- src/libs/musl.zig | 1 - src/libs/wasi_libc.zig | 1 - 4 files changed, 20 insertions(+), 25 deletions(-) delete mode 100644 lib/libc/musl/src/math/frexp.c diff --git a/lib/c/math.zig b/lib/c/math.zig index 5b191c3e93c80d24bfdc4ae823f80b77f1fcb64f..8fe93716dbebc34bf2bc5276d764c4e93ffd8a16 100644 --- a/lib/c/math.zig +++ b/lib/c/math.zig @@ -60,6 +60,7 @@ comptime { symbol(&exp10, "exp10"); symbol(&exp10f, "exp10f"); symbol(&fdim, "fdim"); + symbol(&frexp, "frexp"); symbol(&hypot, "hypot"); symbol(&modf, "modf"); symbol(&pow, "pow"); @@ -161,6 +162,25 @@ fn fdim(x: f64, y: f64) callconv(.c) f64 { return 0; } +fn frexp(x: f64, e: *c_int) callconv(.c) f64 { + // libc expects `*e` to be unspecified in this case; an unspecified C value + // should be a valid value of the relevant type, yet Zig's std + // implementation sets it to `undefined` -- which can even be nonsense + // according to the type (int). Therefore, we're setting it to a valid + // int value in Zig -- a zero. + // + // This mirrors the handling of infinities, where libc also expects + // unspecified for the value `*e` and Zig std sets it to a zero. + if (math.isNan(x)) { + e.* = 0; + return x; + } + + const r = math.frexp(x); + e.* = r.exponent; + return r.significand; +} + fn hypot(x: f64, y: f64) callconv(.c) f64 { return math.hypot(x, y); } diff --git a/lib/libc/musl/src/math/frexp.c b/lib/libc/musl/src/math/frexp.c deleted file mode 100644 index 27b6266ed0c1d7c5dadd06ecc186a994fdcd1c52..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/frexp.c +++ /dev/null @@ -1,23 +0,0 @@ -#include -#include - -double frexp(double x, int *e) -{ - union { double d; uint64_t i; } y = { x }; - int ee = y.i>>52 & 0x7ff; - - if (!ee) { - if (x) { - x = frexp(x*0x1p64, e); - *e -= 64; - } else *e = 0; - return x; - } else if (ee == 0x7ff) { - return x; - } - - *e = ee - 0x3fe; - y.i &= 0x800fffffffffffffull; - y.i |= 0x3fe0000000000000ull; - return y.d; -} diff --git a/src/libs/musl.zig b/src/libs/musl.zig index a3806437d73d7f4734ddd3fde3d69740e8f72e14..e90a08f92c5f62a570728a76f27778160b88f8e4 100644 --- a/src/libs/musl.zig +++ b/src/libs/musl.zig @@ -839,7 +839,6 @@ const src_files = [_][]const u8{ "musl/src/math/__fpclassify.c", "musl/src/math/__fpclassifyf.c", "musl/src/math/__fpclassifyl.c", - "musl/src/math/frexp.c", "musl/src/math/frexpf.c", "musl/src/math/frexpl.c", "musl/src/math/i386/acosl.s", diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig index aaf896806014dc96e21c99701af9427d8a2400d0..d694a0381e385a5ed35b2cb27a0e25e9e144acf3 100644 --- a/src/libs/wasi_libc.zig +++ b/src/libs/wasi_libc.zig @@ -701,7 +701,6 @@ const libc_top_half_src_files = [_][]const u8{ "musl/src/math/finitef.c", "musl/src/math/fma.c", "musl/src/math/fmaf.c", - "musl/src/math/frexp.c", "musl/src/math/frexpf.c", "musl/src/math/frexpl.c", "musl/src/math/ilogb.c", -- 2.54.0 From 521a093334eb92b4fd2430b9018ed94b912857a8 Mon Sep 17 00:00:00 2001 From: mihael Date: Wed, 18 Mar 2026 23:02:53 +0100 Subject: [PATCH 04/20] `libzigc/math`: Implement `frexpf` The `frexp` implementation was generalized so it can be used for `f32` and `c_longdouble` types as well. The changes were tested by running: ``` $ ./build/stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib $ stage4/bin/zig build test-libc -Dlibc-test-path= -Dtest-filter=frexp -fqemu -fwasmtime --summary line Build Summary: 737/737 steps succeeded $ stage4/bin/zig build test-libc -Dlibc-test-path= -Dtest-filter=frexpf -fqemu -fwasmtime --summary line Build Summary: 369/369 steps succeeded ``` --- lib/c/math.zig | 13 +++++++++++-- lib/libc/mingw/math/frexpf.c | 13 ------------- lib/libc/musl/src/math/frexpf.c | 23 ----------------------- src/libs/mingw.zig | 1 - src/libs/musl.zig | 1 - src/libs/wasi_libc.zig | 1 - 6 files changed, 11 insertions(+), 41 deletions(-) delete mode 100644 lib/libc/mingw/math/frexpf.c delete mode 100644 lib/libc/musl/src/math/frexpf.c diff --git a/lib/c/math.zig b/lib/c/math.zig index 8fe93716dbebc34bf2bc5276d764c4e93ffd8a16..1f379792a7c61145815193c33496872e118ea38b 100644 --- a/lib/c/math.zig +++ b/lib/c/math.zig @@ -36,6 +36,7 @@ comptime { if (builtin.target.isMinGW() or builtin.target.isMuslLibC() or builtin.target.isWasiLibC()) { symbol(&coshf, "coshf"); + symbol(&frexpf, "frexpf"); symbol(&hypotf, "hypotf"); symbol(&hypotl, "hypotl"); symbol(&modff, "modff"); @@ -162,7 +163,7 @@ fn fdim(x: f64, y: f64) callconv(.c) f64 { return 0; } -fn frexp(x: f64, e: *c_int) callconv(.c) f64 { +fn frexpGeneric(comptime T: type, x: T, e: *c_int) T { // libc expects `*e` to be unspecified in this case; an unspecified C value // should be a valid value of the relevant type, yet Zig's std // implementation sets it to `undefined` -- which can even be nonsense @@ -170,7 +171,7 @@ fn frexp(x: f64, e: *c_int) callconv(.c) f64 { // int value in Zig -- a zero. // // This mirrors the handling of infinities, where libc also expects - // unspecified for the value `*e` and Zig std sets it to a zero. + // unspecified for the value of `*e` and Zig std sets it to a zero. if (math.isNan(x)) { e.* = 0; return x; @@ -181,6 +182,14 @@ fn frexp(x: f64, e: *c_int) callconv(.c) f64 { return r.significand; } +fn frexp(x: f64, e: *c_int) callconv(.c) f64 { + return frexpGeneric(f64, x, e); +} + +fn frexpf(x: f32, e: *c_int) callconv(.c) f32 { + return frexpGeneric(f32, x, e); +} + fn hypot(x: f64, y: f64) callconv(.c) f64 { return math.hypot(x, y); } diff --git a/lib/libc/mingw/math/frexpf.c b/lib/libc/mingw/math/frexpf.c deleted file mode 100644 index d4b7ab36c20475e51168a19909bfd643e4bac1fa..0000000000000000000000000000000000000000 --- a/lib/libc/mingw/math/frexpf.c +++ /dev/null @@ -1,13 +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. - */ -extern double __cdecl frexp(double _X,int *_Y); - -float frexpf (float, int *); -float frexpf (float x, int *expn) -{ - return (float)frexp(x, expn); -} - diff --git a/lib/libc/musl/src/math/frexpf.c b/lib/libc/musl/src/math/frexpf.c deleted file mode 100644 index 0787097527c22a4be256abb600b9a05061a2ec87..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/frexpf.c +++ /dev/null @@ -1,23 +0,0 @@ -#include -#include - -float frexpf(float x, int *e) -{ - union { float f; uint32_t i; } y = { x }; - int ee = y.i>>23 & 0xff; - - if (!ee) { - if (x) { - x = frexpf(x*0x1p64, e); - *e -= 64; - } else *e = 0; - return x; - } else if (ee == 0xff) { - return x; - } - - *e = ee - 0x7e; - y.i &= 0x807ffffful; - y.i |= 0x3f000000ul; - return y.f; -} diff --git a/src/libs/mingw.zig b/src/libs/mingw.zig index f7b73cc96e8db83c8012be02c561b98c8b082ede..e9417eb9454b22a88ac4a4c9576b7ad6871576df 100644 --- a/src/libs/mingw.zig +++ b/src/libs/mingw.zig @@ -613,7 +613,6 @@ const mingw32_generic_src = [_][]const u8{ "math" ++ path.sep_str ++ "fpclassify.c", "math" ++ path.sep_str ++ "fpclassifyf.c", "math" ++ path.sep_str ++ "fpclassifyl.c", - "math" ++ path.sep_str ++ "frexpf.c", "math" ++ path.sep_str ++ "frexpl.c", "math" ++ path.sep_str ++ "ldexpf.c", "math" ++ path.sep_str ++ "lgamma.c", diff --git a/src/libs/musl.zig b/src/libs/musl.zig index e90a08f92c5f62a570728a76f27778160b88f8e4..86e0995220842f3d6b35442c2bc497df26384522 100644 --- a/src/libs/musl.zig +++ b/src/libs/musl.zig @@ -839,7 +839,6 @@ const src_files = [_][]const u8{ "musl/src/math/__fpclassify.c", "musl/src/math/__fpclassifyf.c", "musl/src/math/__fpclassifyl.c", - "musl/src/math/frexpf.c", "musl/src/math/frexpl.c", "musl/src/math/i386/acosl.s", "musl/src/math/i386/asinf.s", diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig index d694a0381e385a5ed35b2cb27a0e25e9e144acf3..d0743dbe3d10334e7aeafbef8c71bc3ae7043195 100644 --- a/src/libs/wasi_libc.zig +++ b/src/libs/wasi_libc.zig @@ -701,7 +701,6 @@ const libc_top_half_src_files = [_][]const u8{ "musl/src/math/finitef.c", "musl/src/math/fma.c", "musl/src/math/fmaf.c", - "musl/src/math/frexpf.c", "musl/src/math/frexpl.c", "musl/src/math/ilogb.c", "musl/src/math/ilogbf.c", -- 2.54.0 From db1e649575a741527ca0fd7815515c2f97924410 Mon Sep 17 00:00:00 2001 From: mihael Date: Wed, 18 Mar 2026 23:20:46 +0100 Subject: [PATCH 05/20] `libzigc/math`: Implement `frexpl` The changes were tested by running: ``` $ ./build/stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib $ stage4/bin/zig build test-libc -Dlibc-test-path= -Dtest-filter=frexpl -fqemu -fwasmtime --summary line Build Summary: 369/369 steps succeeded ``` --- lib/c/math.zig | 5 +++ lib/libc/mingw/math/frexpl.c | 71 --------------------------------- lib/libc/musl/src/math/frexpl.c | 29 -------------- src/libs/mingw.zig | 1 - src/libs/musl.zig | 1 - src/libs/wasi_libc.zig | 1 - 6 files changed, 5 insertions(+), 103 deletions(-) delete mode 100644 lib/libc/mingw/math/frexpl.c delete mode 100644 lib/libc/musl/src/math/frexpl.c diff --git a/lib/c/math.zig b/lib/c/math.zig index 1f379792a7c61145815193c33496872e118ea38b..1f80e1cc7e5362891884988ecff8423269bf9ce8 100644 --- a/lib/c/math.zig +++ b/lib/c/math.zig @@ -37,6 +37,7 @@ comptime { if (builtin.target.isMinGW() or builtin.target.isMuslLibC() or builtin.target.isWasiLibC()) { symbol(&coshf, "coshf"); symbol(&frexpf, "frexpf"); + symbol(&frexpl, "frexpl"); symbol(&hypotf, "hypotf"); symbol(&hypotl, "hypotl"); symbol(&modff, "modff"); @@ -190,6 +191,10 @@ fn frexpf(x: f32, e: *c_int) callconv(.c) f32 { return frexpGeneric(f32, x, e); } +fn frexpl(x: c_longdouble, e: *c_int) callconv(.c) c_longdouble { + return frexpGeneric(c_longdouble, x, e); +} + fn hypot(x: f64, y: f64) callconv(.c) f64 { return math.hypot(x, y); } diff --git a/lib/libc/mingw/math/frexpl.c b/lib/libc/mingw/math/frexpl.c deleted file mode 100644 index f60b68bcde18c8fafd7451d0774f6b1aaa2aaaa2..0000000000000000000000000000000000000000 --- a/lib/libc/mingw/math/frexpl.c +++ /dev/null @@ -1,71 +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. - */ -long double frexpl(long double value, int* exp); - -#if __SIZEOF_LONG_DOUBLE__ == __SIZEOF_DOUBLE__ - -double frexp(double value, int* exp); - -/* On ARM `long double` is 64 bits. */ -long double frexpl(long double value, int* exp) -{ - return frexp(value, exp); -} - -#elif defined(_AMD64_) || defined(__x86_64__) || defined(_X86_) || defined(__i386__) - -#include - -/* https://en.wikipedia.org/wiki/Extended_precision#x86_extended_precision_format */ -typedef union x87reg_ { - struct __attribute__((__packed__)) { - uint64_t f64; - uint16_t exp : 15; - uint16_t sgn : 1; - }; - long double f; -} x87reg; - -long double frexpl(long double value, int* exp) -{ - int n; - x87reg reg; - reg.f = value; - if(reg.exp == 0x7FFF) { - /* The value is an infinity or NaN. - * Store zero in `*exp`. Return the value as is. */ - *exp = 0; - return reg.f; - } - if(reg.exp != 0) { - /* The value is normalized. - * Extract and zero out the exponent. */ - *exp = reg.exp - 0x3FFE; - reg.exp = 0x3FFE; - return reg.f; - } - if(reg.f64 == 0) { - /* The value is zero. - * Store zero in `*exp`. Return the value as is. - * Note the signness. */ - *exp = 0; - return reg.f; - } - /* The value is denormalized. - * Extract the exponent, normalize the value, then zero out - * the exponent. Note that x87 uses an explicit leading bit. */ - n = __builtin_clzll(reg.f64); - reg.f64 <<= n; - *exp = 1 - 0x3FFE - n; - reg.exp = 0x3FFE; - return reg.f; -} - -#else - -#error Please add `frexpl()` implementation for this platform. - -#endif diff --git a/lib/libc/musl/src/math/frexpl.c b/lib/libc/musl/src/math/frexpl.c deleted file mode 100644 index 3c1b553748eec19a18804362acd508d3877d7fb6..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/frexpl.c +++ /dev/null @@ -1,29 +0,0 @@ -#include "libm.h" - -#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 -long double frexpl(long double x, int *e) -{ - return frexp(x, e); -} -#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 -long double frexpl(long double x, int *e) -{ - union ldshape u = {x}; - int ee = u.i.se & 0x7fff; - - if (!ee) { - if (x) { - x = frexpl(x*0x1p120, e); - *e -= 120; - } else *e = 0; - return x; - } else if (ee == 0x7fff) { - return x; - } - - *e = ee - 0x3ffe; - u.i.se &= 0x8000; - u.i.se |= 0x3ffe; - return u.f; -} -#endif diff --git a/src/libs/mingw.zig b/src/libs/mingw.zig index e9417eb9454b22a88ac4a4c9576b7ad6871576df..254f0e56212c86658317698ce493958291cc823d 100644 --- a/src/libs/mingw.zig +++ b/src/libs/mingw.zig @@ -613,7 +613,6 @@ const mingw32_generic_src = [_][]const u8{ "math" ++ path.sep_str ++ "fpclassify.c", "math" ++ path.sep_str ++ "fpclassifyf.c", "math" ++ path.sep_str ++ "fpclassifyl.c", - "math" ++ path.sep_str ++ "frexpl.c", "math" ++ path.sep_str ++ "ldexpf.c", "math" ++ path.sep_str ++ "lgamma.c", "math" ++ path.sep_str ++ "lgammaf.c", diff --git a/src/libs/musl.zig b/src/libs/musl.zig index 86e0995220842f3d6b35442c2bc497df26384522..73f93a82710d735f2b7998d27811c2d2be2dc0c9 100644 --- a/src/libs/musl.zig +++ b/src/libs/musl.zig @@ -839,7 +839,6 @@ const src_files = [_][]const u8{ "musl/src/math/__fpclassify.c", "musl/src/math/__fpclassifyf.c", "musl/src/math/__fpclassifyl.c", - "musl/src/math/frexpl.c", "musl/src/math/i386/acosl.s", "musl/src/math/i386/asinf.s", "musl/src/math/i386/asinl.s", diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig index d0743dbe3d10334e7aeafbef8c71bc3ae7043195..75840076a50064095e72ae64cdb46013ed756a5e 100644 --- a/src/libs/wasi_libc.zig +++ b/src/libs/wasi_libc.zig @@ -701,7 +701,6 @@ const libc_top_half_src_files = [_][]const u8{ "musl/src/math/finitef.c", "musl/src/math/fma.c", "musl/src/math/fmaf.c", - "musl/src/math/frexpl.c", "musl/src/math/ilogb.c", "musl/src/math/ilogbf.c", "musl/src/math/ilogbl.c", -- 2.54.0 From 2cbd6980e5a06ce4a2828a10b6cdb0efba827f84 Mon Sep 17 00:00:00 2001 From: mihael Date: Fri, 20 Mar 2026 17:44:23 +0100 Subject: [PATCH 06/20] `libzigc/math`: Implement `finite`, `finitef` The behaviour of `libc` `finite` functions is the same as Zig std's, so the function basically delegates to `math.isFinite` and no new tests were added. These functions are obsolete, but still part of `musl`. No results are attached in this case as `libc-test` doesn't have tests for them. --- lib/c/math.zig | 10 ++++++++++ lib/libc/musl/src/math/finite.c | 7 ------- lib/libc/musl/src/math/finitef.c | 7 ------- src/libs/musl.zig | 2 -- src/libs/wasi_libc.zig | 2 -- 5 files changed, 10 insertions(+), 18 deletions(-) delete mode 100644 lib/libc/musl/src/math/finite.c delete mode 100644 lib/libc/musl/src/math/finitef.c diff --git a/lib/c/math.zig b/lib/c/math.zig index 1f80e1cc7e5362891884988ecff8423269bf9ce8..36bdfcfa753e66d3506beefccb7484f8260a3747 100644 --- a/lib/c/math.zig +++ b/lib/c/math.zig @@ -62,6 +62,8 @@ comptime { symbol(&exp10, "exp10"); symbol(&exp10f, "exp10f"); symbol(&fdim, "fdim"); + symbol(&finite, "finite"); + symbol(&finitef, "finitef"); symbol(&frexp, "frexp"); symbol(&hypot, "hypot"); symbol(&modf, "modf"); @@ -164,6 +166,14 @@ fn fdim(x: f64, y: f64) callconv(.c) f64 { return 0; } +fn finite(x: f64) callconv(.c) c_int { + return if (math.isFinite(x)) 1 else 0; +} + +fn finitef(x: f32) callconv(.c) c_int { + return if (math.isFinite(x)) 1 else 0; +} + fn frexpGeneric(comptime T: type, x: T, e: *c_int) T { // libc expects `*e` to be unspecified in this case; an unspecified C value // should be a valid value of the relevant type, yet Zig's std diff --git a/lib/libc/musl/src/math/finite.c b/lib/libc/musl/src/math/finite.c deleted file mode 100644 index 25a0575fb9757d6ddaec5a538db7984d007fdb55..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/finite.c +++ /dev/null @@ -1,7 +0,0 @@ -#define _GNU_SOURCE -#include - -int finite(double x) -{ - return isfinite(x); -} diff --git a/lib/libc/musl/src/math/finitef.c b/lib/libc/musl/src/math/finitef.c deleted file mode 100644 index 2c4c771463a1b244244edf3ea92026f8fe78c67a..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/finitef.c +++ /dev/null @@ -1,7 +0,0 @@ -#define _GNU_SOURCE -#include - -int finitef(float x) -{ - return isfinite(x); -} diff --git a/src/libs/musl.zig b/src/libs/musl.zig index 73f93a82710d735f2b7998d27811c2d2be2dc0c9..c451312bbb9d39620f92364b173df38bc499e795 100644 --- a/src/libs/musl.zig +++ b/src/libs/musl.zig @@ -831,8 +831,6 @@ const src_files = [_][]const u8{ "musl/src/math/__expo2f.c", "musl/src/math/fdimf.c", "musl/src/math/fdiml.c", - "musl/src/math/finite.c", - "musl/src/math/finitef.c", "musl/src/math/fma.c", "musl/src/math/fmaf.c", "musl/src/math/fmal.c", diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig index 75840076a50064095e72ae64cdb46013ed756a5e..fbc5f5a142d7dcec086d2333446a3ea85bd8b17b 100644 --- a/src/libs/wasi_libc.zig +++ b/src/libs/wasi_libc.zig @@ -697,8 +697,6 @@ const libc_top_half_src_files = [_][]const u8{ "musl/src/math/expm1l.c", "musl/src/math/fdimf.c", "musl/src/math/fdiml.c", - "musl/src/math/finite.c", - "musl/src/math/finitef.c", "musl/src/math/fma.c", "musl/src/math/fmaf.c", "musl/src/math/ilogb.c", -- 2.54.0 From e15d8dd3d27047043a371421e59bc11c83a7a1b5 Mon Sep 17 00:00:00 2001 From: mihael Date: Fri, 20 Mar 2026 18:03:36 +0100 Subject: [PATCH 07/20] `libc-test/math`: Enable `rint` tests The `rint` tests were commented out when added along with other functions to the `libc-test` runner, presumably because they were failing. In the meantime, we moved to using a Zig-implemented `rint` and I found that the tests are now passing. To celebrate, I'm putting them back into commission. Having these tests in a working state is nice because functions like `lrint`, `llrint`, and `nearbyint` depend on `rint` implementation being correct. Results: ``` $ ./build/stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib $ stage4/bin/zig build test-libc -Dlibc-test-path= -Dtest-filter=rint -fqemu -fwasmtime --summary line Build Summary: 2209/2209 steps succeeded ``` --- test/libc.zig | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/libc.zig b/test/libc.zig index 423533cc868486599dfaf4b31cf3a3aed0f84fa2..ad6b1af176d5afcc9c7385a45f9964df1ff6af92 100644 --- a/test/libc.zig +++ b/test/libc.zig @@ -293,7 +293,7 @@ pub fn addCases(cases: *tests.LibcContext) void { cases.addLibcTestCase("math/remquo.c", true, .{}); cases.addLibcTestCase("math/remquof.c", true, .{}); cases.addLibcTestCase("math/remquol.c", true, .{}); - // cases.addLibcTestCase("math/rint.c", true, .{}); + cases.addLibcTestCase("math/rint.c", true, .{}); cases.addLibcTestCase("math/rintf.c", true, .{}); // cases.addLibcTestCase("math/rintl.c", true, .{}); cases.addLibcTestCase("math/round.c", true, .{}); -- 2.54.0 From 811bada06d4d88d786119bab752c0e0cff345154 Mon Sep 17 00:00:00 2001 From: mihael Date: Fri, 20 Mar 2026 18:55:37 +0100 Subject: [PATCH 08/20] `libzigc/math`: Implement `lrint` The changes were tested by running: ``` $ ./build/stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib $ stage4/bin/zig build test-libc -Dlibc-test-path= -Dtest-filter=lrint -fqemu -fwasmtime --summary line Build Summary: 737/737 steps succeeded ``` --- lib/c/math.zig | 5 ++ lib/libc/musl/src/math/aarch64/lrint.c | 10 ---- lib/libc/musl/src/math/i386/lrint.c | 8 --- lib/libc/musl/src/math/lrint.c | 72 ------------------------ lib/libc/musl/src/math/powerpc64/lrint.c | 16 ------ lib/libc/musl/src/math/x32/lrint.s | 5 -- lib/libc/musl/src/math/x86_64/lrint.c | 8 --- src/libs/musl.zig | 6 -- src/libs/wasi_libc.zig | 1 - 9 files changed, 5 insertions(+), 126 deletions(-) delete mode 100644 lib/libc/musl/src/math/aarch64/lrint.c delete mode 100644 lib/libc/musl/src/math/i386/lrint.c delete mode 100644 lib/libc/musl/src/math/lrint.c delete mode 100644 lib/libc/musl/src/math/powerpc64/lrint.c delete mode 100644 lib/libc/musl/src/math/x32/lrint.s delete mode 100644 lib/libc/musl/src/math/x86_64/lrint.c diff --git a/lib/c/math.zig b/lib/c/math.zig index 36bdfcfa753e66d3506beefccb7484f8260a3747..d3578d75dd35060a6bdf709b568ebc8f18fbbc46 100644 --- a/lib/c/math.zig +++ b/lib/c/math.zig @@ -66,6 +66,7 @@ comptime { symbol(&finitef, "finitef"); symbol(&frexp, "frexp"); symbol(&hypot, "hypot"); + symbol(&lrint, "lrint"); symbol(&modf, "modf"); symbol(&pow, "pow"); symbol(&pow10, "pow10"); @@ -229,6 +230,10 @@ fn isnanl(x: c_longdouble) callconv(.c) c_int { return if (math.isNan(x)) 1 else 0; } +fn lrint(x: f64) callconv(.c) c_long { + return @intFromFloat(rint(x)); +} + fn modfGeneric(comptime T: type, x: T, iptr: *T) T { if (math.isNegativeInf(x)) { iptr.* = -math.inf(T); diff --git a/lib/libc/musl/src/math/aarch64/lrint.c b/lib/libc/musl/src/math/aarch64/lrint.c deleted file mode 100644 index cb7785ad1e419c076a90064752d561b3dd01b0db..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/aarch64/lrint.c +++ /dev/null @@ -1,10 +0,0 @@ -#include - -long lrint(double x) -{ - long n; - __asm__ ( - "frintx %d1, %d1\n" - "fcvtzs %x0, %d1\n" : "=r"(n), "+w"(x)); - return n; -} diff --git a/lib/libc/musl/src/math/i386/lrint.c b/lib/libc/musl/src/math/i386/lrint.c deleted file mode 100644 index 89563ab26959ede3218de550b97624d7edce556f..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/i386/lrint.c +++ /dev/null @@ -1,8 +0,0 @@ -#include - -long lrint(double x) -{ - long r; - __asm__ ("fistpl %0" : "=m"(r) : "t"(x) : "st"); - return r; -} diff --git a/lib/libc/musl/src/math/lrint.c b/lib/libc/musl/src/math/lrint.c deleted file mode 100644 index ddee7a0d9a689c5e062ad6ef2a59ae6c5cad2683..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/lrint.c +++ /dev/null @@ -1,72 +0,0 @@ -#include -#include -#include -#include "libm.h" - -/* -If the result cannot be represented (overflow, nan), then -lrint raises the invalid exception. - -Otherwise if the input was not an integer then the inexact -exception is raised. - -C99 is a bit vague about whether inexact exception is -allowed to be raised when invalid is raised. -(F.9 explicitly allows spurious inexact exceptions, F.9.6.5 -does not make it clear if that rule applies to lrint, but -IEEE 754r 7.8 seems to forbid spurious inexact exception in -the ineger conversion functions) - -So we try to make sure that no spurious inexact exception is -raised in case of an overflow. - -If the bit size of long > precision of double, then there -cannot be inexact rounding in case the result overflows, -otherwise LONG_MAX and LONG_MIN can be represented exactly -as a double. -*/ - -#if LONG_MAX < 1U<<53 && defined(FE_INEXACT) -#include -#include -#if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1 -#define EPS DBL_EPSILON -#elif FLT_EVAL_METHOD==2 -#define EPS LDBL_EPSILON -#endif -#ifdef __GNUC__ -/* avoid stack frame in lrint */ -__attribute__((noinline)) -#endif -static long lrint_slow(double x) -{ - #pragma STDC FENV_ACCESS ON - int e; - - e = fetestexcept(FE_INEXACT); - x = rint(x); - if (!e && (x > LONG_MAX || x < LONG_MIN)) - feclearexcept(FE_INEXACT); - /* conversion */ - return x; -} - -long lrint(double x) -{ - uint32_t abstop = asuint64(x)>>32 & 0x7fffffff; - uint64_t sign = asuint64(x) & (1ULL << 63); - - if (abstop < 0x41dfffff) { - /* |x| < 0x7ffffc00, no overflow */ - double_t toint = asdouble(asuint64(1/EPS) | sign); - double_t y = x + toint - toint; - return (long)y; - } - return lrint_slow(x); -} -#else -long lrint(double x) -{ - return rint(x); -} -#endif diff --git a/lib/libc/musl/src/math/powerpc64/lrint.c b/lib/libc/musl/src/math/powerpc64/lrint.c deleted file mode 100644 index 4e4b2e00a9f9b69e0522433bacd50b47e7ae07e9..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/powerpc64/lrint.c +++ /dev/null @@ -1,16 +0,0 @@ -#include - -#ifdef _ARCH_PWR5X - -long lrint(double x) -{ - long n; - __asm__ ("fctid %0, %1" : "=d"(n) : "d"(x)); - return n; -} - -#else - -#include "../lrint.c" - -#endif diff --git a/lib/libc/musl/src/math/x32/lrint.s b/lib/libc/musl/src/math/x32/lrint.s deleted file mode 100644 index 15fc2454bc56c721ae3785c1eed138a3d228d66d..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/x32/lrint.s +++ /dev/null @@ -1,5 +0,0 @@ -.global lrint -.type lrint,@function -lrint: - cvtsd2si %xmm0,%rax - ret diff --git a/lib/libc/musl/src/math/x86_64/lrint.c b/lib/libc/musl/src/math/x86_64/lrint.c deleted file mode 100644 index a742fec64c90fa55d876ea7110c9d83a48e613d0..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/x86_64/lrint.c +++ /dev/null @@ -1,8 +0,0 @@ -#include - -long lrint(double x) -{ - long r; - __asm__ ("cvtsd2si %1, %0" : "=r"(r) : "x"(x)); - return r; -} diff --git a/src/libs/musl.zig b/src/libs/musl.zig index c451312bbb9d39620f92364b173df38bc499e795..e7548b16983c19b7f3be9f1cd57b0228a72338d0 100644 --- a/src/libs/musl.zig +++ b/src/libs/musl.zig @@ -787,7 +787,6 @@ const src_files = [_][]const u8{ "musl/src/math/aarch64/llrintf.c", "musl/src/math/aarch64/llround.c", "musl/src/math/aarch64/llroundf.c", - "musl/src/math/aarch64/lrint.c", "musl/src/math/aarch64/lrintf.c", "musl/src/math/aarch64/lround.c", "musl/src/math/aarch64/lroundf.c", @@ -860,7 +859,6 @@ const src_files = [_][]const u8{ "musl/src/math/i386/log1p.s", "musl/src/math/i386/log2l.s", "musl/src/math/i386/logl.s", - "musl/src/math/i386/lrint.c", "musl/src/math/i386/lrintf.c", "musl/src/math/i386/lrintl.c", "musl/src/math/i386/remainder.c", @@ -910,7 +908,6 @@ const src_files = [_][]const u8{ "musl/src/math/logbf.c", "musl/src/math/logbl.c", "musl/src/math/logl.c", - "musl/src/math/lrint.c", "musl/src/math/lrintf.c", "musl/src/math/lrintl.c", "musl/src/math/lround.c", @@ -940,7 +937,6 @@ const src_files = [_][]const u8{ "musl/src/math/pow_data.c", "musl/src/math/powerpc64/fma.c", "musl/src/math/powerpc64/fmaf.c", - "musl/src/math/powerpc64/lrint.c", "musl/src/math/powerpc64/lrintf.c", "musl/src/math/powerpc64/lround.c", "musl/src/math/powerpc64/lroundf.c", @@ -1017,7 +1013,6 @@ const src_files = [_][]const u8{ "musl/src/math/x32/logl.s", "musl/src/math/x32/lrintf.s", "musl/src/math/x32/lrintl.s", - "musl/src/math/x32/lrint.s", "musl/src/math/x32/remainderl.s", "musl/src/math/x32/rintl.s", "musl/src/math/x86_64/acosl.s", @@ -1036,7 +1031,6 @@ const src_files = [_][]const u8{ "musl/src/math/x86_64/log1pl.s", "musl/src/math/x86_64/log2l.s", "musl/src/math/x86_64/logl.s", - "musl/src/math/x86_64/lrint.c", "musl/src/math/x86_64/lrintf.c", "musl/src/math/x86_64/lrintl.c", "musl/src/math/x86_64/remainderl.c", diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig index fbc5f5a142d7dcec086d2333446a3ea85bd8b17b..691deabd088ea8f3831ba2a2561a0e067e43c5ca 100644 --- a/src/libs/wasi_libc.zig +++ b/src/libs/wasi_libc.zig @@ -732,7 +732,6 @@ const libc_top_half_src_files = [_][]const u8{ "musl/src/math/logbf.c", "musl/src/math/logbl.c", "musl/src/math/logl.c", - "musl/src/math/lrint.c", "musl/src/math/lrintf.c", "musl/src/math/lrintl.c", "musl/src/math/lround.c", -- 2.54.0 From 245c160a7b97e1dc8e0fd04e706d7abb36ac9f0c Mon Sep 17 00:00:00 2001 From: mihael Date: Fri, 20 Mar 2026 20:14:25 +0100 Subject: [PATCH 09/20] `libzigc/math`: Implement `rintf` The implementation was ported from `musl` to Zig code, and the `rint` unit tests were generalized so they could be used for `rintf` as well. This was checked both through unit tests and running `libc-test` suite: ``` $ ./build/stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib $ stage4/bin/zig build test-libc -Dlibc-test-path= -Dtest-filter=rintf -fqemu -fwasmtime --summary line Build Summary: 1657/1657 steps succeeded ``` --- lib/c/math.zig | 77 +++++++++++++++++++------- lib/libc/mingw/math/arm64/rintf.c | 12 ---- lib/libc/musl/src/math/aarch64/rintf.c | 7 --- lib/libc/musl/src/math/i386/rintf.c | 7 --- lib/libc/musl/src/math/rintf.c | 30 ---------- lib/libc/musl/src/math/s390x/rintf.c | 15 ----- src/libs/mingw.zig | 1 - src/libs/musl.zig | 4 -- 8 files changed, 56 insertions(+), 97 deletions(-) delete mode 100644 lib/libc/mingw/math/arm64/rintf.c delete mode 100644 lib/libc/musl/src/math/aarch64/rintf.c delete mode 100644 lib/libc/musl/src/math/i386/rintf.c delete mode 100644 lib/libc/musl/src/math/rintf.c delete mode 100644 lib/libc/musl/src/math/s390x/rintf.c diff --git a/lib/c/math.zig b/lib/c/math.zig index d3578d75dd35060a6bdf709b568ebc8f18fbbc46..582c7232c8e05d1ac4fba803d087ed9457df6486 100644 --- a/lib/c/math.zig +++ b/lib/c/math.zig @@ -48,6 +48,10 @@ comptime { symbol(&tanhf, "tanhf"); } + if (builtin.target.isMinGW() or builtin.target.isMuslLibC()) { + symbol(&rintf, "rintf"); + } + if (builtin.target.isMuslLibC() or builtin.target.isWasiLibC()) { symbol(&acos, "acos"); symbol(&acosf, "acosf"); @@ -348,7 +352,7 @@ fn pow10f(x: f32) callconv(.c) f32 { } fn rint(x: f64) callconv(.c) f64 { - const toint: f64 = 1.0 / @as(f64, math.floatEps(f64)); + const toint: f64 = 1.0 / math.floatEps(f64); const a: u64 = @bitCast(x); const e = a >> 52 & 0x7ff; const s = a >> 63; @@ -368,39 +372,70 @@ fn rint(x: f64) callconv(.c) f64 { return y; } -test "rint" { +fn rintf(x: f32) callconv(.c) f32 { + const toint: f32 = 1.0 / math.floatEps(f32); + const a: u32 = @bitCast(x); + const e = a >> 23 & 0xff; + const s = a >> 31; + var y: f32 = undefined; + + if (e >= 0x7f + 23) { + return x; + } + + if (s == 1) { + y = x - toint + toint; + } else { + y = x + toint - toint; + } + + if (y == 0) { + return if (s == 1) -0.0 else 0; + } + return y; +} + +fn testRint(comptime T: type) !void { + const f = switch (T) { + f32 => rintf, + f64 => rint, + else => @compileError("rint not implemented for" ++ @typeName(T)), + }; + // Positive numbers round correctly - try expectEqual(@as(f64, 42.0), rint(42.2)); - try expectEqual(@as(f64, 42.0), rint(41.8)); + try expectEqual(@as(T, 42.0), f(42.2)); + try expectEqual(@as(T, 42.0), f(41.8)); // Negative numbers round correctly - try expectEqual(@as(f64, -6.0), rint(-5.9)); - try expectEqual(@as(f64, -6.0), rint(-6.1)); + try expectEqual(@as(T, -6.0), f(-5.9)); + try expectEqual(@as(T, -6.0), f(-6.1)); // No rounding needed test - try expectEqual(@as(f64, 5.0), rint(5.0)); - try expectEqual(@as(f64, -10.0), rint(-10.0)); - try expectEqual(@as(f64, 0.0), rint(0.0)); + try expectEqual(@as(T, 5.0), f(5.0)); + try expectEqual(@as(T, -10.0), f(-10.0)); + try expectEqual(@as(T, 0.0), f(0.0)); // Very large numbers return unchanged - const large: f64 = 9007199254740992.0; // 2^53 - try expectEqual(large, rint(large)); - try expectEqual(-large, rint(-large)); + const large: T = 9007199254740992.0; // 2^53 + try expectEqual(large, f(large)); + try expectEqual(-large, f(-large)); // Small positive numbers round to zero - const pos_result = rint(0.3); - try expectEqual(@as(f64, 0.0), pos_result); - try expect(@as(u64, @bitCast(pos_result)) == 0); + const pos_result = f(0.3); + try expect(math.isPositiveZero(pos_result)); // Small negative numbers round to negative zero - const neg_result = rint(-0.3); - try expectEqual(@as(f64, 0.0), neg_result); - const bits: u64 = @bitCast(neg_result); - try expect((bits >> 63) == 1); + const neg_result = f(-0.3); + try expect(math.isNegativeZero(neg_result)); // Exact half rounds to nearest even (banker's rounding) - try expectEqual(@as(f64, 2.0), rint(2.5)); - try expectEqual(@as(f64, 4.0), rint(3.5)); + try expectEqual(@as(T, 2.0), f(2.5)); + try expectEqual(@as(T, 4.0), f(3.5)); +} + +test "rint" { + try testRint(f32); + try testRint(f64); } fn tanh(x: f64) callconv(.c) f64 { diff --git a/lib/libc/mingw/math/arm64/rintf.c b/lib/libc/mingw/math/arm64/rintf.c deleted file mode 100644 index 914769a1df063f29b89c8efa9a8a2caa46ae897b..0000000000000000000000000000000000000000 --- a/lib/libc/mingw/math/arm64/rintf.c +++ /dev/null @@ -1,12 +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 - -float rintf (float x) { - float retval = 0.0F; - __asm__ __volatile__ ("frintx %s0, %s1\n\t" : "=w" (retval) : "w" (x)); - return retval; -} diff --git a/lib/libc/musl/src/math/aarch64/rintf.c b/lib/libc/musl/src/math/aarch64/rintf.c deleted file mode 100644 index 1ae7dd254003946c1aa79eafe90b00e218a44f94..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/aarch64/rintf.c +++ /dev/null @@ -1,7 +0,0 @@ -#include - -float rintf(float x) -{ - __asm__ ("frintx %s0, %s1" : "=w"(x) : "w"(x)); - return x; -} diff --git a/lib/libc/musl/src/math/i386/rintf.c b/lib/libc/musl/src/math/i386/rintf.c deleted file mode 100644 index bb4121a4e61ce899fa1cf9d2e879eaf26785f468..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/i386/rintf.c +++ /dev/null @@ -1,7 +0,0 @@ -#include - -float rintf(float x) -{ - __asm__ ("frndint" : "+t"(x)); - return x; -} diff --git a/lib/libc/musl/src/math/rintf.c b/lib/libc/musl/src/math/rintf.c deleted file mode 100644 index 9047688d246a69c2698ce27399135e1db7ac2f65..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/rintf.c +++ /dev/null @@ -1,30 +0,0 @@ -#include -#include -#include - -#if FLT_EVAL_METHOD==0 -#define EPS FLT_EPSILON -#elif FLT_EVAL_METHOD==1 -#define EPS DBL_EPSILON -#elif FLT_EVAL_METHOD==2 -#define EPS LDBL_EPSILON -#endif -static const float_t toint = 1/EPS; - -float rintf(float x) -{ - union {float f; uint32_t i;} u = {x}; - int e = u.i>>23 & 0xff; - int s = u.i>>31; - float_t y; - - if (e >= 0x7f+23) - return x; - if (s) - y = x - toint + toint; - else - y = x + toint - toint; - if (y == 0) - return s ? -0.0f : 0.0f; - return y; -} diff --git a/lib/libc/musl/src/math/s390x/rintf.c b/lib/libc/musl/src/math/s390x/rintf.c deleted file mode 100644 index c9f13130cf9f40164be3b5ee6d38375aca696eec..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/s390x/rintf.c +++ /dev/null @@ -1,15 +0,0 @@ -#include - -#if defined(__HTM__) || __ARCH__ >= 9 - -float rintf(float x) -{ - __asm__ ("fiebr %0, 0, %1" : "=f"(x) : "f"(x)); - return x; -} - -#else - -#include "../rintf.c" - -#endif diff --git a/src/libs/mingw.zig b/src/libs/mingw.zig index 254f0e56212c86658317698ce493958291cc823d..2a3c858030e7b1f1e7f58f71f6262fecf01be14d 100644 --- a/src/libs/mingw.zig +++ b/src/libs/mingw.zig @@ -996,7 +996,6 @@ const mingw32_arm32_src = [_][]const u8{ const mingw32_arm64_src = [_][]const u8{ // mingwex "math" ++ path.sep_str ++ "arm64" ++ path.sep_str ++ "rint.c", - "math" ++ path.sep_str ++ "arm64" ++ path.sep_str ++ "rintf.c", "math" ++ path.sep_str ++ "arm64" ++ path.sep_str ++ "sincos.S", }; diff --git a/src/libs/musl.zig b/src/libs/musl.zig index e7548b16983c19b7f3be9f1cd57b0228a72338d0..e325efd931f61954b68142082cbb96a3d9b7ff24 100644 --- a/src/libs/musl.zig +++ b/src/libs/musl.zig @@ -792,7 +792,6 @@ const src_files = [_][]const u8{ "musl/src/math/aarch64/lroundf.c", "musl/src/math/aarch64/nearbyint.c", "musl/src/math/aarch64/nearbyintf.c", - "musl/src/math/aarch64/rintf.c", "musl/src/math/acosh.c", "musl/src/math/acoshl.c", "musl/src/math/acosl.c", @@ -867,7 +866,6 @@ const src_files = [_][]const u8{ "musl/src/math/i386/remquof.s", "musl/src/math/i386/remquol.s", "musl/src/math/i386/remquo.s", - "musl/src/math/i386/rintf.c", "musl/src/math/i386/rintl.c", "musl/src/math/i386/scalblnf.s", "musl/src/math/i386/scalblnl.s", @@ -955,7 +953,6 @@ const src_files = [_][]const u8{ "musl/src/math/remquo.c", "musl/src/math/remquof.c", "musl/src/math/remquol.c", - "musl/src/math/rintf.c", "musl/src/math/rintl.c", "musl/src/math/riscv32/fma.c", "musl/src/math/riscv32/fmaf.c", @@ -966,7 +963,6 @@ const src_files = [_][]const u8{ "musl/src/math/s390x/nearbyint.c", "musl/src/math/s390x/nearbyintf.c", "musl/src/math/s390x/nearbyintl.c", - "musl/src/math/s390x/rintf.c", "musl/src/math/s390x/rintl.c", "musl/src/math/scalb.c", "musl/src/math/scalbf.c", -- 2.54.0 From bd1dc8948f6093a13f1ac887b71fb71466b24f9e Mon Sep 17 00:00:00 2001 From: mihael Date: Fri, 20 Mar 2026 20:45:22 +0100 Subject: [PATCH 10/20] `libzigc/math`: Remove straggler `rint` file in vendored MinGW Seems like it was missed when implementing `rint`. This was checked by running: ``` $ ./build/stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib $ stage4/bin/zig build test-libc -Dlibc-test-path= -Dtest-filter=rint -fqemu -fwasmtime --summary line Build Summary: 2209/2209 steps succeeded ``` --- lib/c/math.zig | 2 +- lib/libc/mingw/math/arm64/rint.c | 12 ------------ src/libs/mingw.zig | 1 - 3 files changed, 1 insertion(+), 14 deletions(-) delete mode 100644 lib/libc/mingw/math/arm64/rint.c diff --git a/lib/c/math.zig b/lib/c/math.zig index 582c7232c8e05d1ac4fba803d087ed9457df6486..0621b4af9726ec9df87e55a1319c1e9ace74a677 100644 --- a/lib/c/math.zig +++ b/lib/c/math.zig @@ -49,6 +49,7 @@ comptime { } if (builtin.target.isMinGW() or builtin.target.isMuslLibC()) { + symbol(&rint, "rint"); symbol(&rintf, "rintf"); } @@ -81,7 +82,6 @@ comptime { if (builtin.target.isMuslLibC()) { symbol(©sign, "copysign"); symbol(©signf, "copysignf"); - symbol(&rint, "rint"); } symbol(©signl, "copysignl"); diff --git a/lib/libc/mingw/math/arm64/rint.c b/lib/libc/mingw/math/arm64/rint.c deleted file mode 100644 index 03a7ffa0a723b954eedfb5a688291094abf944b4..0000000000000000000000000000000000000000 --- a/lib/libc/mingw/math/arm64/rint.c +++ /dev/null @@ -1,12 +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 - -double rint (double x) { - double retval = 0.0; - __asm__ __volatile__ ("frintx %d0, %d1\n\t" : "=w" (retval) : "w" (x)); - return retval; -} diff --git a/src/libs/mingw.zig b/src/libs/mingw.zig index 2a3c858030e7b1f1e7f58f71f6262fecf01be14d..4a16f4e0a6fd49814374a2e3bd6cf1e9a66b5fb6 100644 --- a/src/libs/mingw.zig +++ b/src/libs/mingw.zig @@ -995,7 +995,6 @@ const mingw32_arm32_src = [_][]const u8{ const mingw32_arm64_src = [_][]const u8{ // mingwex - "math" ++ path.sep_str ++ "arm64" ++ path.sep_str ++ "rint.c", "math" ++ path.sep_str ++ "arm64" ++ path.sep_str ++ "sincos.S", }; -- 2.54.0 From 19334f95c1797af2ead725baf3a67070d90c8586 Mon Sep 17 00:00:00 2001 From: mihael Date: Fri, 20 Mar 2026 22:16:01 +0100 Subject: [PATCH 11/20] `libzigc/math`: Implement `lrintf` This was checked by running: ``` $ ./build/stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib $ stage4/bin/zig build test-libc -Dlibc-test-path= -Dtest-filter=lrintf -fqemu -fwasmtime --summary line Build Summary: 369/369 steps succeeded ``` --- lib/c/math.zig | 5 +++++ lib/libc/musl/src/math/aarch64/lrintf.c | 10 ---------- lib/libc/musl/src/math/i386/lrintf.c | 8 -------- lib/libc/musl/src/math/lrintf.c | 8 -------- lib/libc/musl/src/math/powerpc64/lrintf.c | 16 ---------------- lib/libc/musl/src/math/x32/lrintf.s | 5 ----- lib/libc/musl/src/math/x86_64/lrintf.c | 8 -------- src/libs/musl.zig | 6 ------ src/libs/wasi_libc.zig | 1 - 9 files changed, 5 insertions(+), 62 deletions(-) delete mode 100644 lib/libc/musl/src/math/aarch64/lrintf.c delete mode 100644 lib/libc/musl/src/math/i386/lrintf.c delete mode 100644 lib/libc/musl/src/math/lrintf.c delete mode 100644 lib/libc/musl/src/math/powerpc64/lrintf.c delete mode 100644 lib/libc/musl/src/math/x32/lrintf.s delete mode 100644 lib/libc/musl/src/math/x86_64/lrintf.c diff --git a/lib/c/math.zig b/lib/c/math.zig index 0621b4af9726ec9df87e55a1319c1e9ace74a677..a03a9a6eea894d1392863c81bc6594d496717163 100644 --- a/lib/c/math.zig +++ b/lib/c/math.zig @@ -72,6 +72,7 @@ comptime { symbol(&frexp, "frexp"); symbol(&hypot, "hypot"); symbol(&lrint, "lrint"); + symbol(&lrintf, "lrintf"); symbol(&modf, "modf"); symbol(&pow, "pow"); symbol(&pow10, "pow10"); @@ -238,6 +239,10 @@ fn lrint(x: f64) callconv(.c) c_long { return @intFromFloat(rint(x)); } +fn lrintf(x: f32) callconv(.c) c_long { + return @intFromFloat(rintf(x)); +} + fn modfGeneric(comptime T: type, x: T, iptr: *T) T { if (math.isNegativeInf(x)) { iptr.* = -math.inf(T); diff --git a/lib/libc/musl/src/math/aarch64/lrintf.c b/lib/libc/musl/src/math/aarch64/lrintf.c deleted file mode 100644 index 4d750d699d505f3319edfd765c88a65eb84cdab1..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/aarch64/lrintf.c +++ /dev/null @@ -1,10 +0,0 @@ -#include - -long lrintf(float x) -{ - long n; - __asm__ ( - "frintx %s1, %s1\n" - "fcvtzs %x0, %s1\n" : "=r"(n), "+w"(x)); - return n; -} diff --git a/lib/libc/musl/src/math/i386/lrintf.c b/lib/libc/musl/src/math/i386/lrintf.c deleted file mode 100644 index 0bbf29de06d95dd0435d8ecbc9460c49d5e597a1..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/i386/lrintf.c +++ /dev/null @@ -1,8 +0,0 @@ -#include - -long lrintf(float x) -{ - long r; - __asm__ ("fistpl %0" : "=m"(r) : "t"(x) : "st"); - return r; -} diff --git a/lib/libc/musl/src/math/lrintf.c b/lib/libc/musl/src/math/lrintf.c deleted file mode 100644 index ca0b6a46aaed3cf97d61724923593ea1b9d7aaf3..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/lrintf.c +++ /dev/null @@ -1,8 +0,0 @@ -#include - -/* uses LONG_MAX > 2^24, see comments in lrint.c */ - -long lrintf(float x) -{ - return rintf(x); -} diff --git a/lib/libc/musl/src/math/powerpc64/lrintf.c b/lib/libc/musl/src/math/powerpc64/lrintf.c deleted file mode 100644 index 9070fc03dc41d5801397ca51fd88b0fa0ae4507b..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/powerpc64/lrintf.c +++ /dev/null @@ -1,16 +0,0 @@ -#include - -#ifdef _ARCH_PWR5X - -long lrintf(float x) -{ - long n; - __asm__ ("fctid %0, %1" : "=d"(n) : "f"(x)); - return n; -} - -#else - -#include "../lrintf.c" - -#endif diff --git a/lib/libc/musl/src/math/x32/lrintf.s b/lib/libc/musl/src/math/x32/lrintf.s deleted file mode 100644 index 488423d2173e3fdc0804391337f17971923e8e5b..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/x32/lrintf.s +++ /dev/null @@ -1,5 +0,0 @@ -.global lrintf -.type lrintf,@function -lrintf: - cvtss2si %xmm0,%rax - ret diff --git a/lib/libc/musl/src/math/x86_64/lrintf.c b/lib/libc/musl/src/math/x86_64/lrintf.c deleted file mode 100644 index 2ba5639dc23f57ff6191e9782d481669ae38d449..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/x86_64/lrintf.c +++ /dev/null @@ -1,8 +0,0 @@ -#include - -long lrintf(float x) -{ - long r; - __asm__ ("cvtss2si %1, %0" : "=r"(r) : "x"(x)); - return r; -} diff --git a/src/libs/musl.zig b/src/libs/musl.zig index e325efd931f61954b68142082cbb96a3d9b7ff24..074f9e626eb62772d96e0f58a2b21951a4eb45b7 100644 --- a/src/libs/musl.zig +++ b/src/libs/musl.zig @@ -787,7 +787,6 @@ const src_files = [_][]const u8{ "musl/src/math/aarch64/llrintf.c", "musl/src/math/aarch64/llround.c", "musl/src/math/aarch64/llroundf.c", - "musl/src/math/aarch64/lrintf.c", "musl/src/math/aarch64/lround.c", "musl/src/math/aarch64/lroundf.c", "musl/src/math/aarch64/nearbyint.c", @@ -858,7 +857,6 @@ const src_files = [_][]const u8{ "musl/src/math/i386/log1p.s", "musl/src/math/i386/log2l.s", "musl/src/math/i386/logl.s", - "musl/src/math/i386/lrintf.c", "musl/src/math/i386/lrintl.c", "musl/src/math/i386/remainder.c", "musl/src/math/i386/remainderf.c", @@ -906,7 +904,6 @@ const src_files = [_][]const u8{ "musl/src/math/logbf.c", "musl/src/math/logbl.c", "musl/src/math/logl.c", - "musl/src/math/lrintf.c", "musl/src/math/lrintl.c", "musl/src/math/lround.c", "musl/src/math/lroundf.c", @@ -935,7 +932,6 @@ const src_files = [_][]const u8{ "musl/src/math/pow_data.c", "musl/src/math/powerpc64/fma.c", "musl/src/math/powerpc64/fmaf.c", - "musl/src/math/powerpc64/lrintf.c", "musl/src/math/powerpc64/lround.c", "musl/src/math/powerpc64/lroundf.c", "musl/src/math/powerpc/fma.c", @@ -1007,7 +1003,6 @@ const src_files = [_][]const u8{ "musl/src/math/x32/log1pl.s", "musl/src/math/x32/log2l.s", "musl/src/math/x32/logl.s", - "musl/src/math/x32/lrintf.s", "musl/src/math/x32/lrintl.s", "musl/src/math/x32/remainderl.s", "musl/src/math/x32/rintl.s", @@ -1027,7 +1022,6 @@ const src_files = [_][]const u8{ "musl/src/math/x86_64/log1pl.s", "musl/src/math/x86_64/log2l.s", "musl/src/math/x86_64/logl.s", - "musl/src/math/x86_64/lrintf.c", "musl/src/math/x86_64/lrintl.c", "musl/src/math/x86_64/remainderl.c", "musl/src/math/x86_64/remquol.c", diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig index 691deabd088ea8f3831ba2a2561a0e067e43c5ca..cae605d1241f47d4e4002fc61403c1823476b98a 100644 --- a/src/libs/wasi_libc.zig +++ b/src/libs/wasi_libc.zig @@ -732,7 +732,6 @@ const libc_top_half_src_files = [_][]const u8{ "musl/src/math/logbf.c", "musl/src/math/logbl.c", "musl/src/math/logl.c", - "musl/src/math/lrintf.c", "musl/src/math/lrintl.c", "musl/src/math/lround.c", "musl/src/math/lroundf.c", -- 2.54.0 From c76644caf3051dd2f837a6c40987a64e01b09c54 Mon Sep 17 00:00:00 2001 From: mihael Date: Mon, 23 Mar 2026 22:50:52 +0100 Subject: [PATCH 12/20] `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 From ffd6f6cc6e05e84a4e029b02ff13a3e84cbf1aa4 Mon Sep 17 00:00:00 2001 From: mihael Date: Tue, 24 Mar 2026 01:11:56 +0100 Subject: [PATCH 13/20] `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= -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. --- lib/compiler_rt/tan.zig | 131 +++++++++++++++++++++++++++++---------- lib/compiler_rt/trig.zig | 127 +++++++++++++++++++++++++++++++++++++ 2 files changed, 224 insertions(+), 34 deletions(-) diff --git a/lib/compiler_rt/tan.zig b/lib/compiler_rt/tan.zig index 822ffd4dec0fc21a4fac305ddd995e65ebf8fbc0..6534ca3c273bf8dbf194d1ca8b2a1b7388ed0906 100644 --- a/lib/compiler_rt/tan.zig +++ b/lib/compiler_rt/tan.zig @@ -3,6 +3,7 @@ //! //! https://git.musl-libc.org/cgit/musl/tree/src/math/tanf.c //! https://git.musl-libc.org/cgit/musl/tree/src/math/tan.c +//! https://git.musl-libc.org/cgit/musl/tree/src/math/tanl.c //! https://golang.org/src/math/tan.go const std = @import("std"); @@ -10,10 +11,13 @@ const builtin = @import("builtin"); const math = std.math; const mem = std.mem; const expect = std.testing.expect; +const expectApproxEqAbs = std.testing.expectApproxEqAbs; const kernel = @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 arch = builtin.cpu.arch; const compiler_rt = @import("../compiler_rt.zig"); @@ -116,14 +120,38 @@ pub fn tan(x: f64) callconv(.c) f64 { return kernel.__tan(y[0], y[1], n & 1 != 0); } +fn tanlGeneric(comptime T: type, x: T) T { + if (!(T == f80 or T == f128)) { + @compileError("`tanlGeneric` implemented only for `f80` and `f128`, got: " ++ T); + } + + const se = utils.ldSignExponent(x) & 0x7fff; + if (se == 0x7fff) { + return x - x; + } + + const pi_4 = 0.78539816339744830962; + if (@abs(x) < pi_4) { + if (se < 0x3fff - math.floatMantissaBits(T) / 2) { + if (compiler_rt.want_float_exceptions) { + mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120); + } + return x; + } + return kernel.__tanl(T, x, 0.0, 0); + } + + var y: [2]T = undefined; + const n = rem_pio2l(T, x, &y); + return kernel.__tanl(T, y[0], y[1], n & 1); +} + pub fn __tanx(x: f80) callconv(.c) f80 { - // TODO: more efficient implementation - return @floatCast(tanq(x)); + return tanlGeneric(f80, x); } pub fn tanq(x: f128) callconv(.c) f128 { - // TODO: more correct implementation - return tan(@floatCast(x)); + return tanlGeneric(f128, x); } pub fn tanl(x: c_longdouble) callconv(.c) c_longdouble { @@ -137,45 +165,80 @@ pub fn tanl(x: c_longdouble) callconv(.c) c_longdouble { } } -test "tan" { - try expect(tan(@as(f32, 0.0)) == tanf(0.0)); - try expect(tan(@as(f64, 0.0)) == tan(0.0)); -} - -test "tan32" { +fn testTanNormal(comptime T: type) !void { + const f = switch (T) { + f32 => tanf, + f64 => tan, + else => @compileError("unimplemented"), + }; const epsilon = 0.00001; - try expect(math.approxEqAbs(f32, tanf(0.0), 0.0, epsilon)); - try expect(math.approxEqAbs(f32, tanf(0.2), 0.202710, epsilon)); - try expect(math.approxEqAbs(f32, tanf(0.8923), 1.240422, epsilon)); - try expect(math.approxEqAbs(f32, tanf(1.5), 14.101420, epsilon)); - try expect(math.approxEqAbs(f32, tanf(37.45), -0.254397, epsilon)); - try expect(math.approxEqAbs(f32, tanf(89.123), 2.285852, epsilon)); + try expectApproxEqAbs(@as(T, 0.0), f(0.0), epsilon); + try expectApproxEqAbs(@as(T, 0.202710), f(0.2), epsilon); + try expectApproxEqAbs(@as(T, 1.240422), f(0.8923), epsilon); + try expectApproxEqAbs(@as(T, 14.101420), f(1.5), epsilon); + try expectApproxEqAbs(@as(T, -0.254397), f(37.45), epsilon); + try expectApproxEqAbs(@as(T, 2.285837), f(89.123), epsilon); } -test "tan64" { - const epsilon = 0.000001; +fn testTanSpecial(comptime T: type) !void { + const f = switch (T) { + f32 => tanf, + f64 => tan, + f80 => __tanx, + f128 => tanq, + else => @compileError("unimplemented"), + }; - try expect(math.approxEqAbs(f64, tan(0.0), 0.0, epsilon)); - try expect(math.approxEqAbs(f64, tan(0.2), 0.202710, epsilon)); - try expect(math.approxEqAbs(f64, tan(0.8923), 1.240422, epsilon)); - try expect(math.approxEqAbs(f64, tan(1.5), 14.101420, epsilon)); - try expect(math.approxEqAbs(f64, tan(37.45), -0.254397, epsilon)); - try expect(math.approxEqAbs(f64, tan(89.123), 2.2858376, epsilon)); + try expect(math.isPositiveZero(f(0.0))); + try expect(math.isNegativeZero(f(-0.0))); + try expect(math.isNan(f(math.inf(f32)))); + try expect(math.isNan(f(-math.inf(f32)))); + try expect(math.isNan(f(math.nan(f32)))); +} + +test "tan32.normal" { + try testTanNormal(f32); +} + +test "tan64.normal" { + try testTanNormal(f64); +} + +test "tan80.normal" { + const epsilon = math.floatEps(f80); + + try expectApproxEqAbs(@as(f80, 0.0), __tanx(0.0), epsilon); + try expectApproxEqAbs(@as(f80, 0.2027100355086724833213582716475345), __tanx(0.2), epsilon); + try expectApproxEqAbs(@as(f80, 1.2404217445497097995561220131857544), __tanx(0.8923), epsilon); + try expectApproxEqAbs(@as(f80, 14.10141994717171938764), __tanx(1.5), epsilon); + try expectApproxEqAbs(@as(f80, -0.25439607116885656232), __tanx(37.45), epsilon); + try expectApproxEqAbs(@as(f80, 2.2858376251355320963), __tanx(89.123), epsilon); +} + +test "tan128.normal" { + const epsilon = math.floatEps(f128); + + try expectApproxEqAbs(@as(f128, 0.0), tanq(0.0), epsilon); + try expectApproxEqAbs(@as(f128, 0.2027100355086724833213582716475345), tanq(0.2), epsilon); + try expectApproxEqAbs(@as(f128, 1.2404217445497097995561220131857544), tanq(0.8923), epsilon); + try expectApproxEqAbs(@as(f128, 14.101419947171719387646083651987755), tanq(1.5), epsilon); + try expectApproxEqAbs(@as(f128, -0.2543960711688565630469573224504774), tanq(37.45), epsilon); + try expectApproxEqAbs(@as(f128, 2.2858376251355321074066028114094292), tanq(89.123), epsilon); } test "tan32.special" { - try expect(tanf(0.0) == 0.0); - try expect(tanf(-0.0) == -0.0); - try expect(math.isNan(tanf(math.inf(f32)))); - try expect(math.isNan(tanf(-math.inf(f32)))); - try expect(math.isNan(tanf(math.nan(f32)))); + try testTanSpecial(f32); } test "tan64.special" { - try expect(tan(0.0) == 0.0); - try expect(tan(-0.0) == -0.0); - try expect(math.isNan(tan(math.inf(f64)))); - try expect(math.isNan(tan(-math.inf(f64)))); - try expect(math.isNan(tan(math.nan(f64)))); + try testTanSpecial(f64); +} + +test "tan80.special" { + try testTanSpecial(f80); +} + +test "tan128.special" { + try testTanSpecial(f128); } diff --git a/lib/compiler_rt/trig.zig b/lib/compiler_rt/trig.zig index 616dcf53f216a547bf026130d0cf46c7e6ad99b7..4d1690e4fe032f66bd8e5d68251b8e427f8f483a 100644 --- a/lib/compiler_rt/trig.zig +++ b/lib/compiler_rt/trig.zig @@ -7,6 +7,7 @@ // https://git.musl-libc.org/cgit/musl/tree/src/math/__sindf.c // https://git.musl-libc.org/cgit/musl/tree/src/math/__tand.c // https://git.musl-libc.org/cgit/musl/tree/src/math/__tandf.c +// https://git.musl-libc.org/cgit/musl/tree/src/math/__tanl.c /// kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164 /// Input x is assumed to be bounded by ~pi/4 in magnitude. @@ -271,3 +272,129 @@ pub fn __tandf(x: f64, odd: bool) f32 { const r0 = (x + s * u) + (s * w) * (t + w * r); return @floatCast(if (odd) -1.0 / r0 else r0); } + +pub fn __tanl(comptime T: type, x_: T, y_: T, odd: i32) T { + var x = x_; + var y = y_; + const impl = switch (T) { + f80 => struct { + const pio4: T = 0.785398163397448309628; + const pio4lo: T = -1.25413940316708300586e-20; + + const T3: T = 0.333333333333333333180; + const T5: T = 0.133333333333333372290; + const T7: T = 0.0539682539682504975744; + const T9: f64 = 0.021869488536312216; + const T11: f64 = 0.0088632355256619590; + const T13: f64 = 0.0035921281113786528; + const T15: f64 = 0.0014558334756312418; + const T17: f64 = 0.00059003538700862256; + const T19: f64 = 0.00023907843576635544; + const T21: f64 = 0.000097154625656538905; + const T23: f64 = 0.000038440165747303162; + const T25: f64 = 0.000018082171885432524; + const T27: f64 = 0.0000024196006108814377; + const T29: f64 = 0.0000078293456938132840; + const T31: f64 = -0.0000032609076735050182; + const T33: f64 = 0.0000023261313142559411; + + inline fn rpoly(w: T) T { + return T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + + w * (T25 + w * (T29 + w * T33)))))); + } + + inline fn vpoly(w: T) T { + return T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + + w * (T27 + w * T31))))); + } + }, + f128 => struct { + const pio4: T = 0x1.921fb54442d18469898cc51701b8p-1; + const pio4lo: T = 0x1.cd129024e088a67cc74020bbea60p-116; + + const T3: T = 0x1.5555555555555555555555555553p-2; + const T5: T = 0x1.1111111111111111111111111eb5p-3; + const T7: T = 0x1.ba1ba1ba1ba1ba1ba1ba1b694cd6p-5; + const T9: T = 0x1.664f4882c10f9f32d6bbe09d8bcdp-6; + const T11: T = 0x1.226e355e6c23c8f5b4f5762322eep-7; + const T13: T = 0x1.d6d3d0e157ddfb5fed8e84e27b37p-9; + const T15: T = 0x1.7da36452b75e2b5fce9ee7c2c92ep-10; + const T17: T = 0x1.355824803674477dfcf726649efep-11; + const T19: T = 0x1.f57d7734d1656e0aceb716f614c2p-13; + const T21: T = 0x1.967e18afcb180ed942dfdc518d6cp-14; + const T23: T = 0x1.497d8eea21e95bc7e2aa79b9f2cdp-15; + const T25: T = 0x1.0b132d39f055c81be49eff7afd50p-16; + const T27: T = 0x1.b0f72d33eff7bfa2fbc1059d90b6p-18; + const T29: T = 0x1.5ef2daf21d1113df38d0fbc00267p-19; + const T31: T = 0x1.1c77d6eac0234988cdaa04c96626p-20; + const T33: T = 0x1.cd2a5a292b180e0bdd701057dfe3p-22; + const T35: T = 0x1.75c7357d0298c01a31d0a6f7d518p-23; + const T37: T = 0x1.2f3190f4718a9a520f98f50081fcp-24; + const T39: f64 = 0.000000028443389121318352; + const T41: f64 = 0.000000011981013102001973; + const T43: f64 = 0.0000000038303578044958070; + const T45: f64 = 0.0000000034664378216909893; + const T47: f64 = -0.0000000015090641701997785; + const T49: f64 = 0.0000000029449552300483952; + const T51: f64 = -0.0000000022006995706097711; + const T53: f64 = 0.0000000015468200913196612; + const T55: f64 = -0.00000000061311613386849674; + const T57: f64 = 1.4912469681508012e-10; + + inline fn rpoly(w: T) T { + return T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + + w * (T25 + w * (T29 + w * (T33 + w * (T37 + w * (T41 + + w * (T45 + w * (T49 + w * (T53 + w * T57)))))))))))); + } + + inline fn vpoly(w: T) T { + return T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + + w * (T27 + w * (T31 + w * (T35 + w * (T39 + w * (T43 + + w * (T47 + w * (T51 + w * T55))))))))))); + } + }, + else => @compileError("__tanl supports only f80 and f128, got: " ++ @typeName(T)), + }; + + const big = @abs(x) >= 0.67434; + var sign: i8 = 0; + + if (big) { + if (x < 0) { + sign = -1; + x = -x; + y = -y; + } + x = (impl.pio4 - x) + (impl.pio4lo - y); + y = 0.0; + } + + var z = x * x; + var w = z * z; + var r = impl.rpoly(w); + var v = z * impl.vpoly(w); + var s = z * x; + r = y + z * (s * (r + v) + y) + impl.T3 * s; + w = x + r; + + if (big) { + s = @as(T, @floatFromInt(1 - 2 * odd)); + v = s - 2.0 * (x + (r - w * w / (w + s))); + return if (sign == -1) -v else v; + } + + if (odd == 0) { + return w; + } + + // if allow error up to 2 ulp, simply return + // -1.0 / (x+r) here + // + // compute -1.0 / (x+r) accurately + z = w + 0x1p32 - 0x1p32; + v = r - (z - x); + const a = -1.0 / w; + const t = a + 0x1p32 - 0x1p32; + s = 1.0 + t * z; + return t + a * (s + t * v); +} -- 2.54.0 From b5ec3e597e2cd374853dea5244b2bc2bccccd96a Mon Sep 17 00:00:00 2001 From: mihael Date: Tue, 24 Mar 2026 19:33:55 +0100 Subject: [PATCH 14/20] `libzigc/math`: Implement more precise `sinl` in `compiler_rt` The implementation was ported from `musl`. Unit tests for `f80` and `f128` were also added. The changes were tested by running: ``` $ ./build/stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib $ stage4/bin/zig build test-libc -Dlibc-test-path= -Dtest-filter=sinl -fqemu -fwasmtime --summary line Build Summary: 553/553 steps succeeded ``` --- lib/compiler_rt/math_utils.zig | 2 + lib/compiler_rt/sin.zig | 131 ++++++++++++++++++------ lib/compiler_rt/tan.zig | 3 +- lib/compiler_rt/trig.zig | 100 ++++++++++++++++++ lib/libc/mingw/math/x86/sin.def.h | 65 ------------ lib/libc/mingw/math/x86/sinl.c | 46 --------- lib/libc/mingw/math/x86/sinl_internal.S | 58 ----------- lib/libc/musl/src/math/__sinl.c | 78 -------------- lib/libc/musl/src/math/sinl.c | 41 -------- src/libs/mingw.zig | 2 - src/libs/musl.zig | 2 - src/libs/wasi_libc.zig | 2 - 12 files changed, 201 insertions(+), 329 deletions(-) delete mode 100644 lib/libc/mingw/math/x86/sin.def.h delete mode 100644 lib/libc/mingw/math/x86/sinl.c delete mode 100644 lib/libc/mingw/math/x86/sinl_internal.S delete mode 100644 lib/libc/musl/src/math/__sinl.c delete mode 100644 lib/libc/musl/src/math/sinl.c diff --git a/lib/compiler_rt/math_utils.zig b/lib/compiler_rt/math_utils.zig index 76cbcec468a74602c2a2281d1a07c9a7e1c83039..70e6a4b3b71b252db9699ead62909d65cc688c55 100644 --- a/lib/compiler_rt/math_utils.zig +++ b/lib/compiler_rt/math_utils.zig @@ -1,6 +1,8 @@ const std = @import("std"); pub const U80 = std.meta.Int(.unsigned, 80); +/// pi divided by 4 +pub const pi_4 = 0.78539816339744830962; /// Returns the sign + exponent bits of a `long double` pub fn ldSignExponent(x: anytype) u16 { diff --git a/lib/compiler_rt/sin.zig b/lib/compiler_rt/sin.zig index ef3709bc1882156614a2b69b64e2643195534d84..78a56d9b8c8c1d043f29357e5e870e81d046c772 100644 --- a/lib/compiler_rt/sin.zig +++ b/lib/compiler_rt/sin.zig @@ -3,17 +3,21 @@ //! //! https://git.musl-libc.org/cgit/musl/tree/src/math/sinf.c //! https://git.musl-libc.org/cgit/musl/tree/src/math/sin.c +//! https://git.musl-libc.org/cgit/musl/tree/src/math/sinl.c const std = @import("std"); const math = std.math; const mem = std.mem; const expect = std.testing.expect; +const expectApproxEqAbs = std.testing.expectApproxEqAbs; const compiler_rt = @import("../compiler_rt.zig"); const symbol = @import("../compiler_rt.zig").symbol; 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"); comptime { symbol(&__sinh, "__sinh"); @@ -128,14 +132,39 @@ pub fn sin(x: f64) callconv(.c) f64 { }; } +fn sinlGeneric(comptime T: type, x: T) T { + const se = utils.ldSignExponent(x) & 0x7fff; + if (se == 0x7fff) { + return x - x; + } + + if (@abs(x) < utils.pi_4) { + if (se < 0x3fff - (math.floatMantissaBits(T) / 2)) { + // raise inexact if x!=0 and underflow if subnormal + if (compiler_rt.want_float_exceptions) { + mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120); + } + return x; + } + return trig.__sinl(T, x, 0.0, 0); + } + + var y: [2]T = undefined; + const n = rem_pio2l(T, x, &y); + return switch (n & 3) { + 0 => trig.__sinl(T, y[0], y[1], 1), + 1 => trig.__cosl(T, y[0], y[1]), + 2 => -trig.__sinl(T, y[0], y[1], 1), + else => -trig.__cosl(T, y[0], y[1]), + }; +} + pub fn __sinx(x: f80) callconv(.c) f80 { - // TODO: more efficient implementation - return @floatCast(sinq(x)); + return sinlGeneric(f80, x); } pub fn sinq(x: f128) callconv(.c) f128 { - // TODO: more correct implementation - return sin(@floatCast(x)); + return sinlGeneric(f128, x); } pub fn sinl(x: c_longdouble) callconv(.c) c_longdouble { @@ -149,44 +178,80 @@ pub fn sinl(x: c_longdouble) callconv(.c) c_longdouble { } } -test "sin32" { - const epsilon = 0.00001; +fn testSinSpecial(comptime T: type) !void { + const f = switch (T) { + f32 => sinf, + f64 => sin, + f80 => __sinx, + f128 => sinq, + else => @compileError("unimplemented"), + }; - try expect(math.approxEqAbs(f32, sinf(0.0), 0.0, epsilon)); - try expect(math.approxEqAbs(f32, sinf(0.2), 0.198669, epsilon)); - try expect(math.approxEqAbs(f32, sinf(0.8923), 0.778517, epsilon)); - try expect(math.approxEqAbs(f32, sinf(1.5), 0.997495, epsilon)); - try expect(math.approxEqAbs(f32, sinf(-1.5), -0.997495, epsilon)); - try expect(math.approxEqAbs(f32, sinf(37.45), -0.246544, epsilon)); - try expect(math.approxEqAbs(f32, sinf(89.123), 0.916166, epsilon)); + try expect(math.isPositiveZero(f(0.0))); + try expect(math.isNegativeZero(f(-0.0))); + try expect(math.isNan(f(math.inf(T)))); + try expect(math.isNan(f(-math.inf(T)))); + try expect(math.isNan(f(math.nan(T)))); } -test "sin64" { - const epsilon = 0.000001; - - try expect(math.approxEqAbs(f64, sin(0.0), 0.0, epsilon)); - try expect(math.approxEqAbs(f64, sin(0.2), 0.198669, epsilon)); - try expect(math.approxEqAbs(f64, sin(0.8923), 0.778517, epsilon)); - try expect(math.approxEqAbs(f64, sin(1.5), 0.997495, epsilon)); - try expect(math.approxEqAbs(f64, sin(-1.5), -0.997495, epsilon)); - try expect(math.approxEqAbs(f64, sin(37.45), -0.246543, epsilon)); - try expect(math.approxEqAbs(f64, sin(89.123), 0.916166, epsilon)); +test "sin32.normal" { + const epsilon = math.floatEps(f32); + try expectApproxEqAbs(@as(f32, 0.0), sinf(0.0), epsilon); + try expectApproxEqAbs(@as(f32, 0.19866933), sinf(0.2), epsilon); + try expectApproxEqAbs(@as(f32, 0.77851737), sinf(0.8923), epsilon); + try expectApproxEqAbs(@as(f32, 0.997495), sinf(1.5), epsilon); + try expectApproxEqAbs(@as(f32, -0.997495), sinf(-1.5), epsilon); + try expectApproxEqAbs(@as(f32, -0.24654257), sinf(37.45), epsilon); + try expectApproxEqAbs(@as(f32, 0.9161657), sinf(89.123), epsilon); } test "sin32.special" { - try expect(sinf(0.0) == 0.0); - try expect(sinf(-0.0) == -0.0); - try expect(math.isNan(sinf(math.inf(f32)))); - try expect(math.isNan(sinf(-math.inf(f32)))); - try expect(math.isNan(sinf(math.nan(f32)))); + try testSinSpecial(f32); +} + +test "sin64.normal" { + const epsilon = math.floatEps(f64); + try expectApproxEqAbs(@as(f64, 0.0), sin(0.0), epsilon); + try expectApproxEqAbs(@as(f64, 0.19866933079506122), sin(0.2), epsilon); + try expectApproxEqAbs(@as(f64, 0.7785173385577349), sin(0.8923), epsilon); + try expectApproxEqAbs(@as(f64, 0.9974949866040544), sin(1.5), epsilon); + try expectApproxEqAbs(@as(f64, -0.9974949866040544), sin(-1.5), epsilon); + try expectApproxEqAbs(@as(f64, -0.24654331551411082), sin(37.45), epsilon); + try expectApproxEqAbs(@as(f64, 0.9161652766622714), sin(89.123), epsilon); } test "sin64.special" { - try expect(sin(0.0) == 0.0); - try expect(sin(-0.0) == -0.0); - try expect(math.isNan(sin(math.inf(f64)))); - try expect(math.isNan(sin(-math.inf(f64)))); - try expect(math.isNan(sin(math.nan(f64)))); + try testSinSpecial(f64); +} + +test "sin80.normal" { + const epsilon = math.floatEps(f80); + try expectApproxEqAbs(@as(f80, 0.0), __sinx(0.0), epsilon); + try expectApproxEqAbs(@as(f80, 0.19866933079506121545941262711838975), __sinx(0.2), epsilon); + try expectApproxEqAbs(@as(f80, 0.77851733855773487830689285621486050), __sinx(0.8923), epsilon); + try expectApproxEqAbs(@as(f80, 0.99749498660405443094172337114148732), __sinx(1.5), epsilon); + try expectApproxEqAbs(@as(f80, -0.99749498660405443094172337114148732), __sinx(-1.5), epsilon); + try expectApproxEqAbs(@as(f80, -0.24654331551411356504), __sinx(37.45), epsilon); + try expectApproxEqAbs(@as(f80, 0.91616527666226951006), __sinx(89.123), epsilon); +} + +test "sin80.special" { + try testSinSpecial(f80); +} + +test "sin128.normal" { + const epsilon = math.floatEps(f128); + try expectApproxEqAbs(@as(f128, 0.0), sinq(0.0), epsilon); + try expectApproxEqAbs(@as(f128, 0.19866933079506121545941262711838975), sinq(0.2), epsilon); + try expectApproxEqAbs(@as(f128, 0.77851733855773487830689285621486050), sinq(0.8923), epsilon); + try expectApproxEqAbs(@as(f128, 0.99749498660405443094172337114148732), sinq(1.5), epsilon); + try expectApproxEqAbs(@as(f128, -0.99749498660405443094172337114148732), sinq(-1.5), epsilon); + try expectApproxEqAbs(@as(f128, -0.24654331551411356571238581321661085), sinq(37.45), epsilon); + try expectApproxEqAbs(@as(f128, 0.91616527666226951075019849560482170), sinq(89.123), epsilon); +} + +test "sin128.special" { + try testSinSpecial(f128); } test "sin32 #9901" { diff --git a/lib/compiler_rt/tan.zig b/lib/compiler_rt/tan.zig index 6534ca3c273bf8dbf194d1ca8b2a1b7388ed0906..746982a0d6ac7657bcaffd000e92063c18d5b270 100644 --- a/lib/compiler_rt/tan.zig +++ b/lib/compiler_rt/tan.zig @@ -130,8 +130,7 @@ fn tanlGeneric(comptime T: type, x: T) T { return x - x; } - const pi_4 = 0.78539816339744830962; - if (@abs(x) < pi_4) { + if (@abs(x) < utils.pi_4) { if (se < 0x3fff - math.floatMantissaBits(T) / 2) { if (compiler_rt.want_float_exceptions) { mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120); diff --git a/lib/compiler_rt/trig.zig b/lib/compiler_rt/trig.zig index 4d1690e4fe032f66bd8e5d68251b8e427f8f483a..8a20f2c2ac03c452d0202f825b25471c8f158bab 100644 --- a/lib/compiler_rt/trig.zig +++ b/lib/compiler_rt/trig.zig @@ -7,6 +7,8 @@ // https://git.musl-libc.org/cgit/musl/tree/src/math/__sindf.c // https://git.musl-libc.org/cgit/musl/tree/src/math/__tand.c // https://git.musl-libc.org/cgit/musl/tree/src/math/__tandf.c +// https://git.musl-libc.org/cgit/musl/tree/src/math/__sinl.c +// https://git.musl-libc.org/cgit/musl/tree/src/math/__cosl.c // https://git.musl-libc.org/cgit/musl/tree/src/math/__tanl.c /// kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164 @@ -74,6 +76,52 @@ pub fn __cosdf(x: f64) f32 { return @floatCast(((1.0 + z * C0) + w * C1) + (w * z) * r); } +pub fn __cosl(comptime T: type, x: T, y: T) T { + const impl = switch (T) { + f80 => struct { + const C1: T = 0.0416666666666666666136; + + const C2: f64 = -0.0013888888888888874; + const C3: f64 = 0.000024801587301571716; + const C4: f64 = -0.00000027557319215507120; + const C5: f64 = 0.0000000020876754400407278; + const C6: f64 = -1.1470297442401303e-11; + const C7: f64 = 4.7383039476436467e-14; + + inline fn poly(z: T) T { + return z * (C1 + z * (C2 + z * (C3 + z * (C4 + + z * (C5 + z * (C6 + z * C7)))))); + } + }, + f128 => struct { + const C1: T = 0.04166666666666666666666666666666658424671; + const C2: T = -0.001388888888888888888888888888863490893732; + const C3: T = 0.00002480158730158730158730158600795304914210; + const C4: T = -0.2755731922398589065255474947078934284324e-6; + const C5: T = 0.2087675698786809897659225313136400793948e-8; + const C6: T = -0.1147074559772972315817149986812031204775e-10; + const C7: T = 0.4779477332386808976875457937252120293400e-13; + + const C8: f64 = -0.1561920696721507929516718307820958119868e-15; + const C9: f64 = 0.4110317413744594971475941557607804508039e-18; + const C10: f64 = -0.8896592467191938803288521958313920156409e-21; + const C11: f64 = 0.1601061435794535138244346256065192782581e-23; + + inline fn poly(z: T) T { + return z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 + + z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * C11)))))))))); + } + }, + else => @compileError("__cosl supports only f80 and f128, got: " ++ @typeName(T)), + }; + + const z = x * x; + const r = impl.poly(z); + const hz = 0.5 * z; + const w = 1.0 - hz; + return w + (((1.0 - w) - hz) + (z * r - x * y)); +} + /// kernel sin function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854 /// Input x is assumed to be bounded by ~pi/4 in magnitude. /// Input y is the tail of x. @@ -120,6 +168,58 @@ pub fn __sin(x: f64, y: f64, iy: i32) f64 { } } +pub fn __sinl(comptime T: type, x: T, y: T, iy: i32) T { + const impl = switch (T) { + f80 => struct { + const S1: T = -0.166666666666666666671; + + const S2: f64 = 0.0083333333333333332; + const S3: f64 = -0.00019841269841269427; + const S4: f64 = 0.0000027557319223597490; + const S5: f64 = -0.000000025052108218074604; + const S6: f64 = 1.6059006598854211e-10; + const S7: f64 = -7.6429779983024564e-13; + const S8: f64 = 2.6174587166648325e-15; + + inline fn poly(z: T) T { + return S2 + z * (S3 + z * (S4 + z * (S5 + + z * (S6 + z * (S7 + z * S8))))); + } + }, + f128 => struct { + const S1: T = -0.16666666666666666666666666666666666606732416116558; + const S2: T = 0.0083333333333333333333333333333331135404851288270047; + const S3: T = -0.00019841269841269841269841269839935785325638310428717; + const S4: T = 0.27557319223985890652557316053039946268333231205686e-5; + const S5: T = -0.25052108385441718775048214826384312253862930064745e-7; + const S6: T = 0.16059043836821614596571832194524392581082444805729e-9; + const S7: T = -0.76471637318198151807063387954939213287488216303768e-12; + const S8: T = 0.28114572543451292625024967174638477283187397621303e-14; + + const S9: f64 = -0.82206352458348947812512122163446202498005154296863e-17; + const S10: f64 = 0.19572940011906109418080609928334380560135358385256e-19; + const S11: f64 = -0.38680813379701966970673724299207480965452616911420e-22; + const S12: f64 = 0.64038150078671872796678569586315881020659912139412e-25; + + inline fn poly(z: T) T { + return S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * (S8 + + z * (S9 + z * (S10 + z * (S11 + z * S12))))))))); + } + }, + else => @compileError("__sinl supports only f80 and f128, got: " ++ @typeName(T)), + }; + + const z = x * x; + const v = z * x; + const r = impl.poly(z); + + if (iy == 0) { + return x + v * (impl.S1 + z * r); + } + + return x - ((z * (0.5 * y - v * r) - y) - v * impl.S1); +} + pub fn __sindf(x: f64) f32 { // |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]). const S1 = -0x15555554cbac77.0p-55; // -0.166666666416265235595 diff --git a/lib/libc/mingw/math/x86/sin.def.h b/lib/libc/mingw/math/x86/sin.def.h deleted file mode 100644 index 1673b34be99974755e5ef65ae78215660d39ed10..0000000000000000000000000000000000000000 --- a/lib/libc/mingw/math/x86/sin.def.h +++ /dev/null @@ -1,65 +0,0 @@ -/* - This Software is provided under the Zope Public License (ZPL) Version 2.1. - - Copyright (c) 2009, 2010 by the mingw-w64 project - - See the AUTHORS file for the list of contributors to the mingw-w64 project. - - This license has been certified as open source. It has also been designated - as GPL compatible by the Free Software Foundation (FSF). - - Redistribution and use in source and binary forms, with or without - modification, are permitted provided that the following conditions are met: - - 1. Redistributions in source code must retain the accompanying copyright - notice, this list of conditions, and the following disclaimer. - 2. Redistributions in binary form must reproduce the accompanying - copyright notice, this list of conditions, and the following disclaimer - in the documentation and/or other materials provided with the - distribution. - 3. Names of the copyright holders must not be used to endorse or promote - products derived from this software without prior written permission - from the copyright holders. - 4. The right to distribute this software or to use it for any purpose does - not give you the right to use Servicemarks (sm) or Trademarks (tm) of - the copyright holders. Use of them is covered by separate agreement - with the copyright holders. - 5. If any files are modified, you must cause the modified files to carry - prominent notices stating that you changed the files and the date of - any change. - - Disclaimer - - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY EXPRESSED - OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES - OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO - EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY DIRECT, INDIRECT, - INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT - LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, - OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF - LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING - NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, - EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -*/ - -#include "../complex/complex_internal.h" -#include - -extern long double __sinl_internal (long double); - -__FLT_TYPE -__FLT_ABI(sin) (__FLT_TYPE x) -{ - int x_class = fpclassify (x); - if (x_class == FP_NAN) - { - __FLT_RPT_DOMAIN ("sin", x, 0.0, x); - return x; - } - else if (x_class == FP_INFINITE) - { - __FLT_RPT_DOMAIN ("sin", x, 0.0, __FLT_NAN); - return __FLT_NAN; - } - return (__FLT_TYPE) __sinl_internal ((long double) x); -} diff --git a/lib/libc/mingw/math/x86/sinl.c b/lib/libc/mingw/math/x86/sinl.c deleted file mode 100644 index 0bbb71d3792f8c8144675c486ce0f0c343af2058..0000000000000000000000000000000000000000 --- a/lib/libc/mingw/math/x86/sinl.c +++ /dev/null @@ -1,46 +0,0 @@ -/* - This Software is provided under the Zope Public License (ZPL) Version 2.1. - - Copyright (c) 2009, 2010 by the mingw-w64 project - - See the AUTHORS file for the list of contributors to the mingw-w64 project. - - This license has been certified as open source. It has also been designated - as GPL compatible by the Free Software Foundation (FSF). - - Redistribution and use in source and binary forms, with or without - modification, are permitted provided that the following conditions are met: - - 1. Redistributions in source code must retain the accompanying copyright - notice, this list of conditions, and the following disclaimer. - 2. Redistributions in binary form must reproduce the accompanying - copyright notice, this list of conditions, and the following disclaimer - in the documentation and/or other materials provided with the - distribution. - 3. Names of the copyright holders must not be used to endorse or promote - products derived from this software without prior written permission - from the copyright holders. - 4. The right to distribute this software or to use it for any purpose does - not give you the right to use Servicemarks (sm) or Trademarks (tm) of - the copyright holders. Use of them is covered by separate agreement - with the copyright holders. - 5. If any files are modified, you must cause the modified files to carry - prominent notices stating that you changed the files and the date of - any change. - - Disclaimer - - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY EXPRESSED - OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES - OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO - EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY DIRECT, INDIRECT, - INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT - LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, - OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF - LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING - NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, - EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -*/ - -#define _NEW_COMPLEX_LDOUBLE 1 -#include "sin.def.h" diff --git a/lib/libc/mingw/math/x86/sinl_internal.S b/lib/libc/mingw/math/x86/sinl_internal.S deleted file mode 100644 index 6d766b098d52b3ddb85bfaaba94c7e18b7ab69a5..0000000000000000000000000000000000000000 --- a/lib/libc/mingw/math/x86/sinl_internal.S +++ /dev/null @@ -1,58 +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 "sinl_internal.S" - .text -#ifdef __x86_64__ - .align 8 -#else - .align 4 -#endif -.globl __MINGW_USYMBOL(__sinl_internal) - .def __MINGW_USYMBOL(__sinl_internal); .scl 2; .type 32; .endef -__MINGW_USYMBOL(__sinl_internal): -#ifdef __x86_64__ - fldt (%rdx) - fsin - fnstsw %ax - testl $0x400,%eax - jnz 1f - movq %rcx,%rax - movq $0,8(%rcx) - fstpt (%rcx) - ret -1: fldpi - fadd %st(0) - fxch %st(1) -2: fprem1 - fnstsw %ax - testl $0x400,%eax - jnz 2b - fstp %st(1) - fsin - movq %rcx,%rax - movq $0,8(%rcx) - fstpt (%rcx) - ret -#else - fldt 4(%esp) - fsin - fnstsw %ax - testl $0x400,%eax - jnz 1f - ret -1: fldpi - fadd %st(0) - fxch %st(1) -2: fprem1 - fnstsw %ax - testl $0x400,%eax - jnz 2b - fstp %st(1) - fsin - ret -#endif diff --git a/lib/libc/musl/src/math/__sinl.c b/lib/libc/musl/src/math/__sinl.c deleted file mode 100644 index 2525bbe86689080b4457ba7d5b1ea53793b90d9c..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/__sinl.c +++ /dev/null @@ -1,78 +0,0 @@ -/* origin: FreeBSD /usr/src/lib/msun/ld80/k_sinl.c */ -/* origin: FreeBSD /usr/src/lib/msun/ld128/k_sinl.c */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#include "libm.h" - -#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 -#if LDBL_MANT_DIG == 64 -/* - * ld80 version of __sin.c. See __sin.c for most comments. - */ -/* - * Domain [-0.7854, 0.7854], range ~[-1.89e-22, 1.915e-22] - * |sin(x)/x - s(x)| < 2**-72.1 - * - * See __cosl.c for more details about the polynomial. - */ -static const long double -S1 = -0.166666666666666666671L; /* -0xaaaaaaaaaaaaaaab.0p-66 */ -static const double -S2 = 0.0083333333333333332, /* 0x11111111111111.0p-59 */ -S3 = -0.00019841269841269427, /* -0x1a01a01a019f81.0p-65 */ -S4 = 0.0000027557319223597490, /* 0x171de3a55560f7.0p-71 */ -S5 = -0.000000025052108218074604, /* -0x1ae64564f16cad.0p-78 */ -S6 = 1.6059006598854211e-10, /* 0x161242b90243b5.0p-85 */ -S7 = -7.6429779983024564e-13, /* -0x1ae42ebd1b2e00.0p-93 */ -S8 = 2.6174587166648325e-15; /* 0x179372ea0b3f64.0p-101 */ -#define POLY(z) (S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*S8)))))) -#elif LDBL_MANT_DIG == 113 -/* - * ld128 version of __sin.c. See __sin.c for most comments. - */ -/* - * Domain [-0.7854, 0.7854], range ~[-1.53e-37, 1.659e-37] - * |sin(x)/x - s(x)| < 2**-122.1 - * - * See __cosl.c for more details about the polynomial. - */ -static const long double -S1 = -0.16666666666666666666666666666666666606732416116558L, -S2 = 0.0083333333333333333333333333333331135404851288270047L, -S3 = -0.00019841269841269841269841269839935785325638310428717L, -S4 = 0.27557319223985890652557316053039946268333231205686e-5L, -S5 = -0.25052108385441718775048214826384312253862930064745e-7L, -S6 = 0.16059043836821614596571832194524392581082444805729e-9L, -S7 = -0.76471637318198151807063387954939213287488216303768e-12L, -S8 = 0.28114572543451292625024967174638477283187397621303e-14L; -static const double -S9 = -0.82206352458348947812512122163446202498005154296863e-17, -S10 = 0.19572940011906109418080609928334380560135358385256e-19, -S11 = -0.38680813379701966970673724299207480965452616911420e-22, -S12 = 0.64038150078671872796678569586315881020659912139412e-25; -#define POLY(z) (S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*(S8+ \ - z*(S9+z*(S10+z*(S11+z*S12)))))))))) -#endif - -long double __sinl(long double x, long double y, int iy) -{ - long double z,r,v; - - z = x*x; - v = z*x; - r = POLY(z); - if (iy == 0) - return x+v*(S1+z*r); - return x-((z*(0.5*y-v*r)-y)-v*S1); -} -#endif diff --git a/lib/libc/musl/src/math/sinl.c b/lib/libc/musl/src/math/sinl.c deleted file mode 100644 index 9c0b16eed3492f1e843b9425f11b4ba57c980c7f..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/sinl.c +++ /dev/null @@ -1,41 +0,0 @@ -#include "libm.h" - -#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 -long double sinl(long double x) -{ - return sin(x); -} -#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 -long double sinl(long double x) -{ - union ldshape u = {x}; - unsigned n; - long double y[2], hi, lo; - - u.i.se &= 0x7fff; - if (u.i.se == 0x7fff) - return x - x; - if (u.f < M_PI_4) { - if (u.i.se < 0x3fff - LDBL_MANT_DIG/2) { - /* raise inexact if x!=0 and underflow if subnormal */ - FORCE_EVAL(u.i.se == 0 ? x*0x1p-120f : x+0x1p120f); - return x; - } - return __sinl(x, 0.0, 0); - } - n = __rem_pio2l(x, y); - hi = y[0]; - lo = y[1]; - switch (n & 3) { - case 0: - return __sinl(hi, lo, 1); - case 1: - return __cosl(hi, lo); - case 2: - return -__sinl(hi, lo, 1); - case 3: - default: - return -__cosl(hi, lo); - } -} -#endif diff --git a/src/libs/mingw.zig b/src/libs/mingw.zig index 4a16f4e0a6fd49814374a2e3bd6cf1e9a66b5fb6..f8070ccbf8ff494f2dadacca1923ba9337dfc22f 100644 --- a/src/libs/mingw.zig +++ b/src/libs/mingw.zig @@ -963,8 +963,6 @@ const mingw32_x86_src = [_][]const u8{ "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "scalbn.S", "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "scalbnf.S", "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "scalbnl.S", - "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "sinl.c", - "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "sinl_internal.S", // ucrtbase "math" ++ path.sep_str ++ "nextafterl.c", "math" ++ path.sep_str ++ "nexttoward.c", diff --git a/src/libs/musl.zig b/src/libs/musl.zig index 074f9e626eb62772d96e0f58a2b21951a4eb45b7..d0375fb03efc42eec6dd22552e90aaa62512463a 100644 --- a/src/libs/musl.zig +++ b/src/libs/musl.zig @@ -979,8 +979,6 @@ const src_files = [_][]const u8{ "musl/src/math/sinh.c", "musl/src/math/sinhf.c", "musl/src/math/sinhl.c", - "musl/src/math/__sinl.c", - "musl/src/math/sinl.c", "musl/src/math/__tan.c", "musl/src/math/__tandf.c", "musl/src/math/tanhl.c", diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig index cae605d1241f47d4e4002fc61403c1823476b98a..09897473d9bb3aa2ef89fc8c459990a0f7430ed1 100644 --- a/src/libs/wasi_libc.zig +++ b/src/libs/wasi_libc.zig @@ -780,8 +780,6 @@ const libc_top_half_src_files = [_][]const u8{ "musl/src/math/__sin.c", "musl/src/math/__sindf.c", "musl/src/math/sinhl.c", - "musl/src/math/__sinl.c", - "musl/src/math/sinl.c", "musl/src/math/__tan.c", "musl/src/math/__tandf.c", "musl/src/math/tanhl.c", -- 2.54.0 From a6d6af947d09fbd52662c96c6f0a12510ef15f70 Mon Sep 17 00:00:00 2001 From: mihael Date: Wed, 25 Mar 2026 00:20:21 +0100 Subject: [PATCH 15/20] `libzigc/math`: Implement more precise `cosl` in `compiler_rt` The implementation was ported from `musl`. Unit tests for `f80` and `f128` were also added. `__cosl` was already implemented in `trig.zig` while working on `sinl`. The changes were tested by running: ``` $ ./build/stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib $ stage4/bin/zig build test-libc -Dlibc-test-path= -Dtest-filter='math.cosl' -fqemu -fwasmtime --summary line Build Summary: 553/553 steps succeeded ``` --- lib/compiler_rt/cos.zig | 134 ++++++++++++++++++------ lib/libc/mingw/math/x86/cos.def.h | 65 ------------ lib/libc/mingw/math/x86/cosl.c | 46 -------- lib/libc/mingw/math/x86/cosl_internal.S | 55 ---------- lib/libc/musl/src/math/__cosl.c | 96 ----------------- lib/libc/musl/src/math/cosl.c | 39 ------- src/libs/mingw.zig | 2 - src/libs/musl.zig | 2 - src/libs/wasi_libc.zig | 2 - 9 files changed, 103 insertions(+), 338 deletions(-) delete mode 100644 lib/libc/mingw/math/x86/cos.def.h delete mode 100644 lib/libc/mingw/math/x86/cosl.c delete mode 100644 lib/libc/mingw/math/x86/cosl_internal.S delete mode 100644 lib/libc/musl/src/math/__cosl.c delete mode 100644 lib/libc/musl/src/math/cosl.c diff --git a/lib/compiler_rt/cos.zig b/lib/compiler_rt/cos.zig index 3534e7837cc2c144c167c38698115e907e6469b7..891133711467c84d655e08f09ee414b92b23ebc7 100644 --- a/lib/compiler_rt/cos.zig +++ b/lib/compiler_rt/cos.zig @@ -1,13 +1,23 @@ +//! 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/cosf.c +//! https://git.musl-libc.org/cgit/musl/tree/src/math/cos.c +//! https://git.musl-libc.org/cgit/musl/tree/src/math/cosl.c + const std = @import("std"); const math = std.math; const mem = std.mem; const expect = std.testing.expect; +const expectApproxEqAbs = std.testing.expectApproxEqAbs; const compiler_rt = @import("../compiler_rt.zig"); const symbol = @import("../compiler_rt.zig").symbol; 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"); comptime { symbol(&__cosh, "__cosh"); @@ -112,14 +122,36 @@ pub fn cos(x: f64) callconv(.c) f64 { }; } -pub fn __cosx(a: f80) callconv(.c) f80 { - // TODO: more efficient implementation - return @floatCast(cosq(a)); +fn coslGeneric(comptime T: type, x: T) T { + const se = utils.ldSignExponent(x) & 0x7fff; + if (se == 0x7fff) { + return x - x; + } + + if (@abs(x) < utils.pi_4) { + if (se < 0x3fff - math.floatMantissaBits(T)) { + // raise inexact if x!=0 + return 1.0 + x; + } + return trig.__cosl(T, x, 0.0); + } + + var y: [2]T = undefined; + const n = rem_pio2l(T, x, &y); + return switch (n & 3) { + 0 => trig.__cosl(T, y[0], y[1]), + 1 => -trig.__sinl(T, y[0], y[1], 1), + 2 => -trig.__cosl(T, y[0], y[1]), + else => trig.__sinl(T, y[0], y[1], 1), + }; +} + +pub fn __cosx(x: f80) callconv(.c) f80 { + return coslGeneric(f80, x); } -pub fn cosq(a: f128) callconv(.c) f128 { - // TODO: more correct implementation - return cos(@floatCast(a)); +pub fn cosq(x: f128) callconv(.c) f128 { + return coslGeneric(f128, x); } pub fn cosl(x: c_longdouble) callconv(.c) c_longdouble { @@ -133,38 +165,78 @@ pub fn cosl(x: c_longdouble) callconv(.c) c_longdouble { } } -test "cos32" { - const epsilon = 0.00001; +fn testCosSpecial(comptime T: type) !void { + const f = switch (T) { + f32 => cosf, + f64 => cos, + f80 => __cosx, + f128 => cosq, + else => @compileError("unimplemented"), + }; - try expect(math.approxEqAbs(f32, cosf(0.0), 1.0, epsilon)); - try expect(math.approxEqAbs(f32, cosf(0.2), 0.980067, epsilon)); - try expect(math.approxEqAbs(f32, cosf(0.8923), 0.627623, epsilon)); - try expect(math.approxEqAbs(f32, cosf(1.5), 0.070737, epsilon)); - try expect(math.approxEqAbs(f32, cosf(-1.5), 0.070737, epsilon)); - try expect(math.approxEqAbs(f32, cosf(37.45), 0.969132, epsilon)); - try expect(math.approxEqAbs(f32, cosf(89.123), 0.400798, epsilon)); + try expect(f(0.0) == 1.0); + try expect(f(-0.0) == 1.0); + try expect(math.isNan(f(math.inf(T)))); + try expect(math.isNan(f(-math.inf(T)))); + try expect(math.isNan(f(math.nan(T)))); } -test "cos64" { - const epsilon = 0.000001; - - try expect(math.approxEqAbs(f64, cos(0.0), 1.0, epsilon)); - try expect(math.approxEqAbs(f64, cos(0.2), 0.980067, epsilon)); - try expect(math.approxEqAbs(f64, cos(0.8923), 0.627623, epsilon)); - try expect(math.approxEqAbs(f64, cos(1.5), 0.070737, epsilon)); - try expect(math.approxEqAbs(f64, cos(-1.5), 0.070737, epsilon)); - try expect(math.approxEqAbs(f64, cos(37.45), 0.969132, epsilon)); - try expect(math.approxEqAbs(f64, cos(89.123), 0.40080, epsilon)); +test "cos32.normal" { + const epsilon = math.floatEps(f32); + try expectApproxEqAbs(@as(f32, 1.0), cosf(0.0), epsilon); + try expectApproxEqAbs(@as(f32, 0.9800666), cosf(0.2), epsilon); + try expectApproxEqAbs(@as(f32, 0.6276231), cosf(0.8923), epsilon); + try expectApproxEqAbs(@as(f32, 0.0707372), cosf(1.5), epsilon); + try expectApproxEqAbs(@as(f32, 0.0707372), cosf(-1.5), epsilon); + try expectApproxEqAbs(@as(f32, 0.96913195), cosf(37.45), epsilon); + try expectApproxEqAbs(@as(f32, 0.40079966), cosf(89.123), epsilon); } test "cos32.special" { - try expect(math.isNan(cosf(math.inf(f32)))); - try expect(math.isNan(cosf(-math.inf(f32)))); - try expect(math.isNan(cosf(math.nan(f32)))); + try testCosSpecial(f32); +} + +test "cos64.normal" { + const epsilon = math.floatEps(f64); + try expectApproxEqAbs(@as(f64, 1.0), cos(0.0), epsilon); + try expectApproxEqAbs(@as(f64, 0.9800665778412416), cos(0.2), epsilon); + try expectApproxEqAbs(@as(f64, 0.6276230983360804), cos(0.8923), epsilon); + try expectApproxEqAbs(@as(f64, 0.0707372016677029), cos(1.5), epsilon); + try expectApproxEqAbs(@as(f64, 0.0707372016677029), cos(-1.5), epsilon); + try expectApproxEqAbs(@as(f64, 0.9691317730707778), cos(37.45), epsilon); + try expectApproxEqAbs(@as(f64, 0.4008006809354791), cos(89.123), epsilon); } test "cos64.special" { - try expect(math.isNan(cos(math.inf(f64)))); - try expect(math.isNan(cos(-math.inf(f64)))); - try expect(math.isNan(cos(math.nan(f64)))); + try testCosSpecial(f64); +} + +test "cos80.normal" { + const epsilon = math.floatEps(f80); + try expectApproxEqAbs(@as(f80, 1.0), __cosx(0.0), epsilon); + try expectApproxEqAbs(@as(f80, 0.98006657784124163112419651674816888), __cosx(0.2), epsilon); + try expectApproxEqAbs(@as(f80, 0.62762309833608037003563995939286067), __cosx(0.8923), epsilon); + try expectApproxEqAbs(@as(f80, 0.070737201667702910088189851434268747), __cosx(1.5), epsilon); + try expectApproxEqAbs(@as(f80, 0.070737201667702910088189851434268747), __cosx(-1.5), epsilon); + try expectApproxEqAbs(@as(f80, 0.9691317730707771246), __cosx(37.45), epsilon); + try expectApproxEqAbs(@as(f80, 0.4008006809354834001), __cosx(89.123), epsilon); +} + +test "cos80.special" { + try testCosSpecial(f80); +} + +test "cos128.normal" { + const epsilon = math.floatEps(f128); + try expectApproxEqAbs(@as(f128, 1.0), cosq(0.0), epsilon); + try expectApproxEqAbs(@as(f128, 0.98006657784124163112419651674816888), cosq(0.2), epsilon); + try expectApproxEqAbs(@as(f128, 0.62762309833608037003563995939286067), cosq(0.8923), epsilon); + try expectApproxEqAbs(@as(f128, 0.070737201667702910088189851434268747), cosq(1.5), epsilon); + try expectApproxEqAbs(@as(f128, 0.070737201667702910088189851434268747), cosq(-1.5), epsilon); + try expectApproxEqAbs(@as(f128, 0.96913177307077712443149563847233230), cosq(37.45), epsilon); + try expectApproxEqAbs(@as(f128, 0.40080068093548339848199454493704702), cosq(89.123), epsilon); +} + +test "cos128.special" { + try testCosSpecial(f128); } diff --git a/lib/libc/mingw/math/x86/cos.def.h b/lib/libc/mingw/math/x86/cos.def.h deleted file mode 100644 index f082f3740845a0351dbef31c308aab6fa18756b1..0000000000000000000000000000000000000000 --- a/lib/libc/mingw/math/x86/cos.def.h +++ /dev/null @@ -1,65 +0,0 @@ -/* - This Software is provided under the Zope Public License (ZPL) Version 2.1. - - Copyright (c) 2009, 2010 by the mingw-w64 project - - See the AUTHORS file for the list of contributors to the mingw-w64 project. - - This license has been certified as open source. It has also been designated - as GPL compatible by the Free Software Foundation (FSF). - - Redistribution and use in source and binary forms, with or without - modification, are permitted provided that the following conditions are met: - - 1. Redistributions in source code must retain the accompanying copyright - notice, this list of conditions, and the following disclaimer. - 2. Redistributions in binary form must reproduce the accompanying - copyright notice, this list of conditions, and the following disclaimer - in the documentation and/or other materials provided with the - distribution. - 3. Names of the copyright holders must not be used to endorse or promote - products derived from this software without prior written permission - from the copyright holders. - 4. The right to distribute this software or to use it for any purpose does - not give you the right to use Servicemarks (sm) or Trademarks (tm) of - the copyright holders. Use of them is covered by separate agreement - with the copyright holders. - 5. If any files are modified, you must cause the modified files to carry - prominent notices stating that you changed the files and the date of - any change. - - Disclaimer - - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY EXPRESSED - OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES - OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO - EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY DIRECT, INDIRECT, - INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT - LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, - OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF - LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING - NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, - EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -*/ - -#include "../complex/complex_internal.h" -#include - -extern long double __cosl_internal (long double); - -__FLT_TYPE -__FLT_ABI(cos) (__FLT_TYPE x) -{ - int x_class = fpclassify (x); - if (x_class == FP_NAN) - { - __FLT_RPT_DOMAIN ("cos", x, 0.0, x); - return x; - } - else if (x_class == FP_INFINITE) - { - __FLT_RPT_DOMAIN ("cos", x, 0.0, __FLT_NAN); - return __FLT_NAN; - } - return (__FLT_TYPE) __cosl_internal ((long double) x); -} diff --git a/lib/libc/mingw/math/x86/cosl.c b/lib/libc/mingw/math/x86/cosl.c deleted file mode 100644 index f798862d3d934dc25870bea9288fdd18b42d554d..0000000000000000000000000000000000000000 --- a/lib/libc/mingw/math/x86/cosl.c +++ /dev/null @@ -1,46 +0,0 @@ -/* - This Software is provided under the Zope Public License (ZPL) Version 2.1. - - Copyright (c) 2009, 2010 by the mingw-w64 project - - See the AUTHORS file for the list of contributors to the mingw-w64 project. - - This license has been certified as open source. It has also been designated - as GPL compatible by the Free Software Foundation (FSF). - - Redistribution and use in source and binary forms, with or without - modification, are permitted provided that the following conditions are met: - - 1. Redistributions in source code must retain the accompanying copyright - notice, this list of conditions, and the following disclaimer. - 2. Redistributions in binary form must reproduce the accompanying - copyright notice, this list of conditions, and the following disclaimer - in the documentation and/or other materials provided with the - distribution. - 3. Names of the copyright holders must not be used to endorse or promote - products derived from this software without prior written permission - from the copyright holders. - 4. The right to distribute this software or to use it for any purpose does - not give you the right to use Servicemarks (sm) or Trademarks (tm) of - the copyright holders. Use of them is covered by separate agreement - with the copyright holders. - 5. If any files are modified, you must cause the modified files to carry - prominent notices stating that you changed the files and the date of - any change. - - Disclaimer - - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY EXPRESSED - OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES - OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO - EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY DIRECT, INDIRECT, - INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT - LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, - OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF - LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING - NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, - EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -*/ - -#define _NEW_COMPLEX_LDOUBLE 1 -#include "cos.def.h" diff --git a/lib/libc/mingw/math/x86/cosl_internal.S b/lib/libc/mingw/math/x86/cosl_internal.S deleted file mode 100644 index 3c8f60d143e1cde3abc31488bf288aa6d1c64876..0000000000000000000000000000000000000000 --- a/lib/libc/mingw/math/x86/cosl_internal.S +++ /dev/null @@ -1,55 +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 "cosl_internal.S" - .text -#ifdef __x86_64__ - .align 8 -#else - .align 4 -#endif -.globl __MINGW_USYMBOL(__cosl_internal) - .def __MINGW_USYMBOL(__cosl_internal); .scl 2; .type 32; .endef -__MINGW_USYMBOL(__cosl_internal): -#ifdef __x86_64__ - fldt (%rdx) - fcos - fnstsw %ax - testl $0x400,%eax - jz 1f - fldpi - fadd %st(0) - fxch %st(1) -2: fprem1 - fnstsw %ax - testl $0x400,%eax - jnz 2b - fstp %st(1) - fcos -1: movq %rcx,%rax - movq $0,8(%rcx) - fstpt (%rcx) - ret -#else - fldt 4(%esp) - fcos - fnstsw %ax - testl $0x400,%eax - jnz 1f - ret -1: fldpi - fadd %st(0) - fxch %st(1) -2: fprem1 - fnstsw %ax - testl $0x400,%eax - jnz 2b - fstp %st(1) - fcos - ret -#endif - diff --git a/lib/libc/musl/src/math/__cosl.c b/lib/libc/musl/src/math/__cosl.c deleted file mode 100644 index fa522ddd7ae263b02a0dfb00b04486aa4ca90af5..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/__cosl.c +++ /dev/null @@ -1,96 +0,0 @@ -/* origin: FreeBSD /usr/src/lib/msun/ld80/k_cosl.c */ -/* origin: FreeBSD /usr/src/lib/msun/ld128/k_cosl.c */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - - -#include "libm.h" - -#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 -#if LDBL_MANT_DIG == 64 -/* - * ld80 version of __cos.c. See __cos.c for most comments. - */ -/* - * Domain [-0.7854, 0.7854], range ~[-2.43e-23, 2.425e-23]: - * |cos(x) - c(x)| < 2**-75.1 - * - * The coefficients of c(x) were generated by a pari-gp script using - * a Remez algorithm that searches for the best higher coefficients - * after rounding leading coefficients to a specified precision. - * - * Simpler methods like Chebyshev or basic Remez barely suffice for - * cos() in 64-bit precision, because we want the coefficient of x^2 - * to be precisely -0.5 so that multiplying by it is exact, and plain - * rounding of the coefficients of a good polynomial approximation only - * gives this up to about 64-bit precision. Plain rounding also gives - * a mediocre approximation for the coefficient of x^4, but a rounding - * error of 0.5 ulps for this coefficient would only contribute ~0.01 - * ulps to the final error, so this is unimportant. Rounding errors in - * higher coefficients are even less important. - * - * In fact, coefficients above the x^4 one only need to have 53-bit - * precision, and this is more efficient. We get this optimization - * almost for free from the complications needed to search for the best - * higher coefficients. - */ -static const long double -C1 = 0.0416666666666666666136L; /* 0xaaaaaaaaaaaaaa9b.0p-68 */ -static const double -C2 = -0.0013888888888888874, /* -0x16c16c16c16c10.0p-62 */ -C3 = 0.000024801587301571716, /* 0x1a01a01a018e22.0p-68 */ -C4 = -0.00000027557319215507120, /* -0x127e4fb7602f22.0p-74 */ -C5 = 0.0000000020876754400407278, /* 0x11eed8caaeccf1.0p-81 */ -C6 = -1.1470297442401303e-11, /* -0x19393412bd1529.0p-89 */ -C7 = 4.7383039476436467e-14; /* 0x1aac9d9af5c43e.0p-97 */ -#define POLY(z) (z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*C7))))))) -#elif LDBL_MANT_DIG == 113 -/* - * ld128 version of __cos.c. See __cos.c for most comments. - */ -/* - * Domain [-0.7854, 0.7854], range ~[-1.80e-37, 1.79e-37]: - * |cos(x) - c(x))| < 2**-122.0 - * - * 113-bit precision requires more care than 64-bit precision, since - * simple methods give a minimax polynomial with coefficient for x^2 - * that is 1 ulp below 0.5, but we want it to be precisely 0.5. See - * above for more details. - */ -static const long double -C1 = 0.04166666666666666666666666666666658424671L, -C2 = -0.001388888888888888888888888888863490893732L, -C3 = 0.00002480158730158730158730158600795304914210L, -C4 = -0.2755731922398589065255474947078934284324e-6L, -C5 = 0.2087675698786809897659225313136400793948e-8L, -C6 = -0.1147074559772972315817149986812031204775e-10L, -C7 = 0.4779477332386808976875457937252120293400e-13L; -static const double -C8 = -0.1561920696721507929516718307820958119868e-15, -C9 = 0.4110317413744594971475941557607804508039e-18, -C10 = -0.8896592467191938803288521958313920156409e-21, -C11 = 0.1601061435794535138244346256065192782581e-23; -#define POLY(z) (z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*(C7+ \ - z*(C8+z*(C9+z*(C10+z*C11))))))))))) -#endif - -long double __cosl(long double x, long double y) -{ - long double hz,z,r,w; - - z = x*x; - r = POLY(z); - hz = 0.5*z; - w = 1.0-hz; - return w + (((1.0-w)-hz) + (z*r-x*y)); -} -#endif diff --git a/lib/libc/musl/src/math/cosl.c b/lib/libc/musl/src/math/cosl.c deleted file mode 100644 index 79c41c77fa6a0addc57e33c62f31d879392b36e6..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/cosl.c +++ /dev/null @@ -1,39 +0,0 @@ -#include "libm.h" - -#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 -long double cosl(long double x) { - return cos(x); -} -#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 -long double cosl(long double x) -{ - union ldshape u = {x}; - unsigned n; - long double y[2], hi, lo; - - u.i.se &= 0x7fff; - if (u.i.se == 0x7fff) - return x - x; - x = u.f; - if (x < M_PI_4) { - if (u.i.se < 0x3fff - LDBL_MANT_DIG) - /* raise inexact if x!=0 */ - return 1.0 + x; - return __cosl(x, 0); - } - n = __rem_pio2l(x, y); - hi = y[0]; - lo = y[1]; - switch (n & 3) { - case 0: - return __cosl(hi, lo); - case 1: - return -__sinl(hi, lo, 1); - case 2: - return -__cosl(hi, lo); - case 3: - default: - return __sinl(hi, lo, 1); - } -} -#endif diff --git a/src/libs/mingw.zig b/src/libs/mingw.zig index f8070ccbf8ff494f2dadacca1923ba9337dfc22f..f0a44ca4be1082811abd85af6026567ed8a80d7c 100644 --- a/src/libs/mingw.zig +++ b/src/libs/mingw.zig @@ -940,8 +940,6 @@ const mingw32_x86_src = [_][]const u8{ "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "atan2l.c", "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "atanhl.c", "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "atanl.c", - "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "cosl.c", - "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "cosl_internal.S", "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "cossinl.c", "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "exp2l.S", "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "expl.c", diff --git a/src/libs/musl.zig b/src/libs/musl.zig index d0375fb03efc42eec6dd22552e90aaa62512463a..f460fabc91b4fcccb6f8adab1378329826f31a48 100644 --- a/src/libs/musl.zig +++ b/src/libs/musl.zig @@ -811,8 +811,6 @@ const src_files = [_][]const u8{ "musl/src/math/__cos.c", "musl/src/math/__cosdf.c", "musl/src/math/coshl.c", - "musl/src/math/__cosl.c", - "musl/src/math/cosl.c", "musl/src/math/erf.c", "musl/src/math/erff.c", "musl/src/math/erfl.c", diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig index 09897473d9bb3aa2ef89fc8c459990a0f7430ed1..ae092cb60b6c17fb61f417278541f490ec8a00b1 100644 --- a/src/libs/wasi_libc.zig +++ b/src/libs/wasi_libc.zig @@ -682,8 +682,6 @@ const libc_top_half_src_files = [_][]const u8{ "musl/src/math/__cos.c", "musl/src/math/__cosdf.c", "musl/src/math/coshl.c", - "musl/src/math/__cosl.c", - "musl/src/math/cosl.c", "musl/src/math/erf.c", "musl/src/math/erff.c", "musl/src/math/erfl.c", -- 2.54.0 From 9a53be5fb444873ee561eb54a0ea1aaeaaad844d Mon Sep 17 00:00:00 2001 From: mihael Date: Thu, 26 Mar 2026 21:43:55 +0100 Subject: [PATCH 16/20] `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 From f6ac1a6e05cbd46f646cf09e036921ca282694cf Mon Sep 17 00:00:00 2001 From: mihael Date: Thu, 26 Mar 2026 23:38:26 +0100 Subject: [PATCH 17/20] `compiler_rt`: Symbolically link `cosl` & `sinl` impl to `__cosl` & `__sinl` Some of the functions (e.g., `lgammal`) depended on the `__sinl` and `__cosl` functions being present. --- lib/compiler_rt/cos.zig | 1 + lib/compiler_rt/sin.zig | 1 + 2 files changed, 2 insertions(+) diff --git a/lib/compiler_rt/cos.zig b/lib/compiler_rt/cos.zig index 891133711467c84d655e08f09ee414b92b23ebc7..3c4bc4adb43e05584c08d0842217249967c3d75a 100644 --- a/lib/compiler_rt/cos.zig +++ b/lib/compiler_rt/cos.zig @@ -21,6 +21,7 @@ const utils = @import("math_utils.zig"); comptime { symbol(&__cosh, "__cosh"); + symbol(&cosl, "__cosl"); symbol(&cosf, "cosf"); symbol(&cos, "cos"); symbol(&__cosx, "__cosx"); diff --git a/lib/compiler_rt/sin.zig b/lib/compiler_rt/sin.zig index 78a56d9b8c8c1d043f29357e5e870e81d046c772..ca8c9943bfcad2d047d2dea27f4d9dba8a17243b 100644 --- a/lib/compiler_rt/sin.zig +++ b/lib/compiler_rt/sin.zig @@ -21,6 +21,7 @@ const utils = @import("math_utils.zig"); comptime { symbol(&__sinh, "__sinh"); + symbol(&sinl, "__sinl"); symbol(&sinf, "sinf"); symbol(&sin, "sin"); symbol(&__sinx, "__sinx"); -- 2.54.0 From 61161132b63d641657df3f099511e2fa5062d0d1 Mon Sep 17 00:00:00 2001 From: mihael Date: Fri, 27 Mar 2026 22:45:56 +0100 Subject: [PATCH 18/20] `compiler_rt`: Dissolve `math_utils.zig` into `long_double.zig` & `trig.zig` The closest namespace the pi/4 constant could belong to is `trig.zig` since it's used across trig function implementations. On the other hand, chucking `long double` bit slicing functions into `trig.zig` seems a little more awkward, so they're put into their own namespace. --- lib/compiler_rt/cos.zig | 6 +++--- .../{math_utils.zig => long_double.zig} | 12 ++++++------ lib/compiler_rt/rem_pio2l.zig | 16 ++++++++-------- lib/compiler_rt/sin.zig | 6 +++--- lib/compiler_rt/sincos.zig | 6 +++--- lib/compiler_rt/tan.zig | 6 +++--- lib/compiler_rt/trig.zig | 3 +++ 7 files changed, 29 insertions(+), 26 deletions(-) rename lib/compiler_rt/{math_utils.zig => long_double.zig} (68%) diff --git a/lib/compiler_rt/cos.zig b/lib/compiler_rt/cos.zig index 3c4bc4adb43e05584c08d0842217249967c3d75a..582de530f5186b9d435081ad117ed6cf2d7e7f40 100644 --- a/lib/compiler_rt/cos.zig +++ b/lib/compiler_rt/cos.zig @@ -17,7 +17,7 @@ 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 ld = @import("long_double.zig"); comptime { symbol(&__cosh, "__cosh"); @@ -124,12 +124,12 @@ pub fn cos(x: f64) callconv(.c) f64 { } fn coslGeneric(comptime T: type, x: T) T { - const se = utils.ldSignExponent(x) & 0x7fff; + const se = ld.signExponent(x) & 0x7fff; if (se == 0x7fff) { return x - x; } - if (@abs(x) < utils.pi_4) { + if (@abs(x) < trig.pi_4) { if (se < 0x3fff - math.floatMantissaBits(T)) { // raise inexact if x!=0 return 1.0 + x; diff --git a/lib/compiler_rt/math_utils.zig b/lib/compiler_rt/long_double.zig similarity index 68% rename from lib/compiler_rt/math_utils.zig rename to lib/compiler_rt/long_double.zig index 70e6a4b3b71b252db9699ead62909d65cc688c55..8c018b529fe90c4815f575c5b3c0bf447c1d04b2 100644 --- a/lib/compiler_rt/math_utils.zig +++ b/lib/compiler_rt/long_double.zig @@ -1,11 +1,11 @@ +//! Utilities for dealing with the `long double` type (`f80` or `f128`) + const std = @import("std"); pub const U80 = std.meta.Int(.unsigned, 80); -/// pi divided by 4 -pub const pi_4 = 0.78539816339744830962; /// Returns the sign + exponent bits of a `long double` -pub fn ldSignExponent(x: anytype) u16 { +pub fn signExponent(x: anytype) u16 { const T = @TypeOf(x); switch (T) { f80 => { @@ -16,12 +16,12 @@ pub fn ldSignExponent(x: anytype) u16 { const bits: u128 = @bitCast(x); return @intCast(bits >> 112); }, - else => @compileError("`ldSignExponent` supports only `f80` and `f128`, got: " ++ @typeName(T)), + else => @compileError("`signExponent` 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 { +pub fn mantissaTop(x: anytype) u16 { const T = @TypeOf(x); switch (T) { f80 => { @@ -32,6 +32,6 @@ pub fn ldMantissaTop(x: anytype) u16 { const bits: u128 = @bitCast(x); return @intCast((bits >> 96) & 0xFFFF); }, - else => @compileError("`ldMantissaTop` supports only `f80` and `f128`, got: " ++ @typeName(T)), + else => @compileError("`mantissaTop` supports only `f80` and `f128`, got: " ++ @typeName(T)), } } diff --git a/lib/compiler_rt/rem_pio2l.zig b/lib/compiler_rt/rem_pio2l.zig index 1ca9d170bcdcbbfbc065f23367581d9afc63e980..fe878f116d91204b000f6fba1d2e5280b15723da 100644 --- a/lib/compiler_rt/rem_pio2l.zig +++ b/lib/compiler_rt/rem_pio2l.zig @@ -6,7 +6,7 @@ const std = @import("std"); const math = std.math; -const utils = @import("math_utils.zig"); +const ld = @import("long_double.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 { @@ -34,8 +34,8 @@ pub fn rem_pio2l(comptime T: type, x: T, y: *[2]T) i32 { 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 se = ld.signExponent(x_val); + const top = ld.mantissaTop(x_val); const lhs = (@as(u32, se & 0x7fff) << 16) | top; const rhs: u32 = ((0x3fff + 25) << 16) | 0x921f >> 1 | 0x8000; return lhs < rhs; @@ -62,8 +62,8 @@ pub fn rem_pio2l(comptime T: type, x: T, y: *[2]T) i32 { 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 se = ld.signExponent(x_val); + const top = ld.mantissaTop(x_val); const lhs = (@as(u32, se & 0x7fff) << 16) | top; const rhs: u32 = ((0x3fff + 45) << 16) | 0x921f; return lhs < rhs; @@ -77,7 +77,7 @@ pub fn rem_pio2l(comptime T: type, x: T, y: *[2]T) i32 { else => @compileError("rem_pio2l supports only f80 and f128, got: " ++ @typeName(T)), }; - const x_se = utils.ldSignExponent(x); + const x_se = ld.signExponent(x); const ex: i32 = @intCast(x_se & 0x7fff); if (impl.small(x)) { @@ -105,14 +105,14 @@ pub fn rem_pio2l(comptime T: type, x: T, y: *[2]T) i32 { y[0] = r - w; - const ey: i32 = @intCast(utils.ldSignExponent(y[0]) & 0x7fff); + const ey: i32 = @intCast(ld.signExponent(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); + const ey2: i32 = @intCast(ld.signExponent(y[0]) & 0x7fff); if (ex - ey2 > impl.round2) { t = r; w = fn_ * impl.pio2_3; diff --git a/lib/compiler_rt/sin.zig b/lib/compiler_rt/sin.zig index ca8c9943bfcad2d047d2dea27f4d9dba8a17243b..0cff43997e06e31606cd0be4ffb66d3cb36203a2 100644 --- a/lib/compiler_rt/sin.zig +++ b/lib/compiler_rt/sin.zig @@ -17,7 +17,7 @@ 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 ld = @import("long_double.zig"); comptime { symbol(&__sinh, "__sinh"); @@ -134,12 +134,12 @@ pub fn sin(x: f64) callconv(.c) f64 { } fn sinlGeneric(comptime T: type, x: T) T { - const se = utils.ldSignExponent(x) & 0x7fff; + const se = ld.signExponent(x) & 0x7fff; if (se == 0x7fff) { return x - x; } - if (@abs(x) < utils.pi_4) { + if (@abs(x) < trig.pi_4) { if (se < 0x3fff - (math.floatMantissaBits(T) / 2)) { // raise inexact if x!=0 and underflow if subnormal if (compiler_rt.want_float_exceptions) { diff --git a/lib/compiler_rt/sincos.zig b/lib/compiler_rt/sincos.zig index bc43525fbe73143978bbb172aa1b7408d3fa6470..1616cc5556459503d5ac6fd9eeab045a9ba52c07 100644 --- a/lib/compiler_rt/sincos.zig +++ b/lib/compiler_rt/sincos.zig @@ -9,7 +9,7 @@ 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 ld = @import("long_double.zig"); const compiler_rt = @import("../compiler_rt.zig"); const symbol = compiler_rt.symbol; @@ -197,7 +197,7 @@ fn sincoslGeneric(comptime T: type, x: T, r_sin: *T, r_cos: *T) void { @compileError("`sincoslGeneric` implemented only for `f80` and `f128`, got: " ++ @typeName(T)); } - const se = utils.ldSignExponent(x) & 0x7fff; + const se = ld.signExponent(x) & 0x7fff; if (se == 0x7fff) { const result = x - x; r_sin.* = result; @@ -205,7 +205,7 @@ fn sincoslGeneric(comptime T: type, x: T, r_sin: *T, r_cos: *T) void { return; } - if (@abs(x) < utils.pi_4) { + if (@abs(x) < trig.pi_4) { if (se < 0x3fff - math.floatMantissaBits(T)) { // raise underflow if subnormal if (compiler_rt.want_float_exceptions and se == 0) { diff --git a/lib/compiler_rt/tan.zig b/lib/compiler_rt/tan.zig index 746982a0d6ac7657bcaffd000e92063c18d5b270..1f925e3e8c6c6abac73d65d33f39ff2680377b03 100644 --- a/lib/compiler_rt/tan.zig +++ b/lib/compiler_rt/tan.zig @@ -17,7 +17,7 @@ const kernel = @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 ld = @import("long_double.zig"); const arch = builtin.cpu.arch; const compiler_rt = @import("../compiler_rt.zig"); @@ -125,12 +125,12 @@ fn tanlGeneric(comptime T: type, x: T) T { @compileError("`tanlGeneric` implemented only for `f80` and `f128`, got: " ++ T); } - const se = utils.ldSignExponent(x) & 0x7fff; + const se = ld.signExponent(x) & 0x7fff; if (se == 0x7fff) { return x - x; } - if (@abs(x) < utils.pi_4) { + if (@abs(x) < kernel.pi_4) { if (se < 0x3fff - math.floatMantissaBits(T) / 2) { if (compiler_rt.want_float_exceptions) { mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120); diff --git a/lib/compiler_rt/trig.zig b/lib/compiler_rt/trig.zig index 8a20f2c2ac03c452d0202f825b25471c8f158bab..e68089baeab1b6a9ff4cb818e81db5ca7e14952a 100644 --- a/lib/compiler_rt/trig.zig +++ b/lib/compiler_rt/trig.zig @@ -11,6 +11,9 @@ // https://git.musl-libc.org/cgit/musl/tree/src/math/__cosl.c // https://git.musl-libc.org/cgit/musl/tree/src/math/__tanl.c +/// pi divided by 4 +pub const pi_4 = 0.78539816339744830962; + /// kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164 /// Input x is assumed to be bounded by ~pi/4 in magnitude. /// Input y is the tail of x. -- 2.54.0 From ad10c7600793547393848528e7d4a608d9fde59e Mon Sep 17 00:00:00 2001 From: mihael Date: Thu, 2 Apr 2026 21:51:02 +0200 Subject: [PATCH 19/20] `compiler_rt`: Stop using dunders for trig functions names Additionally, use `std.math.pi / 4.0` for `pi/4` in `trig.zig` instead of hardcoding the value. --- lib/compiler_rt/cos.zig | 70 +++++++++++++++++------------------ lib/compiler_rt/sin.zig | 70 +++++++++++++++++------------------ lib/compiler_rt/sincos.zig | 76 +++++++++++++++++++------------------- lib/compiler_rt/tan.zig | 46 +++++++++++------------ lib/compiler_rt/trig.zig | 29 ++++++++------- 5 files changed, 146 insertions(+), 145 deletions(-) diff --git a/lib/compiler_rt/cos.zig b/lib/compiler_rt/cos.zig index 582de530f5186b9d435081ad117ed6cf2d7e7f40..aa9282625fba5a3b95aa0b55cb7fc23a6df090ed 100644 --- a/lib/compiler_rt/cos.zig +++ b/lib/compiler_rt/cos.zig @@ -20,11 +20,11 @@ const rem_pio2l = @import("rem_pio2l.zig").rem_pio2l; const ld = @import("long_double.zig"); comptime { - symbol(&__cosh, "__cosh"); + symbol(&cosh, "__cosh"); symbol(&cosl, "__cosl"); symbol(&cosf, "cosf"); symbol(&cos, "cos"); - symbol(&__cosx, "__cosx"); + symbol(&cosx, "__cosx"); if (compiler_rt.want_ppc_abi) { symbol(&cosq, "cosf128"); } @@ -32,7 +32,7 @@ comptime { symbol(&cosl, "cosl"); } -pub fn __cosh(a: f16) callconv(.c) f16 { +pub fn cosh(a: f16) callconv(.c) f16 { // TODO: more efficient implementation return @floatCast(cosf(a)); } @@ -54,27 +54,27 @@ pub fn cosf(x: f32) callconv(.c) f32 { if (compiler_rt.want_float_exceptions) mem.doNotOptimizeAway(x + 0x1p120); return 1.0; } - return trig.__cosdf(x); + return trig.cosdf(x); } if (ix <= 0x407b53d1) { // |x| ~<= 5*pi/4 if (ix > 0x4016cbe3) { // |x| ~> 3*pi/4 - return -trig.__cosdf(if (sign) x + c2pio2 else x - c2pio2); + return -trig.cosdf(if (sign) x + c2pio2 else x - c2pio2); } else { if (sign) { - return trig.__sindf(x + c1pio2); + return trig.sindf(x + c1pio2); } else { - return trig.__sindf(c1pio2 - x); + return trig.sindf(c1pio2 - x); } } } if (ix <= 0x40e231d5) { // |x| ~<= 9*pi/4 if (ix > 0x40afeddf) { // |x| ~> 7*pi/4 - return trig.__cosdf(if (sign) x + c4pio2 else x - c4pio2); + return trig.cosdf(if (sign) x + c4pio2 else x - c4pio2); } else { if (sign) { - return trig.__sindf(-x - c3pio2); + return trig.sindf(-x - c3pio2); } else { - return trig.__sindf(x - c3pio2); + return trig.sindf(x - c3pio2); } } } @@ -87,10 +87,10 @@ pub fn cosf(x: f32) callconv(.c) f32 { var y: f64 = undefined; const n = rem_pio2f(x, &y); return switch (n & 3) { - 0 => trig.__cosdf(y), - 1 => trig.__sindf(-y), - 2 => -trig.__cosdf(y), - else => trig.__sindf(y), + 0 => trig.cosdf(y), + 1 => trig.sindf(-y), + 2 => -trig.cosdf(y), + else => trig.sindf(y), }; } @@ -105,7 +105,7 @@ pub fn cos(x: f64) callconv(.c) f64 { if (compiler_rt.want_float_exceptions) mem.doNotOptimizeAway(x + 0x1p120); return 1.0; } - return trig.__cos(x, 0); + return trig.cos(x, 0); } // cos(Inf or NaN) is NaN @@ -116,10 +116,10 @@ pub fn cos(x: f64) callconv(.c) f64 { var y: [2]f64 = undefined; const n = rem_pio2(x, &y); return switch (n & 3) { - 0 => trig.__cos(y[0], y[1]), - 1 => -trig.__sin(y[0], y[1], 1), - 2 => -trig.__cos(y[0], y[1]), - else => trig.__sin(y[0], y[1], 1), + 0 => trig.cos(y[0], y[1]), + 1 => -trig.sin(y[0], y[1], 1), + 2 => -trig.cos(y[0], y[1]), + else => trig.sin(y[0], y[1], 1), }; } @@ -134,20 +134,20 @@ fn coslGeneric(comptime T: type, x: T) T { // raise inexact if x!=0 return 1.0 + x; } - return trig.__cosl(T, x, 0.0); + return trig.cosl(T, x, 0.0); } var y: [2]T = undefined; const n = rem_pio2l(T, x, &y); return switch (n & 3) { - 0 => trig.__cosl(T, y[0], y[1]), - 1 => -trig.__sinl(T, y[0], y[1], 1), - 2 => -trig.__cosl(T, y[0], y[1]), - else => trig.__sinl(T, y[0], y[1], 1), + 0 => trig.cosl(T, y[0], y[1]), + 1 => -trig.sinl(T, y[0], y[1], 1), + 2 => -trig.cosl(T, y[0], y[1]), + else => trig.sinl(T, y[0], y[1], 1), }; } -pub fn __cosx(x: f80) callconv(.c) f80 { +pub fn cosx(x: f80) callconv(.c) f80 { return coslGeneric(f80, x); } @@ -157,10 +157,10 @@ pub fn cosq(x: f128) callconv(.c) f128 { pub fn cosl(x: c_longdouble) callconv(.c) c_longdouble { switch (@typeInfo(c_longdouble).float.bits) { - 16 => return __cosh(x), + 16 => return cosh(x), 32 => return cosf(x), 64 => return cos(x), - 80 => return __cosx(x), + 80 => return cosx(x), 128 => return cosq(x), else => @compileError("unreachable"), } @@ -170,7 +170,7 @@ fn testCosSpecial(comptime T: type) !void { const f = switch (T) { f32 => cosf, f64 => cos, - f80 => __cosx, + f80 => cosx, f128 => cosq, else => @compileError("unimplemented"), }; @@ -214,13 +214,13 @@ test "cos64.special" { test "cos80.normal" { const epsilon = math.floatEps(f80); - try expectApproxEqAbs(@as(f80, 1.0), __cosx(0.0), epsilon); - try expectApproxEqAbs(@as(f80, 0.98006657784124163112419651674816888), __cosx(0.2), epsilon); - try expectApproxEqAbs(@as(f80, 0.62762309833608037003563995939286067), __cosx(0.8923), epsilon); - try expectApproxEqAbs(@as(f80, 0.070737201667702910088189851434268747), __cosx(1.5), epsilon); - try expectApproxEqAbs(@as(f80, 0.070737201667702910088189851434268747), __cosx(-1.5), epsilon); - try expectApproxEqAbs(@as(f80, 0.9691317730707771246), __cosx(37.45), epsilon); - try expectApproxEqAbs(@as(f80, 0.4008006809354834001), __cosx(89.123), epsilon); + try expectApproxEqAbs(@as(f80, 1.0), cosx(0.0), epsilon); + try expectApproxEqAbs(@as(f80, 0.98006657784124163112419651674816888), cosx(0.2), epsilon); + try expectApproxEqAbs(@as(f80, 0.62762309833608037003563995939286067), cosx(0.8923), epsilon); + try expectApproxEqAbs(@as(f80, 0.070737201667702910088189851434268747), cosx(1.5), epsilon); + try expectApproxEqAbs(@as(f80, 0.070737201667702910088189851434268747), cosx(-1.5), epsilon); + try expectApproxEqAbs(@as(f80, 0.9691317730707771246), cosx(37.45), epsilon); + try expectApproxEqAbs(@as(f80, 0.4008006809354834001), cosx(89.123), epsilon); } test "cos80.special" { diff --git a/lib/compiler_rt/sin.zig b/lib/compiler_rt/sin.zig index 0cff43997e06e31606cd0be4ffb66d3cb36203a2..84f970cff5432a672715597c8a60bed9fd5cd46d 100644 --- a/lib/compiler_rt/sin.zig +++ b/lib/compiler_rt/sin.zig @@ -20,11 +20,11 @@ const rem_pio2l = @import("rem_pio2l.zig").rem_pio2l; const ld = @import("long_double.zig"); comptime { - symbol(&__sinh, "__sinh"); + symbol(&sinh, "__sinh"); symbol(&sinl, "__sinl"); symbol(&sinf, "sinf"); symbol(&sin, "sin"); - symbol(&__sinx, "__sinx"); + symbol(&sinx, "__sinx"); if (compiler_rt.want_ppc_abi) { symbol(&sinq, "sinf128"); } @@ -32,7 +32,7 @@ comptime { symbol(&sinl, "sinl"); } -pub fn __sinh(x: f16) callconv(.c) f16 { +pub fn sinh(x: f16) callconv(.c) f16 { // TODO: more efficient implementation return @floatCast(sinf(x)); } @@ -60,27 +60,27 @@ pub fn sinf(x: f32) callconv(.c) f32 { } return x; } - return trig.__sindf(x); + return trig.sindf(x); } if (ix <= 0x407b53d1) { // |x| ~<= 5*pi/4 if (ix <= 0x4016cbe3) { // |x| ~<= 3pi/4 if (sign) { - return -trig.__cosdf(x + s1pio2); + return -trig.cosdf(x + s1pio2); } else { - return trig.__cosdf(x - s1pio2); + return trig.cosdf(x - s1pio2); } } - return trig.__sindf(if (sign) -(x + s2pio2) else -(x - s2pio2)); + return trig.sindf(if (sign) -(x + s2pio2) else -(x - s2pio2)); } if (ix <= 0x40e231d5) { // |x| ~<= 9*pi/4 if (ix <= 0x40afeddf) { // |x| ~<= 7*pi/4 if (sign) { - return trig.__cosdf(x + s3pio2); + return trig.cosdf(x + s3pio2); } else { - return -trig.__cosdf(x - s3pio2); + return -trig.cosdf(x - s3pio2); } } - return trig.__sindf(if (sign) x + s4pio2 else x - s4pio2); + return trig.sindf(if (sign) x + s4pio2 else x - s4pio2); } // sin(Inf or NaN) is NaN @@ -91,10 +91,10 @@ pub fn sinf(x: f32) callconv(.c) f32 { var y: f64 = undefined; const n = rem_pio2f(x, &y); return switch (n & 3) { - 0 => trig.__sindf(y), - 1 => trig.__cosdf(y), - 2 => trig.__sindf(-y), - else => -trig.__cosdf(y), + 0 => trig.sindf(y), + 1 => trig.cosdf(y), + 2 => trig.sindf(-y), + else => -trig.cosdf(y), }; } @@ -115,7 +115,7 @@ pub fn sin(x: f64) callconv(.c) f64 { } return x; } - return trig.__sin(x, 0.0, 0); + return trig.sin(x, 0.0, 0); } // sin(Inf or NaN) is NaN @@ -126,10 +126,10 @@ pub fn sin(x: f64) callconv(.c) f64 { var y: [2]f64 = undefined; const n = rem_pio2(x, &y); return switch (n & 3) { - 0 => trig.__sin(y[0], y[1], 1), - 1 => trig.__cos(y[0], y[1]), - 2 => -trig.__sin(y[0], y[1], 1), - else => -trig.__cos(y[0], y[1]), + 0 => trig.sin(y[0], y[1], 1), + 1 => trig.cos(y[0], y[1]), + 2 => -trig.sin(y[0], y[1], 1), + else => -trig.cos(y[0], y[1]), }; } @@ -147,20 +147,20 @@ fn sinlGeneric(comptime T: type, x: T) T { } return x; } - return trig.__sinl(T, x, 0.0, 0); + return trig.sinl(T, x, 0.0, 0); } var y: [2]T = undefined; const n = rem_pio2l(T, x, &y); return switch (n & 3) { - 0 => trig.__sinl(T, y[0], y[1], 1), - 1 => trig.__cosl(T, y[0], y[1]), - 2 => -trig.__sinl(T, y[0], y[1], 1), - else => -trig.__cosl(T, y[0], y[1]), + 0 => trig.sinl(T, y[0], y[1], 1), + 1 => trig.cosl(T, y[0], y[1]), + 2 => -trig.sinl(T, y[0], y[1], 1), + else => -trig.cosl(T, y[0], y[1]), }; } -pub fn __sinx(x: f80) callconv(.c) f80 { +pub fn sinx(x: f80) callconv(.c) f80 { return sinlGeneric(f80, x); } @@ -170,10 +170,10 @@ pub fn sinq(x: f128) callconv(.c) f128 { pub fn sinl(x: c_longdouble) callconv(.c) c_longdouble { switch (@typeInfo(c_longdouble).float.bits) { - 16 => return __sinh(x), + 16 => return sinh(x), 32 => return sinf(x), 64 => return sin(x), - 80 => return __sinx(x), + 80 => return sinx(x), 128 => return sinq(x), else => @compileError("unreachable"), } @@ -183,7 +183,7 @@ fn testSinSpecial(comptime T: type) !void { const f = switch (T) { f32 => sinf, f64 => sin, - f80 => __sinx, + f80 => sinx, f128 => sinq, else => @compileError("unimplemented"), }; @@ -227,13 +227,13 @@ test "sin64.special" { test "sin80.normal" { const epsilon = math.floatEps(f80); - try expectApproxEqAbs(@as(f80, 0.0), __sinx(0.0), epsilon); - try expectApproxEqAbs(@as(f80, 0.19866933079506121545941262711838975), __sinx(0.2), epsilon); - try expectApproxEqAbs(@as(f80, 0.77851733855773487830689285621486050), __sinx(0.8923), epsilon); - try expectApproxEqAbs(@as(f80, 0.99749498660405443094172337114148732), __sinx(1.5), epsilon); - try expectApproxEqAbs(@as(f80, -0.99749498660405443094172337114148732), __sinx(-1.5), epsilon); - try expectApproxEqAbs(@as(f80, -0.24654331551411356504), __sinx(37.45), epsilon); - try expectApproxEqAbs(@as(f80, 0.91616527666226951006), __sinx(89.123), epsilon); + try expectApproxEqAbs(@as(f80, 0.0), sinx(0.0), epsilon); + try expectApproxEqAbs(@as(f80, 0.19866933079506121545941262711838975), sinx(0.2), epsilon); + try expectApproxEqAbs(@as(f80, 0.77851733855773487830689285621486050), sinx(0.8923), epsilon); + try expectApproxEqAbs(@as(f80, 0.99749498660405443094172337114148732), sinx(1.5), epsilon); + try expectApproxEqAbs(@as(f80, -0.99749498660405443094172337114148732), sinx(-1.5), epsilon); + try expectApproxEqAbs(@as(f80, -0.24654331551411356504), sinx(37.45), epsilon); + try expectApproxEqAbs(@as(f80, 0.91616527666226951006), sinx(89.123), epsilon); } test "sin80.special" { diff --git a/lib/compiler_rt/sincos.zig b/lib/compiler_rt/sincos.zig index 1616cc5556459503d5ac6fd9eeab045a9ba52c07..3f67627405dc3f3a8b2bd5b789eec15a253cb91a 100644 --- a/lib/compiler_rt/sincos.zig +++ b/lib/compiler_rt/sincos.zig @@ -14,10 +14,10 @@ const compiler_rt = @import("../compiler_rt.zig"); const symbol = compiler_rt.symbol; comptime { - symbol(&__sincosh, "__sincosh"); + symbol(&sincosh, "__sincosh"); symbol(&sincosf, "sincosf"); symbol(&sincos, "sincos"); - symbol(&__sincosx, "__sincosx"); + symbol(&sincosx, "__sincosx"); if (compiler_rt.want_ppc_abi) { symbol(&sincosq, "sincosf128"); } @@ -25,7 +25,7 @@ comptime { symbol(&sincosl, "sincosl"); } -pub fn __sincosh(x: f16, r_sin: *f16, r_cos: *f16) callconv(.c) void { +pub fn sincosh(x: f16, r_sin: *f16, r_cos: *f16) callconv(.c) void { // TODO: more efficient implementation var big_sin: f32 = undefined; var big_cos: f32 = undefined; @@ -60,8 +60,8 @@ pub fn sincosf(x: f32, r_sin: *f32, r_cos: *f32) callconv(.c) void { r_cos.* = 1.0; return; } - r_sin.* = trig.__sindf(x); - r_cos.* = trig.__cosdf(x); + r_sin.* = trig.sindf(x); + r_cos.* = trig.cosdf(x); return; } @@ -70,17 +70,17 @@ pub fn sincosf(x: f32, r_sin: *f32, r_cos: *f32) callconv(.c) void { // |x| ~<= 3pi/4 if (ix <= 0x4016cbe3) { if (sign) { - r_sin.* = -trig.__cosdf(x + sc1pio2); - r_cos.* = trig.__sindf(x + sc1pio2); + r_sin.* = -trig.cosdf(x + sc1pio2); + r_cos.* = trig.sindf(x + sc1pio2); } else { - r_sin.* = trig.__cosdf(sc1pio2 - x); - r_cos.* = trig.__sindf(sc1pio2 - x); + r_sin.* = trig.cosdf(sc1pio2 - x); + r_cos.* = trig.sindf(sc1pio2 - x); } return; } // -sin(x+c) is not correct if x+c could be 0: -0 vs +0 - r_sin.* = -trig.__sindf(if (sign) x + sc2pio2 else x - sc2pio2); - r_cos.* = -trig.__cosdf(if (sign) x + sc2pio2 else x - sc2pio2); + r_sin.* = -trig.sindf(if (sign) x + sc2pio2 else x - sc2pio2); + r_cos.* = -trig.cosdf(if (sign) x + sc2pio2 else x - sc2pio2); return; } @@ -89,16 +89,16 @@ pub fn sincosf(x: f32, r_sin: *f32, r_cos: *f32) callconv(.c) void { // |x| ~<= 7*pi/4 if (ix <= 0x40afeddf) { if (sign) { - r_sin.* = trig.__cosdf(x + sc3pio2); - r_cos.* = -trig.__sindf(x + sc3pio2); + r_sin.* = trig.cosdf(x + sc3pio2); + r_cos.* = -trig.sindf(x + sc3pio2); } else { - r_sin.* = -trig.__cosdf(x - sc3pio2); - r_cos.* = trig.__sindf(x - sc3pio2); + r_sin.* = -trig.cosdf(x - sc3pio2); + r_cos.* = trig.sindf(x - sc3pio2); } return; } - r_sin.* = trig.__sindf(if (sign) x + sc4pio2 else x - sc4pio2); - r_cos.* = trig.__cosdf(if (sign) x + sc4pio2 else x - sc4pio2); + r_sin.* = trig.sindf(if (sign) x + sc4pio2 else x - sc4pio2); + r_cos.* = trig.cosdf(if (sign) x + sc4pio2 else x - sc4pio2); return; } @@ -113,8 +113,8 @@ pub fn sincosf(x: f32, r_sin: *f32, r_cos: *f32) callconv(.c) void { // general argument reduction needed var y: f64 = undefined; const n = rem_pio2f(x, &y); - const s = trig.__sindf(y); - const c = trig.__cosdf(y); + const s = trig.sindf(y); + const c = trig.cosdf(y); switch (n & 3) { 0 => { r_sin.* = s; @@ -154,8 +154,8 @@ pub fn sincos(x: f64, r_sin: *f64, r_cos: *f64) callconv(.c) void { r_cos.* = 1.0; return; } - r_sin.* = trig.__sin(x, 0.0, 0); - r_cos.* = trig.__cos(x, 0.0); + r_sin.* = trig.sin(x, 0.0, 0); + r_cos.* = trig.cos(x, 0.0); return; } @@ -170,8 +170,8 @@ pub fn sincos(x: f64, r_sin: *f64, r_cos: *f64) callconv(.c) void { // argument reduction needed var y: [2]f64 = undefined; const n = rem_pio2(x, &y); - const s = trig.__sin(y[0], y[1], 1); - const c = trig.__cos(y[0], y[1]); + const s = trig.sin(y[0], y[1], 1); + const c = trig.cos(y[0], y[1]); switch (n & 3) { 0 => { r_sin.* = s; @@ -216,15 +216,15 @@ fn sincoslGeneric(comptime T: type, x: T, r_sin: *T, r_cos: *T) void { r_cos.* = 1.0 + x; return; } - r_sin.* = trig.__sinl(T, x, 0.0, 0); - r_cos.* = trig.__cosl(T, x, 0.0); + 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]); + 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; @@ -245,7 +245,7 @@ fn sincoslGeneric(comptime T: type, x: T, r_sin: *T, r_cos: *T) void { } } -pub fn __sincosx(x: f80, r_sin: *f80, r_cos: *f80) callconv(.c) void { +pub fn sincosx(x: f80, r_sin: *f80, r_cos: *f80) callconv(.c) void { return sincoslGeneric(f80, x, r_sin, r_cos); } @@ -255,10 +255,10 @@ pub fn sincosq(x: f128, r_sin: *f128, r_cos: *f128) callconv(.c) void { 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), + 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), + 80 => return sincosx(x, r_sin, r_cos), 128 => return sincosq(x, r_sin, r_cos), else => @compileError("unreachable"), } @@ -268,7 +268,7 @@ fn testSincosSpecial(comptime T: type) !void { const f = switch (T) { f32 => sincosf, f64 => sincos, - f80 => __sincosx, + f80 => sincosx, f128 => sincosq, else => @compileError("unimplemented"), }; @@ -378,31 +378,31 @@ test "sincos80.normal" { var s: f80 = undefined; var c: f80 = undefined; - __sincosx(0.0, &s, &c); + 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); + 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); + 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); + 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); + 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); + 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); + sincosx(89.123, &s, &c); try expectApproxEqAbs(@as(f80, 0.91616527666226951006), s, epsilon); try expectApproxEqAbs(@as(f80, 0.4008006809354834001), c, epsilon); } diff --git a/lib/compiler_rt/tan.zig b/lib/compiler_rt/tan.zig index 1f925e3e8c6c6abac73d65d33f39ff2680377b03..5d5ba3196d52ee1f6f5161a8e0e36a587626cb30 100644 --- a/lib/compiler_rt/tan.zig +++ b/lib/compiler_rt/tan.zig @@ -24,10 +24,10 @@ const compiler_rt = @import("../compiler_rt.zig"); const symbol = @import("../compiler_rt.zig").symbol; comptime { - symbol(&__tanh, "__tanh"); + symbol(&tanh, "__tanh"); symbol(&tanf, "tanf"); symbol(&tan, "tan"); - symbol(&__tanx, "__tanx"); + symbol(&tanx, "__tanx"); if (compiler_rt.want_ppc_abi) { symbol(&tanq, "tanf128"); } @@ -35,7 +35,7 @@ comptime { symbol(&tanl, "tanl"); } -pub fn __tanh(x: f16) callconv(.c) f16 { +pub fn tanh(x: f16) callconv(.c) f16 { // TODO: more efficient implementation return @floatCast(tanf(x)); } @@ -63,20 +63,20 @@ pub fn tanf(x: f32) callconv(.c) f32 { } return x; } - return kernel.__tandf(x, false); + return kernel.tandf(x, false); } if (ix <= 0x407b53d1) { // |x| ~<= 5*pi/4 if (ix <= 0x4016cbe3) { // |x| ~<= 3pi/4 - return kernel.__tandf((if (sign) x + t1pio2 else x - t1pio2), true); + return kernel.tandf((if (sign) x + t1pio2 else x - t1pio2), true); } else { - return kernel.__tandf((if (sign) x + t2pio2 else x - t2pio2), false); + return kernel.tandf((if (sign) x + t2pio2 else x - t2pio2), false); } } if (ix <= 0x40e231d5) { // |x| ~<= 9*pi/4 if (ix <= 0x40afeddf) { // |x| ~<= 7*pi/4 - return kernel.__tandf((if (sign) x + t3pio2 else x - t3pio2), true); + return kernel.tandf((if (sign) x + t3pio2 else x - t3pio2), true); } else { - return kernel.__tandf((if (sign) x + t4pio2 else x - t4pio2), false); + return kernel.tandf((if (sign) x + t4pio2 else x - t4pio2), false); } } @@ -87,7 +87,7 @@ pub fn tanf(x: f32) callconv(.c) f32 { var y: f64 = undefined; const n = rem_pio2f(x, &y); - return kernel.__tandf(y, n & 1 != 0); + return kernel.tandf(y, n & 1 != 0); } pub fn tan(x: f64) callconv(.c) f64 { @@ -107,7 +107,7 @@ pub fn tan(x: f64) callconv(.c) f64 { } return x; } - return kernel.__tan(x, 0.0, false); + return kernel.tan(x, 0.0, false); } // tan(Inf or NaN) is NaN @@ -117,7 +117,7 @@ pub fn tan(x: f64) callconv(.c) f64 { var y: [2]f64 = undefined; const n = rem_pio2(x, &y); - return kernel.__tan(y[0], y[1], n & 1 != 0); + return kernel.tan(y[0], y[1], n & 1 != 0); } fn tanlGeneric(comptime T: type, x: T) T { @@ -137,15 +137,15 @@ fn tanlGeneric(comptime T: type, x: T) T { } return x; } - return kernel.__tanl(T, x, 0.0, 0); + return kernel.tanl(T, x, 0.0, 0); } var y: [2]T = undefined; const n = rem_pio2l(T, x, &y); - return kernel.__tanl(T, y[0], y[1], n & 1); + return kernel.tanl(T, y[0], y[1], n & 1); } -pub fn __tanx(x: f80) callconv(.c) f80 { +pub fn tanx(x: f80) callconv(.c) f80 { return tanlGeneric(f80, x); } @@ -155,10 +155,10 @@ pub fn tanq(x: f128) callconv(.c) f128 { pub fn tanl(x: c_longdouble) callconv(.c) c_longdouble { switch (@typeInfo(c_longdouble).float.bits) { - 16 => return __tanh(x), + 16 => return tanh(x), 32 => return tanf(x), 64 => return tan(x), - 80 => return __tanx(x), + 80 => return tanx(x), 128 => return tanq(x), else => @compileError("unreachable"), } @@ -184,7 +184,7 @@ fn testTanSpecial(comptime T: type) !void { const f = switch (T) { f32 => tanf, f64 => tan, - f80 => __tanx, + f80 => tanx, f128 => tanq, else => @compileError("unimplemented"), }; @@ -207,12 +207,12 @@ test "tan64.normal" { test "tan80.normal" { const epsilon = math.floatEps(f80); - try expectApproxEqAbs(@as(f80, 0.0), __tanx(0.0), epsilon); - try expectApproxEqAbs(@as(f80, 0.2027100355086724833213582716475345), __tanx(0.2), epsilon); - try expectApproxEqAbs(@as(f80, 1.2404217445497097995561220131857544), __tanx(0.8923), epsilon); - try expectApproxEqAbs(@as(f80, 14.10141994717171938764), __tanx(1.5), epsilon); - try expectApproxEqAbs(@as(f80, -0.25439607116885656232), __tanx(37.45), epsilon); - try expectApproxEqAbs(@as(f80, 2.2858376251355320963), __tanx(89.123), epsilon); + try expectApproxEqAbs(@as(f80, 0.0), tanx(0.0), epsilon); + try expectApproxEqAbs(@as(f80, 0.2027100355086724833213582716475345), tanx(0.2), epsilon); + try expectApproxEqAbs(@as(f80, 1.2404217445497097995561220131857544), tanx(0.8923), epsilon); + try expectApproxEqAbs(@as(f80, 14.10141994717171938764), tanx(1.5), epsilon); + try expectApproxEqAbs(@as(f80, -0.25439607116885656232), tanx(37.45), epsilon); + try expectApproxEqAbs(@as(f80, 2.2858376251355320963), tanx(89.123), epsilon); } test "tan128.normal" { diff --git a/lib/compiler_rt/trig.zig b/lib/compiler_rt/trig.zig index e68089baeab1b6a9ff4cb818e81db5ca7e14952a..7c106100ec3cecffb6c3b1a2c9423d84855c04b9 100644 --- a/lib/compiler_rt/trig.zig +++ b/lib/compiler_rt/trig.zig @@ -11,8 +11,9 @@ // https://git.musl-libc.org/cgit/musl/tree/src/math/__cosl.c // https://git.musl-libc.org/cgit/musl/tree/src/math/__tanl.c -/// pi divided by 4 -pub const pi_4 = 0.78539816339744830962; +const std = @import("std"); + +pub const pi_4 = std.math.pi / 4.0; /// kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164 /// Input x is assumed to be bounded by ~pi/4 in magnitude. @@ -49,7 +50,7 @@ pub const pi_4 = 0.78539816339744830962; /// expression for cos(). Retention happens in all cases tested /// under FreeBSD, so don't pessimize things by forcibly clipping /// any extra precision in w. -pub fn __cos(x: f64, y: f64) f64 { +pub fn cos(x: f64, y: f64) f64 { const C1 = 4.16666666666666019037e-02; // 0x3FA55555, 0x5555554C const C2 = -1.38888888888741095749e-03; // 0xBF56C16C, 0x16C15177 const C3 = 2.48015872894767294178e-05; // 0x3EFA01A0, 0x19CB1590 @@ -65,7 +66,7 @@ pub fn __cos(x: f64, y: f64) f64 { return w + (((1.0 - w) - hz) + (z * r - x * y)); } -pub fn __cosdf(x: f64) f32 { +pub fn cosdf(x: f64) f32 { // |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). const C0 = -0x1ffffffd0c5e81.0p-54; // -0.499999997251031003120 const C1 = 0x155553e1053a42.0p-57; // 0.0416666233237390631894 @@ -79,7 +80,7 @@ pub fn __cosdf(x: f64) f32 { return @floatCast(((1.0 + z * C0) + w * C1) + (w * z) * r); } -pub fn __cosl(comptime T: type, x: T, y: T) T { +pub fn cosl(comptime T: type, x: T, y: T) T { const impl = switch (T) { f80 => struct { const C1: T = 0.0416666666666666666136; @@ -115,7 +116,7 @@ pub fn __cosl(comptime T: type, x: T, y: T) T { z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * C11)))))))))); } }, - else => @compileError("__cosl supports only f80 and f128, got: " ++ @typeName(T)), + else => @compileError("cosl supports only f80 and f128, got: " ++ @typeName(T)), }; const z = x * x; @@ -152,7 +153,7 @@ pub fn __cosl(comptime T: type, x: T, y: T) T { /// r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6)))) /// then 3 2 /// sin(x) = x + (S1*x + (x *(r-y/2)+y)) -pub fn __sin(x: f64, y: f64, iy: i32) f64 { +pub fn sin(x: f64, y: f64, iy: i32) f64 { const S1 = -1.66666666666666324348e-01; // 0xBFC55555, 0x55555549 const S2 = 8.33333333332248946124e-03; // 0x3F811111, 0x1110F8A6 const S3 = -1.98412698298579493134e-04; // 0xBF2A01A0, 0x19C161D5 @@ -171,7 +172,7 @@ pub fn __sin(x: f64, y: f64, iy: i32) f64 { } } -pub fn __sinl(comptime T: type, x: T, y: T, iy: i32) T { +pub fn sinl(comptime T: type, x: T, y: T, iy: i32) T { const impl = switch (T) { f80 => struct { const S1: T = -0.166666666666666666671; @@ -209,7 +210,7 @@ pub fn __sinl(comptime T: type, x: T, y: T, iy: i32) T { z * (S9 + z * (S10 + z * (S11 + z * S12))))))))); } }, - else => @compileError("__sinl supports only f80 and f128, got: " ++ @typeName(T)), + else => @compileError("sinl supports only f80 and f128, got: " ++ @typeName(T)), }; const z = x * x; @@ -223,7 +224,7 @@ pub fn __sinl(comptime T: type, x: T, y: T, iy: i32) T { return x - ((z * (0.5 * y - v * r) - y) - v * impl.S1); } -pub fn __sindf(x: f64) f32 { +pub fn sindf(x: f64) f32 { // |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]). const S1 = -0x15555554cbac77.0p-55; // -0.166666666416265235595 const S2 = 0x111110896efbb2.0p-59; // 0.0083333293858894631756 @@ -270,7 +271,7 @@ pub fn __sindf(x: f64) f32 { /// 4. For x in [0.67434,pi/4], let y = pi/4 - x, then /// tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y)) /// = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y))) -pub fn __tan(x_: f64, y_: f64, odd: bool) f64 { +pub fn tan(x_: f64, y_: f64, odd: bool) f64 { var x = x_; var y = y_; @@ -343,7 +344,7 @@ pub fn __tan(x_: f64, y_: f64, odd: bool) f64 { return a0 + a * (1.0 + a0 * w0 + a0 * v); } -pub fn __tandf(x: f64, odd: bool) f32 { +pub fn tandf(x: f64, odd: bool) f32 { // |tan(x)/x - t(x)| < 2**-25.5 (~[-2e-08, 2e-08]). const T = [_]f64{ 0x15554d3418c99f.0p-54, // 0.333331395030791399758 @@ -376,7 +377,7 @@ pub fn __tandf(x: f64, odd: bool) f32 { return @floatCast(if (odd) -1.0 / r0 else r0); } -pub fn __tanl(comptime T: type, x_: T, y_: T, odd: i32) T { +pub fn tanl(comptime T: type, x_: T, y_: T, odd: i32) T { var x = x_; var y = y_; const impl = switch (T) { @@ -456,7 +457,7 @@ pub fn __tanl(comptime T: type, x_: T, y_: T, odd: i32) T { w * (T47 + w * (T51 + w * T55))))))))))); } }, - else => @compileError("__tanl supports only f80 and f128, got: " ++ @typeName(T)), + else => @compileError("tanl supports only f80 and f128, got: " ++ @typeName(T)), }; const big = @abs(x) >= 0.67434; -- 2.54.0 From 4ccac1de416bff8846352ef3bf2fb592c5a419a9 Mon Sep 17 00:00:00 2001 From: mihael Date: Thu, 2 Apr 2026 23:53:13 +0200 Subject: [PATCH 20/20] `compiler_rt`: Make `long double` trig impls non-generic --- lib/compiler_rt/cos.zig | 68 ++++--- lib/compiler_rt/sin.zig | 47 +++-- lib/compiler_rt/sincos.zig | 147 +++++++++------ lib/compiler_rt/tan.zig | 63 ++++--- lib/compiler_rt/trig.zig | 365 +++++++++++++++++++------------------ 5 files changed, 395 insertions(+), 295 deletions(-) diff --git a/lib/compiler_rt/cos.zig b/lib/compiler_rt/cos.zig index aa9282625fba5a3b95aa0b55cb7fc23a6df090ed..3620e31835b20461f80a096c4ccc7a59365fd812 100644 --- a/lib/compiler_rt/cos.zig +++ b/lib/compiler_rt/cos.zig @@ -123,36 +123,52 @@ pub fn cos(x: f64) callconv(.c) f64 { }; } -fn coslGeneric(comptime T: type, x: T) T { - const se = ld.signExponent(x) & 0x7fff; - if (se == 0x7fff) { - return x - x; - } - - if (@abs(x) < trig.pi_4) { - if (se < 0x3fff - math.floatMantissaBits(T)) { - // raise inexact if x!=0 - return 1.0 + x; - } - return trig.cosl(T, x, 0.0); - } - - var y: [2]T = undefined; - const n = rem_pio2l(T, x, &y); - return switch (n & 3) { - 0 => trig.cosl(T, y[0], y[1]), - 1 => -trig.sinl(T, y[0], y[1], 1), - 2 => -trig.cosl(T, y[0], y[1]), - else => trig.sinl(T, y[0], y[1], 1), - }; -} - pub fn cosx(x: f80) callconv(.c) f80 { - return coslGeneric(f80, x); + const se = ld.signExponent(x) & 0x7fff; + if (se == 0x7fff) { + return x - x; + } + + if (@abs(x) < trig.pi_4) { + if (se < 0x3fff - math.floatMantissaBits(f80)) { + // raise inexact if x!=0 + return 1.0 + x; + } + return trig.cosx(x, 0.0); + } + + var y: [2]f80 = undefined; + const n = rem_pio2l(f80, x, &y); + return switch (n & 3) { + 0 => trig.cosx(y[0], y[1]), + 1 => -trig.sinx(y[0], y[1], 1), + 2 => -trig.cosx(y[0], y[1]), + else => trig.sinx(y[0], y[1], 1), + }; } pub fn cosq(x: f128) callconv(.c) f128 { - return coslGeneric(f128, x); + const se = ld.signExponent(x) & 0x7fff; + if (se == 0x7fff) { + return x - x; + } + + if (@abs(x) < trig.pi_4) { + if (se < 0x3fff - math.floatMantissaBits(f128)) { + // raise inexact if x!=0 + return 1.0 + x; + } + return trig.cosq(x, 0.0); + } + + var y: [2]f128 = undefined; + const n = rem_pio2l(f128, x, &y); + return switch (n & 3) { + 0 => trig.cosq(y[0], y[1]), + 1 => -trig.sinq(y[0], y[1], 1), + 2 => -trig.cosq(y[0], y[1]), + else => trig.sinq(y[0], y[1], 1), + }; } pub fn cosl(x: c_longdouble) callconv(.c) c_longdouble { diff --git a/lib/compiler_rt/sin.zig b/lib/compiler_rt/sin.zig index 84f970cff5432a672715597c8a60bed9fd5cd46d..3f1509383bacbef90c8f81debc5a6eaeedbfdf7d 100644 --- a/lib/compiler_rt/sin.zig +++ b/lib/compiler_rt/sin.zig @@ -133,39 +133,58 @@ pub fn sin(x: f64) callconv(.c) f64 { }; } -fn sinlGeneric(comptime T: type, x: T) T { +fn sinx(x: f80) callconv(.c) f80 { const se = ld.signExponent(x) & 0x7fff; if (se == 0x7fff) { return x - x; } if (@abs(x) < trig.pi_4) { - if (se < 0x3fff - (math.floatMantissaBits(T) / 2)) { + if (se < 0x3fff - (math.floatMantissaBits(f80) / 2)) { // raise inexact if x!=0 and underflow if subnormal if (compiler_rt.want_float_exceptions) { mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120); } return x; } - return trig.sinl(T, x, 0.0, 0); + return trig.sinx(x, 0.0, 0); } - var y: [2]T = undefined; - const n = rem_pio2l(T, x, &y); + var y: [2]f80 = undefined; + const n = rem_pio2l(f80, x, &y); return switch (n & 3) { - 0 => trig.sinl(T, y[0], y[1], 1), - 1 => trig.cosl(T, y[0], y[1]), - 2 => -trig.sinl(T, y[0], y[1], 1), - else => -trig.cosl(T, y[0], y[1]), + 0 => trig.sinx(y[0], y[1], 1), + 1 => trig.cosx(y[0], y[1]), + 2 => -trig.sinx(y[0], y[1], 1), + else => -trig.cosx(y[0], y[1]), }; } -pub fn sinx(x: f80) callconv(.c) f80 { - return sinlGeneric(f80, x); -} - pub fn sinq(x: f128) callconv(.c) f128 { - return sinlGeneric(f128, x); + const se = ld.signExponent(x) & 0x7fff; + if (se == 0x7fff) { + return x - x; + } + + if (@abs(x) < trig.pi_4) { + if (se < 0x3fff - (math.floatMantissaBits(f128) / 2)) { + // raise inexact if x!=0 and underflow if subnormal + if (compiler_rt.want_float_exceptions) { + mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120); + } + return x; + } + return trig.sinq(x, 0.0, 0); + } + + var y: [2]f128 = undefined; + const n = rem_pio2l(f128, x, &y); + return switch (n & 3) { + 0 => trig.sinq(y[0], y[1], 1), + 1 => trig.cosq(y[0], y[1]), + 2 => -trig.sinq(y[0], y[1], 1), + else => -trig.cosq(y[0], y[1]), + }; } pub fn sinl(x: c_longdouble) callconv(.c) c_longdouble { diff --git a/lib/compiler_rt/sincos.zig b/lib/compiler_rt/sincos.zig index 3f67627405dc3f3a8b2bd5b789eec15a253cb91a..7208325c2b86f65dd6736437572eea3c09235ff3 100644 --- a/lib/compiler_rt/sincos.zig +++ b/lib/compiler_rt/sincos.zig @@ -192,65 +192,102 @@ 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 = ld.signExponent(x) & 0x7fff; - if (se == 0x7fff) { - const result = x - x; - r_sin.* = result; - r_cos.* = result; - return; - } - - if (@abs(x) < trig.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 { - return sincoslGeneric(f80, x, r_sin, r_cos); + const se = ld.signExponent(x) & 0x7fff; + if (se == 0x7fff) { + const result = x - x; + r_sin.* = result; + r_cos.* = result; + return; + } + + if (@abs(x) < trig.pi_4) { + if (se < 0x3fff - math.floatMantissaBits(f80)) { + // 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.sinx(x, 0.0, 0); + r_cos.* = trig.cosx(x, 0.0); + return; + } + + var y: [2]f80 = undefined; + const n = rem_pio2l(f80, x, &y); + const s = trig.sinx(y[0], y[1], 1); + const c = trig.cosx(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 sincosq(x: f128, r_sin: *f128, r_cos: *f128) callconv(.c) void { - return sincoslGeneric(f128, x, r_sin, r_cos); + const se = ld.signExponent(x) & 0x7fff; + if (se == 0x7fff) { + const result = x - x; + r_sin.* = result; + r_cos.* = result; + return; + } + + if (@abs(x) < trig.pi_4) { + if (se < 0x3fff - math.floatMantissaBits(f128)) { + // 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.sinq(x, 0.0, 0); + r_cos.* = trig.cosq(x, 0.0); + return; + } + + var y: [2]f128 = undefined; + const n = rem_pio2l(f128, x, &y); + const s = trig.sinq(y[0], y[1], 1); + const c = trig.cosq(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 sincosl(x: c_longdouble, r_sin: *c_longdouble, r_cos: *c_longdouble) callconv(.c) void { diff --git a/lib/compiler_rt/tan.zig b/lib/compiler_rt/tan.zig index 5d5ba3196d52ee1f6f5161a8e0e36a587626cb30..378491aad6d25de04a4cbeb15747eb20e118b25b 100644 --- a/lib/compiler_rt/tan.zig +++ b/lib/compiler_rt/tan.zig @@ -120,37 +120,46 @@ pub fn tan(x: f64) callconv(.c) f64 { return kernel.tan(y[0], y[1], n & 1 != 0); } -fn tanlGeneric(comptime T: type, x: T) T { - if (!(T == f80 or T == f128)) { - @compileError("`tanlGeneric` implemented only for `f80` and `f128`, got: " ++ T); - } - - const se = ld.signExponent(x) & 0x7fff; - if (se == 0x7fff) { - return x - x; - } - - if (@abs(x) < kernel.pi_4) { - if (se < 0x3fff - math.floatMantissaBits(T) / 2) { - if (compiler_rt.want_float_exceptions) { - mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120); - } - return x; - } - return kernel.tanl(T, x, 0.0, 0); - } - - var y: [2]T = undefined; - const n = rem_pio2l(T, x, &y); - return kernel.tanl(T, y[0], y[1], n & 1); -} - pub fn tanx(x: f80) callconv(.c) f80 { - return tanlGeneric(f80, x); + const se = ld.signExponent(x) & 0x7fff; + if (se == 0x7fff) { + return x - x; + } + + if (@abs(x) < kernel.pi_4) { + if (se < 0x3fff - math.floatMantissaBits(f80) / 2) { + if (compiler_rt.want_float_exceptions) { + mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120); + } + return x; + } + return kernel.tanx(x, 0.0, 0); + } + + var y: [2]f80 = undefined; + const n = rem_pio2l(f80, x, &y); + return kernel.tanx(y[0], y[1], n & 1); } pub fn tanq(x: f128) callconv(.c) f128 { - return tanlGeneric(f128, x); + const se = ld.signExponent(x) & 0x7fff; + if (se == 0x7fff) { + return x - x; + } + + if (@abs(x) < kernel.pi_4) { + if (se < 0x3fff - math.floatMantissaBits(f128) / 2) { + if (compiler_rt.want_float_exceptions) { + mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120); + } + return x; + } + return kernel.tanq(x, 0.0, 0); + } + + var y: [2]f128 = undefined; + const n = rem_pio2l(f128, x, &y); + return kernel.tanq(y[0], y[1], n & 1); } pub fn tanl(x: c_longdouble) callconv(.c) c_longdouble { diff --git a/lib/compiler_rt/trig.zig b/lib/compiler_rt/trig.zig index 7c106100ec3cecffb6c3b1a2c9423d84855c04b9..e357613948625ca83779ddaa51a6754cc5870218 100644 --- a/lib/compiler_rt/trig.zig +++ b/lib/compiler_rt/trig.zig @@ -80,49 +80,43 @@ pub fn cosdf(x: f64) f32 { return @floatCast(((1.0 + z * C0) + w * C1) + (w * z) * r); } -pub fn cosl(comptime T: type, x: T, y: T) T { - const impl = switch (T) { - f80 => struct { - const C1: T = 0.0416666666666666666136; +pub fn cosx(x: f80, y: f80) f80 { + const C1: f80 = 0.0416666666666666666136; + const C2: f64 = -0.0013888888888888874; + const C3: f64 = 0.000024801587301571716; + const C4: f64 = -0.00000027557319215507120; + const C5: f64 = 0.0000000020876754400407278; + const C6: f64 = -1.1470297442401303e-11; + const C7: f64 = 4.7383039476436467e-14; - const C2: f64 = -0.0013888888888888874; - const C3: f64 = 0.000024801587301571716; - const C4: f64 = -0.00000027557319215507120; - const C5: f64 = 0.0000000020876754400407278; - const C6: f64 = -1.1470297442401303e-11; - const C7: f64 = 4.7383039476436467e-14; - - inline fn poly(z: T) T { - return z * (C1 + z * (C2 + z * (C3 + z * (C4 + - z * (C5 + z * (C6 + z * C7)))))); - } - }, - f128 => struct { - const C1: T = 0.04166666666666666666666666666666658424671; - const C2: T = -0.001388888888888888888888888888863490893732; - const C3: T = 0.00002480158730158730158730158600795304914210; - const C4: T = -0.2755731922398589065255474947078934284324e-6; - const C5: T = 0.2087675698786809897659225313136400793948e-8; - const C6: T = -0.1147074559772972315817149986812031204775e-10; - const C7: T = 0.4779477332386808976875457937252120293400e-13; + const z = x * x; + const r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + + z * (C5 + z * (C6 + z * C7)))))); + const hz = 0.5 * z; + const w = 1.0 - hz; - const C8: f64 = -0.1561920696721507929516718307820958119868e-15; - const C9: f64 = 0.4110317413744594971475941557607804508039e-18; - const C10: f64 = -0.8896592467191938803288521958313920156409e-21; - const C11: f64 = 0.1601061435794535138244346256065192782581e-23; + return w + (((1.0 - w) - hz) + (z * r - x * y)); +} - inline fn poly(z: T) T { - return z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 + - z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * C11)))))))))); - } - }, - else => @compileError("cosl supports only f80 and f128, got: " ++ @typeName(T)), - }; +pub fn cosq(x: f128, y: f128) f128 { + const C1: f128 = 0.04166666666666666666666666666666658424671; + const C2: f128 = -0.001388888888888888888888888888863490893732; + const C3: f128 = 0.00002480158730158730158730158600795304914210; + const C4: f128 = -0.2755731922398589065255474947078934284324e-6; + const C5: f128 = 0.2087675698786809897659225313136400793948e-8; + const C6: f128 = -0.1147074559772972315817149986812031204775e-10; + const C7: f128 = 0.4779477332386808976875457937252120293400e-13; + const C8: f64 = -0.1561920696721507929516718307820958119868e-15; + const C9: f64 = 0.4110317413744594971475941557607804508039e-18; + const C10: f64 = -0.8896592467191938803288521958313920156409e-21; + const C11: f64 = 0.1601061435794535138244346256065192782581e-23; const z = x * x; - const r = impl.poly(z); + const r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 + + z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * C11)))))))))); const hz = 0.5 * z; const w = 1.0 - hz; + return w + (((1.0 - w) - hz) + (z * r - x * y)); } @@ -172,58 +166,6 @@ pub fn sin(x: f64, y: f64, iy: i32) f64 { } } -pub fn sinl(comptime T: type, x: T, y: T, iy: i32) T { - const impl = switch (T) { - f80 => struct { - const S1: T = -0.166666666666666666671; - - const S2: f64 = 0.0083333333333333332; - const S3: f64 = -0.00019841269841269427; - const S4: f64 = 0.0000027557319223597490; - const S5: f64 = -0.000000025052108218074604; - const S6: f64 = 1.6059006598854211e-10; - const S7: f64 = -7.6429779983024564e-13; - const S8: f64 = 2.6174587166648325e-15; - - inline fn poly(z: T) T { - return S2 + z * (S3 + z * (S4 + z * (S5 + - z * (S6 + z * (S7 + z * S8))))); - } - }, - f128 => struct { - const S1: T = -0.16666666666666666666666666666666666606732416116558; - const S2: T = 0.0083333333333333333333333333333331135404851288270047; - const S3: T = -0.00019841269841269841269841269839935785325638310428717; - const S4: T = 0.27557319223985890652557316053039946268333231205686e-5; - const S5: T = -0.25052108385441718775048214826384312253862930064745e-7; - const S6: T = 0.16059043836821614596571832194524392581082444805729e-9; - const S7: T = -0.76471637318198151807063387954939213287488216303768e-12; - const S8: T = 0.28114572543451292625024967174638477283187397621303e-14; - - const S9: f64 = -0.82206352458348947812512122163446202498005154296863e-17; - const S10: f64 = 0.19572940011906109418080609928334380560135358385256e-19; - const S11: f64 = -0.38680813379701966970673724299207480965452616911420e-22; - const S12: f64 = 0.64038150078671872796678569586315881020659912139412e-25; - - inline fn poly(z: T) T { - return S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * (S8 + - z * (S9 + z * (S10 + z * (S11 + z * S12))))))))); - } - }, - else => @compileError("sinl supports only f80 and f128, got: " ++ @typeName(T)), - }; - - const z = x * x; - const v = z * x; - const r = impl.poly(z); - - if (iy == 0) { - return x + v * (impl.S1 + z * r); - } - - return x - ((z * (0.5 * y - v * r) - y) - v * impl.S1); -} - pub fn sindf(x: f64) f32 { // |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]). const S1 = -0x15555554cbac77.0p-55; // -0.166666666416265235595 @@ -239,6 +181,52 @@ pub fn sindf(x: f64) f32 { return @floatCast((x + s * (S1 + z * S2)) + s * w * r); } +pub fn sinx(x: f80, y: f80, iy: i32) f80 { + const S1: f80 = -0.166666666666666666671; + const S2: f64 = 0.0083333333333333332; + const S3: f64 = -0.00019841269841269427; + const S4: f64 = 0.0000027557319223597490; + const S5: f64 = -0.000000025052108218074604; + const S6: f64 = 1.6059006598854211e-10; + const S7: f64 = -7.6429779983024564e-13; + const S8: f64 = 2.6174587166648325e-15; + + const z = x * x; + const v = z * x; + const r = S2 + z * (S3 + z * (S4 + z * (S5 + + z * (S6 + z * (S7 + z * S8))))); + + if (iy == 0) + return x + v * (S1 + z * r); + + return x - ((z * (0.5 * y - v * r) - y) - v * S1); +} + +pub fn sinq(x: f128, y: f128, iy: i32) f128 { + const S1: f128 = -0.16666666666666666666666666666666666606732416116558; + const S2: f128 = 0.0083333333333333333333333333333331135404851288270047; + const S3: f128 = -0.00019841269841269841269841269839935785325638310428717; + const S4: f128 = 0.27557319223985890652557316053039946268333231205686e-5; + const S5: f128 = -0.25052108385441718775048214826384312253862930064745e-7; + const S6: f128 = 0.16059043836821614596571832194524392581082444805729e-9; + const S7: f128 = -0.76471637318198151807063387954939213287488216303768e-12; + const S8: f128 = 0.28114572543451292625024967174638477283187397621303e-14; + const S9: f64 = -0.82206352458348947812512122163446202498005154296863e-17; + const S10: f64 = 0.19572940011906109418080609928334380560135358385256e-19; + const S11: f64 = -0.38680813379701966970673724299207480965452616911420e-22; + const S12: f64 = 0.64038150078671872796678569586315881020659912139412e-25; + + const z = x * x; + const v = z * x; + const r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * (S8 + + z * (S9 + z * (S10 + z * (S11 + z * S12))))))))); + + if (iy == 0) + return x + v * (S1 + z * r); + + return x - ((z * (0.5 * y - v * r) - y) - v * S1); +} + /// kernel tan function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854 /// Input x is assumed to be bounded by ~pi/4 in magnitude. /// Input y is the tail of x. @@ -377,88 +365,112 @@ pub fn tandf(x: f64, odd: bool) f32 { return @floatCast(if (odd) -1.0 / r0 else r0); } -pub fn tanl(comptime T: type, x_: T, y_: T, odd: i32) T { +pub fn tanx(x_: f80, y_: f80, odd: i32) f80 { + const pio4: f80 = 0.785398163397448309628; + const pio4lo: f80 = -1.25413940316708300586e-20; + + const T3: f80 = 0.333333333333333333180; + const T5: f80 = 0.133333333333333372290; + const T7: f80 = 0.0539682539682504975744; + const T9: f64 = 0.021869488536312216; + const T11: f64 = 0.0088632355256619590; + const T13: f64 = 0.0035921281113786528; + const T15: f64 = 0.0014558334756312418; + const T17: f64 = 0.00059003538700862256; + const T19: f64 = 0.00023907843576635544; + const T21: f64 = 0.000097154625656538905; + const T23: f64 = 0.000038440165747303162; + const T25: f64 = 0.000018082171885432524; + const T27: f64 = 0.0000024196006108814377; + const T29: f64 = 0.0000078293456938132840; + const T31: f64 = -0.0000032609076735050182; + const T33: f64 = 0.0000023261313142559411; + + var x = x_; + var y = y_; + const big = @abs(x) >= 0.67434; + var sign: i8 = 0; + + if (big) { + if (x < 0) { + sign = -1; + x = -x; + y = -y; + } + x = (pio4 - x) + (pio4lo - y); + y = 0.0; + } + + var z = x * x; + var w = z * z; + + var r = T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + + w * (T25 + w * (T29 + w * T33)))))); + + var v = z * (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + + w * (T27 + w * T31)))))); + + var s = z * x; + r = y + z * (s * (r + v) + y) + T3 * s; + w = x + r; + + if (big) { + s = @as(f80, @floatFromInt(1 - 2 * odd)); + v = s - 2.0 * (x + (r - w * w / (w + s))); + return if (sign == -1) -v else v; + } + + if (odd == 0) { + return w; + } + + // if allow error up to 2 ulp, simply return + // -1.0 / (x+r) here + // + // compute -1.0 / (x+r) accurately + z = w + 0x1p32 - 0x1p32; + v = r - (z - x); + const a = -1.0 / w; + const t = a + 0x1p32 - 0x1p32; + s = 1.0 + t * z; + return t + a * (s + t * v); +} + +pub fn tanq(x_: f128, y_: f128, odd: i32) f128 { + const pio4: f128 = 0x1.921fb54442d18469898cc51701b8p-1; + const pio4lo: f128 = 0x1.cd129024e088a67cc74020bbea60p-116; + + const T3: f128 = 0x1.5555555555555555555555555553p-2; + const T5: f128 = 0x1.1111111111111111111111111eb5p-3; + const T7: f128 = 0x1.ba1ba1ba1ba1ba1ba1ba1b694cd6p-5; + const T9: f128 = 0x1.664f4882c10f9f32d6bbe09d8bcdp-6; + const T11: f128 = 0x1.226e355e6c23c8f5b4f5762322eep-7; + const T13: f128 = 0x1.d6d3d0e157ddfb5fed8e84e27b37p-9; + const T15: f128 = 0x1.7da36452b75e2b5fce9ee7c2c92ep-10; + const T17: f128 = 0x1.355824803674477dfcf726649efep-11; + const T19: f128 = 0x1.f57d7734d1656e0aceb716f614c2p-13; + const T21: f128 = 0x1.967e18afcb180ed942dfdc518d6cp-14; + const T23: f128 = 0x1.497d8eea21e95bc7e2aa79b9f2cdp-15; + const T25: f128 = 0x1.0b132d39f055c81be49eff7afd50p-16; + const T27: f128 = 0x1.b0f72d33eff7bfa2fbc1059d90b6p-18; + const T29: f128 = 0x1.5ef2daf21d1113df38d0fbc00267p-19; + const T31: f128 = 0x1.1c77d6eac0234988cdaa04c96626p-20; + const T33: f128 = 0x1.cd2a5a292b180e0bdd701057dfe3p-22; + const T35: f128 = 0x1.75c7357d0298c01a31d0a6f7d518p-23; + const T37: f128 = 0x1.2f3190f4718a9a520f98f50081fcp-24; + const T39: f64 = 0.000000028443389121318352; + const T41: f64 = 0.000000011981013102001973; + const T43: f64 = 0.0000000038303578044958070; + const T45: f64 = 0.0000000034664378216909893; + const T47: f64 = -0.0000000015090641701997785; + const T49: f64 = 0.0000000029449552300483952; + const T51: f64 = -0.0000000022006995706097711; + const T53: f64 = 0.0000000015468200913196612; + const T55: f64 = -0.00000000061311613386849674; + const T57: f64 = 1.4912469681508012e-10; + var x = x_; var y = y_; - const impl = switch (T) { - f80 => struct { - const pio4: T = 0.785398163397448309628; - const pio4lo: T = -1.25413940316708300586e-20; - - const T3: T = 0.333333333333333333180; - const T5: T = 0.133333333333333372290; - const T7: T = 0.0539682539682504975744; - const T9: f64 = 0.021869488536312216; - const T11: f64 = 0.0088632355256619590; - const T13: f64 = 0.0035921281113786528; - const T15: f64 = 0.0014558334756312418; - const T17: f64 = 0.00059003538700862256; - const T19: f64 = 0.00023907843576635544; - const T21: f64 = 0.000097154625656538905; - const T23: f64 = 0.000038440165747303162; - const T25: f64 = 0.000018082171885432524; - const T27: f64 = 0.0000024196006108814377; - const T29: f64 = 0.0000078293456938132840; - const T31: f64 = -0.0000032609076735050182; - const T33: f64 = 0.0000023261313142559411; - - inline fn rpoly(w: T) T { - return T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + - w * (T25 + w * (T29 + w * T33)))))); - } - - inline fn vpoly(w: T) T { - return T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + - w * (T27 + w * T31))))); - } - }, - f128 => struct { - const pio4: T = 0x1.921fb54442d18469898cc51701b8p-1; - const pio4lo: T = 0x1.cd129024e088a67cc74020bbea60p-116; - - const T3: T = 0x1.5555555555555555555555555553p-2; - const T5: T = 0x1.1111111111111111111111111eb5p-3; - const T7: T = 0x1.ba1ba1ba1ba1ba1ba1ba1b694cd6p-5; - const T9: T = 0x1.664f4882c10f9f32d6bbe09d8bcdp-6; - const T11: T = 0x1.226e355e6c23c8f5b4f5762322eep-7; - const T13: T = 0x1.d6d3d0e157ddfb5fed8e84e27b37p-9; - const T15: T = 0x1.7da36452b75e2b5fce9ee7c2c92ep-10; - const T17: T = 0x1.355824803674477dfcf726649efep-11; - const T19: T = 0x1.f57d7734d1656e0aceb716f614c2p-13; - const T21: T = 0x1.967e18afcb180ed942dfdc518d6cp-14; - const T23: T = 0x1.497d8eea21e95bc7e2aa79b9f2cdp-15; - const T25: T = 0x1.0b132d39f055c81be49eff7afd50p-16; - const T27: T = 0x1.b0f72d33eff7bfa2fbc1059d90b6p-18; - const T29: T = 0x1.5ef2daf21d1113df38d0fbc00267p-19; - const T31: T = 0x1.1c77d6eac0234988cdaa04c96626p-20; - const T33: T = 0x1.cd2a5a292b180e0bdd701057dfe3p-22; - const T35: T = 0x1.75c7357d0298c01a31d0a6f7d518p-23; - const T37: T = 0x1.2f3190f4718a9a520f98f50081fcp-24; - const T39: f64 = 0.000000028443389121318352; - const T41: f64 = 0.000000011981013102001973; - const T43: f64 = 0.0000000038303578044958070; - const T45: f64 = 0.0000000034664378216909893; - const T47: f64 = -0.0000000015090641701997785; - const T49: f64 = 0.0000000029449552300483952; - const T51: f64 = -0.0000000022006995706097711; - const T53: f64 = 0.0000000015468200913196612; - const T55: f64 = -0.00000000061311613386849674; - const T57: f64 = 1.4912469681508012e-10; - - inline fn rpoly(w: T) T { - return T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + - w * (T25 + w * (T29 + w * (T33 + w * (T37 + w * (T41 + - w * (T45 + w * (T49 + w * (T53 + w * T57)))))))))))); - } - - inline fn vpoly(w: T) T { - return T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + - w * (T27 + w * (T31 + w * (T35 + w * (T39 + w * (T43 + - w * (T47 + w * (T51 + w * T55))))))))))); - } - }, - else => @compileError("tanl supports only f80 and f128, got: " ++ @typeName(T)), - }; const big = @abs(x) >= 0.67434; var sign: i8 = 0; @@ -469,20 +481,27 @@ pub fn tanl(comptime T: type, x_: T, y_: T, odd: i32) T { x = -x; y = -y; } - x = (impl.pio4 - x) + (impl.pio4lo - y); + x = (pio4 - x) + (pio4lo - y); y = 0.0; } var z = x * x; var w = z * z; - var r = impl.rpoly(w); - var v = z * impl.vpoly(w); + + var r = T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + + w * (T25 + w * (T29 + w * (T33 + w * (T37 + w * (T41 + + w * (T45 + w * (T49 + w * (T53 + w * T57)))))))))))); + + var v = z * (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + + w * (T27 + w * (T31 + w * (T35 + w * (T39 + w * (T43 + + w * (T47 + w * (T51 + w * T55)))))))))))); + var s = z * x; - r = y + z * (s * (r + v) + y) + impl.T3 * s; + r = y + z * (s * (r + v) + y) + T3 * s; w = x + r; if (big) { - s = @as(T, @floatFromInt(1 - 2 * odd)); + s = @as(f128, @floatFromInt(1 - 2 * odd)); v = s - 2.0 * (x + (r - w * w / (w + s))); return if (sign == -1) -v else v; } -- 2.54.0