1//! Ported from musl, which is MIT licensed.
2//! https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//!
4//! https://git.musl-libc.org/cgit/musl/tree/src/math/sqrtf.c
5//! https://git.musl-libc.org/cgit/musl/tree/src/math/sqrt.c
6//! https://git.musl-libc.org/cgit/musl/tree/src/math/sqrtl.c
7
8const std = @import("std");
9const builtin = @import("builtin");
10const arch = builtin.cpu.arch;
11const math = std.math;
12const compiler_rt = @import("../compiler_rt.zig");
13const symbol = compiler_rt.symbol;
14
15comptime {
16 symbol(&__sqrth, "__sqrth");
17 symbol(&sqrtf, "sqrtf");
18 symbol(&sqrt, "sqrt");
19 symbol(&__sqrtx, "__sqrtx");
20 symbol(&sqrtq, "sqrtf128");
21 if (compiler_rt.want_sparc64_abi) {
22 symbol(&_Qp_sqrt, "_Qp_sqrt");
23 } else if (compiler_rt.want_sparc32_abi) {
24 symbol(&sqrtq, "_Q_sqrt");
25 }
26 symbol(&sqrtl, "sqrtl");
27}
28
29fn __sqrth(x: compiler_rt.f16.Abi) callconv(.c) compiler_rt.f16.Abi {
30 return compiler_rt.f16.toAbi(sqrt_f16(compiler_rt.f16.fromAbi(x)));
31}
32pub fn sqrt_f16(x: f16) f16 {
33 var ix: u16 = @bitCast(x);
34 var top = ix >> 10;
35
36 // special case handling.
37 if (top -% 0x01 >= 0x1F - 0x01) {
38 @branchHint(.unlikely);
39 // x < 0x1p-14 or inf or nan.
40 if (ix & 0x7FFF == 0) return x;
41 if (ix == 0x7C00) return x;
42 if (ix > 0x7C00) return math.nan(f16);
43 // x is subnormal, normalize it.
44 ix = @bitCast(x * 0x1p10);
45 top = (ix >> 10) -% 10;
46 }
47
48 // argument reduction:
49 // x = 4^e m; with integer e, and m in [1, 4)
50 // m: fixed point representation [2.14]
51 // 2^e is the exponent part of the result.
52 const even = (top & 1) != 0;
53 const m = if (even) (ix << 4) & 0x7FFF else (ix << 5) | 0x8000;
54 top = (top +% 0x0F) >> 1;
55
56 // approximate r ~ 1/sqrt(m) and s ~ sqrt(m) when m in [1,4)
57 // the fixed point representations are
58 // m: 2.14 r: 0.16, s: 2.14, d: 2.14, u: 2.14, three: 2.14
59 const three: u16 = 0xC000;
60 const i: usize = @intCast((ix >> 4) & 0x7F);
61 const r = rsqrt_tab[i];
62 // |r*sqrt(m) - 1| < 0x1p-8
63 var s = mul16(m, r);
64 // |s/sqrt(m) - 1| < 0x1p-8
65 const d = mul16(s, r);
66 const u = three - d;
67 s = mul16(s, u); // repr: 3.13
68 // -0x1.20p-13 < s/sqrt(m) - 1 < 0x7Dp-16
69 s = (s - 1) >> 3; // repr: 6.10
70 // s < sqrt(m) < s + 0x1.24p-10
71
72 // compute nearest rounded result:
73 // the nearest result to 10 bits is either s or s+0x1p-10,
74 // we can decide by comparing (2^10 s + 0.5)^2 to 2^20 m.
75 const d0 = (m << 6) -% s *% s;
76 const d1 = s -% d0;
77 const d2 = d1 +% s +% 1;
78 s += d1 >> 15;
79 s &= 0x03FF;
80 s |= top << 10;
81 const y: f16 = @bitCast(s);
82
83 // handle rounding modes and inexact exception:
84 // only (s+1)^2 == 2^6 m case is exact otherwise
85 // add a tiny value to cause the fenv effects.
86 if (d2 != 0) {
87 @branchHint(.likely);
88 var tiny: u16 = 0x0001;
89 tiny |= (d1 ^ d2) & 0x8000;
90 const t: f16 = @bitCast(tiny);
91 return y + t;
92 }
93
94 return y;
95}
96
97fn sqrtf(x: compiler_rt.f32.Abi) callconv(.c) compiler_rt.f32.Abi {
98 return compiler_rt.f32.toAbi(sqrt_f32(compiler_rt.f32.fromAbi(x)));
99}
100pub fn sqrt_f32(x: f32) f32 {
101 var ix: u32 = @bitCast(x);
102
103 if (ix < @as(u32, @bitCast(@as(f32, 0x1p-126))) or @as(u32, @bitCast(std.math.inf(f32))) <= ix) {
104 @branchHint(.unlikely);
105
106 if (ix & 0x7fffffff == 0)
107 return x;
108
109 if (ix == @as(u32, @bitCast(std.math.inf(f32))))
110 return x;
111
112 if (ix > @as(u32, @bitCast(std.math.inf(f32))))
113 return if (compiler_rt.want_float_exceptions) (x - x) / 0.0 else math.nan(f32);
114
115 ix = @as(u32, @bitCast(@as(i32, @bitCast(x * 0x1p23)) - (23 << 23)));
116 }
117
118 const m: u32 = if (ix & 0x00800000 != 0)
119 (ix << 7) & 0x7fffffff
120 else
121 (ix << 8) | 0x80000000;
122
123 const ey = ((ix >> 1) + (0x3f800000 >> 1)) & 0x7f800000;
124 // const ey = ((ix + 0x3f800000) & 0xff000000) >> 1;
125
126 const three = 0xc0000000;
127 const i = (ix >> 17) & 0x7f;
128 var r = @as(u32, rsqrt_tab[i]) << 16;
129
130 var s = mul32(m, r);
131 var d = mul32(s, r);
132 var u = three - d;
133 r = mul32(r, u) << 1;
134 s = mul32(s, u) << 1;
135 d = mul32(s, r);
136 u = three - d;
137 s = mul32(s, u);
138 s = (s - 1) >> 6;
139
140 const d0 = (m << 16) -% s *% s;
141 const d1 = s -% d0;
142 const d2 = d1 +% s +% 1;
143 const y: f32 = @bitCast(((s + (d1 >> 31)) & 0x007fffff) | ey);
144
145 const tiny: u32 = if (d2 == 0) blk: {
146 @branchHint(.unlikely);
147 break :blk 0;
148 } else 0x01000000;
149 const t: f32 = @bitCast(tiny | ((d1 ^ d2) & 0x80000000));
150
151 return y + t;
152}
153
154fn sqrt(x: compiler_rt.f64.Abi) callconv(.c) compiler_rt.f64.Abi {
155 return compiler_rt.f64.toAbi(sqrt_f64(compiler_rt.f64.fromAbi(x)));
156}
157pub fn sqrt_f64(x: f64) f64 {
158 var ix: u64 = @bitCast(x);
159 var top = ix >> 52;
160
161 // special case handling.
162 if (top -% 0x001 >= 0x7FF - 0x001) {
163 @branchHint(.unlikely);
164 // x < 0x1p-1022 or inf or nan.
165 if (ix & 0x7FFF_FFFF_FFFF_FFFF == 0) return x;
166 if (ix == 0x7FF0_0000_0000_0000) return x;
167 if (ix > 0x7FF0_0000_0000_0000) return if (compiler_rt.want_float_exceptions) (x - x) / 0.0 else math.nan(f64);
168 // x is subnormal, normalize it.
169 ix = @bitCast(x * 0x1p52);
170 top = (ix >> 52) -% 52;
171 }
172
173 // argument reduction:
174 // x = 4^e m; with integer e, and m in [1, 4)
175 // m: fixed point representation [2.62]
176 // 2^e is the exponent part of the result.
177 const even = (top & 1) != 0;
178 const m = if (even) (ix << 10) & 0x7FFF_FFFF_FFFF_FFFF else (ix << 11) | 0x8000_0000_0000_0000;
179 top = (top +% 0x3FF) >> 1;
180
181 // approximate r ~ 1/sqrt(m) and s ~ sqrt(m) when m in [1,4)
182 //
183 // initial estimate:
184 // 7bit table lookup (1bit exponent and 6bit significand).
185 //
186 // iterative approximation:
187 // using 2 goldschmidt iterations with 32bit int arithmetics
188 // and a final iteration with 64bit int arithmetics.
189 //
190 // details:
191 //
192 // the relative error (e = r0 sqrt(m)-1) of a linear estimate
193 // (r0 = a m + b) is |e| < 0.085955 ~ 0x1.6p-4 at best,
194 // a table lookup is faster and needs one less iteration
195 // 6 bit lookup table (128b) gives |e| < 0x1.f9p-8
196 // 7 bit lookup table (256b) gives |e| < 0x1.fdp-9
197 // for single and double prec 6bit is enough but for quad
198 // prec 7bit is needed (or modified iterations). to avoid
199 // one more iteration >=13bit table would be needed (16k).
200 //
201 // a newton-raphson iteration for r is
202 // w = r*r
203 // u = 3 - m*w
204 // r = r*u/2
205 // can use a goldschmidt iteration for s at the end or
206 // s = m*r
207 //
208 // first goldschmidt iteration is
209 // s = m*r
210 // u = 3 - s*r
211 // r = r*u/2
212 // s = s*u/2
213 // next goldschmidt iteration is
214 // u = 3 - s*r
215 // r = r*u/2
216 // s = s*u/2
217 // and at the end r is not computed only s.
218 //
219 // they use the same amount of operations and converge at the
220 // same quadratic rate, i.e. if
221 // r1 sqrt(m) - 1 = e, then
222 // r2 sqrt(m) - 1 = -3/2 e^2 - 1/2 e^3
223 // the advantage of goldschmidt is that the mul for s and r
224 // are independent (computed in parallel), however it is not
225 // "self synchronizing": it only uses the input m in the
226 // first iteration so rounding errors accumulate. at the end
227 // or when switching to larger precision arithmetics rounding
228 // errors dominate so the first iteration should be used.
229 //
230 // the fixed point representations are
231 // m: 2.30 r: 0.32, s: 2.30, d: 2.30, u: 2.30, three: 2.30
232 // and after switching to 64 bit
233 // m: 2.62 r: 0.64, s: 2.62, d: 2.62, u: 2.62, three: 2.62
234 const three: struct { u32, u64 } = .{
235 0xC000_0000,
236 0xC000_0000_0000_0000,
237 };
238 var r: struct { u32, u64 } = undefined;
239 var s: struct { u32, u64 } = undefined;
240 var d: struct { u32, u64 } = undefined;
241 var u: struct { u32, u64 } = undefined;
242 const i: usize = @intCast((ix >> 46) & 0x7F);
243 r[0] = @intCast(rsqrt_tab[i]);
244 r[0] <<= 16;
245 // |r sqrt(m) - 1| < 0x1.fdp-9
246 s[0] = mul32(@intCast(m >> 32), r[0]);
247 // |s/sqrt(m) - 1| < 0x1.fdp-9
248 d[0] = mul32(s[0], r[0]);
249 u[0] = three[0] - d[0];
250 r[0] = mul32(r[0], u[0]) << 1;
251 // |r sqrt(m) - 1| < 0x1.7bp-16
252 s[0] = mul32(s[0], u[0]) << 1;
253 // |s/sqrt(m) - 1| < 0x1.7bp-16
254 d[0] = mul32(s[0], r[0]);
255 u[0] = three[0] - d[0];
256 r[0] = mul32(r[0], u[0]) << 1;
257 // |r sqrt(m) - 1| < 0x1.3704p-29 (measured worst-case)
258 r[1] = @intCast(r[0]);
259 r[1] <<= 32;
260 s[1] = mul64(m, r[1]);
261 d[1] = mul64(s[1], r[1]);
262 u[1] = three[1] - d[1];
263 s[1] = mul64(s[1], u[1]); // repr: 3.61
264 // -0x1p-57 < s - sqrt(m) < 0x1.8001p-61
265 s[1] = (s[1] - 2) >> 9; // repr: 12.52
266 // -0x1.09p-52 < s - sqrt(m) < -0x1.fffcp-63
267
268 // s < sqrt(m) < s + 0x1.09p-52
269 // compute nearest rounded result:
270 // the nearest result to 52 bits is either s or s+0x1p-52,
271 // we can decide by comparing (2^52 s + 0.5)^2 to 2^104 m.
272 const d0 = (m << 42) -% s[1] *% s[1];
273 const d1 = s[1] -% d0;
274 const d2 = d1 +% s[1] +% 1;
275 s[1] += d1 >> 63;
276 s[1] &= 0x000F_FFFF_FFFF_FFFF;
277 s[1] |= top << 52;
278 const y: f64 = @bitCast(s[1]);
279
280 // handle rounding modes and inexact exception:
281 // only (s+1)^2 == 2^42 m case is exact otherwise
282 // add a tiny value to cause the fenv effects.
283 if (d2 != 0) {
284 @branchHint(.likely);
285 var tiny: u64 = 0x0010_0000_0000_0000;
286 tiny |= (d1 ^ d2) & 0x8000_0000_0000_0000;
287 const t: f64 = @bitCast(tiny);
288 return y + t;
289 }
290
291 return y;
292}
293
294fn __sqrtx(x: compiler_rt.f80.Abi) callconv(.c) compiler_rt.f80.Abi {
295 return compiler_rt.f80.toAbi(sqrt_f80(compiler_rt.f80.fromAbi(x)));
296}
297pub fn sqrt_f80(x: f80) f80 {
298 var ix: u80 = @bitCast(x);
299 var top = ix >> 64;
300
301 // special case handling.
302 if (top -% 0x0001 >= 0x7FFF - 0x0001) {
303 @branchHint(.unlikely);
304 // x < 0x1p-16382 or inf or nan.
305 if (ix & 0x7FFF_FFFF_FFFF_FFFF_FFFF == 0) return x;
306 if (ix == 0x7FFF_8000_0000_0000_0000) return x;
307 if (ix > 0x7FFF_8000_0000_0000_0000) return if (compiler_rt.want_float_exceptions) (x - x) / 0.0 else math.nan(f80);
308 // x is subnormal, normalize it.
309 ix = @bitCast(x * 0x1p63);
310 top = (ix >> 64) -% 63;
311 }
312
313 // argument reduction:
314 // x = 4^e m; with integer e, and m in [1, 4)
315 // m: fixed point representation [2.78]
316 // 2^e is the exponent part of the result.
317 const even = (top & 1) != 0;
318 const m = if (even) (ix << 15) & 0x7FFF_FFFF_FFFF_FFFF_FFFF else ix << 16;
319 top = (top +% 0x3FFF) >> 1;
320
321 // approximate r ~ 1/sqrt(m) and s ~ sqrt(m) when m in [1,4)
322 // the fixed point representations are
323 // m: 2.30 r: 0.32, s: 2.30, d: 2.30, u: 2.30, three: 2.30
324 // and after switching to 64 bit
325 // m: 2.62 r: 0.64, s: 2.62, d: 2.62, u: 2.62, three: 2.62
326 // and after switching to 80 bit
327 // m: 2.78 r: 0.80, s: 2.78, d: 2.78, u: 2.78, three: 2.78
328 const three: struct { u32, u64, u80 } = .{
329 0xC000_0000,
330 0xC000_0000_0000_0000,
331 0xC000_0000_0000_0000_0000,
332 };
333 var r: struct { u32, u64, u80 } = undefined;
334 var s: struct { u32, u64, u80 } = undefined;
335 var d: struct { u32, u64, u80 } = undefined;
336 var u: struct { u32, u64, u80 } = undefined;
337 var i: usize = @intCast((ix >> 57) & 0x3F);
338 if (even) i += 64;
339 r[0] = @intCast(rsqrt_tab[i]);
340 r[0] <<= 16;
341 // |r sqrt(m) - 1| < 0x1p-8
342 s[0] = mul32(@intCast(m >> 48), r[0]);
343 d[0] = mul32(s[0], r[0]);
344 u[0] = three[0] - d[0];
345 r[0] = mul32(u[0], r[0]) << 1;
346 // |r sqrt(m) - 1| < 0x1.7bp-16, switch to 64bit
347 r[1] = @intCast(r[0]);
348 r[1] <<= 32;
349 s[1] = mul64(@intCast(m >> 16), r[1]);
350 d[1] = mul64(s[1], r[1]);
351 u[1] = three[1] - d[1];
352 r[1] = mul64(u[1], r[1]) << 1;
353 // |r sqrt(m) - 1| < 0x1.a5p-31
354 s[1] = mul64(u[1], s[1]) << 1;
355 d[1] = mul64(s[1], r[1]);
356 u[1] = three[1] - d[1];
357 r[1] = mul64(u[1], r[1]) << 1;
358 // |r sqrt(m) - 1| < 0x1.c001p-59, switch to 80bit
359 r[2] = @intCast(r[1]);
360 r[2] <<= 16;
361 s[2] = mul80(m, r[2]);
362 d[2] = mul80(s[2], r[2]);
363 u[2] = three[2] - d[2];
364 s[2] = mul80(u[2], s[2]); // repr: 3.77
365 s[2] = (s[2] - 4) >> 14; // repr: 17.63
366 // s < sqrt(m) < s + 1 ULP + tiny
367
368 // compute nearest rounded result:
369 // the nearest result to 63 bits is either s or s+0x1p-63,
370 // we can decide by comparing (2^63 s + 0.5)^2 to 2^126 m
371 const d0 = (m << 48) -% mul80_tail(s[2], s[2]);
372 const d1 = s[2] -% d0;
373 const d2 = d1 +% s[2] +% 1;
374 s[2] += d1 >> 79;
375 s[2] &= 0x0000_7FFF_FFFF_FFFF_FFFF;
376 s[2] |= 0x0000_8000_0000_0000_0000;
377 s[2] |= top << 64;
378 const y: f80 = @bitCast(s[2]);
379
380 // handle rounding modes and inexact exception:
381 // only (s+1)^2 == 2^48 m case is exact otherwise
382 // add a tiny value to cause the fenv effects.
383 if (d2 != 0) {
384 @branchHint(.likely);
385 var tiny: u80 = 0x0001_8000_0000_0000_0000;
386 tiny |= (d1 ^ d2) & 0x8000_0000_0000_0000_0000;
387 const t: f80 = @bitCast(tiny);
388 return y + t;
389 }
390
391 return y;
392}
393
394fn sqrtq(x: compiler_rt.f128.Abi) callconv(.c) compiler_rt.f128.Abi {
395 return compiler_rt.f128.toAbi(sqrt_f128(compiler_rt.f128.fromAbi(x)));
396}
397pub fn sqrt_f128(x: f128) f128 {
398 var ix: u128 = @bitCast(x);
399 var top = ix >> 112;
400
401 // special case handling.
402 if (top -% 0x0001 >= 0x7FFF - 0x0001) {
403 @branchHint(.unlikely);
404 // x < 0x1p-16382 or inf or nan.
405 if (ix & 0x7FFF_FFFF_FFFF_FFFF_FFFF_FFFF_FFFF_FFFF == 0) return x;
406 if (ix == 0x7FFF_0000_0000_0000_0000_0000_0000_0000) return x;
407 if (ix > 0x7FFF_0000_0000_0000_0000_0000_0000_0000) return math.nan(f128);
408 // x is subnormal, normalize it.
409 ix = @bitCast(x * 0x1p112);
410 top = (ix >> 112) -% 112;
411 }
412
413 // argument reduction:
414 // x = 4^e m; with integer e, and m in [1, 4)
415 // m: fixed point representation [2.126]
416 // 2^e is the exponent part of the result.
417 const even = (top & 1) != 0;
418 const m = if (even) (ix << 14) & 0x7FFF_FFFF_FFFF_FFFF_FFFF_FFFF_FFFF_FFFF else (ix << 15) | 0x8000_0000_0000_0000_0000_0000_0000_0000;
419 top = (top +% 0x3FFF) >> 1;
420
421 // approximate r ~ 1/sqrt(m) and s ~ sqrt(m) when m in [1,4)
422 // the fixed point representations are
423 // m: 2.30 r: 0.32, s: 2.30, d: 2.30, u: 2.30, three: 2.30
424 // and after switching to 64 bit
425 // m: 2.62 r: 0.64, s: 2.62, d: 2.62, u: 2.62, three: 2.62
426 // and after switching to 128 bit
427 // m: 2.126 r: 0.128, s: 2.126, d: 2.126, u: 2.126, three: 2.126
428 const three: struct { u32, u64, u128 } = .{
429 0xC000_0000,
430 0xC000_0000_0000_0000,
431 0xC000_0000_0000_0000_0000_0000_0000_0000,
432 };
433 var r: struct { u32, u64, u128 } = undefined;
434 var s: struct { u32, u64, u128 } = undefined;
435 var d: struct { u32, u64, u128 } = undefined;
436 var u: struct { u32, u64, u128 } = undefined;
437 const i: usize = @intCast((ix >> 106) & 0x7F);
438 r[0] = @intCast(rsqrt_tab[i]);
439 r[0] <<= 16;
440 // |r sqrt(m) - 1| < 0x1p-8
441 s[0] = mul32(@intCast(m >> 96), r[0]);
442 d[0] = mul32(s[0], r[0]);
443 u[0] = three[0] - d[0];
444 r[0] = mul32(u[0], r[0]) << 1;
445 // |r sqrt(m) - 1| < 0x1.7bp-16, switch to 64bit
446 r[1] = @intCast(r[0]);
447 r[1] <<= 32;
448 s[1] = mul64(@intCast(m >> 64), r[1]);
449 d[1] = mul64(s[1], r[1]);
450 u[1] = three[1] - d[1];
451 r[1] = mul64(u[1], r[1]) << 1;
452 // |r sqrt(m) - 1| < 0x1.a5p-31
453 s[1] = mul64(u[1], s[1]) << 1;
454 d[1] = mul64(s[1], r[1]);
455 u[1] = three[1] - d[1];
456 r[1] = mul64(u[1], r[1]) << 1;
457 // |r sqrt(m) - 1| < 0x1.c001p-59, switch to 128bit
458 r[2] = @intCast(r[1]);
459 r[2] <<= 64;
460 s[2] = mul128(m, r[2]);
461 d[2] = mul128(s[2], r[2]);
462 u[2] = three[2] - d[2];
463 s[2] = mul128(u[2], s[2]); // repr: 3.125
464 // -0x1p-116 < s - sqrt(m) < 0x3.8001p-125
465 s[2] = (s[2] - 4) >> 13; // repr: 16.122
466 // s < sqrt(m) < s + 1 ULP + tiny
467
468 // compute nearest rounded result:
469 // the nearest result to 122 bits is either s or s+0x1p-122,
470 // we can decide by comparing (2^122 s + 0.5)^2 to 2^244 m
471 const d0 = (m << 98) -% s[2] *% s[2];
472 const d1 = s[2] -% d0;
473 const d2 = d1 +% s[2] +% 1;
474 s[2] += d1 >> 127;
475 s[2] &= 0x0000_FFFF_FFFF_FFFF_FFFF_FFFF_FFFF_FFFF;
476 s[2] |= top << 112;
477 const y: f128 = @bitCast(s[2]);
478
479 // handle rounding modes and inexact exception:
480 // only (s+1)^2 == 2^98 m case is exact otherwise
481 // add a tiny value to cause the fenv effects.
482 if (d2 != 0) {
483 @branchHint(.likely);
484 var tiny: u128 = 0x0001_0000_0000_0000_0000_0000_0000_0000;
485 tiny |= (d1 ^ d2) & 0x8000_0000_0000_0000_0000_0000_0000_0000;
486 const t: f128 = @bitCast(tiny);
487 return y + t;
488 }
489
490 return y;
491}
492
493fn _Qp_sqrt(c: *f128, a: *f128) callconv(.c) void {
494 c.* = sqrt(@floatCast(a.*));
495}
496
497pub fn sqrtl(x: c_longdouble) callconv(.c) c_longdouble {
498 switch (@typeInfo(c_longdouble).float.bits) {
499 64 => return sqrt_f64(x),
500 80 => return sqrt_f80(x),
501 128 => return sqrt_f128(x),
502 else => comptime unreachable,
503 }
504}
505
506const rsqrt_tab: [128]u16 = .{
507 0xB451, 0xB2F0, 0xB196, 0xB044, 0xAEF9, 0xADB6, 0xAC79, 0xAB43,
508 0xAA14, 0xA8EB, 0xA7C8, 0xA6AA, 0xA592, 0xA480, 0xA373, 0xA26B,
509 0xA168, 0xA06A, 0x9F70, 0x9E7B, 0x9D8A, 0x9C9D, 0x9BB5, 0x9AD1,
510 0x99F0, 0x9913, 0x983A, 0x9765, 0x9693, 0x95C4, 0x94F8, 0x9430,
511 0x936B, 0x92A9, 0x91EA, 0x912E, 0x9075, 0x8FBE, 0x8F0A, 0x8E59,
512 0x8DAA, 0x8CFE, 0x8C54, 0x8BAC, 0x8B07, 0x8A64, 0x89C4, 0x8925,
513 0x8889, 0x87EE, 0x8756, 0x86C0, 0x862B, 0x8599, 0x8508, 0x8479,
514 0x83EC, 0x8361, 0x82D8, 0x8250, 0x81C9, 0x8145, 0x80C2, 0x8040,
515 0xFF02, 0xFD0E, 0xFB25, 0xF947, 0xF773, 0xF5AA, 0xF3EA, 0xF234,
516 0xF087, 0xEEE3, 0xED47, 0xEBB3, 0xEA27, 0xE8A3, 0xE727, 0xE5B2,
517 0xE443, 0xE2DC, 0xE17A, 0xE020, 0xDECB, 0xDD7D, 0xDC34, 0xDAF1,
518 0xD9B3, 0xD87B, 0xD748, 0xD61A, 0xD4F1, 0xD3CD, 0xD2AD, 0xD192,
519 0xD07B, 0xCF69, 0xCE5B, 0xCD51, 0xCC4A, 0xCB48, 0xCA4A, 0xC94F,
520 0xC858, 0xC764, 0xC674, 0xC587, 0xC49D, 0xC3B7, 0xC2D4, 0xC1F4,
521 0xC116, 0xC03C, 0xBF65, 0xBE90, 0xBDBE, 0xBCEF, 0xBC23, 0xBB59,
522 0xBA91, 0xB9CC, 0xB90A, 0xB84A, 0xB78C, 0xB6D0, 0xB617, 0xB560,
523};
524
525inline fn mul16(a: u16, b: u16) u16 {
526 return @intCast(@as(u32, a) * b >> 16);
527}
528
529inline fn mul32(a: u32, b: u32) u32 {
530 return @intCast(@as(u64, a) * b >> 32);
531}
532
533inline fn mul64(a: u64, b: u64) u64 {
534 return @intCast(@as(u128, a) * b >> 64);
535}
536
537inline fn mul80(a: u80, b: u80) u80 {
538 const ahi = a >> 40;
539 const alo = a & 0xFF_FFFF_FFFF;
540 const bhi = b >> 40;
541 const blo = b & 0xFF_FFFF_FFFF;
542 return ahi * bhi + (ahi * blo >> 40) + (alo * bhi >> 40);
543}
544
545inline fn mul128(a: u128, b: u128) u128 {
546 const ahi = a >> 64;
547 const alo = a & 0xFFFF_FFFF_FFFF_FFFF;
548 const bhi = b >> 64;
549 const blo = b & 0xFFFF_FFFF_FFFF_FFFF;
550 return ahi * bhi + (ahi * blo >> 64) + (alo * bhi >> 64);
551}
552
553inline fn mul80_tail(a: u80, b: u80) u80 {
554 const ahi = a >> 40;
555 const alo = a & 0xFF_FFFF_FFFF;
556 const bhi = b >> 40;
557 const blo = b & 0xFF_FFFF_FFFF;
558 return alo * blo +% ((ahi * blo) << 40) +% ((alo * bhi) << 40);
559}
560
561test "sqrt_f16" {
562 // sqrt(±0) is ±0
563 try std.testing.expectEqual(sqrt_f16(0x0.0p0), 0x0.0p0);
564 try std.testing.expectEqual(sqrt_f16(-0x0.0p0), -0x0.0p0);
565 // sqrt(+max) is finite
566 try std.testing.expectEqual(sqrt_f16(0x1.FFCp15), 0x1.FFCp7);
567 // sqrt(4)=2
568 try std.testing.expectEqual(sqrt_f16(0x1p2), 0x1p1);
569 // sqrt(x) for x=1, 1±ulp
570 try std.testing.expectEqual(sqrt_f16(0x1p0), 0x1p0);
571 try std.testing.expectEqual(sqrt_f16(0x1.004p0), 0x1p0);
572 try std.testing.expectEqual(sqrt_f16(0x1.FF8p-1), 0x1.FFCp-1);
573 // sqrt(+min) is non-zero
574 try std.testing.expectEqual(sqrt_f16(0x1p-14), 0x1p-7);
575 // sqrt(min subnormal) is non-zero
576 try std.testing.expectEqual(sqrt_f16(0x0.004p-14), 0x1p-12);
577 // sqrt(inf) is inf
578 try std.testing.expect(math.isInf(sqrt_f16(math.inf(f16))));
579 // sqrt(nan) is nan
580 try std.testing.expect(math.isNan(sqrt_f16(math.nan(f16))));
581 // sqrt(-ve) is nan
582 try std.testing.expect(math.isNan(sqrt_f16(-0x1p-14)));
583 try std.testing.expect(math.isNan(sqrt_f16(-0x1p+0)));
584 try std.testing.expect(math.isNan(sqrt_f16(-math.inf(f16))));
585 // random arguments
586 try std.testing.expectEqual(sqrt_f16(0x1.1p14), 0x1.08p7);
587 try std.testing.expectEqual(sqrt_f16(0x1.C9p-12), 0x1.56p-6);
588 try std.testing.expectEqual(sqrt_f16(0x1.CE8p-7), 0x1.E68p-4);
589 try std.testing.expectEqual(sqrt_f16(0x1.134p-7), 0x1.778p-4);
590 try std.testing.expectEqual(sqrt_f16(0x1.E9Cp-10), 0x1.62p-5);
591 try std.testing.expectEqual(sqrt_f16(0x1.3Dp9), 0x1.92Cp4);
592 try std.testing.expectEqual(sqrt_f16(0x1.AA4p8), 0x1.4A4p4);
593 try std.testing.expectEqual(sqrt_f16(0x1.8A8p4), 0x1.3DCp2);
594 try std.testing.expectEqual(sqrt_f16(0x1.8Fp-7), 0x1.C4p-4);
595 try std.testing.expectEqual(sqrt_f16(0x1.584p-11), 0x1.A3Cp-6);
596}
597
598test "sqrt_f32" {
599 // sqrt(±0) is ±0
600 try std.testing.expectEqual(sqrt_f32(0x0.0p0), 0x0.0p0);
601 try std.testing.expectEqual(sqrt_f32(-0x0.0p0), -0x0.0p0);
602 // sqrt(+max) is finite
603 try std.testing.expectEqual(sqrt_f32(0x1.FFFFFEp127), 0x1.FFFFFEp63);
604 // sqrt(4)=2
605 try std.testing.expectEqual(sqrt_f32(0x1p2), 0x1p1);
606 // sqrt(x) for x=1, 1±ulp
607 try std.testing.expectEqual(sqrt_f32(0x1p0), 0x1p0);
608 try std.testing.expectEqual(sqrt_f32(0x1.000002p0), 0x1p0);
609 try std.testing.expectEqual(sqrt_f32(0x1.FFFFFEp-1), 0x1.FFFFFEp-1);
610 // sqrt(+min) is non-zero
611 try std.testing.expectEqual(sqrt_f32(0x1p-126), 0x1p-63);
612 // sqrt(min subnormal) is non-zero
613 try std.testing.expectEqual(sqrt_f32(0x0.000002p-126), 0x1.6a09e6p-75);
614 // sqrt(inf) is inf
615 try std.testing.expect(math.isInf(sqrt_f32(math.inf(f32))));
616 // sqrt(nan) is nan
617 try std.testing.expect(math.isNan(sqrt_f32(math.nan(f32))));
618 // sqrt(-ve) is nan
619 try std.testing.expect(math.isNan(sqrt_f32(-0x1p-149)));
620 try std.testing.expect(math.isNan(sqrt_f32(-0x1p0)));
621 try std.testing.expect(math.isNan(sqrt_f32(-math.inf(f32))));
622 // random arguments
623 try std.testing.expectEqual(sqrt_f32(0x1.4DD57Ep77), 0x1.9D6DA8p38);
624 try std.testing.expectEqual(sqrt_f32(0x1.871848p102), 0x1.3C6AFAp51);
625 try std.testing.expectEqual(sqrt_f32(0x1.A1D748p-112), 0x1.470EFCp-56);
626 try std.testing.expectEqual(sqrt_f32(0x1.E626C2p18), 0x1.60C80Ep9);
627 try std.testing.expectEqual(sqrt_f32(0x1.E80E66p-29), 0x1.F3E282p-15);
628 try std.testing.expectEqual(sqrt_f32(0x1.B47204p89), 0x1.D8B732p44);
629 try std.testing.expectEqual(sqrt_f32(0x1.77F45p15), 0x1.B6BC3Ap7);
630 try std.testing.expectEqual(sqrt_f32(0x1.AD5F5p-48), 0x1.4B8A72p-24);
631 try std.testing.expectEqual(sqrt_f32(0x1.91A39p-76), 0x1.40A7A8p-38);
632 try std.testing.expectEqual(sqrt_f32(0x1.DAE088p79), 0x1.ED16DCp39);
633}
634
635test "sqrt_f64" {
636 // sqrt(±0) is ±0
637 try std.testing.expectEqual(sqrt_f64(0x0.0p0), 0x0.0p0);
638 try std.testing.expectEqual(sqrt_f64(-0x0.0p0), -0x0.0p0);
639 // sqrt(+max) is finite
640 try std.testing.expectEqual(sqrt_f64(math.floatMax(f64)), 0x1.FFFFFFFFFFFFFp511);
641 // sqrt(4)=2
642 try std.testing.expectEqual(sqrt_f64(0x1p2), 0x1p1);
643 // sqrt(x) for x=1, 1±ulp
644 try std.testing.expectEqual(sqrt_f64(0x1p0), 0x1p0);
645 try std.testing.expectEqual(sqrt_f64(0x1p0 + math.floatEps(f64)), 0x1p0);
646 try std.testing.expectEqual(sqrt_f64(0x1p0 - math.floatEps(f64)), 0x1.FFFFFFFFFFFFFp-1);
647 // sqrt(+min) is non-zero
648 try std.testing.expectEqual(sqrt_f64(math.floatMin(f64)), 0x1p-511);
649 // sqrt(min subnormal) is non-zero
650 try std.testing.expectEqual(sqrt_f64(math.floatTrueMin(f64)), 0x1p-537);
651 // sqrt(inf) is inf
652 try std.testing.expect(math.isInf(sqrt_f64(math.inf(f64))));
653 // sqrt(nan) is nan
654 try std.testing.expect(math.isNan(sqrt_f64(math.nan(f64))));
655 // sqrt(-ve) is nan
656 try std.testing.expect(math.isNan(sqrt_f64(-0x1p-1074)));
657 try std.testing.expect(math.isNan(sqrt_f64(-0x1p0)));
658 try std.testing.expect(math.isNan(sqrt_f64(-math.inf(f64))));
659 // random arguments
660 try std.testing.expectEqual(sqrt_f64(0x1.27D3510D4789Bp471), 0x1.852E97E58CFB7p235);
661 try std.testing.expectEqual(sqrt_f64(0x1.8C4FCD5A07846p791), 0x1.C27504E56D938p395);
662 try std.testing.expectEqual(sqrt_f64(0x1.B1B69324F96E7p-137), 0x1.D73BD0414D8BFp-69);
663 try std.testing.expectEqual(sqrt_f64(0x1.1CBD179A811FEp278), 0x1.0DFCB9A114A61p139);
664 try std.testing.expectEqual(sqrt_f64(0x1.1D0C7EFB04A56p917), 0x1.7E0708A25DDCDp458);
665 try std.testing.expectEqual(sqrt_f64(0x1.21B355DA8C94Bp-249), 0x1.8121CBE2608E3p-125);
666 try std.testing.expectEqual(sqrt_f64(0x1.63024D4C5E987p487), 0x1.AA56AEA589DCDp243);
667 try std.testing.expectEqual(sqrt_f64(0x1.45AC3BE941F6Ep339), 0x1.9857F3F453E2Dp169);
668 try std.testing.expectEqual(sqrt_f64(0x1.3B719C733AA24p267), 0x1.91E12E3AC8F71p133);
669 try std.testing.expectEqual(sqrt_f64(0x1.0B150433A2275p357), 0x1.71CAB87F8277Cp178);
670}
671
672test "__sqrtx" {
673 // sqrt(±0) is ±0
674 try std.testing.expectEqual(sqrt_f80(0x0.0p0), 0x0.0p0);
675 try std.testing.expectEqual(sqrt_f80(-0x0.0p0), -0x0.0p0);
676 // sqrt(+max) is finite
677 try std.testing.expectEqual(sqrt_f80(math.floatMax(f80)), 0x1.FFFFFFFFFFFFFFFEp8191);
678 // sqrt(4)=2
679 try std.testing.expectEqual(sqrt_f80(0x1p2), 0x1p1);
680 // sqrt(x) for x=1, 1±ulp
681 try std.testing.expectEqual(sqrt_f80(0x1p0), 0x1p0);
682 try std.testing.expectEqual(sqrt_f80(0x1p0 + math.floatEps(f80)), 0x1p0);
683 try std.testing.expectEqual(sqrt_f80(0x1p0 - math.floatEps(f80)), 0x1.FFFFFFFFFFFFFFFEp-1);
684 // sqrt(+min) is non-zero
685 try std.testing.expectEqual(sqrt_f80(math.floatMin(f80)), 0x1p-8191);
686 // sqrt(min subnormal) is non-zero
687 try std.testing.expectEqual(sqrt_f80(math.floatTrueMin(f80)), 0x1.6A09E667F3BCC908p-8223);
688 // sqrt(inf) is inf
689 try std.testing.expect(math.isInf(sqrt_f80(math.inf(f80))));
690 // sqrt(nan) is nan
691 try std.testing.expect(math.isNan(sqrt_f80(math.nan(f80))));
692 // sqrt(-ve) is nan
693 try std.testing.expect(math.isNan(sqrt_f80(-0x1p-16442)));
694 try std.testing.expect(math.isNan(sqrt_f80(-0x1p0)));
695 try std.testing.expect(math.isNan(sqrt_f80(-math.inf(f80))));
696 // random arguments
697 try std.testing.expectEqual(sqrt_f80(0x1.087F3953486918A4p15482), 0x1.0436BBE03D02F32p7741);
698 try std.testing.expectEqual(sqrt_f80(0x1.530CF9E2AE84D8Fp-6330), 0x1.269CFEF51933BE58p-3165);
699 try std.testing.expectEqual(sqrt_f80(0x1.3F971515EADD574Ap5713), 0x1.9483232AB780B006p2856);
700 try std.testing.expectEqual(sqrt_f80(0x1.4CC0DC7379222954p864), 0x1.23DD4D0A4758C2Cp432);
701 try std.testing.expectEqual(sqrt_f80(0x1.920E5649559A839Ep-3181), 0x1.C5B5BC0F98DD83D2p-1591);
702 try std.testing.expectEqual(sqrt_f80(0x1.2E59726F87CD1746p-629), 0x1.8973327E95CB350Cp-315);
703 try std.testing.expectEqual(sqrt_f80(0x1.D3A16391F57B4D64p-9034), 0x1.59FF08B7DEEF5DB2p-4517);
704 try std.testing.expectEqual(sqrt_f80(0x1.E7053D8DAA49BCEEp-11411), 0x1.F35AA3EA5E18E344p-5706);
705 try std.testing.expectEqual(sqrt_f80(0x1.797ED0B05DD4A984p7521), 0x1.B7A22E40C6A7867Ap3760);
706 try std.testing.expectEqual(sqrt_f80(0x1.FC50806445C7226Ap15371), 0x1.FE2766142653F5BEp7685);
707}
708
709test "sqrt_f128" {
710 // sqrt(±0) is ±0
711 try std.testing.expectEqual(sqrt_f128(0x0.0p0), 0x0.0p0);
712 try std.testing.expectEqual(sqrt_f128(-0x0.0p0), -0x0.0p0);
713 // sqrt(+max) is finite
714 try std.testing.expectEqual(sqrt_f128(math.floatMax(f128)), 0x1.FFFFFFFFFFFFFFFFFFFFFFFFFFFFp8191);
715 // sqrt(4)=2
716 try std.testing.expectEqual(sqrt_f128(0x1p2), 0x1p1);
717 // sqrt(x) for x=1, 1±ulp
718 try std.testing.expectEqual(sqrt_f128(0x1p0), 0x1p0);
719 try std.testing.expectEqual(sqrt_f128(0x1p0 + math.floatEps(f128)), 0x1p0);
720 try std.testing.expectEqual(sqrt_f128(0x1p0 - math.floatEps(f128)), 0x1.FFFFFFFFFFFFFFFFFFFFFFFFFFFFp-1);
721 // sqrt(+min) is non-zero
722 try std.testing.expectEqual(sqrt_f128(math.floatMin(f128)), 0x1p-8191);
723 // sqrt(min subnormal) is non-zero
724 try std.testing.expectEqual(sqrt_f128(math.floatTrueMin(f128)), 0x1p-8247);
725 // sqrt(inf) is inf
726 try std.testing.expect(math.isInf(sqrt_f128(math.inf(f128))));
727 // sqrt(nan) is nan
728 try std.testing.expect(math.isNan(sqrt_f128(math.nan(f128))));
729 // sqrt(-ve) is nan
730 try std.testing.expect(math.isNan(sqrt_f128(-0x1p-16442)));
731 try std.testing.expect(math.isNan(sqrt_f128(-0x1p0)));
732 try std.testing.expect(math.isNan(sqrt_f128(-math.inf(f128))));
733 // random arguments
734 try std.testing.expectEqual(sqrt_f128(0x1.B6942D29A331751600C9F3AF7E5Fp3363), 0x1.D9DE9AFEF0F2D25586A50CA39D4Dp1681);
735 try std.testing.expectEqual(sqrt_f128(0x1.5E65C405F84D471A8070ADD7A42Dp11765), 0x1.A78F7F9452B4D9EC2403C81D9D42p5882);
736 try std.testing.expectEqual(sqrt_f128(0x1.B42334D68F8016D8AE6F5E22B044p-5624), 0x1.4E247A7F2FF2A325E9377BB09C8p-2812);
737 try std.testing.expectEqual(sqrt_f128(0x1.E61715047F80F2E0B9382B38E06Bp10062), 0x1.60C25D9DFDC0116B78EF5AFDE0E9p5031);
738 try std.testing.expectEqual(sqrt_f128(0x1.2ED0B53B494CB55A7B04E653D40Ep-1026), 0x1.166CE78D658D2453D700B04C5748p-513);
739 try std.testing.expectEqual(sqrt_f128(0x1.1BA756B9790E78A4E6F0B083AA89p1835), 0x1.7D1767EA3303DB7A46940033988p917);
740 try std.testing.expectEqual(sqrt_f128(0x1.5B6C574319C1120335C8E1609704p4512), 0x1.2A3A8A415BB1648C548FBA2A4182p2256);
741 try std.testing.expectEqual(sqrt_f128(0x1.FF91E8CDEE1552A2B74E77B602Ep14953), 0x1.FFC8F171267D4FE75CBE7AB4D851p7476);
742 try std.testing.expectEqual(sqrt_f128(0x1.9B1837CFC629A1B6B1BB97099E7Dp2892), 0x1.4468511B909EAF8641BD59105A6Bp1446);
743 try std.testing.expectEqual(sqrt_f128(0x1.0E2115475E64A92340914E7F7B37p-13951), 0x1.73E536F82F414134012F55BA5368p-6976);
744}