const-frac 0.0.4

use core::mem;
use gcd;

pub const fn reduce(a: u128, b: u128) -> (u128, u128) {
    let gcd = gcd::binary_u128(a, b);

    (a / gcd, b / gcd)
}

pub const fn min(a: u8, b: u8) -> u8 {
    if a > b { b } else { a }
}

pub const fn bitlen(i: u128) -> u8 {
    (u128::BITS - i.leading_zeros()) as u8
}

pub const fn bitlen64(i: u64) -> u8 {
    (u64::BITS - i.leading_zeros()) as u8
}

const EXP_LEN: u8 = (u64::BITS - f64::MANTISSA_DIGITS) as u8;
pub const BIAS: u16 = (1 << EXP_LEN as u16 - 1) - 2;
//const EXP_NAN: u16 = !(!0 << EXP_LEN);

pub const fn encode_to_f64(neg: bool, mantissa: u64, exp2: i16) -> f64 {
    let sign = if neg { 1 << u64::BITS - 1 } else { 0 };
    let exp2 = exp2 + BIAS as i16;
    let exp2 = (exp2 as u64) << f64::MANTISSA_DIGITS >> 1;
    let mantissa = mantissa & !0 >> EXP_LEN + 1;

    unsafe { mem::transmute(sign | exp2 | mantissa) }
}

pub const fn decode_f64(d: f64) -> (bool, u64, i16) {
    let bits: u64 = unsafe { mem::transmute(d) };
    let neg = bits & 1 << u64::BITS - 1 != 0;
    let exp2 = (bits & !0 << f64::MANTISSA_DIGITS >> 1) >> f64::MANTISSA_DIGITS - 1;
    let mantissa = bits & !0 >> EXP_LEN + 1;

    (
        neg,
        if exp2 == 0 {
            mantissa
        } else {
            mantissa | 1 << f64::MANTISSA_DIGITS
        },
        exp2 as i16 - BIAS as i16
    )
}

pub const fn const_div(numer: u128, denom: u128, neg: bool, exp2: i16) -> f64 {
    match (numer, denom) {
        (0, 0) => f64::NAN,
        (_, 0) => if neg { f64::NEG_INFINITY } else { f64::INFINITY },
        (0, _) => 0f64,
        _ => div_impl(numer, denom, neg, exp2),
    }
}

const MLEN: u8 = f64::MANTISSA_DIGITS as u8 + 1;

/// return (mantissa, remain, rest bits, fpos)
const fn first_div(numer: u128, denom: u128) -> (u64, u128, u8, i16) {
    match numer / denom {
        0 => {
            let nbits = bitlen(numer);
            let dbits = bitlen(denom);
            let n = dbits - nbits; // >= 0
            let numer = numer << n;
            let (r, fpos) = match numer.checked_sub(denom) {
                Some(r) => (r, MLEN + n - 1),
                None => {
                    let numer = numer << 1;

                    match numer.overflowing_sub(denom) {
                        (d, true) => (-(d as i128) as u128, MLEN - 1),
                        (d, false) => (d, MLEN + n),
                    }
                },
            };

            (1u64 << MLEN - 1, r, MLEN - 1, fpos as i16)
        },
        q => {
            let r = numer - q * denom;

            match MLEN.overflowing_sub(bitlen(q)) {
                (d, true) => {
                    let exceeded = -(d as i8) as i16; // > 0
                    ((q >> exceeded) as u64, r, 0, -exceeded)
                },
                (d, false) => ((q << d) as u64, r, d, d as i16),
            }
        },
    }
}

const fn round(mut mantissa: u64) -> (u64, bool) {
    let mut overflow = false;

    if mantissa & 1 == 1 {
        mantissa += 1;
        if bitlen64(mantissa) > MLEN {
            mantissa >>= 1;
            overflow = true;
        }
    }
    (mantissa >> 1, overflow)
}

const fn fill_qbits(mut r: u128, denom: u128, mut mantissa: u64, mut qbits: u8) -> u64 {
    assert!(r < denom);
    while qbits > 0 {
        let _hi = r >> u128::BITS / 2 ;
        let n0 = r.leading_zeros() as u8;
        let n = min(qbits, n0);
        let numer = r << n;
        let q = match numer / denom {
            0 => {
                if n == qbits {
                    break;
                }
                r = match (numer << 1).overflowing_sub(denom) {
                    (diff, true) => diff,
                    (_diff, false) => unreachable!(),
                };
                qbits -= if n == 0 { 2 } else { n + 1 };
                1
            },
            q => {
                assert!(n != 0);
                r = numer - q * denom;
                qbits -= n;
                q
            },
        };

        mantissa |= (q << qbits) as u64;
    }
    mantissa
}

/// numer, denom have to be already reduced.
/// numer/denom must be binary repeating decimal.
const fn div_impl(numer: u128, denom: u128, neg: bool, mut exp2: i16) -> f64 {
    let (mantissa, r, qbits, fpos) = first_div(numer, denom);
    let mantissa =  if r == 0 {
        mantissa
    } else {
        fill_qbits(r, denom, mantissa, qbits)
    };
    let (mantissa, overflow) = round(mantissa);
    let fpos = fpos - 1 - overflow as i16;

    exp2 += MLEN as i16 - 1 - fpos ;
    if exp2 > f64::MAX_EXP as i16 {
        if neg {
            f64::NEG_INFINITY
        } else {
            f64::INFINITY
        }
    } else if exp2 < f64::MIN_EXP as i16 {
        let n = f64::MIN_EXP as i16 - exp2;

        encode_to_f64(neg, mantissa >> n, 0)
    } else {
        encode_to_f64(neg, mantissa as u64, exp2)
    }
}

#[cfg(test)]
mod tests {
    extern crate std;
    use super::*;
    use num::Float;

    const NUMER: u128 = 1 << 63;
    const DENOM: u128 = 1;
    const QUOT: f64 = const_div(NUMER, DENOM, false, 0);

    #[test]
    fn test_const() {
        let left = QUOT.integer_decode();
        let right = (NUMER as f64).integer_decode();

        assert_eq!(left, right, "\n{:b}\n{:b}\n", left.0, right.0);
    }

    #[test]
    fn test_simpl_div() {
        let q = const_div(5, 9, false, 0);
        let left = q.integer_decode();
        let right = (5f64 / 9f64).integer_decode();

        assert_eq!(left, right, "\n{:b}\n{:b}\n", left.0, right.0);
    }

    const PAI: u128 = 3_14159_26535_89793_23846_26433_83279_50288;
    const EX: u128 = 1_00000_00000_00000_00000_00000_00000_00000;

    #[test]
    fn test_pi() {
        let q = const_div(EX, PAI, false, 0);
        let left = q.integer_decode();
        let right = (1f64 / std::f64::consts::PI).integer_decode();

        assert_eq!(left, right, "\n{:b}\n{:b}\n", left.0, right.0);
    }

    #[test]
    #[ignore]
    fn test_inv_pi() {
        let q = const_div(
            5u128.pow(38),
            0xEC58_DFA7_4641_AF52_AD0D_16E7_7D57_6623,
            false,
            38
        ); // x10-38
        let left = q.integer_decode();
        let right = (1f64 / std::f64::consts::PI).integer_decode();

        assert_eq!(left, right, "\n{:b}, {}\n{:b}, {}\n", left.0, left.1, right.0, right.1);
    }

    #[test]
    fn test_small() {
        let q = const_div(1, EX, false, 0);
        //let q = const_div(1, 5421010862427522u128 << 64, false, 0);
        let left = q.integer_decode();
        let right = (1f64 / EX as f64).integer_decode();

        assert_eq!(left, right,
            "\n{:b}\n{:b}\n{:b}\n",
            left.0,
            right.0,
            (EX as f64).integer_decode().0
        );
    }

    #[test]
    fn test_big_div() {
        const DENOM: u128 = !0;
        const NUMER: u128 = DENOM - 2;

        let q = const_div(NUMER, DENOM, false, 0);
        let left = q.integer_decode();
        let right = (NUMER as f64 / DENOM as f64).integer_decode();

        assert_eq!(left, right, "\n{:b}\n{:b}\n", left.0, right.0);
    }
}