authorgravatar for hi@mihaelm.commihael <hi@mihaelm.com> 2026-03-24 19:33:55+01:00
committergravatar for hi@mihaelm.commihael <hi@mihaelm.com> 2026-04-02 23:54:20+02:00
logb5ec3e597e2cd374853dea5244b2bc2bccccd96a
tree199c5273b0b637e3b12d5543875e5ca176ea1f3e
parentffd6f6cc6e05e84a4e029b02ff13a3e84cbf1aa4
signaturebadge-check Signed by SSH key SHA256:aoFoShdYLdrqMichqKXSSieTKUfACUIDJHsKc4V2tQg

`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=<LIBC-TEST-PATH> -Dtest-filter=sinl -fqemu -fwasmtime --summary line Build Summary: 553/553 steps succeeded ```

12 files changed, 201 insertions(+), 329 deletions(-)

lib/compiler_rt/math_utils.zig+2
...@@ -1,6 +1,8 @@...@@ -1,6 +1,8 @@
1const std = @import("std");1const std = @import("std");
22
3pub const U80 = std.meta.Int(.unsigned, 80);3pub const U80 = std.meta.Int(.unsigned, 80);
4/// pi divided by 4
5pub const pi_4 = 0.78539816339744830962;
46
5/// Returns the sign + exponent bits of a `long double`7/// Returns the sign + exponent bits of a `long double`
6pub fn ldSignExponent(x: anytype) u16 {8pub fn ldSignExponent(x: anytype) u16 {
lib/compiler_rt/sin.zig+98-33
...@@ -3,17 +3,21 @@...@@ -3,17 +3,21 @@
3//!3//!
4//! https://git.musl-libc.org/cgit/musl/tree/src/math/sinf.c4//! https://git.musl-libc.org/cgit/musl/tree/src/math/sinf.c
5//! https://git.musl-libc.org/cgit/musl/tree/src/math/sin.c5//! https://git.musl-libc.org/cgit/musl/tree/src/math/sin.c
6//! https://git.musl-libc.org/cgit/musl/tree/src/math/sinl.c
67
7const std = @import("std");8const std = @import("std");
8const math = std.math;9const math = std.math;
9const mem = std.mem;10const mem = std.mem;
10const expect = std.testing.expect;11const expect = std.testing.expect;
12const expectApproxEqAbs = std.testing.expectApproxEqAbs;
1113
12const compiler_rt = @import("../compiler_rt.zig");14const compiler_rt = @import("../compiler_rt.zig");
13const symbol = @import("../compiler_rt.zig").symbol;15const symbol = @import("../compiler_rt.zig").symbol;
14const trig = @import("trig.zig");16const trig = @import("trig.zig");
15const rem_pio2 = @import("rem_pio2.zig").rem_pio2;17const rem_pio2 = @import("rem_pio2.zig").rem_pio2;
16const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f;18const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f;
19const rem_pio2l = @import("rem_pio2l.zig").rem_pio2l;
20const utils = @import("math_utils.zig");
1721
18comptime {22comptime {
19 symbol(&__sinh, "__sinh");23 symbol(&__sinh, "__sinh");
...@@ -128,14 +132,39 @@ pub fn sin(x: f64) callconv(.c) f64 {...@@ -128,14 +132,39 @@ pub fn sin(x: f64) callconv(.c) f64 {
128 };132 };
129}133}
130134
135fn sinlGeneric(comptime T: type, x: T) T {
136 const se = utils.ldSignExponent(x) & 0x7fff;
137 if (se == 0x7fff) {
138 return x - x;
139 }
140
141 if (@abs(x) < utils.pi_4) {
142 if (se < 0x3fff - (math.floatMantissaBits(T) / 2)) {
143 // raise inexact if x!=0 and underflow if subnormal
144 if (compiler_rt.want_float_exceptions) {
145 mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120);
146 }
147 return x;
148 }
149 return trig.__sinl(T, x, 0.0, 0);
150 }
151
152 var y: [2]T = undefined;
153 const n = rem_pio2l(T, x, &y);
154 return switch (n & 3) {
155 0 => trig.__sinl(T, y[0], y[1], 1),
156 1 => trig.__cosl(T, y[0], y[1]),
157 2 => -trig.__sinl(T, y[0], y[1], 1),
158 else => -trig.__cosl(T, y[0], y[1]),
159 };
160}
161
131pub fn __sinx(x: f80) callconv(.c) f80 {162pub fn __sinx(x: f80) callconv(.c) f80 {
132 // TODO: more efficient implementation163 return sinlGeneric(f80, x);
133 return @floatCast(sinq(x));
134}164}
135165
136pub fn sinq(x: f128) callconv(.c) f128 {166pub fn sinq(x: f128) callconv(.c) f128 {
137 // TODO: more correct implementation167 return sinlGeneric(f128, x);
138 return sin(@floatCast(x));
139}168}
140169
141pub fn sinl(x: c_longdouble) callconv(.c) c_longdouble {170pub fn sinl(x: c_longdouble) callconv(.c) c_longdouble {
...@@ -149,44 +178,80 @@ pub fn sinl(x: c_longdouble) callconv(.c) c_longdouble {...@@ -149,44 +178,80 @@ pub fn sinl(x: c_longdouble) callconv(.c) c_longdouble {
149 }178 }
150}179}
151180
152test "sin32" {181fn testSinSpecial(comptime T: type) !void {
153 const epsilon = 0.00001;182 const f = switch (T) {
183 f32 => sinf,
184 f64 => sin,
185 f80 => __sinx,
186 f128 => sinq,
187 else => @compileError("unimplemented"),
188 };
154189
155 try expect(math.approxEqAbs(f32, sinf(0.0), 0.0, epsilon));190 try expect(math.isPositiveZero(f(0.0)));
156 try expect(math.approxEqAbs(f32, sinf(0.2), 0.198669, epsilon));191 try expect(math.isNegativeZero(f(-0.0)));
157 try expect(math.approxEqAbs(f32, sinf(0.8923), 0.778517, epsilon));192 try expect(math.isNan(f(math.inf(T))));
158 try expect(math.approxEqAbs(f32, sinf(1.5), 0.997495, epsilon));193 try expect(math.isNan(f(-math.inf(T))));
159 try expect(math.approxEqAbs(f32, sinf(-1.5), -0.997495, epsilon));194 try expect(math.isNan(f(math.nan(T))));
160 try expect(math.approxEqAbs(f32, sinf(37.45), -0.246544, epsilon));
161 try expect(math.approxEqAbs(f32, sinf(89.123), 0.916166, epsilon));
162}195}
163196
164test "sin64" {197test "sin32.normal" {
165 const epsilon = 0.000001;198 const epsilon = math.floatEps(f32);
166199 try expectApproxEqAbs(@as(f32, 0.0), sinf(0.0), epsilon);
167 try expect(math.approxEqAbs(f64, sin(0.0), 0.0, epsilon));200 try expectApproxEqAbs(@as(f32, 0.19866933), sinf(0.2), epsilon);
168 try expect(math.approxEqAbs(f64, sin(0.2), 0.198669, epsilon));201 try expectApproxEqAbs(@as(f32, 0.77851737), sinf(0.8923), epsilon);
169 try expect(math.approxEqAbs(f64, sin(0.8923), 0.778517, epsilon));202 try expectApproxEqAbs(@as(f32, 0.997495), sinf(1.5), epsilon);
170 try expect(math.approxEqAbs(f64, sin(1.5), 0.997495, epsilon));203 try expectApproxEqAbs(@as(f32, -0.997495), sinf(-1.5), epsilon);
171 try expect(math.approxEqAbs(f64, sin(-1.5), -0.997495, epsilon));204 try expectApproxEqAbs(@as(f32, -0.24654257), sinf(37.45), epsilon);
172 try expect(math.approxEqAbs(f64, sin(37.45), -0.246543, epsilon));205 try expectApproxEqAbs(@as(f32, 0.9161657), sinf(89.123), epsilon);
173 try expect(math.approxEqAbs(f64, sin(89.123), 0.916166, epsilon));
174}206}
175207
176test "sin32.special" {208test "sin32.special" {
177 try expect(sinf(0.0) == 0.0);209 try testSinSpecial(f32);
178 try expect(sinf(-0.0) == -0.0);210}
179 try expect(math.isNan(sinf(math.inf(f32))));211
180 try expect(math.isNan(sinf(-math.inf(f32))));212test "sin64.normal" {
181 try expect(math.isNan(sinf(math.nan(f32))));213 const epsilon = math.floatEps(f64);
214 try expectApproxEqAbs(@as(f64, 0.0), sin(0.0), epsilon);
215 try expectApproxEqAbs(@as(f64, 0.19866933079506122), sin(0.2), epsilon);
216 try expectApproxEqAbs(@as(f64, 0.7785173385577349), sin(0.8923), epsilon);
217 try expectApproxEqAbs(@as(f64, 0.9974949866040544), sin(1.5), epsilon);
218 try expectApproxEqAbs(@as(f64, -0.9974949866040544), sin(-1.5), epsilon);
219 try expectApproxEqAbs(@as(f64, -0.24654331551411082), sin(37.45), epsilon);
220 try expectApproxEqAbs(@as(f64, 0.9161652766622714), sin(89.123), epsilon);
182}221}
183222
184test "sin64.special" {223test "sin64.special" {
185 try expect(sin(0.0) == 0.0);224 try testSinSpecial(f64);
186 try expect(sin(-0.0) == -0.0);225}
187 try expect(math.isNan(sin(math.inf(f64))));226
188 try expect(math.isNan(sin(-math.inf(f64))));227test "sin80.normal" {
189 try expect(math.isNan(sin(math.nan(f64))));228 const epsilon = math.floatEps(f80);
229 try expectApproxEqAbs(@as(f80, 0.0), __sinx(0.0), epsilon);
230 try expectApproxEqAbs(@as(f80, 0.19866933079506121545941262711838975), __sinx(0.2), epsilon);
231 try expectApproxEqAbs(@as(f80, 0.77851733855773487830689285621486050), __sinx(0.8923), epsilon);
232 try expectApproxEqAbs(@as(f80, 0.99749498660405443094172337114148732), __sinx(1.5), epsilon);
233 try expectApproxEqAbs(@as(f80, -0.99749498660405443094172337114148732), __sinx(-1.5), epsilon);
234 try expectApproxEqAbs(@as(f80, -0.24654331551411356504), __sinx(37.45), epsilon);
235 try expectApproxEqAbs(@as(f80, 0.91616527666226951006), __sinx(89.123), epsilon);
236}
237
238test "sin80.special" {
239 try testSinSpecial(f80);
240}
241
242test "sin128.normal" {
243 const epsilon = math.floatEps(f128);
244 try expectApproxEqAbs(@as(f128, 0.0), sinq(0.0), epsilon);
245 try expectApproxEqAbs(@as(f128, 0.19866933079506121545941262711838975), sinq(0.2), epsilon);
246 try expectApproxEqAbs(@as(f128, 0.77851733855773487830689285621486050), sinq(0.8923), epsilon);
247 try expectApproxEqAbs(@as(f128, 0.99749498660405443094172337114148732), sinq(1.5), epsilon);
248 try expectApproxEqAbs(@as(f128, -0.99749498660405443094172337114148732), sinq(-1.5), epsilon);
249 try expectApproxEqAbs(@as(f128, -0.24654331551411356571238581321661085), sinq(37.45), epsilon);
250 try expectApproxEqAbs(@as(f128, 0.91616527666226951075019849560482170), sinq(89.123), epsilon);
251}
252
253test "sin128.special" {
254 try testSinSpecial(f128);
190}255}
191256
192test "sin32 #9901" {257test "sin32 #9901" {
lib/compiler_rt/tan.zig+1-2
...@@ -130,8 +130,7 @@ fn tanlGeneric(comptime T: type, x: T) T {...@@ -130,8 +130,7 @@ fn tanlGeneric(comptime T: type, x: T) T {
130 return x - x;130 return x - x;
131 }131 }
132132
133 const pi_4 = 0.78539816339744830962;133 if (@abs(x) < utils.pi_4) {
134 if (@abs(x) < pi_4) {
135 if (se < 0x3fff - math.floatMantissaBits(T) / 2) {134 if (se < 0x3fff - math.floatMantissaBits(T) / 2) {
136 if (compiler_rt.want_float_exceptions) {135 if (compiler_rt.want_float_exceptions) {
137 mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120);136 mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120);
lib/compiler_rt/trig.zig+100
...@@ -7,6 +7,8 @@...@@ -7,6 +7,8 @@
7// https://git.musl-libc.org/cgit/musl/tree/src/math/__sindf.c7// https://git.musl-libc.org/cgit/musl/tree/src/math/__sindf.c
8// https://git.musl-libc.org/cgit/musl/tree/src/math/__tand.c8// https://git.musl-libc.org/cgit/musl/tree/src/math/__tand.c
9// https://git.musl-libc.org/cgit/musl/tree/src/math/__tandf.c9// https://git.musl-libc.org/cgit/musl/tree/src/math/__tandf.c
10// https://git.musl-libc.org/cgit/musl/tree/src/math/__sinl.c
11// https://git.musl-libc.org/cgit/musl/tree/src/math/__cosl.c
10// https://git.musl-libc.org/cgit/musl/tree/src/math/__tanl.c12// https://git.musl-libc.org/cgit/musl/tree/src/math/__tanl.c
1113
12/// kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.78539816414/// kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
...@@ -74,6 +76,52 @@ pub fn __cosdf(x: f64) f32 {...@@ -74,6 +76,52 @@ pub fn __cosdf(x: f64) f32 {
74 return @floatCast(((1.0 + z * C0) + w * C1) + (w * z) * r);76 return @floatCast(((1.0 + z * C0) + w * C1) + (w * z) * r);
75}77}
7678
79pub fn __cosl(comptime T: type, x: T, y: T) T {
80 const impl = switch (T) {
81 f80 => struct {
82 const C1: T = 0.0416666666666666666136;
83
84 const C2: f64 = -0.0013888888888888874;
85 const C3: f64 = 0.000024801587301571716;
86 const C4: f64 = -0.00000027557319215507120;
87 const C5: f64 = 0.0000000020876754400407278;
88 const C6: f64 = -1.1470297442401303e-11;
89 const C7: f64 = 4.7383039476436467e-14;
90
91 inline fn poly(z: T) T {
92 return z * (C1 + z * (C2 + z * (C3 + z * (C4 +
93 z * (C5 + z * (C6 + z * C7))))));
94 }
95 },
96 f128 => struct {
97 const C1: T = 0.04166666666666666666666666666666658424671;
98 const C2: T = -0.001388888888888888888888888888863490893732;
99 const C3: T = 0.00002480158730158730158730158600795304914210;
100 const C4: T = -0.2755731922398589065255474947078934284324e-6;
101 const C5: T = 0.2087675698786809897659225313136400793948e-8;
102 const C6: T = -0.1147074559772972315817149986812031204775e-10;
103 const C7: T = 0.4779477332386808976875457937252120293400e-13;
104
105 const C8: f64 = -0.1561920696721507929516718307820958119868e-15;
106 const C9: f64 = 0.4110317413744594971475941557607804508039e-18;
107 const C10: f64 = -0.8896592467191938803288521958313920156409e-21;
108 const C11: f64 = 0.1601061435794535138244346256065192782581e-23;
109
110 inline fn poly(z: T) T {
111 return z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 +
112 z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * C11))))))))));
113 }
114 },
115 else => @compileError("__cosl supports only f80 and f128, got: " ++ @typeName(T)),
116 };
117
118 const z = x * x;
119 const r = impl.poly(z);
120 const hz = 0.5 * z;
121 const w = 1.0 - hz;
122 return w + (((1.0 - w) - hz) + (z * r - x * y));
123}
124
77/// kernel sin function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854125/// kernel sin function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
78/// Input x is assumed to be bounded by ~pi/4 in magnitude.126/// Input x is assumed to be bounded by ~pi/4 in magnitude.
79/// Input y is the tail of x.127/// Input y is the tail of x.
...@@ -120,6 +168,58 @@ pub fn __sin(x: f64, y: f64, iy: i32) f64 {...@@ -120,6 +168,58 @@ pub fn __sin(x: f64, y: f64, iy: i32) f64 {
120 }168 }
121}169}
122170
171pub fn __sinl(comptime T: type, x: T, y: T, iy: i32) T {
172 const impl = switch (T) {
173 f80 => struct {
174 const S1: T = -0.166666666666666666671;
175
176 const S2: f64 = 0.0083333333333333332;
177 const S3: f64 = -0.00019841269841269427;
178 const S4: f64 = 0.0000027557319223597490;
179 const S5: f64 = -0.000000025052108218074604;
180 const S6: f64 = 1.6059006598854211e-10;
181 const S7: f64 = -7.6429779983024564e-13;
182 const S8: f64 = 2.6174587166648325e-15;
183
184 inline fn poly(z: T) T {
185 return S2 + z * (S3 + z * (S4 + z * (S5 +
186 z * (S6 + z * (S7 + z * S8)))));
187 }
188 },
189 f128 => struct {
190 const S1: T = -0.16666666666666666666666666666666666606732416116558;
191 const S2: T = 0.0083333333333333333333333333333331135404851288270047;
192 const S3: T = -0.00019841269841269841269841269839935785325638310428717;
193 const S4: T = 0.27557319223985890652557316053039946268333231205686e-5;
194 const S5: T = -0.25052108385441718775048214826384312253862930064745e-7;
195 const S6: T = 0.16059043836821614596571832194524392581082444805729e-9;
196 const S7: T = -0.76471637318198151807063387954939213287488216303768e-12;
197 const S8: T = 0.28114572543451292625024967174638477283187397621303e-14;
198
199 const S9: f64 = -0.82206352458348947812512122163446202498005154296863e-17;
200 const S10: f64 = 0.19572940011906109418080609928334380560135358385256e-19;
201 const S11: f64 = -0.38680813379701966970673724299207480965452616911420e-22;
202 const S12: f64 = 0.64038150078671872796678569586315881020659912139412e-25;
203
204 inline fn poly(z: T) T {
205 return S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * (S8 +
206 z * (S9 + z * (S10 + z * (S11 + z * S12)))))))));
207 }
208 },
209 else => @compileError("__sinl supports only f80 and f128, got: " ++ @typeName(T)),
210 };
211
212 const z = x * x;
213 const v = z * x;
214 const r = impl.poly(z);
215
216 if (iy == 0) {
217 return x + v * (impl.S1 + z * r);
218 }
219
220 return x - ((z * (0.5 * y - v * r) - y) - v * impl.S1);
221}
222
123pub fn __sindf(x: f64) f32 {223pub fn __sindf(x: f64) f32 {
124 // |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]).224 // |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]).
125 const S1 = -0x15555554cbac77.0p-55; // -0.166666666416265235595225 const S1 = -0x15555554cbac77.0p-55; // -0.166666666416265235595
lib/libc/mingw/math/x86/sin.def.h deleted-65
...@@ -1,65 +0,0 @@
1/*
2 This Software is provided under the Zope Public License (ZPL) Version 2.1.
3
4 Copyright (c) 2009, 2010 by the mingw-w64 project
5
6 See the AUTHORS file for the list of contributors to the mingw-w64 project.
7
8 This license has been certified as open source. It has also been designated
9 as GPL compatible by the Free Software Foundation (FSF).
10
11 Redistribution and use in source and binary forms, with or without
12 modification, are permitted provided that the following conditions are met:
13
14 1. Redistributions in source code must retain the accompanying copyright
15 notice, this list of conditions, and the following disclaimer.
16 2. Redistributions in binary form must reproduce the accompanying
17 copyright notice, this list of conditions, and the following disclaimer
18 in the documentation and/or other materials provided with the
19 distribution.
20 3. Names of the copyright holders must not be used to endorse or promote
21 products derived from this software without prior written permission
22 from the copyright holders.
23 4. The right to distribute this software or to use it for any purpose does
24 not give you the right to use Servicemarks (sm) or Trademarks (tm) of
25 the copyright holders. Use of them is covered by separate agreement
26 with the copyright holders.
27 5. If any files are modified, you must cause the modified files to carry
28 prominent notices stating that you changed the files and the date of
29 any change.
30
31 Disclaimer
32
33 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY EXPRESSED
34 OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
35 OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
36 EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY DIRECT, INDIRECT,
37 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
38 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
39 OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
40 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
41 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
42 EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
43*/
44
45#include "../complex/complex_internal.h"
46#include <errno.h>
47
48extern long double __sinl_internal (long double);
49
50__FLT_TYPE
51__FLT_ABI(sin) (__FLT_TYPE x)
52{
53 int x_class = fpclassify (x);
54 if (x_class == FP_NAN)
55 {
56 __FLT_RPT_DOMAIN ("sin", x, 0.0, x);
57 return x;
58 }
59 else if (x_class == FP_INFINITE)
60 {
61 __FLT_RPT_DOMAIN ("sin", x, 0.0, __FLT_NAN);
62 return __FLT_NAN;
63 }
64 return (__FLT_TYPE) __sinl_internal ((long double) x);
65}
lib/libc/mingw/math/x86/sinl.c deleted-46
...@@ -1,46 +0,0 @@
1/*
2 This Software is provided under the Zope Public License (ZPL) Version 2.1.
3
4 Copyright (c) 2009, 2010 by the mingw-w64 project
5
6 See the AUTHORS file for the list of contributors to the mingw-w64 project.
7
8 This license has been certified as open source. It has also been designated
9 as GPL compatible by the Free Software Foundation (FSF).
10
11 Redistribution and use in source and binary forms, with or without
12 modification, are permitted provided that the following conditions are met:
13
14 1. Redistributions in source code must retain the accompanying copyright
15 notice, this list of conditions, and the following disclaimer.
16 2. Redistributions in binary form must reproduce the accompanying
17 copyright notice, this list of conditions, and the following disclaimer
18 in the documentation and/or other materials provided with the
19 distribution.
20 3. Names of the copyright holders must not be used to endorse or promote
21 products derived from this software without prior written permission
22 from the copyright holders.
23 4. The right to distribute this software or to use it for any purpose does
24 not give you the right to use Servicemarks (sm) or Trademarks (tm) of
25 the copyright holders. Use of them is covered by separate agreement
26 with the copyright holders.
27 5. If any files are modified, you must cause the modified files to carry
28 prominent notices stating that you changed the files and the date of
29 any change.
30
31 Disclaimer
32
33 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY EXPRESSED
34 OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
35 OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
36 EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY DIRECT, INDIRECT,
37 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
38 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
39 OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
40 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
41 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
42 EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
43*/
44
45#define _NEW_COMPLEX_LDOUBLE 1
46#include "sin.def.h"
lib/libc/mingw/math/x86/sinl_internal.S deleted-58
...@@ -1,58 +0,0 @@
1/**
2 * This file has no copyright assigned and is placed in the Public Domain.
3 * This file is part of the mingw-w64 runtime package.
4 * No warranty is given; refer to the file DISCLAIMER.PD within this package.
5 */
6#include <_mingw_mac.h>
7
8 .file "sinl_internal.S"
9 .text
10#ifdef __x86_64__
11 .align 8
12#else
13 .align 4
14#endif
15.globl __MINGW_USYMBOL(__sinl_internal)
16 .def __MINGW_USYMBOL(__sinl_internal); .scl 2; .type 32; .endef
17__MINGW_USYMBOL(__sinl_internal):
18#ifdef __x86_64__
19 fldt (%rdx)
20 fsin
21 fnstsw %ax
22 testl $0x400,%eax
23 jnz 1f
24 movq %rcx,%rax
25 movq $0,8(%rcx)
26 fstpt (%rcx)
27 ret
281: fldpi
29 fadd %st(0)
30 fxch %st(1)
312: fprem1
32 fnstsw %ax
33 testl $0x400,%eax
34 jnz 2b
35 fstp %st(1)
36 fsin
37 movq %rcx,%rax
38 movq $0,8(%rcx)
39 fstpt (%rcx)
40 ret
41#else
42 fldt 4(%esp)
43 fsin
44 fnstsw %ax
45 testl $0x400,%eax
46 jnz 1f
47 ret
481: fldpi
49 fadd %st(0)
50 fxch %st(1)
512: fprem1
52 fnstsw %ax
53 testl $0x400,%eax
54 jnz 2b
55 fstp %st(1)
56 fsin
57 ret
58#endif
lib/libc/musl/src/math/__sinl.c deleted-78
...@@ -1,78 +0,0 @@
1/* origin: FreeBSD /usr/src/lib/msun/ld80/k_sinl.c */
2/* origin: FreeBSD /usr/src/lib/msun/ld128/k_sinl.c */
3/*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
7 *
8 * Developed at SunSoft, a Sun Microsystems, Inc. business.
9 * Permission to use, copy, modify, and distribute this
10 * software is freely granted, provided that this notice
11 * is preserved.
12 * ====================================================
13 */
14
15#include "libm.h"
16
17#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
18#if LDBL_MANT_DIG == 64
19/*
20 * ld80 version of __sin.c. See __sin.c for most comments.
21 */
22/*
23 * Domain [-0.7854, 0.7854], range ~[-1.89e-22, 1.915e-22]
24 * |sin(x)/x - s(x)| < 2**-72.1
25 *
26 * See __cosl.c for more details about the polynomial.
27 */
28static const long double
29S1 = -0.166666666666666666671L; /* -0xaaaaaaaaaaaaaaab.0p-66 */
30static const double
31S2 = 0.0083333333333333332, /* 0x11111111111111.0p-59 */
32S3 = -0.00019841269841269427, /* -0x1a01a01a019f81.0p-65 */
33S4 = 0.0000027557319223597490, /* 0x171de3a55560f7.0p-71 */
34S5 = -0.000000025052108218074604, /* -0x1ae64564f16cad.0p-78 */
35S6 = 1.6059006598854211e-10, /* 0x161242b90243b5.0p-85 */
36S7 = -7.6429779983024564e-13, /* -0x1ae42ebd1b2e00.0p-93 */
37S8 = 2.6174587166648325e-15; /* 0x179372ea0b3f64.0p-101 */
38#define POLY(z) (S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*S8))))))
39#elif LDBL_MANT_DIG == 113
40/*
41 * ld128 version of __sin.c. See __sin.c for most comments.
42 */
43/*
44 * Domain [-0.7854, 0.7854], range ~[-1.53e-37, 1.659e-37]
45 * |sin(x)/x - s(x)| < 2**-122.1
46 *
47 * See __cosl.c for more details about the polynomial.
48 */
49static const long double
50S1 = -0.16666666666666666666666666666666666606732416116558L,
51S2 = 0.0083333333333333333333333333333331135404851288270047L,
52S3 = -0.00019841269841269841269841269839935785325638310428717L,
53S4 = 0.27557319223985890652557316053039946268333231205686e-5L,
54S5 = -0.25052108385441718775048214826384312253862930064745e-7L,
55S6 = 0.16059043836821614596571832194524392581082444805729e-9L,
56S7 = -0.76471637318198151807063387954939213287488216303768e-12L,
57S8 = 0.28114572543451292625024967174638477283187397621303e-14L;
58static const double
59S9 = -0.82206352458348947812512122163446202498005154296863e-17,
60S10 = 0.19572940011906109418080609928334380560135358385256e-19,
61S11 = -0.38680813379701966970673724299207480965452616911420e-22,
62S12 = 0.64038150078671872796678569586315881020659912139412e-25;
63#define POLY(z) (S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*(S8+ \
64 z*(S9+z*(S10+z*(S11+z*S12))))))))))
65#endif
66
67long double __sinl(long double x, long double y, int iy)
68{
69 long double z,r,v;
70
71 z = x*x;
72 v = z*x;
73 r = POLY(z);
74 if (iy == 0)
75 return x+v*(S1+z*r);
76 return x-((z*(0.5*y-v*r)-y)-v*S1);
77}
78#endif
lib/libc/musl/src/math/sinl.c deleted-41
...@@ -1,41 +0,0 @@
1#include "libm.h"
2
3#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
4long double sinl(long double x)
5{
6 return sin(x);
7}
8#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
9long double sinl(long double x)
10{
11 union ldshape u = {x};
12 unsigned n;
13 long double y[2], hi, lo;
14
15 u.i.se &= 0x7fff;
16 if (u.i.se == 0x7fff)
17 return x - x;
18 if (u.f < M_PI_4) {
19 if (u.i.se < 0x3fff - LDBL_MANT_DIG/2) {
20 /* raise inexact if x!=0 and underflow if subnormal */
21 FORCE_EVAL(u.i.se == 0 ? x*0x1p-120f : x+0x1p120f);
22 return x;
23 }
24 return __sinl(x, 0.0, 0);
25 }
26 n = __rem_pio2l(x, y);
27 hi = y[0];
28 lo = y[1];
29 switch (n & 3) {
30 case 0:
31 return __sinl(hi, lo, 1);
32 case 1:
33 return __cosl(hi, lo);
34 case 2:
35 return -__sinl(hi, lo, 1);
36 case 3:
37 default:
38 return -__cosl(hi, lo);
39 }
40}
41#endif
src/libs/mingw.zig-2
...@@ -963,8 +963,6 @@ const mingw32_x86_src = [_][]const u8{...@@ -963,8 +963,6 @@ const mingw32_x86_src = [_][]const u8{
963 "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "scalbn.S",963 "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "scalbn.S",
964 "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "scalbnf.S",964 "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "scalbnf.S",
965 "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "scalbnl.S",965 "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "scalbnl.S",
966 "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "sinl.c",
967 "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "sinl_internal.S",
968 // ucrtbase966 // ucrtbase
969 "math" ++ path.sep_str ++ "nextafterl.c",967 "math" ++ path.sep_str ++ "nextafterl.c",
970 "math" ++ path.sep_str ++ "nexttoward.c",968 "math" ++ path.sep_str ++ "nexttoward.c",
src/libs/musl.zig-2
...@@ -979,8 +979,6 @@ const src_files = [_][]const u8{...@@ -979,8 +979,6 @@ const src_files = [_][]const u8{
979 "musl/src/math/sinh.c",979 "musl/src/math/sinh.c",
980 "musl/src/math/sinhf.c",980 "musl/src/math/sinhf.c",
981 "musl/src/math/sinhl.c",981 "musl/src/math/sinhl.c",
982 "musl/src/math/__sinl.c",
983 "musl/src/math/sinl.c",
984 "musl/src/math/__tan.c",982 "musl/src/math/__tan.c",
985 "musl/src/math/__tandf.c",983 "musl/src/math/__tandf.c",
986 "musl/src/math/tanhl.c",984 "musl/src/math/tanhl.c",
src/libs/wasi_libc.zig-2
...@@ -780,8 +780,6 @@ const libc_top_half_src_files = [_][]const u8{...@@ -780,8 +780,6 @@ const libc_top_half_src_files = [_][]const u8{
780 "musl/src/math/__sin.c",780 "musl/src/math/__sin.c",
781 "musl/src/math/__sindf.c",781 "musl/src/math/__sindf.c",
782 "musl/src/math/sinhl.c",782 "musl/src/math/sinhl.c",
783 "musl/src/math/__sinl.c",
784 "musl/src/math/sinl.c",
785 "musl/src/math/__tan.c",783 "musl/src/math/__tan.c",
786 "musl/src/math/__tandf.c",784 "musl/src/math/__tandf.c",
787 "musl/src/math/tanhl.c",785 "musl/src/math/tanhl.c",