blob: da458ea1d05f16831c6b90367445c92a3c4ba0e2 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
|
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)
}
|
