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: 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)] |
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 | |