authorgravatar for rpkak@noreply.codeberg.orgrpkak <rpkak@noreply.codeberg.org> 2026-01-09 07:00:23+01:00
committergravatar for rpkak@noreply.codeberg.orgrpkak <rpkak@noreply.codeberg.org> 2026-01-10 00:09:54+01:00
log08d37d6e14886ad6fe1f5680226f92e9c2634388
tree5ad0f01a9e99e00413d3fcc85ee8c4488f07191d
parent813ae89208ae39181ecfbd6b76a0e8fd9a38e97b
signaturebadge-check Signed by SSH key SHA256:A2g6ATENrgLKQcUdthsn72P0njmQxnHSimzwnyJsZzo

compiler_rt: fix sqrtf implementation; libc: remove sqrt implementations

39 files changed, 42 insertions(+), 1119 deletions(-)

lib/compiler_rt/sqrt.zig+42-60
...@@ -54,7 +54,7 @@ pub fn __sqrth(x: f16) callconv(.c) f16 {...@@ -54,7 +54,7 @@ pub fn __sqrth(x: f16) callconv(.c) f16 {
54 // m: 2.14 r: 0.16, s: 2.14, d: 2.14, u: 2.14, three: 2.1454 // m: 2.14 r: 0.16, s: 2.14, d: 2.14, u: 2.14, three: 2.14
55 const three: u16 = 0xC000;55 const three: u16 = 0xC000;
56 const i: usize = @intCast((ix >> 4) & 0x7F);56 const i: usize = @intCast((ix >> 4) & 0x7F);
57 const r = __rsqrt_tab[i];57 const r = rsqrt_tab[i];
58 // |r*sqrt(m) - 1| < 0x1p-858 // |r*sqrt(m) - 1| < 0x1p-8
59 var s = mul16(m, r);59 var s = mul16(m, r);
60 // |s/sqrt(m) - 1| < 0x1p-860 // |s/sqrt(m) - 1| < 0x1p-8
...@@ -92,74 +92,56 @@ pub fn __sqrth(x: f16) callconv(.c) f16 {...@@ -92,74 +92,56 @@ pub fn __sqrth(x: f16) callconv(.c) f16 {
9292
93pub fn sqrtf(x: f32) callconv(.c) f32 {93pub fn sqrtf(x: f32) callconv(.c) f32 {
94 var ix: u32 = @bitCast(x);94 var ix: u32 = @bitCast(x);
95 var top = ix >> 23;
9695
97 // special case handling.96 if (ix < @as(u32, @bitCast(@as(f32, 0x1p-126))) or @as(u32, @bitCast(std.math.inf(f32))) <= ix) {
98 if (top -% 0x01 >= 0xFF - 0x01) {
99 @branchHint(.unlikely);97 @branchHint(.unlikely);
100 // x < 0x1p-126 or inf or nan.98
101 if (ix & 0x7FFF_FFFF == 0) return x;99 if (ix & 0x7fffffff == 0)
102 if (ix == 0x7F80_0000) return x;100 return x;
103 if (ix > 0x7F80_0000) return math.nan(f32);101
104 // x is subnormal, normalize it.102 if (ix == @as(u32, @bitCast(std.math.inf(f32))))
105 ix = @bitCast(x * 0x1p23);103 return x;
106 top = (ix >> 23) -% 23;104
105 if (ix > @as(u32, @bitCast(std.math.inf(f32))))
106 return if (common.want_float_exceptions) (x - x) / 0.0 else math.nan(f32);
107
108 ix = @as(u32, @bitCast(@as(i32, @bitCast(x * 0x1p23)) - (23 << 23)));
107 }109 }
108110
109 // argument reduction:111 const m: u32 = if (ix & 0x00800000 != 0)
110 // x = 4^e m; with integer e, and m in [1, 4)112 (ix << 7) & 0x7fffffff
111 // m: fixed point representation [2.30]113 else
112 // 2^e is the exponent part of the result.114 (ix << 8) | 0x80000000;
113 const even = (top & 1) != 0;115
114 const m = if (even) (ix << 7) & 0x7FFF_FFFF else (ix << 8) | 0x8000_0000;116 const ey = ((ix >> 1) + (0x3f800000 >> 1)) & 0x7f800000;
115 top = (top +% 0x7F) >> 1;117 // const ey = ((ix + 0x3f800000) & 0xff000000) >> 1;
118
119 const three = 0xc0000000;
120 const i = (ix >> 17) & 0x7f;
121 var r = @as(u32, rsqrt_tab[i]) << 16;
116122
117 // approximate r ~ 1/sqrt(m) and s ~ sqrt(m) when m in [1,4)
118 // the fixed point representations are
119 // m: 2.30 r: 0.32, s: 2.30, d: 2.30, u: 2.30, three: 2.30
120 const three: u32 = 0xC000_0000;
121 var i: usize = @intCast((ix >> 17) & 0x3F);
122 if (even) i += 64;
123 var r = @as(u32, @intCast(__rsqrt_tab[i])) << 16;
124 // |r*sqrt(m) - 1| < 0x1p-8
125 var s = mul32(m, r);123 var s = mul32(m, r);
126 // |s/sqrt(m) - 1| < 0x1p-8
127 var d = mul32(s, r);124 var d = mul32(s, r);
128 var u = three - d;125 var u = three - d;
129 r = mul32(r, u) << 1;126 r = mul32(r, u) << 1;
130 // |r*sqrt(m) - 1| < 0x1.7bp-16
131 s = mul32(s, u) << 1;127 s = mul32(s, u) << 1;
132 // |s/sqrt(m) - 1| < 0x1.7bp-16
133 d = mul32(s, r);128 d = mul32(s, r);
134 u = three - d;129 u = three - d;
135 s = mul32(s, u); // repr: 3.29130 s = mul32(s, u);
136 // -0x1.03p-28 < s/sqrt(m) - 1 < 0x1.fp-31131 s = (s - 1) >> 6;
137 s = (s - 1) >> 6; // repr: 9.23
138 // s < sqrt(m) < s + 0x1.08p-23
139132
140 // compute nearest rounded result:
141 // the nearest result to 23 bits is either s or s+0x1p-23,
142 // we can decide by comparing (2^23 s + 0.5)^2 to 2^46 m.
143 const d0 = (m << 16) -% s *% s;133 const d0 = (m << 16) -% s *% s;
144 const d1 = s -% d0;134 const d1 = s -% d0;
145 const d2 = d1 +% s +% 1;135 const d2 = d1 +% s +% 1;
146 s += d1 >> 31;136 const y: f32 = @bitCast(((s + (d1 >> 31)) & 0x007fffff) | ey);
147 s &= 0x007F_FFFF;
148 s |= top << 23;
149 const y: f32 = @bitCast(s);
150137
151 // handle rounding modes and inexact exception:138 const tiny: u32 = if (d2 == 0) blk: {
152 // only (s+1)^2 == 2^16 m case is exact otherwise139 @branchHint(.unlikely);
153 // add a tiny value to cause the fenv effects.140 break :blk 0;
154 if (d2 != 0) {141 } else 0x01000000;
155 @branchHint(.likely);142 const t: f32 = @bitCast(tiny | ((d1 ^ d2) & 0x80000000));
156 var tiny: u32 = 0x0100_0000;
157 tiny |= (d1 ^ d2) & 0x8000_0000;
158 const t: f32 = @bitCast(tiny);
159 return y + t;
160 }
161143
162 return y;144 return y + t;
163}145}
164146
165pub fn sqrt(x: f64) callconv(.c) f64 {147pub fn sqrt(x: f64) callconv(.c) f64 {
...@@ -172,7 +154,7 @@ pub fn sqrt(x: f64) callconv(.c) f64 {...@@ -172,7 +154,7 @@ pub fn sqrt(x: f64) callconv(.c) f64 {
172 // x < 0x1p-1022 or inf or nan.154 // x < 0x1p-1022 or inf or nan.
173 if (ix & 0x7FFF_FFFF_FFFF_FFFF == 0) return x;155 if (ix & 0x7FFF_FFFF_FFFF_FFFF == 0) return x;
174 if (ix == 0x7FF0_0000_0000_0000) return x;156 if (ix == 0x7FF0_0000_0000_0000) return x;
175 if (ix > 0x7FF0_0000_0000_0000) return math.nan(f64);157 if (ix > 0x7FF0_0000_0000_0000) return if (common.want_float_exceptions) (x - x) / 0.0 else math.nan(f64);
176 // x is subnormal, normalize it.158 // x is subnormal, normalize it.
177 ix = @bitCast(x * 0x1p52);159 ix = @bitCast(x * 0x1p52);
178 top = (ix >> 52) -% 52;160 top = (ix >> 52) -% 52;
...@@ -248,7 +230,7 @@ pub fn sqrt(x: f64) callconv(.c) f64 {...@@ -248,7 +230,7 @@ pub fn sqrt(x: f64) callconv(.c) f64 {
248 var d: struct { u32, u64 } = undefined;230 var d: struct { u32, u64 } = undefined;
249 var u: struct { u32, u64 } = undefined;231 var u: struct { u32, u64 } = undefined;
250 const i: usize = @intCast((ix >> 46) & 0x7F);232 const i: usize = @intCast((ix >> 46) & 0x7F);
251 r[0] = @intCast(__rsqrt_tab[i]);233 r[0] = @intCast(rsqrt_tab[i]);
252 r[0] <<= 16;234 r[0] <<= 16;
253 // |r sqrt(m) - 1| < 0x1.fdp-9235 // |r sqrt(m) - 1| < 0x1.fdp-9
254 s[0] = mul32(@intCast(m >> 32), r[0]);236 s[0] = mul32(@intCast(m >> 32), r[0]);
...@@ -309,7 +291,7 @@ pub fn __sqrtx(x: f80) callconv(.c) f80 {...@@ -309,7 +291,7 @@ pub fn __sqrtx(x: f80) callconv(.c) f80 {
309 // x < 0x1p-16382 or inf or nan.291 // x < 0x1p-16382 or inf or nan.
310 if (ix & 0x7FFF_FFFF_FFFF_FFFF_FFFF == 0) return x;292 if (ix & 0x7FFF_FFFF_FFFF_FFFF_FFFF == 0) return x;
311 if (ix == 0x7FFF_8000_0000_0000_0000) return x;293 if (ix == 0x7FFF_8000_0000_0000_0000) return x;
312 if (ix > 0x7FFF_8000_0000_0000_0000) return math.nan(f80);294 if (ix > 0x7FFF_8000_0000_0000_0000) return if (common.want_float_exceptions) (x - x) / 0.0 else math.nan(f80);
313 // x is subnormal, normalize it.295 // x is subnormal, normalize it.
314 ix = @bitCast(x * 0x1p63);296 ix = @bitCast(x * 0x1p63);
315 top = (ix >> 64) -% 63;297 top = (ix >> 64) -% 63;
...@@ -341,7 +323,7 @@ pub fn __sqrtx(x: f80) callconv(.c) f80 {...@@ -341,7 +323,7 @@ pub fn __sqrtx(x: f80) callconv(.c) f80 {
341 var u: struct { u32, u64, u80 } = undefined;323 var u: struct { u32, u64, u80 } = undefined;
342 var i: usize = @intCast((ix >> 57) & 0x3F);324 var i: usize = @intCast((ix >> 57) & 0x3F);
343 if (even) i += 64;325 if (even) i += 64;
344 r[0] = @intCast(__rsqrt_tab[i]);326 r[0] = @intCast(rsqrt_tab[i]);
345 r[0] <<= 16;327 r[0] <<= 16;
346 // |r sqrt(m) - 1| < 0x1p-8328 // |r sqrt(m) - 1| < 0x1p-8
347 s[0] = mul32(@intCast(m >> 48), r[0]);329 s[0] = mul32(@intCast(m >> 48), r[0]);
...@@ -437,7 +419,7 @@ pub fn sqrtq(x: f128) callconv(.c) f128 {...@@ -437,7 +419,7 @@ pub fn sqrtq(x: f128) callconv(.c) f128 {
437 var d: struct { u32, u64, u128 } = undefined;419 var d: struct { u32, u64, u128 } = undefined;
438 var u: struct { u32, u64, u128 } = undefined;420 var u: struct { u32, u64, u128 } = undefined;
439 const i: usize = @intCast((ix >> 106) & 0x7F);421 const i: usize = @intCast((ix >> 106) & 0x7F);
440 r[0] = @intCast(__rsqrt_tab[i]);422 r[0] = @intCast(rsqrt_tab[i]);
441 r[0] <<= 16;423 r[0] <<= 16;
442 // |r sqrt(m) - 1| < 0x1p-8424 // |r sqrt(m) - 1| < 0x1p-8
443 s[0] = mul32(@intCast(m >> 96), r[0]);425 s[0] = mul32(@intCast(m >> 96), r[0]);
...@@ -507,7 +489,7 @@ pub fn sqrtl(x: c_longdouble) callconv(.c) c_longdouble {...@@ -507,7 +489,7 @@ pub fn sqrtl(x: c_longdouble) callconv(.c) c_longdouble {
507 }489 }
508}490}
509491
510const __rsqrt_tab: [128]u16 = .{492const rsqrt_tab: [128]u16 = .{
511 0xB451, 0xB2F0, 0xB196, 0xB044, 0xAEF9, 0xADB6, 0xAC79, 0xAB43,493 0xB451, 0xB2F0, 0xB196, 0xB044, 0xAEF9, 0xADB6, 0xAC79, 0xAB43,
512 0xAA14, 0xA8EB, 0xA7C8, 0xA6AA, 0xA592, 0xA480, 0xA373, 0xA26B,494 0xAA14, 0xA8EB, 0xA7C8, 0xA6AA, 0xA592, 0xA480, 0xA373, 0xA26B,
513 0xA168, 0xA06A, 0x9F70, 0x9E7B, 0x9D8A, 0x9C9D, 0x9BB5, 0x9AD1,495 0xA168, 0xA06A, 0x9F70, 0x9E7B, 0x9D8A, 0x9C9D, 0x9BB5, 0x9AD1,
...@@ -527,15 +509,15 @@ const __rsqrt_tab: [128]u16 = .{...@@ -527,15 +509,15 @@ const __rsqrt_tab: [128]u16 = .{
527};509};
528510
529inline fn mul16(a: u16, b: u16) u16 {511inline fn mul16(a: u16, b: u16) u16 {
530 return @intCast(@as(u32, @intCast(a)) * @as(u32, @intCast(b)) >> 16);512 return @intCast(@as(u32, a) * b >> 16);
531}513}
532514
533inline fn mul32(a: u32, b: u32) u32 {515inline fn mul32(a: u32, b: u32) u32 {
534 return @intCast(@as(u64, @intCast(a)) * @as(u64, @intCast(b)) >> 32);516 return @intCast(@as(u64, a) * b >> 32);
535}517}
536518
537inline fn mul64(a: u64, b: u64) u64 {519inline fn mul64(a: u64, b: u64) u64 {
538 return @intCast(@as(u128, @intCast(a)) * @as(u128, @intCast(b)) >> 64);520 return @intCast(@as(u128, a) * b >> 64);
539}521}
540522
541inline fn mul80(a: u80, b: u80) u80 {523inline fn mul80(a: u80, b: u80) u80 {
lib/libc/mingw/math/sqrt.def.h deleted-92
...@@ -1,92 +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
48__FLT_TYPE
49__FLT_ABI (sqrt) (__FLT_TYPE x)
50{
51 __FLT_TYPE res = __FLT_CST (0.0);
52 int x_class = fpclassify (x);
53 if (x_class == FP_NAN || signbit (x))
54 {
55 if (x_class == FP_ZERO)
56 return __FLT_CST (-0.0);
57
58 if (x_class == FP_NAN)
59 {
60 __FLT_RPT_DOMAIN ("sqrt", x, 0.0, x);
61 return x;
62 }
63
64 res = -__FLT_NAN;
65 __FLT_RPT_DOMAIN ("sqrt", x, 0.0, res);
66 return res;
67 }
68 else if (x_class == FP_ZERO)
69 return __FLT_CST (0.0);
70 else if (x_class == FP_INFINITE)
71 return __FLT_HUGE_VAL;
72 else if (x == __FLT_CST (1.0))
73 return __FLT_CST (1.0);
74#if defined(__arm__) || defined(_ARM_)
75#if _NEW_COMPLEX_FLOAT
76 asm volatile ("fsqrts %[dst], %[src];\n" : [dst] "=t" (res) : [src] "t" (x));
77#else
78 asm volatile ("fsqrtd %[dst], %[src];\n" : [dst] "=w" (res) : [src] "w" (x));
79#endif
80#elif defined(__aarch64__) || defined(_ARM64_) || defined(__arm64ec__) || defined(_ARM64EC_)
81#if _NEW_COMPLEX_FLOAT
82 asm volatile ("fsqrt %s[dst], %s[src]\n" : [dst] "=w" (res) : [src] "w" (x));
83#else
84 asm volatile ("fsqrt %d[dst], %d[src]\n" : [dst] "=w" (res) : [src] "w" (x));
85#endif
86#elif defined(_X86_) || defined(__i386__) || defined(_AMD64_) || defined(__x86_64__)
87 asm volatile ("fsqrt" : "=t" (res) : "0" (x));
88#else
89#error Not supported on your platform yet
90#endif
91 return res;
92}
lib/libc/mingw/math/sqrtf.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_FLOAT 1
46#include "sqrt.def.h"
lib/libc/mingw/math/sqrtl.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 "sqrt.def.h"
lib/libc/musl/src/math/aarch64/sqrt.c deleted-7
...@@ -1,7 +0,0 @@
1#include <math.h>
2
3double sqrt(double x)
4{
5 __asm__ ("fsqrt %d0, %d1" : "=w"(x) : "w"(x));
6 return x;
7}
lib/libc/musl/src/math/aarch64/sqrtf.c deleted-7
...@@ -1,7 +0,0 @@
1#include <math.h>
2
3float sqrtf(float x)
4{
5 __asm__ ("fsqrt %s0, %s1" : "=w"(x) : "w"(x));
6 return x;
7}
lib/libc/musl/src/math/arm/sqrt.c deleted-15
...@@ -1,15 +0,0 @@
1#include <math.h>
2
3#if (__ARM_PCS_VFP || (__VFP_FP__ && !__SOFTFP__)) && (__ARM_FP&8)
4
5double sqrt(double x)
6{
7 __asm__ ("vsqrt.f64 %P0, %P1" : "=w"(x) : "w"(x));
8 return x;
9}
10
11#else
12
13#include "../sqrt.c"
14
15#endif
lib/libc/musl/src/math/arm/sqrtf.c deleted-15
...@@ -1,15 +0,0 @@
1#include <math.h>
2
3#if (__ARM_PCS_VFP || (__VFP_FP__ && !__SOFTFP__)) && !BROKEN_VFP_ASM
4
5float sqrtf(float x)
6{
7 __asm__ ("vsqrt.f32 %0, %1" : "=t"(x) : "t"(x));
8 return x;
9}
10
11#else
12
13#include "../sqrtf.c"
14
15#endif
lib/libc/musl/src/math/i386/sqrt.c deleted-15
...@@ -1,15 +0,0 @@
1#include "libm.h"
2
3double sqrt(double x)
4{
5 union ldshape ux;
6 unsigned fpsr;
7 __asm__ ("fsqrt; fnstsw %%ax": "=t"(ux.f), "=a"(fpsr) : "0"(x));
8 if ((ux.i.m & 0x7ff) != 0x400)
9 return (double)ux.f;
10 /* Rounding to double would have encountered an exact halfway case.
11 Adjust mantissa downwards if fsqrt rounded up, else upwards.
12 (result of fsqrt could not have been exact) */
13 ux.i.m ^= (fpsr & 0x200) + 0x300;
14 return (double)ux.f;
15}
lib/libc/musl/src/math/i386/sqrtf.c deleted-12
...@@ -1,12 +0,0 @@
1#include <math.h>
2
3float sqrtf(float x)
4{
5 long double t;
6 /* The long double result has sufficient precision so that
7 * second rounding to float still keeps the returned value
8 * correctly rounded, see Pierre Roux, "Innocuous Double
9 * Rounding of Basic Arithmetic Operations". */
10 __asm__ ("fsqrt" : "=t"(t) : "0"(x));
11 return (float)t;
12}
lib/libc/musl/src/math/i386/sqrtl.c deleted-7
...@@ -1,7 +0,0 @@
1#include <math.h>
2
3long double sqrtl(long double x)
4{
5 __asm__ ("fsqrt" : "+t"(x));
6 return x;
7}
lib/libc/musl/src/math/m68k/sqrtl.c deleted-15
...@@ -1,15 +0,0 @@
1#include <math.h>
2
3#if __HAVE_68881__
4
5long double sqrtl(long double x)
6{
7 __asm__ ("fsqrt.x %1,%0" : "=f"(x) : "fm"(x));
8 return x;
9}
10
11#else
12
13#include "../sqrtl.c"
14
15#endif
lib/libc/musl/src/math/mips/sqrt.c deleted-16
...@@ -1,16 +0,0 @@
1#if !defined(__mips_soft_float) && __mips >= 3
2
3#include <math.h>
4
5double sqrt(double x)
6{
7 double r;
8 __asm__("sqrt.d %0,%1" : "=f"(r) : "f"(x));
9 return r;
10}
11
12#else
13
14#include "../sqrt.c"
15
16#endif
lib/libc/musl/src/math/mips/sqrtf.c deleted-16
...@@ -1,16 +0,0 @@
1#if !defined(__mips_soft_float) && __mips >= 2
2
3#include <math.h>
4
5float sqrtf(float x)
6{
7 float r;
8 __asm__("sqrt.s %0,%1" : "=f"(r) : "f"(x));
9 return r;
10}
11
12#else
13
14#include "../sqrtf.c"
15
16#endif
lib/libc/musl/src/math/powerpc/sqrt.c deleted-15
...@@ -1,15 +0,0 @@
1#include <math.h>
2
3#if !defined _SOFT_FLOAT && defined _ARCH_PPCSQ
4
5double sqrt(double x)
6{
7 __asm__ ("fsqrt %0, %1\n" : "=d" (x) : "d" (x));
8 return x;
9}
10
11#else
12
13#include "../sqrt.c"
14
15#endif
lib/libc/musl/src/math/powerpc/sqrtf.c deleted-15
...@@ -1,15 +0,0 @@
1#include <math.h>
2
3#if !defined _SOFT_FLOAT && defined _ARCH_PPCSQ
4
5float sqrtf(float x)
6{
7 __asm__ ("fsqrts %0, %1\n" : "=f" (x) : "f" (x));
8 return x;
9}
10
11#else
12
13#include "../sqrtf.c"
14
15#endif
lib/libc/musl/src/math/powerpc64/sqrt.c deleted-7
...@@ -1,7 +0,0 @@
1#include <math.h>
2
3double sqrt(double x)
4{
5 __asm__ ("fsqrt %0, %1" : "=d"(x) : "d"(x));
6 return x;
7}
lib/libc/musl/src/math/powerpc64/sqrtf.c deleted-7
...@@ -1,7 +0,0 @@
1#include <math.h>
2
3float sqrtf(float x)
4{
5 __asm__ ("fsqrts %0, %1" : "=f"(x) : "f"(x));
6 return x;
7}
lib/libc/musl/src/math/riscv32/sqrt.c deleted-15
...@@ -1,15 +0,0 @@
1#include <math.h>
2
3#if __riscv_flen >= 64
4
5double sqrt(double x)
6{
7 __asm__ ("fsqrt.d %0, %1" : "=f"(x) : "f"(x));
8 return x;
9}
10
11#else
12
13#include "../sqrt.c"
14
15#endif
lib/libc/musl/src/math/riscv32/sqrtf.c deleted-15
...@@ -1,15 +0,0 @@
1#include <math.h>
2
3#if __riscv_flen >= 32
4
5float sqrtf(float x)
6{
7 __asm__ ("fsqrt.s %0, %1" : "=f"(x) : "f"(x));
8 return x;
9}
10
11#else
12
13#include "../sqrtf.c"
14
15#endif
lib/libc/musl/src/math/riscv64/sqrt.c deleted-15
...@@ -1,15 +0,0 @@
1#include <math.h>
2
3#if __riscv_flen >= 64
4
5double sqrt(double x)
6{
7 __asm__ ("fsqrt.d %0, %1" : "=f"(x) : "f"(x));
8 return x;
9}
10
11#else
12
13#include "../sqrt.c"
14
15#endif
lib/libc/musl/src/math/riscv64/sqrtf.c deleted-15
...@@ -1,15 +0,0 @@
1#include <math.h>
2
3#if __riscv_flen >= 32
4
5float sqrtf(float x)
6{
7 __asm__ ("fsqrt.s %0, %1" : "=f"(x) : "f"(x));
8 return x;
9}
10
11#else
12
13#include "../sqrtf.c"
14
15#endif
lib/libc/musl/src/math/s390x/sqrt.c deleted-15
...@@ -1,15 +0,0 @@
1#include <math.h>
2
3#if defined(__HTM__) || __ARCH__ >= 9
4
5double sqrt(double x)
6{
7 __asm__ ("sqdbr %0, %1" : "=f"(x) : "f"(x));
8 return x;
9}
10
11#else
12
13#include "../sqrt.c"
14
15#endif
lib/libc/musl/src/math/s390x/sqrtf.c deleted-15
...@@ -1,15 +0,0 @@
1#include <math.h>
2
3#if defined(__HTM__) || __ARCH__ >= 9
4
5float sqrtf(float x)
6{
7 __asm__ ("sqebr %0, %1" : "=f"(x) : "f"(x));
8 return x;
9}
10
11#else
12
13#include "../sqrtf.c"
14
15#endif
lib/libc/musl/src/math/s390x/sqrtl.c deleted-15
...@@ -1,15 +0,0 @@
1#include <math.h>
2
3#if defined(__HTM__) || __ARCH__ >= 9
4
5long double sqrtl(long double x)
6{
7 __asm__ ("sqxbr %0, %1" : "=f"(x) : "f"(x));
8 return x;
9}
10
11#else
12
13#include "../sqrtl.c"
14
15#endif
lib/libc/musl/src/math/sqrt.c deleted-158
...@@ -1,158 +0,0 @@
1#include <stdint.h>
2#include <math.h>
3#include "libm.h"
4#include "sqrt_data.h"
5
6#define FENV_SUPPORT 1
7
8/* returns a*b*2^-32 - e, with error 0 <= e < 1. */
9static inline uint32_t mul32(uint32_t a, uint32_t b)
10{
11 return (uint64_t)a*b >> 32;
12}
13
14/* returns a*b*2^-64 - e, with error 0 <= e < 3. */
15static inline uint64_t mul64(uint64_t a, uint64_t b)
16{
17 uint64_t ahi = a>>32;
18 uint64_t alo = a&0xffffffff;
19 uint64_t bhi = b>>32;
20 uint64_t blo = b&0xffffffff;
21 return ahi*bhi + (ahi*blo >> 32) + (alo*bhi >> 32);
22}
23
24double sqrt(double x)
25{
26 uint64_t ix, top, m;
27
28 /* special case handling. */
29 ix = asuint64(x);
30 top = ix >> 52;
31 if (predict_false(top - 0x001 >= 0x7ff - 0x001)) {
32 /* x < 0x1p-1022 or inf or nan. */
33 if (ix * 2 == 0)
34 return x;
35 if (ix == 0x7ff0000000000000)
36 return x;
37 if (ix > 0x7ff0000000000000)
38 return __math_invalid(x);
39 /* x is subnormal, normalize it. */
40 ix = asuint64(x * 0x1p52);
41 top = ix >> 52;
42 top -= 52;
43 }
44
45 /* argument reduction:
46 x = 4^e m; with integer e, and m in [1, 4)
47 m: fixed point representation [2.62]
48 2^e is the exponent part of the result. */
49 int even = top & 1;
50 m = (ix << 11) | 0x8000000000000000;
51 if (even) m >>= 1;
52 top = (top + 0x3ff) >> 1;
53
54 /* approximate r ~ 1/sqrt(m) and s ~ sqrt(m) when m in [1,4)
55
56 initial estimate:
57 7bit table lookup (1bit exponent and 6bit significand).
58
59 iterative approximation:
60 using 2 goldschmidt iterations with 32bit int arithmetics
61 and a final iteration with 64bit int arithmetics.
62
63 details:
64
65 the relative error (e = r0 sqrt(m)-1) of a linear estimate
66 (r0 = a m + b) is |e| < 0.085955 ~ 0x1.6p-4 at best,
67 a table lookup is faster and needs one less iteration
68 6 bit lookup table (128b) gives |e| < 0x1.f9p-8
69 7 bit lookup table (256b) gives |e| < 0x1.fdp-9
70 for single and double prec 6bit is enough but for quad
71 prec 7bit is needed (or modified iterations). to avoid
72 one more iteration >=13bit table would be needed (16k).
73
74 a newton-raphson iteration for r is
75 w = r*r
76 u = 3 - m*w
77 r = r*u/2
78 can use a goldschmidt iteration for s at the end or
79 s = m*r
80
81 first goldschmidt iteration is
82 s = m*r
83 u = 3 - s*r
84 r = r*u/2
85 s = s*u/2
86 next goldschmidt iteration is
87 u = 3 - s*r
88 r = r*u/2
89 s = s*u/2
90 and at the end r is not computed only s.
91
92 they use the same amount of operations and converge at the
93 same quadratic rate, i.e. if
94 r1 sqrt(m) - 1 = e, then
95 r2 sqrt(m) - 1 = -3/2 e^2 - 1/2 e^3
96 the advantage of goldschmidt is that the mul for s and r
97 are independent (computed in parallel), however it is not
98 "self synchronizing": it only uses the input m in the
99 first iteration so rounding errors accumulate. at the end
100 or when switching to larger precision arithmetics rounding
101 errors dominate so the first iteration should be used.
102
103 the fixed point representations are
104 m: 2.30 r: 0.32, s: 2.30, d: 2.30, u: 2.30, three: 2.30
105 and after switching to 64 bit
106 m: 2.62 r: 0.64, s: 2.62, d: 2.62, u: 2.62, three: 2.62 */
107
108 static const uint64_t three = 0xc0000000;
109 uint64_t r, s, d, u, i;
110
111 i = (ix >> 46) % 128;
112 r = (uint32_t)__rsqrt_tab[i] << 16;
113 /* |r sqrt(m) - 1| < 0x1.fdp-9 */
114 s = mul32(m>>32, r);
115 /* |s/sqrt(m) - 1| < 0x1.fdp-9 */
116 d = mul32(s, r);
117 u = three - d;
118 r = mul32(r, u) << 1;
119 /* |r sqrt(m) - 1| < 0x1.7bp-16 */
120 s = mul32(s, u) << 1;
121 /* |s/sqrt(m) - 1| < 0x1.7bp-16 */
122 d = mul32(s, r);
123 u = three - d;
124 r = mul32(r, u) << 1;
125 /* |r sqrt(m) - 1| < 0x1.3704p-29 (measured worst-case) */
126 r = r << 32;
127 s = mul64(m, r);
128 d = mul64(s, r);
129 u = (three<<32) - d;
130 s = mul64(s, u); /* repr: 3.61 */
131 /* -0x1p-57 < s - sqrt(m) < 0x1.8001p-61 */
132 s = (s - 2) >> 9; /* repr: 12.52 */
133 /* -0x1.09p-52 < s - sqrt(m) < -0x1.fffcp-63 */
134
135 /* s < sqrt(m) < s + 0x1.09p-52,
136 compute nearest rounded result:
137 the nearest result to 52 bits is either s or s+0x1p-52,
138 we can decide by comparing (2^52 s + 0.5)^2 to 2^104 m. */
139 uint64_t d0, d1, d2;
140 double y, t;
141 d0 = (m << 42) - s*s;
142 d1 = s - d0;
143 d2 = d1 + s + 1;
144 s += d1 >> 63;
145 s &= 0x000fffffffffffff;
146 s |= top << 52;
147 y = asdouble(s);
148 if (FENV_SUPPORT) {
149 /* handle rounding modes and inexact exception:
150 only (s+1)^2 == 2^42 m case is exact otherwise
151 add a tiny value to cause the fenv effects. */
152 uint64_t tiny = predict_false(d2==0) ? 0 : 0x0010000000000000;
153 tiny |= (d1^d2) & 0x8000000000000000;
154 t = asdouble(tiny);
155 y = eval_as_double(y + t);
156 }
157 return y;
158}
lib/libc/musl/src/math/sqrt_data.c deleted-19
...@@ -1,19 +0,0 @@
1#include "sqrt_data.h"
2const uint16_t __rsqrt_tab[128] = {
30xb451,0xb2f0,0xb196,0xb044,0xaef9,0xadb6,0xac79,0xab43,
40xaa14,0xa8eb,0xa7c8,0xa6aa,0xa592,0xa480,0xa373,0xa26b,
50xa168,0xa06a,0x9f70,0x9e7b,0x9d8a,0x9c9d,0x9bb5,0x9ad1,
60x99f0,0x9913,0x983a,0x9765,0x9693,0x95c4,0x94f8,0x9430,
70x936b,0x92a9,0x91ea,0x912e,0x9075,0x8fbe,0x8f0a,0x8e59,
80x8daa,0x8cfe,0x8c54,0x8bac,0x8b07,0x8a64,0x89c4,0x8925,
90x8889,0x87ee,0x8756,0x86c0,0x862b,0x8599,0x8508,0x8479,
100x83ec,0x8361,0x82d8,0x8250,0x81c9,0x8145,0x80c2,0x8040,
110xff02,0xfd0e,0xfb25,0xf947,0xf773,0xf5aa,0xf3ea,0xf234,
120xf087,0xeee3,0xed47,0xebb3,0xea27,0xe8a3,0xe727,0xe5b2,
130xe443,0xe2dc,0xe17a,0xe020,0xdecb,0xdd7d,0xdc34,0xdaf1,
140xd9b3,0xd87b,0xd748,0xd61a,0xd4f1,0xd3cd,0xd2ad,0xd192,
150xd07b,0xcf69,0xce5b,0xcd51,0xcc4a,0xcb48,0xca4a,0xc94f,
160xc858,0xc764,0xc674,0xc587,0xc49d,0xc3b7,0xc2d4,0xc1f4,
170xc116,0xc03c,0xbf65,0xbe90,0xbdbe,0xbcef,0xbc23,0xbb59,
180xba91,0xb9cc,0xb90a,0xb84a,0xb78c,0xb6d0,0xb617,0xb560,
19};
lib/libc/musl/src/math/sqrt_data.h deleted-13
...@@ -1,13 +0,0 @@
1#ifndef _SQRT_DATA_H
2#define _SQRT_DATA_H
3
4#include <features.h>
5#include <stdint.h>
6
7/* if x in [1,2): i = (int)(64*x);
8 if x in [2,4): i = (int)(32*x-64);
9 __rsqrt_tab[i]*2^-16 is estimating 1/sqrt(x) with small relative error:
10 |__rsqrt_tab[i]*0x1p-16*sqrt(x) - 1| < -0x1.fdp-9 < 2^-8 */
11extern hidden const uint16_t __rsqrt_tab[128];
12
13#endif
lib/libc/musl/src/math/sqrtf.c deleted-83
...@@ -1,83 +0,0 @@
1#include <stdint.h>
2#include <math.h>
3#include "libm.h"
4#include "sqrt_data.h"
5
6#define FENV_SUPPORT 1
7
8static inline uint32_t mul32(uint32_t a, uint32_t b)
9{
10 return (uint64_t)a*b >> 32;
11}
12
13/* see sqrt.c for more detailed comments. */
14
15float sqrtf(float x)
16{
17 uint32_t ix, m, m1, m0, even, ey;
18
19 ix = asuint(x);
20 if (predict_false(ix - 0x00800000 >= 0x7f800000 - 0x00800000)) {
21 /* x < 0x1p-126 or inf or nan. */
22 if (ix * 2 == 0)
23 return x;
24 if (ix == 0x7f800000)
25 return x;
26 if (ix > 0x7f800000)
27 return __math_invalidf(x);
28 /* x is subnormal, normalize it. */
29 ix = asuint(x * 0x1p23f);
30 ix -= 23 << 23;
31 }
32
33 /* x = 4^e m; with int e and m in [1, 4). */
34 even = ix & 0x00800000;
35 m1 = (ix << 8) | 0x80000000;
36 m0 = (ix << 7) & 0x7fffffff;
37 m = even ? m0 : m1;
38
39 /* 2^e is the exponent part of the return value. */
40 ey = ix >> 1;
41 ey += 0x3f800000 >> 1;
42 ey &= 0x7f800000;
43
44 /* compute r ~ 1/sqrt(m), s ~ sqrt(m) with 2 goldschmidt iterations. */
45 static const uint32_t three = 0xc0000000;
46 uint32_t r, s, d, u, i;
47 i = (ix >> 17) % 128;
48 r = (uint32_t)__rsqrt_tab[i] << 16;
49 /* |r*sqrt(m) - 1| < 0x1p-8 */
50 s = mul32(m, r);
51 /* |s/sqrt(m) - 1| < 0x1p-8 */
52 d = mul32(s, r);
53 u = three - d;
54 r = mul32(r, u) << 1;
55 /* |r*sqrt(m) - 1| < 0x1.7bp-16 */
56 s = mul32(s, u) << 1;
57 /* |s/sqrt(m) - 1| < 0x1.7bp-16 */
58 d = mul32(s, r);
59 u = three - d;
60 s = mul32(s, u);
61 /* -0x1.03p-28 < s/sqrt(m) - 1 < 0x1.fp-31 */
62 s = (s - 1)>>6;
63 /* s < sqrt(m) < s + 0x1.08p-23 */
64
65 /* compute nearest rounded result. */
66 uint32_t d0, d1, d2;
67 float y, t;
68 d0 = (m << 16) - s*s;
69 d1 = s - d0;
70 d2 = d1 + s + 1;
71 s += d1 >> 31;
72 s &= 0x007fffff;
73 s |= ey;
74 y = asfloat(s);
75 if (FENV_SUPPORT) {
76 /* handle rounding and inexact exception. */
77 uint32_t tiny = predict_false(d2==0) ? 0 : 0x01000000;
78 tiny |= (d1^d2) & 0x80000000;
79 t = asfloat(tiny);
80 y = eval_as_float(y + t);
81 }
82 return y;
83}
lib/libc/musl/src/math/sqrtl.c deleted-259
...@@ -1,259 +0,0 @@
1#include <stdint.h>
2#include <math.h>
3#include <float.h>
4#include "libm.h"
5
6#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
7long double sqrtl(long double x)
8{
9 return sqrt(x);
10}
11#elif (LDBL_MANT_DIG == 113 || LDBL_MANT_DIG == 64) && LDBL_MAX_EXP == 16384
12#include "sqrt_data.h"
13
14#define FENV_SUPPORT 1
15
16typedef struct {
17 uint64_t hi;
18 uint64_t lo;
19} u128;
20
21/* top: 16 bit sign+exponent, x: significand. */
22static inline long double mkldbl(uint64_t top, u128 x)
23{
24 union ldshape u;
25#if LDBL_MANT_DIG == 113
26 u.i2.hi = x.hi;
27 u.i2.lo = x.lo;
28 u.i2.hi &= 0x0000ffffffffffff;
29 u.i2.hi |= top << 48;
30#elif LDBL_MANT_DIG == 64
31 u.i.se = top;
32 u.i.m = x.lo;
33 /* force the top bit on non-zero (and non-subnormal) results. */
34 if (top & 0x7fff)
35 u.i.m |= 0x8000000000000000;
36#endif
37 return u.f;
38}
39
40/* return: top 16 bit is sign+exp and following bits are the significand. */
41static inline u128 asu128(long double x)
42{
43 union ldshape u = {.f=x};
44 u128 r;
45#if LDBL_MANT_DIG == 113
46 r.hi = u.i2.hi;
47 r.lo = u.i2.lo;
48#elif LDBL_MANT_DIG == 64
49 r.lo = u.i.m<<49;
50 /* ignore the top bit: pseudo numbers are not handled. */
51 r.hi = u.i.m>>15;
52 r.hi &= 0x0000ffffffffffff;
53 r.hi |= (uint64_t)u.i.se << 48;
54#endif
55 return r;
56}
57
58/* returns a*b*2^-32 - e, with error 0 <= e < 1. */
59static inline uint32_t mul32(uint32_t a, uint32_t b)
60{
61 return (uint64_t)a*b >> 32;
62}
63
64/* returns a*b*2^-64 - e, with error 0 <= e < 3. */
65static inline uint64_t mul64(uint64_t a, uint64_t b)
66{
67 uint64_t ahi = a>>32;
68 uint64_t alo = a&0xffffffff;
69 uint64_t bhi = b>>32;
70 uint64_t blo = b&0xffffffff;
71 return ahi*bhi + (ahi*blo >> 32) + (alo*bhi >> 32);
72}
73
74static inline u128 add64(u128 a, uint64_t b)
75{
76 u128 r;
77 r.lo = a.lo + b;
78 r.hi = a.hi;
79 if (r.lo < a.lo)
80 r.hi++;
81 return r;
82}
83
84static inline u128 add128(u128 a, u128 b)
85{
86 u128 r;
87 r.lo = a.lo + b.lo;
88 r.hi = a.hi + b.hi;
89 if (r.lo < a.lo)
90 r.hi++;
91 return r;
92}
93
94static inline u128 sub64(u128 a, uint64_t b)
95{
96 u128 r;
97 r.lo = a.lo - b;
98 r.hi = a.hi;
99 if (a.lo < b)
100 r.hi--;
101 return r;
102}
103
104static inline u128 sub128(u128 a, u128 b)
105{
106 u128 r;
107 r.lo = a.lo - b.lo;
108 r.hi = a.hi - b.hi;
109 if (a.lo < b.lo)
110 r.hi--;
111 return r;
112}
113
114/* a<<n, 0 <= n <= 127 */
115static inline u128 lsh(u128 a, int n)
116{
117 if (n == 0)
118 return a;
119 if (n >= 64) {
120 a.hi = a.lo<<(n-64);
121 a.lo = 0;
122 } else {
123 a.hi = (a.hi<<n) | (a.lo>>(64-n));
124 a.lo = a.lo<<n;
125 }
126 return a;
127}
128
129/* a>>n, 0 <= n <= 127 */
130static inline u128 rsh(u128 a, int n)
131{
132 if (n == 0)
133 return a;
134 if (n >= 64) {
135 a.lo = a.hi>>(n-64);
136 a.hi = 0;
137 } else {
138 a.lo = (a.lo>>n) | (a.hi<<(64-n));
139 a.hi = a.hi>>n;
140 }
141 return a;
142}
143
144/* returns a*b exactly. */
145static inline u128 mul64_128(uint64_t a, uint64_t b)
146{
147 u128 r;
148 uint64_t ahi = a>>32;
149 uint64_t alo = a&0xffffffff;
150 uint64_t bhi = b>>32;
151 uint64_t blo = b&0xffffffff;
152 uint64_t lo1 = ((ahi*blo)&0xffffffff) + ((alo*bhi)&0xffffffff) + (alo*blo>>32);
153 uint64_t lo2 = (alo*blo)&0xffffffff;
154 r.hi = ahi*bhi + (ahi*blo>>32) + (alo*bhi>>32) + (lo1>>32);
155 r.lo = (lo1<<32) + lo2;
156 return r;
157}
158
159/* returns a*b*2^-128 - e, with error 0 <= e < 7. */
160static inline u128 mul128(u128 a, u128 b)
161{
162 u128 hi = mul64_128(a.hi, b.hi);
163 uint64_t m1 = mul64(a.hi, b.lo);
164 uint64_t m2 = mul64(a.lo, b.hi);
165 return add64(add64(hi, m1), m2);
166}
167
168/* returns a*b % 2^128. */
169static inline u128 mul128_tail(u128 a, u128 b)
170{
171 u128 lo = mul64_128(a.lo, b.lo);
172 lo.hi += a.hi*b.lo + a.lo*b.hi;
173 return lo;
174}
175
176
177/* see sqrt.c for detailed comments. */
178
179long double sqrtl(long double x)
180{
181 u128 ix, ml;
182 uint64_t top;
183
184 ix = asu128(x);
185 top = ix.hi >> 48;
186 if (predict_false(top - 0x0001 >= 0x7fff - 0x0001)) {
187 /* x < 0x1p-16382 or inf or nan. */
188 if (2*ix.hi == 0 && ix.lo == 0)
189 return x;
190 if (ix.hi == 0x7fff000000000000 && ix.lo == 0)
191 return x;
192 if (top >= 0x7fff)
193 return __math_invalidl(x);
194 /* x is subnormal, normalize it. */
195 ix = asu128(x * 0x1p112);
196 top = ix.hi >> 48;
197 top -= 112;
198 }
199
200 /* x = 4^e m; with int e and m in [1, 4) */
201 int even = top & 1;
202 ml = lsh(ix, 15);
203 ml.hi |= 0x8000000000000000;
204 if (even) ml = rsh(ml, 1);
205 top = (top + 0x3fff) >> 1;
206
207 /* r ~ 1/sqrt(m) */
208 const uint64_t three = 0xc0000000;
209 uint64_t r, s, d, u, i;
210 i = (ix.hi >> 42) % 128;
211 r = (uint32_t)__rsqrt_tab[i] << 16;
212 /* |r sqrt(m) - 1| < 0x1p-8 */
213 s = mul32(ml.hi>>32, r);
214 d = mul32(s, r);
215 u = three - d;
216 r = mul32(u, r) << 1;
217 /* |r sqrt(m) - 1| < 0x1.7bp-16, switch to 64bit */
218 r = r<<32;
219 s = mul64(ml.hi, r);
220 d = mul64(s, r);
221 u = (three<<32) - d;
222 r = mul64(u, r) << 1;
223 /* |r sqrt(m) - 1| < 0x1.a5p-31 */
224 s = mul64(u, s) << 1;
225 d = mul64(s, r);
226 u = (three<<32) - d;
227 r = mul64(u, r) << 1;
228 /* |r sqrt(m) - 1| < 0x1.c001p-59, switch to 128bit */
229
230 const u128 threel = {.hi=three<<32, .lo=0};
231 u128 rl, sl, dl, ul;
232 rl.hi = r;
233 rl.lo = 0;
234 sl = mul128(ml, rl);
235 dl = mul128(sl, rl);
236 ul = sub128(threel, dl);
237 sl = mul128(ul, sl); /* repr: 3.125 */
238 /* -0x1p-116 < s - sqrt(m) < 0x3.8001p-125 */
239 sl = rsh(sub64(sl, 4), 125-(LDBL_MANT_DIG-1));
240 /* s < sqrt(m) < s + 1 ULP + tiny */
241
242 long double y;
243 u128 d2, d1, d0;
244 d0 = sub128(lsh(ml, 2*(LDBL_MANT_DIG-1)-126), mul128_tail(sl,sl));
245 d1 = sub128(sl, d0);
246 d2 = add128(add64(sl, 1), d1);
247 sl = add64(sl, d1.hi >> 63);
248 y = mkldbl(top, sl);
249 if (FENV_SUPPORT) {
250 /* handle rounding modes and inexact exception. */
251 top = predict_false((d2.hi|d2.lo)==0) ? 0 : 1;
252 top |= ((d1.hi^d2.hi)&0x8000000000000000) >> 48;
253 y += mkldbl(top, (u128){0});
254 }
255 return y;
256}
257#else
258#error unsupported long double format
259#endif
lib/libc/musl/src/math/x32/sqrt.s deleted-4
...@@ -1,4 +0,0 @@
1.global sqrt
2.type sqrt,@function
3sqrt: sqrtsd %xmm0, %xmm0
4 ret
lib/libc/musl/src/math/x32/sqrtf.s deleted-4
...@@ -1,4 +0,0 @@
1.global sqrtf
2.type sqrtf,@function
3sqrtf: sqrtss %xmm0, %xmm0
4 ret
lib/libc/musl/src/math/x32/sqrtl.s deleted-5
...@@ -1,5 +0,0 @@
1.global sqrtl
2.type sqrtl,@function
3sqrtl: fldt 8(%esp)
4 fsqrt
5 ret
lib/libc/musl/src/math/x86_64/sqrt.c deleted-7
...@@ -1,7 +0,0 @@
1#include <math.h>
2
3double sqrt(double x)
4{
5 __asm__ ("sqrtsd %1, %0" : "=x"(x) : "x"(x));
6 return x;
7}
lib/libc/musl/src/math/x86_64/sqrtf.c deleted-7
...@@ -1,7 +0,0 @@
1#include <math.h>
2
3float sqrtf(float x)
4{
5 __asm__ ("sqrtss %1, %0" : "=x"(x) : "x"(x));
6 return x;
7}
lib/libc/musl/src/math/x86_64/sqrtl.c deleted-7
...@@ -1,7 +0,0 @@
1#include <math.h>
2
3long double sqrtl(long double x)
4{
5 __asm__ ("fsqrt" : "+t"(x));
6 return x;
7}
src/libs/mingw.zig-2
...@@ -630,7 +630,6 @@ const mingw32_generic_src = [_][]const u8{...@@ -630,7 +630,6 @@ const mingw32_generic_src = [_][]const u8{
630 "math" ++ path.sep_str ++ "signbitl.c",630 "math" ++ path.sep_str ++ "signbitl.c",
631 "math" ++ path.sep_str ++ "signgam.c",631 "math" ++ path.sep_str ++ "signgam.c",
632 "math" ++ path.sep_str ++ "sinhl.c",632 "math" ++ path.sep_str ++ "sinhl.c",
633 "math" ++ path.sep_str ++ "sqrtl.c",
634 "math" ++ path.sep_str ++ "tanhl.c",633 "math" ++ path.sep_str ++ "tanhl.c",
635 "misc" ++ path.sep_str ++ "alarm.c",634 "misc" ++ path.sep_str ++ "alarm.c",
636 "misc" ++ path.sep_str ++ "btowc.c",635 "misc" ++ path.sep_str ++ "btowc.c",
...@@ -995,7 +994,6 @@ const mingw32_x86_32_src = [_][]const u8{...@@ -995,7 +994,6 @@ const mingw32_x86_32_src = [_][]const u8{
995 "math" ++ path.sep_str ++ "modff.c",994 "math" ++ path.sep_str ++ "modff.c",
996 "math" ++ path.sep_str ++ "powf.c",995 "math" ++ path.sep_str ++ "powf.c",
997 "math" ++ path.sep_str ++ "sinhf.c",996 "math" ++ path.sep_str ++ "sinhf.c",
998 "math" ++ path.sep_str ++ "sqrtf.c",
999 "math" ++ path.sep_str ++ "tanhf.c",997 "math" ++ path.sep_str ++ "tanhf.c",
1000 "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "acosf.c",998 "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "acosf.c",
1001 "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "asinf.c",999 "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "asinf.c",
src/libs/musl.zig-31
...@@ -826,8 +826,6 @@ const src_files = [_][]const u8{...@@ -826,8 +826,6 @@ const src_files = [_][]const u8{
826 "musl/src/math/aarch64/rintf.c",826 "musl/src/math/aarch64/rintf.c",
827 "musl/src/math/aarch64/round.c",827 "musl/src/math/aarch64/round.c",
828 "musl/src/math/aarch64/roundf.c",828 "musl/src/math/aarch64/roundf.c",
829 "musl/src/math/aarch64/sqrt.c",
830 "musl/src/math/aarch64/sqrtf.c",
831 "musl/src/math/acos.c",829 "musl/src/math/acos.c",
832 "musl/src/math/acosf.c",830 "musl/src/math/acosf.c",
833 "musl/src/math/acosh.c",831 "musl/src/math/acosh.c",
...@@ -836,8 +834,6 @@ const src_files = [_][]const u8{...@@ -836,8 +834,6 @@ const src_files = [_][]const u8{
836 "musl/src/math/acosl.c",834 "musl/src/math/acosl.c",
837 "musl/src/math/arm/fma.c",835 "musl/src/math/arm/fma.c",
838 "musl/src/math/arm/fmaf.c",836 "musl/src/math/arm/fmaf.c",
839 "musl/src/math/arm/sqrt.c",
840 "musl/src/math/arm/sqrtf.c",
841 "musl/src/math/asin.c",837 "musl/src/math/asin.c",
842 "musl/src/math/asinf.c",838 "musl/src/math/asinf.c",
843 "musl/src/math/asinh.c",839 "musl/src/math/asinh.c",
...@@ -944,9 +940,6 @@ const src_files = [_][]const u8{...@@ -944,9 +940,6 @@ const src_files = [_][]const u8{
944 "musl/src/math/i386/scalbnf.s",940 "musl/src/math/i386/scalbnf.s",
945 "musl/src/math/i386/scalbnl.s",941 "musl/src/math/i386/scalbnl.s",
946 "musl/src/math/i386/scalbn.s",942 "musl/src/math/i386/scalbn.s",
947 "musl/src/math/i386/sqrt.c",
948 "musl/src/math/i386/sqrtf.c",
949 "musl/src/math/i386/sqrtl.c",
950 "musl/src/math/ilogb.c",943 "musl/src/math/ilogb.c",
951 "musl/src/math/ilogbf.c",944 "musl/src/math/ilogbf.c",
952 "musl/src/math/ilogbl.c",945 "musl/src/math/ilogbl.c",
...@@ -986,7 +979,6 @@ const src_files = [_][]const u8{...@@ -986,7 +979,6 @@ const src_files = [_][]const u8{
986 "musl/src/math/lround.c",979 "musl/src/math/lround.c",
987 "musl/src/math/lroundf.c",980 "musl/src/math/lroundf.c",
988 "musl/src/math/lroundl.c",981 "musl/src/math/lroundl.c",
989 "musl/src/math/m68k/sqrtl.c",
990 "musl/src/math/__math_divzero.c",982 "musl/src/math/__math_divzero.c",
991 "musl/src/math/__math_divzerof.c",983 "musl/src/math/__math_divzerof.c",
992 "musl/src/math/__math_invalid.c",984 "musl/src/math/__math_invalid.c",
...@@ -998,8 +990,6 @@ const src_files = [_][]const u8{...@@ -998,8 +990,6 @@ const src_files = [_][]const u8{
998 "musl/src/math/__math_uflowf.c",990 "musl/src/math/__math_uflowf.c",
999 "musl/src/math/__math_xflow.c",991 "musl/src/math/__math_xflow.c",
1000 "musl/src/math/__math_xflowf.c",992 "musl/src/math/__math_xflowf.c",
1001 "musl/src/math/mips/sqrt.c",
1002 "musl/src/math/mips/sqrtf.c",
1003 "musl/src/math/modf.c",993 "musl/src/math/modf.c",
1004 "musl/src/math/modff.c",994 "musl/src/math/modff.c",
1005 "musl/src/math/modfl.c",995 "musl/src/math/modfl.c",
...@@ -1023,12 +1013,8 @@ const src_files = [_][]const u8{...@@ -1023,12 +1013,8 @@ const src_files = [_][]const u8{
1023 "musl/src/math/powerpc64/lroundf.c",1013 "musl/src/math/powerpc64/lroundf.c",
1024 "musl/src/math/powerpc64/round.c",1014 "musl/src/math/powerpc64/round.c",
1025 "musl/src/math/powerpc64/roundf.c",1015 "musl/src/math/powerpc64/roundf.c",
1026 "musl/src/math/powerpc64/sqrt.c",
1027 "musl/src/math/powerpc64/sqrtf.c",
1028 "musl/src/math/powerpc/fma.c",1016 "musl/src/math/powerpc/fma.c",
1029 "musl/src/math/powerpc/fmaf.c",1017 "musl/src/math/powerpc/fmaf.c",
1030 "musl/src/math/powerpc/sqrt.c",
1031 "musl/src/math/powerpc/sqrtf.c",
1032 "musl/src/math/powf.c",1018 "musl/src/math/powf.c",
1033 "musl/src/math/powf_data.c",1019 "musl/src/math/powf_data.c",
1034 "musl/src/math/powl.c",1020 "musl/src/math/powl.c",
...@@ -1047,12 +1033,8 @@ const src_files = [_][]const u8{...@@ -1047,12 +1033,8 @@ const src_files = [_][]const u8{
1047 "musl/src/math/rintl.c",1033 "musl/src/math/rintl.c",
1048 "musl/src/math/riscv32/fma.c",1034 "musl/src/math/riscv32/fma.c",
1049 "musl/src/math/riscv32/fmaf.c",1035 "musl/src/math/riscv32/fmaf.c",
1050 "musl/src/math/riscv32/sqrt.c",
1051 "musl/src/math/riscv32/sqrtf.c",
1052 "musl/src/math/riscv64/fma.c",1036 "musl/src/math/riscv64/fma.c",
1053 "musl/src/math/riscv64/fmaf.c",1037 "musl/src/math/riscv64/fmaf.c",
1054 "musl/src/math/riscv64/sqrt.c",
1055 "musl/src/math/riscv64/sqrtf.c",
1056 "musl/src/math/round.c",1038 "musl/src/math/round.c",
1057 "musl/src/math/roundf.c",1039 "musl/src/math/roundf.c",
1058 "musl/src/math/roundl.c",1040 "musl/src/math/roundl.c",
...@@ -1067,9 +1049,6 @@ const src_files = [_][]const u8{...@@ -1067,9 +1049,6 @@ const src_files = [_][]const u8{
1067 "musl/src/math/s390x/round.c",1049 "musl/src/math/s390x/round.c",
1068 "musl/src/math/s390x/roundf.c",1050 "musl/src/math/s390x/roundf.c",
1069 "musl/src/math/s390x/roundl.c",1051 "musl/src/math/s390x/roundl.c",
1070 "musl/src/math/s390x/sqrt.c",
1071 "musl/src/math/s390x/sqrtf.c",
1072 "musl/src/math/s390x/sqrtl.c",
1073 "musl/src/math/scalb.c",1052 "musl/src/math/scalb.c",
1074 "musl/src/math/scalbf.c",1053 "musl/src/math/scalbf.c",
1075 "musl/src/math/scalbln.c",1054 "musl/src/math/scalbln.c",
...@@ -1092,10 +1071,6 @@ const src_files = [_][]const u8{...@@ -1092,10 +1071,6 @@ const src_files = [_][]const u8{
1092 "musl/src/math/sinhl.c",1071 "musl/src/math/sinhl.c",
1093 "musl/src/math/__sinl.c",1072 "musl/src/math/__sinl.c",
1094 "musl/src/math/sinl.c",1073 "musl/src/math/sinl.c",
1095 "musl/src/math/sqrt.c",
1096 "musl/src/math/sqrt_data.c",
1097 "musl/src/math/sqrtf.c",
1098 "musl/src/math/sqrtl.c",
1099 "musl/src/math/__tan.c",1074 "musl/src/math/__tan.c",
1100 "musl/src/math/__tandf.c",1075 "musl/src/math/__tandf.c",
1101 "musl/src/math/tanh.c",1076 "musl/src/math/tanh.c",
...@@ -1128,9 +1103,6 @@ const src_files = [_][]const u8{...@@ -1128,9 +1103,6 @@ const src_files = [_][]const u8{
1128 "musl/src/math/x32/lrint.s",1103 "musl/src/math/x32/lrint.s",
1129 "musl/src/math/x32/remainderl.s",1104 "musl/src/math/x32/remainderl.s",
1130 "musl/src/math/x32/rintl.s",1105 "musl/src/math/x32/rintl.s",
1131 "musl/src/math/x32/sqrtf.s",
1132 "musl/src/math/x32/sqrtl.s",
1133 "musl/src/math/x32/sqrt.s",
1134 "musl/src/math/x86_64/acosl.s",1106 "musl/src/math/x86_64/acosl.s",
1135 "musl/src/math/x86_64/asinl.s",1107 "musl/src/math/x86_64/asinl.s",
1136 "musl/src/math/x86_64/atan2l.s",1108 "musl/src/math/x86_64/atan2l.s",
...@@ -1154,9 +1126,6 @@ const src_files = [_][]const u8{...@@ -1154,9 +1126,6 @@ const src_files = [_][]const u8{
1154 "musl/src/math/x86_64/remainderl.c",1126 "musl/src/math/x86_64/remainderl.c",
1155 "musl/src/math/x86_64/remquol.c",1127 "musl/src/math/x86_64/remquol.c",
1156 "musl/src/math/x86_64/rintl.c",1128 "musl/src/math/x86_64/rintl.c",
1157 "musl/src/math/x86_64/sqrt.c",
1158 "musl/src/math/x86_64/sqrtf.c",
1159 "musl/src/math/x86_64/sqrtl.c",
1160 "musl/src/misc/a64l.c",1129 "musl/src/misc/a64l.c",
1161 "musl/src/misc/basename.c",1130 "musl/src/misc/basename.c",
1162 "musl/src/misc/dirname.c",1131 "musl/src/misc/dirname.c",
src/libs/wasi_libc.zig-2
...@@ -834,8 +834,6 @@ const libc_top_half_src_files = [_][]const u8{...@@ -834,8 +834,6 @@ const libc_top_half_src_files = [_][]const u8{
834 "musl/src/math/sinhl.c",834 "musl/src/math/sinhl.c",
835 "musl/src/math/__sinl.c",835 "musl/src/math/__sinl.c",
836 "musl/src/math/sinl.c",836 "musl/src/math/sinl.c",
837 "musl/src/math/sqrt_data.c",
838 "musl/src/math/sqrtl.c",
839 "musl/src/math/__tan.c",837 "musl/src/math/__tan.c",
840 "musl/src/math/__tandf.c",838 "musl/src/math/__tandf.c",
841 "musl/src/math/tanh.c",839 "musl/src/math/tanh.c",