libm/math/
hypot.rs

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

libm/math/
hypot.rs
