authorgravatar for pivoc@protonmail.comPivok <pivoc@protonmail.com> 2026-02-05 21:57:32+01:00
committergravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2026-02-05 21:57:32+01:00
logd0b39c7f2b95522c262a92b97bf5654721ace1c7
tree5e2b0a6739be02cb6c61d9b5f23723c298d085a1
parentc38f9336a3157ae863a795b005db670c6cc04906

libzigc: hypot (#31104)

First time contribution. Implements hypot for libzigc #30978. Commands i run: ``` $ stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib $ stage4/bin/zig build test-libc -Dlibc-test-path=../../libc-test -Dtest-filter=hypot --summary line -fqemu -fwasmtime Build Summary: 725/737 steps succeeded (12 skipped) ``` I also changed std.math.hypot becuase some libc-tests raised fp exceptions. Example: ``` ../../libc-test/src/math/special/hypot.h:8: bad fp exception: RN hypot(0x1p-1074,0x0p+0)=0x1p-1074, want 0 got INEXACT|UNDERFLOW ../../libc-test/src/math/special/hypot.h:9: bad fp exception: RN hypot(0x1p-1074,-0x0p+0)=0x1p-1074, want 0 got INEXACT|UNDERFLOW ``` I also run this command as a quick sanity check: ``` $ stage4/bin/zig build test-std -Dtest-filter=hypot -Dtest-target-filter=x86_64-linux-musl --summary line Build Summary: 5/5 steps succeeded; 136/136 tests passed ``` Reviewed-on: https://codeberg.org/ziglang/zig/pulls/31104 Reviewed-by: Andrew Kelley <andrew@ziglang.org> Co-authored-by: Pivok <pivoc@protonmail.com> Co-committed-by: Pivok <pivoc@protonmail.com>

6 files changed, 9 insertions(+), 119 deletions(-)

lib/c/math.zig+5
......@@ -41,6 +41,7 @@ comptime {
4141 @export(&atanl, .{ .name = "atanl", .linkage = common.linkage, .visibility = common.visibility });
4242 @export(&cbrt, .{ .name = "cbrt", .linkage = common.linkage, .visibility = common.visibility });
4343 @export(&cbrtf, .{ .name = "cbrtf", .linkage = common.linkage, .visibility = common.visibility });
44 @export(&hypot, .{ .name = "hypot", .linkage = common.linkage, .visibility = common.visibility });
4445 @export(&pow, .{ .name = "pow", .linkage = common.linkage, .visibility = common.visibility });
4546 }
4647
......@@ -118,6 +119,10 @@ fn cbrtf(x: f32) callconv(.c) f32 {
118119 return math.cbrt(x);
119120}
120121
122fn hypot(x: f64, y: f64) callconv(.c) f64 {
123 return math.hypot(x, y);
124}
125
121126fn pow(x: f64, y: f64) callconv(.c) f64 {
122127 return math.pow(f64, x, y);
123128}
lib/libc/musl/src/math/hypot.c deleted-67
......@@ -1,67 +0,0 @@
1#include <math.h>
2#include <stdint.h>
3#include <float.h>
4
5#if FLT_EVAL_METHOD > 1U && LDBL_MANT_DIG == 64
6#define SPLIT (0x1p32 + 1)
7#else
8#define SPLIT (0x1p27 + 1)
9#endif
10
11static void sq(double_t *hi, double_t *lo, double x)
12{
13 double_t xh, xl, xc;
14
15 xc = (double_t)x*SPLIT;
16 xh = x - xc + xc;
17 xl = x - xh;
18 *hi = (double_t)x*x;
19 *lo = xh*xh - *hi + 2*xh*xl + xl*xl;
20}
21
22double hypot(double x, double y)
23{
24 union {double f; uint64_t i;} ux = {x}, uy = {y}, ut;
25 int ex, ey;
26 double_t hx, lx, hy, ly, z;
27
28 /* arrange |x| >= |y| */
29 ux.i &= -1ULL>>1;
30 uy.i &= -1ULL>>1;
31 if (ux.i < uy.i) {
32 ut = ux;
33 ux = uy;
34 uy = ut;
35 }
36
37 /* special cases */
38 ex = ux.i>>52;
39 ey = uy.i>>52;
40 x = ux.f;
41 y = uy.f;
42 /* note: hypot(inf,nan) == inf */
43 if (ey == 0x7ff)
44 return y;
45 if (ex == 0x7ff || uy.i == 0)
46 return x;
47 /* note: hypot(x,y) ~= x + y*y/x/2 with inexact for small y/x */
48 /* 64 difference is enough for ld80 double_t */
49 if (ex - ey > 64)
50 return x + y;
51
52 /* precise sqrt argument in nearest rounding mode without overflow */
53 /* xh*xh must not overflow and xl*xl must not underflow in sq */
54 z = 1;
55 if (ex > 0x3ff+510) {
56 z = 0x1p700;
57 x *= 0x1p-700;
58 y *= 0x1p-700;
59 } else if (ey < 0x3ff-450) {
60 z = 0x1p-700;
61 x *= 0x1p700;
62 y *= 0x1p700;
63 }
64 sq(&hx, &lx, x);
65 sq(&hy, &ly, y);
66 return z*sqrt(ly+lx+hy+hx);
67}
lib/libc/musl/src/math/i386/hypot.s deleted-45
......@@ -1,45 +0,0 @@
1.global hypot
2.type hypot,@function
3hypot:
4 mov 8(%esp),%eax
5 mov 16(%esp),%ecx
6 add %eax,%eax
7 add %ecx,%ecx
8 and %eax,%ecx
9 cmp $0xffe00000,%ecx
10 jae 2f
11 or 4(%esp),%eax
12 jnz 1f
13 fldl 12(%esp)
14 fabs
15 ret
161: mov 16(%esp),%eax
17 add %eax,%eax
18 or 12(%esp),%eax
19 jnz 1f
20 fldl 4(%esp)
21 fabs
22 ret
231: fldl 4(%esp)
24 fld %st(0)
25 fmulp
26 fldl 12(%esp)
27 fld %st(0)
28 fmulp
29 faddp
30 fsqrt
31 ret
322: sub $0xffe00000,%eax
33 or 4(%esp),%eax
34 jnz 1f
35 fldl 4(%esp)
36 fabs
37 ret
381: mov 16(%esp),%eax
39 add %eax,%eax
40 sub $0xffe00000,%eax
41 or 12(%esp),%eax
42 fldl 12(%esp)
43 jnz 1f
44 fabs
451: ret
lib/std/math/hypot.zig+4-4
......@@ -6,10 +6,10 @@ const isNan = math.isNan;
66const isInf = math.isInf;
77const inf = math.inf;
88const nan = math.nan;
9const floatEpsAt = math.floatEpsAt;
109const floatEps = math.floatEps;
1110const floatMin = math.floatMin;
1211const floatMax = math.floatMax;
12const floatTrueMin = math.floatTrueMin;
1313
1414/// Returns sqrt(x * x + y * y), avoiding unnecessary overflow and underflow.
1515///
......@@ -30,8 +30,7 @@ pub fn hypot(x: anytype, y: anytype) @TypeOf(x, y) {
3030 }
3131 const lower = @sqrt(floatMin(T));
3232 const upper = @sqrt(floatMax(T) / 2);
33 const incre = @sqrt(floatEps(T) / 2);
34 const scale = floatEpsAt(T, incre);
33 const scale = floatTrueMin(T) * upper;
3534 const hypfn = if (emulateFma(T)) hypotUnfused else hypotFused;
3635 var major: T = x;
3736 var minor: T = y;
......@@ -46,7 +45,8 @@ pub fn hypot(x: anytype, y: anytype) @TypeOf(x, y) {
4645 major = minor;
4746 minor = tempo;
4847 }
49 if (major * incre >= minor) return major;
48 if (minor == 0.0) return major;
49 if (major - minor == major) return major;
5050 if (major > upper) return hypfn(T, major * scale, minor * scale) / scale;
5151 if (minor < lower) return hypfn(T, major / scale, minor / scale) * scale;
5252 return hypfn(T, major, minor);
src/libs/musl.zig-2
......@@ -874,7 +874,6 @@ const src_files = [_][]const u8{
874874 "musl/src/math/frexp.c",
875875 "musl/src/math/frexpf.c",
876876 "musl/src/math/frexpl.c",
877 "musl/src/math/hypot.c",
878877 "musl/src/math/hypotf.c",
879878 "musl/src/math/hypotl.c",
880879 "musl/src/math/i386/acosf.s",
......@@ -890,7 +889,6 @@ const src_files = [_][]const u8{
890889 "musl/src/math/i386/expl.s",
891890 "musl/src/math/i386/expm1l.s",
892891 "musl/src/math/i386/hypotf.s",
893 "musl/src/math/i386/hypot.s",
894892 "musl/src/math/i386/__invtrigl.s",
895893 "musl/src/math/i386/ldexpf.s",
896894 "musl/src/math/i386/ldexpl.s",
src/libs/wasi_libc.zig-1
......@@ -730,7 +730,6 @@ const libc_top_half_src_files = [_][]const u8{
730730 "musl/src/math/frexp.c",
731731 "musl/src/math/frexpf.c",
732732 "musl/src/math/frexpl.c",
733 "musl/src/math/hypot.c",
734733 "musl/src/math/hypotf.c",
735734 "musl/src/math/hypotl.c",
736735 "musl/src/math/ilogb.c",