//! Eisel-Lemire 算法的实现。

use crate::num::dec2flt::common::BiasedFp;
use crate::num::dec2flt::float::RawFloat;
use crate::num::dec2flt::table::{
    LARGEST_POWER_OF_FIVE, POWER_OF_FIVE_128, SMALLEST_POWER_OF_FIVE,
};

/// 使用扩展精度表示计算浮点数。
///
/// 将浮点数的有效数字和十进制指数快速转换为具有二进制浮点数的扩展表示。
/// 该算法将准确解析绝大多数情况，并使用 128 位表示 (后备 192 位表示)。
///
/// 该算法使用预先计算的 5 次幂按十进制指数缩放指数，并计算表示是否可以明确四舍五入到最近的机器浮点数。
/// 这里不处理接近一半的情况，并由负的、有偏差的二进制指数表示。
///
/// 该算法在第 5 节 "Fast Algorithm" 和第 6 节 "Exact Numbers And Ties" 中的 "Daniel Lemire, Number Parsing at a Gigabyte per Second" 中有详细描述，可在线获取: <https://arxiv.org/abs/2101.11408.pdf>.
///
///
///
///
///
///
///
///
///
pub fn compute_float<F: RawFloat>(q: i64, mut w: u64) -> BiasedFp {
    let fp_zero = BiasedFp::zero_pow2(0);
    let fp_inf = BiasedFp::zero_pow2(F::INFINITE_POWER);
    let fp_error = BiasedFp::zero_pow2(-1);

    // 短路如果值只能是一个字面量 0 或无穷大。
    if w == 0 || q < F::SMALLEST_POWER_OF_TEN as i64 {
        return fp_zero;
    } else if q > F::LARGEST_POWER_OF_TEN as i64 {
        return fp_inf;
    }
    // 规范化我们的有效数字，因此设置了最高有效位。
    let lz = w.leading_zeros();
    w <<= lz;
    let (lo, hi) = compute_product_approx(q, w, F::MANTISSA_EXPLICIT_BITS + 3);
    if lo == 0xFFFF_FFFF_FFFF_FFFF {
        // 如果我们未能用我们的 128 位值近似 w x 5^-q。
        // 由于加 1 可能导致溢出，然后在中点处向上舍入，这可能导致浮点数舍入不当。
        //
        // 然而，这只能在 q ∈ [-27, 55] 时发生。
        // q 的上限是 55，因为 5^55 < 2^128，但是，这只能在 5^q > 2^64 时发生，否则乘积可以用 64 位表示，从而产生精确的结果。
        //
        // 对于负指数，仅当 5^-q < 2^64 时才会发生舍入到偶数。
        //
        // 有关负指数舍入的详细说明，请参见 <https://arxiv.org/pdf/2101.11408.pdf#section.9.1>。
        // 有关正指数舍入的详细说明，请参见 <https://arxiv.org/pdf/2101.11408.pdf#section.8>。
        //
        //
        //
        //
        //
        let inside_safe_exponent = (q >= -27) && (q <= 55);
        if !inside_safe_exponent {
            return fp_error;
        }
    }
    let upperbit = (hi >> 63) as i32;
    let mut mantissa = hi >> (upperbit + 64 - F::MANTISSA_EXPLICIT_BITS as i32 - 3);
    let mut power2 = power(q as i32) + upperbit - lz as i32 - F::MINIMUM_EXPONENT;
    if power2 <= 0 {
        if -power2 + 1 >= 64 {
            // 低于最小指数超过 64 位，必须为 0。
            return fp_zero;
        }
        // 有一个 subnormal 值。
        mantissa >>= -power2 + 1;
        mantissa += mantissa & 1;
        mantissa >>= 1;
        power2 = (mantissa >= (1_u64 << F::MANTISSA_EXPLICIT_BITS)) as i32;
        return BiasedFp { f: mantissa, e: power2 };
    }
    // 需要处理四舍五入的关系。
    // 通常情况下，我们需要四舍五入，但如果我们正好介于两者之间并且我们有一个偶数基础，我们就需要四舍五入。
    //
    //
    // 这只会在以下情况下发生:
    //  1. 128 位表示的低 64 位为 0。
    //      IE, 5^q 适合单个 64 位字。
    //  2. 截断尾数之前的最低有效位是奇数。
    //  3. 移至尾数位 + 1 时截断的所有位均为 0。
    //
    // 或者，我们可能介于两个浮动之间: 我们正好在一半。
    if lo <= 1
        && q >= F::MIN_EXPONENT_ROUND_TO_EVEN as i64
        && q <= F::MAX_EXPONENT_ROUND_TO_EVEN as i64
        && mantissa & 3 == 1
        && (mantissa << (upperbit + 64 - F::MANTISSA_EXPLICIT_BITS as i32 - 3)) == hi
    {
        // 将最低位清零，所以我们不四舍五入。
        mantissa &= !1_u64;
    }
    // 舍入到偶数，然后将有效数字移到位。
    mantissa += mantissa & 1;
    mantissa >>= 1;
    if mantissa >= (2_u64 << F::MANTISSA_EXPLICIT_BITS) {
        // 向上舍入溢出，因此设置进位位。
        // 将尾数设置为 1 (仅设置隐式隐藏位) 并增加指数。
        //
        mantissa = 1_u64 << F::MANTISSA_EXPLICIT_BITS;
        power2 += 1;
    }
    // 将隐藏位清零。
    mantissa &= !(1_u64 << F::MANTISSA_EXPLICIT_BITS);
    if power2 >= F::INFINITE_POWER {
        // 指数高于最大正常值，必须是无穷大。
        return fp_inf;
    }
    BiasedFp { f: mantissa, e: power2 }
}

/// 从十进制指数计算基数为 2 的指数。
/// 这对 log2(10) 使用预先计算的整数近似值，其中 217706/2^16 对于整个非有限十进制指数范围是准确的。
///
///
fn power(q: i32) -> i32 {
    (q.wrapping_mul(152_170 + 65536) >> 16) + 63
}

fn full_multiplication(a: u64, b: u64) -> (u64, u64) {
    let r = (a as u128) * (b as u128);
    (r as u64, (r >> 64) as u64)
}

// 这将计算或近似 w * 5**q 并返回一对近似结果的 64 位字，其中 "high" 部分对应于最高有效位，低部分对应于最低有效位。
//
//
fn compute_product_approx(q: i64, w: u64, precision: usize) -> (u64, u64) {
    debug_assert!(q >= SMALLEST_POWER_OF_FIVE as i64);
    debug_assert!(q <= LARGEST_POWER_OF_FIVE as i64);
    debug_assert!(precision <= 64);

    let mask = if precision < 64 {
        0xFFFF_FFFF_FFFF_FFFF_u64 >> precision
    } else {
        0xFFFF_FFFF_FFFF_FFFF_u64
    };

    // 5^q < 2^64，那么乘法总是提供一个精确的值。
    // 这意味着每当我们需要将关系四舍五入到偶数时，我们总是有一个精确的值。
    //
    let index = (q - SMALLEST_POWER_OF_FIVE as i64) as usize;
    let (lo5, hi5) = POWER_OF_FIVE_128[index];
    // 只要有 1 个零，但在显式尾数位中，只需要一次乘法，+1 表示隐藏位，+1 确定舍入方向，+1 表示计算的乘积是否有前导零。
    //
    //
    //
    let (mut first_lo, mut first_hi) = full_multiplication(w, lo5);
    if first_hi & mask == mask {
        // 需要进行第二次乘法以获得较低乘积的更好精度。这在 q < 55 时总是准确的，因为 5^55 < 2^128。
        //
        // 如果这样包装，那么我们需要将 hi 产品四舍五入。
        //
        let (_, second_hi) = full_multiplication(w, hi5);
        first_lo = first_lo.wrapping_add(second_hi);
        if second_hi > first_lo {
            first_hi += 1;
        }
    }
    (first_lo, first_hi)
}