c/
math_hypot.rs

1// Copyright (c) 2026 Meng Zhuo
2// SPDX-License-Identifier: BSD-3-Clause
3
4// Translated from musl: hypot.c, hypotf.c
5// Uses only integer/bit ops and helper functions from math_helpers.rs.
6// SPLIT constant for double_t (FLT_EVAL_METHOD == 0, standard x86-64).
7
8const HYPOT_SPLIT: f64 = asdouble(0x41A0000002000000); // 2^27 + 1
9
10// ============================================================
11// hypotf (float)
12// ============================================================
13
14const HYPOTF_SCALE_UP: f32 = asfloat(0x6C800000); // 0x1p90f
15const HYPOTF_SCALE_DN: f32 = asfloat(0x12800000); // 0x1p-90f
16
17#[no_mangle]
18pub extern "C" fn hypotf(x: f32, y: f32) -> f32 {
19    let mut ux = asuint(x) & 0x7FFFFFFF;
20    let mut uy = asuint(y) & 0x7FFFFFFF;
21
22    if ux < uy {
23        let ut = ux;
24        ux = uy;
25        uy = ut;
26    }
27
28    let mut x = asfloat(ux);
29    let mut y = asfloat(uy);
30
31    // hypot(inf, nan) == inf
32    if uy == 0xFF << 23 {
33        return y;
34    }
35    if ux >= 0xFF << 23 || uy == 0 || ux.wrapping_sub(uy) >= 25 << 23 {
36        return x + y;
37    }
38
39    let mut z: f32 = 1.0;
40    if ux >= (0x7F + 60) << 23 {
41        z = HYPOTF_SCALE_UP;
42        x *= HYPOTF_SCALE_DN;
43        y *= HYPOTF_SCALE_DN;
44    } else if uy < (0x7F - 60) << 23 {
45        z = HYPOTF_SCALE_DN;
46        x *= HYPOTF_SCALE_UP;
47        y *= HYPOTF_SCALE_UP;
48    }
49    z * sqrtf(((x as f64) * (x as f64) + (y as f64) * (y as f64)) as f32)
50}
51
52// ============================================================
53// hypot (double)
54// ============================================================
55
56const HYPOT_SCALE_UP: f64 = asdouble(0x6BB0000000000000); // 0x1p700
57const HYPOT_SCALE_DN: f64 = asdouble(0x1430000000000000); // 0x1p-700
58
59/// Extended precision square: computes x*x as (hi, lo) where hi+lo ≈ x*x
60/// with much higher precision using Dekker splitting.
61#[inline]
62fn hypot_sq(x: f64) -> (f64, f64) {
63    let xc = x * HYPOT_SPLIT;
64    let xh = x - xc + xc;
65    let xl = x - xh;
66    let hi = x * x;
67    let lo = xh * xh - hi + 2.0 * xh * xl + xl * xl;
68    (hi, lo)
69}
70
71#[no_mangle]
72pub extern "C" fn hypot(x: f64, y: f64) -> f64 {
73    let mut ux = asuint64(x) & 0x7FFFFFFFFFFFFFFF;
74    let mut uy = asuint64(y) & 0x7FFFFFFFFFFFFFFF;
75
76    // arrange |x| >= |y|
77    if ux < uy {
78        let ut = ux;
79        ux = uy;
80        uy = ut;
81    }
82
83    let ex = (ux >> 52) as i32;
84    let ey = (uy >> 52) as i32;
85    let x = asdouble(ux);
86    let y = asdouble(uy);
87
88    // hypot(inf, nan) == inf
89    if ey == 0x7FF {
90        return y;
91    }
92    if ex == 0x7FF || uy == 0 {
93        return x;
94    }
95    // 64 difference is enough for ld80 double_t
96    if ex - ey > 64 {
97        return x + y;
98    }
99
100    // precise sqrt argument in nearest rounding mode without overflow
101    // xh*xh must not overflow and xl*xl must not underflow in sq
102    let mut z: f64 = 1.0;
103    let mut x = x;
104    let mut y = y;
105    if ex > 0x3FF + 510 {
106        z = HYPOT_SCALE_UP;
107        x *= HYPOT_SCALE_DN;
108        y *= HYPOT_SCALE_DN;
109    } else if ey < 0x3FF - 450 {
110        z = HYPOT_SCALE_DN;
111        x *= HYPOT_SCALE_UP;
112        y *= HYPOT_SCALE_UP;
113    }
114
115    let (hx, lx) = hypot_sq(x);
116    let (hy, ly) = hypot_sq(y);
117    z * sqrt(ly + lx + hy + hx)
118}
c/math_hypot.rs

c/
math_hypot.rs
