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 floatEpsAt = math.floatEpsAt;
  9const floatEps = math.floatEps;
 10const floatMin = math.floatMin;
 11const floatMax = math.floatMax;
 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 incre = @sqrt(floatEps(T) / 2);
 33    const scale = floatEpsAt(T, incre);
 34    const hypfn = if (emulateFma(T)) hypotUnfused else hypotFused;
 35    var major: T = x;
 36    var minor: T = y;
 37    if (isInf(major) or isInf(minor)) return inf(T);
 38    if (isNan(major) or isNan(minor)) return nan(T);
 39    if (T == f16) return @floatCast(@sqrt(@mulAdd(f32, x, x, @as(f32, y) * y)));
 40    if (T == f32) return @floatCast(@sqrt(@mulAdd(f64, x, x, @as(f64, y) * y)));
 41    major = @abs(major);
 42    minor = @abs(minor);
 43    if (minor > major) {
 44        const tempo = major;
 45        major = minor;
 46        minor = tempo;
 47    }
 48    if (major * incre >= minor) 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}