]>
Commit | Line | Data |
---|---|---|
8faf50e0 XL |
1 | use core::f64; |
2 | ||
3 | use super::sqrt; | |
4 | ||
5 | const SPLIT: f64 = 134217728. + 1.; // 0x1p27 + 1 === (2 ^ 27) + 1 | |
6 | ||
8faf50e0 XL |
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 | ||
48663c56 | 20 | #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] |
8faf50e0 XL |
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); | |
48663c56 | 73 | z * sqrt(ly + lx + hy + hx) |
8faf50e0 | 74 | } |