hypot.rs - source
Description: Source of the Rust file `/Users/erlendbasso/.cargo/registry/src/github.com-1ecc6299db9ec823/libm-0.2.6/src/math/hypot.rs`.
Keywords: rust, rustlang, rust-lang

    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)
    }