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: f64 = x * x;
16 let lo: f64 = 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}
75