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)

}