1 use core::f64;
2 
3 use super::sqrt;
4 
5 const SPLIT: f64 = 134217728. + 1.; // 0x1p27 + 1 === (2 ^ 27) + 1
6 
7 fn 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)]
21 pub 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 }
75