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