libm/math/
hypot.rs

1use core::f64;
2
3use super::sqrt;
4
5const SPLIT: f64 = 134217728. + 1.; // 0x1p27 + 1 === (2 ^ 27) + 1
6
7fn sq(x: f64) -> (f64, f64) {
8    let xh: f64;
9    let xl: f64;
10    let xc: f64;
11
12    xc = x * SPLIT;
13    xh = x - xc + xc;
14    xl = x - xh;
15    let hi = x * x;
16    let lo = xh * xh - hi + 2. * xh * xl + xl * xl;
17    (hi, lo)
18}
19
20#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
21pub fn hypot(mut x: f64, mut y: f64) -> f64 {
22    let x1p700 = f64::from_bits(0x6bb0000000000000); // 0x1p700 === 2 ^ 700
23    let x1p_700 = f64::from_bits(0x1430000000000000); // 0x1p-700 === 2 ^ -700
24
25    let mut uxi = x.to_bits();
26    let mut uyi = y.to_bits();
27    let uti;
28    let ex: i64;
29    let ey: i64;
30    let mut z: f64;
31
32    /* arrange |x| >= |y| */
33    uxi &= -1i64 as u64 >> 1;
34    uyi &= -1i64 as u64 >> 1;
35    if uxi < uyi {
36        uti = uxi;
37        uxi = uyi;
38        uyi = uti;
39    }
40
41    /* special cases */
42    ex = (uxi >> 52) as i64;
43    ey = (uyi >> 52) as i64;
44    x = f64::from_bits(uxi);
45    y = f64::from_bits(uyi);
46    /* note: hypot(inf,nan) == inf */
47    if ey == 0x7ff {
48        return y;
49    }
50    if ex == 0x7ff || uyi == 0 {
51        return x;
52    }
53    /* note: hypot(x,y) ~= x + y*y/x/2 with inexact for small y/x */
54    /* 64 difference is enough for ld80 double_t */
55    if ex - ey > 64 {
56        return x + y;
57    }
58
59    /* precise sqrt argument in nearest rounding mode without overflow */
60    /* xh*xh must not overflow and xl*xl must not underflow in sq */
61    z = 1.;
62    if ex > 0x3ff + 510 {
63        z = x1p700;
64        x *= x1p_700;
65        y *= x1p_700;
66    } else if ey < 0x3ff - 450 {
67        z = x1p_700;
68        x *= x1p700;
69        y *= x1p700;
70    }
71    let (hx, lx) = sq(x);
72    let (hy, ly) = sq(y);
73    z * sqrt(ly + lx + hy + hx)
74}
libm/math/hypot.rs

libm/math/
hypot.rs
