| use core::f64; |
| |
| use super::sqrt; |
| |
| const SPLIT: f64 = 134217728. + 1.; // 0x1p27 + 1 === (2 ^ 27) + 1 |
| |
| fn sq(x: f64) -> (f64, f64) { |
| let xh: f64; |
| let xl: f64; |
| let xc: f64; |
| |
| xc = x * SPLIT; |
| xh = x - xc + xc; |
| xl = x - xh; |
| let hi = x * x; |
| let lo = xh * xh - hi + 2. * xh * xl + xl * xl; |
| (hi, lo) |
| } |
| |
| #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] |
| pub fn hypot(mut x: f64, mut y: f64) -> f64 { |
| let x1p700 = f64::from_bits(0x6bb0000000000000); // 0x1p700 === 2 ^ 700 |
| let x1p_700 = f64::from_bits(0x1430000000000000); // 0x1p-700 === 2 ^ -700 |
| |
| let mut uxi = x.to_bits(); |
| let mut uyi = y.to_bits(); |
| let uti; |
| let ex: i64; |
| let ey: i64; |
| let mut z: f64; |
| |
| /* arrange |x| >= |y| */ |
| uxi &= -1i64 as u64 >> 1; |
| uyi &= -1i64 as u64 >> 1; |
| if uxi < uyi { |
| uti = uxi; |
| uxi = uyi; |
| uyi = uti; |
| } |
| |
| /* special cases */ |
| ex = (uxi >> 52) as i64; |
| ey = (uyi >> 52) as i64; |
| x = f64::from_bits(uxi); |
| y = f64::from_bits(uyi); |
| /* note: hypot(inf,nan) == inf */ |
| if ey == 0x7ff { |
| return y; |
| } |
| if ex == 0x7ff || uyi == 0 { |
| return x; |
| } |
| /* note: hypot(x,y) ~= x + y*y/x/2 with inexact for small y/x */ |
| /* 64 difference is enough for ld80 double_t */ |
| if ex - ey > 64 { |
| return x + y; |
| } |
| |
| /* precise sqrt argument in nearest rounding mode without overflow */ |
| /* xh*xh must not overflow and xl*xl must not underflow in sq */ |
| z = 1.; |
| if ex > 0x3ff + 510 { |
| z = x1p700; |
| x *= x1p_700; |
| y *= x1p_700; |
| } else if ey < 0x3ff - 450 { |
| z = x1p_700; |
| x *= x1p700; |
| y *= x1p700; |
| } |
| let (hx, lx) = sq(x); |
| let (hy, ly) = sq(y); |
| z * sqrt(ly + lx + hy + hx) |
| } |
