1const std = @import("../std.zig");
2const math = std.math;
3const expect = std.testing.expect;
4const isNan = math.isNan;
5const isInf = math.isInf;
6const inf = math.inf;
7const nan = math.nan;
8const floatEps = math.floatEps;
9const floatMin = math.floatMin;
10const floatMax = math.floatMax;
11const floatTrueMin = math.floatTrueMin;
12
13/// Returns sqrt(x * x + y * y), avoiding unnecessary overflow and underflow.
14///
15/// Special Cases:
16///
17/// | x | y | hypot |
18/// |-------|-------|-------|
19/// | +-inf | any | +inf |
20/// | any | +-inf | +inf |
21/// | nan | fin | nan |
22/// | fin | nan | nan |
23pub fn hypot(x: anytype, y: anytype) @TypeOf(x, y) {
24 const T = @TypeOf(x, y);
25 switch (@typeInfo(T)) {
26 .float => {},
27 .comptime_float => return @sqrt(x * x + y * y),
28 else => @compileError("hypot not implemented for " ++ @typeName(T)),
29 }
30 const lower = @sqrt(floatMin(T));
31 const upper = @sqrt(floatMax(T) / 2);
32 const scale = floatTrueMin(T) * upper;
33 const hypfn = if (emulateFma(T)) hypotUnfused else hypotFused;
34 var major: T = x;
35 var minor: T = y;
36 if (isInf(major) or isInf(minor)) return inf(T);
37 if (isNan(major) or isNan(minor)) return nan(T);
38 if (T == f16) return @floatCast(@sqrt(@mulAdd(f32, x, x, @as(f32, y) * y)));
39 if (T == f32) return @floatCast(@sqrt(@mulAdd(f64, x, x, @as(f64, y) * y)));
40 major = @abs(major);
41 minor = @abs(minor);
42 if (minor > major) {
43 const tempo = major;
44 major = minor;
45 minor = tempo;
46 }
47 if (minor == 0.0) return major;
48 if (major - minor == major) return major;
49 if (major > upper) return hypfn(T, major * scale, minor * scale) / scale;
50 if (minor < lower) return hypfn(T, major / scale, minor / scale) * scale;
51 return hypfn(T, major, minor);
52}
53
54inline fn emulateFma(comptime T: type) bool {
55 // If @mulAdd lowers to the software implementation,
56 // hypotUnfused should be used in place of hypotFused.
57 // This takes an educated guess, but ideally we should
58 // properly detect at comptime when that fallback will
59 // occur.
60 return (T == f128 or T == f80);
61}
62
63inline fn hypotFused(comptime F: type, x: F, y: F) F {
64 const r = @sqrt(@mulAdd(F, x, x, y * y));
65 const rr = r * r;
66 const xx = x * x;
67 const z = @mulAdd(F, -y, y, rr - xx) + @mulAdd(F, r, r, -rr) - @mulAdd(F, x, x, -xx);
68 return r - z / (2 * r);
69}
70
71inline fn hypotUnfused(comptime F: type, x: F, y: F) F {
72 const r = @sqrt(x * x + y * y);
73 if (r <= 2 * y) { // 30deg or steeper
74 const dx = r - y;
75 const z = x * (2 * dx - x) + (dx - 2 * (x - y)) * dx;
76 return r - z / (2 * r);
77 } else { // shallower than 30 deg
78 const dy = r - x;
79 const z = 2 * dy * (x - 2 * y) + (4 * dy - y) * y + dy * dy;
80 return r - z / (2 * r);
81 }
82}
83
84const hypot_test_cases = .{
85 .{ 0.0, -1.2, 1.2 },
86 .{ 0.2, -0.34, 0.3944616584663203993612799816649560759946493601889826495362 },
87 .{ 0.8923, 2.636890, 2.7837722899152509525110650481670176852603253522923737962880 },
88 .{ 1.5, 5.25, 5.4600824169603887033229768686452745953332522619323580787836 },
89 .{ 37.45, 159.835, 164.16372840856167640478217141034363907565754072954443805164 },
90 .{ 89.123, 382.028905, 392.28687638576315875933966414927490685367196874260165618371 },
91 .{ 123123.234375, 529428.707813, 543556.88524707706887251269205923830745438413088753096759371 },
92};
93
94test hypot {
95 try expect(hypot(0.3, 0.4) == 0.5);
96}
97
98test "hypot.correct" {
99 inline for (.{ f16, f32, f64, f128 }) |T| {
100 inline for (hypot_test_cases) |v| {
101 const a: T, const b: T, const c: T = v;
102 try expect(math.approxEqRel(T, hypot(a, b), c, @sqrt(floatEps(T))));
103 }
104 }
105}
106
107test "hypot.precise" {
108 inline for (.{ f16, f32, f64 }) |T| { // f128 seems to be 5 ulp
109 inline for (hypot_test_cases) |v| {
110 const a: T, const b: T, const c: T = v;
111 try expect(math.approxEqRel(T, hypot(a, b), c, floatEps(T)));
112 }
113 }
114}
115
116test "hypot.special" {
117 @setEvalBranchQuota(2000);
118 inline for (.{ f16, f32, f64, f128 }) |T| {
119 try expect(math.isNan(hypot(nan(T), 0.0)));
120 try expect(math.isNan(hypot(0.0, nan(T))));
121
122 try expect(math.isPositiveInf(hypot(inf(T), 0.0)));
123 try expect(math.isPositiveInf(hypot(0.0, inf(T))));
124 try expect(math.isPositiveInf(hypot(inf(T), nan(T))));
125 try expect(math.isPositiveInf(hypot(nan(T), inf(T))));
126
127 try expect(math.isPositiveInf(hypot(-inf(T), 0.0)));
128 try expect(math.isPositiveInf(hypot(0.0, -inf(T))));
129 try expect(math.isPositiveInf(hypot(-inf(T), nan(T))));
130 try expect(math.isPositiveInf(hypot(nan(T), -inf(T))));
131 }
132}