diff --git a/lib/compiler_rt/math_utils.zig b/lib/compiler_rt/math_utils.zig index 76cbcec468..70e6a4b3b7 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 ef3709bc18..78a56d9b8c 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 6534ca3c27..746982a0d6 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 4d1690e4fe..8a20f2c2ac 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 1673b34be9..0000000000 --- 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 0bbb71d379..0000000000 --- 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 6d766b098d..0000000000 --- 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 2525bbe866..0000000000 --- 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 9c0b16eed3..0000000000 --- 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 4a16f4e0a6..f8070ccbf8 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 074f9e626e..d0375fb03e 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 cae605d124..09897473d9 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",