dashu-float 0.6.1

use dashu_base::{BitTest, DivRem, Sign};
use dashu_int::{DoubleWord, IBig, UBig, Word};

#[inline]
pub const fn base_as_ibig<const B: Word>() -> IBig {
    IBig::from_parts_const(Sign::Positive, B as DoubleWord)
}

/// Calculate the number of digits in base `B`.
///
/// Returns the integer `k` such that `B^(k-1) <= value < B^k`.
/// If value is `0`, then `k = 0` is returned.
#[inline]
pub fn digit_len<const B: Word>(value: &IBig) -> usize {
    if value.is_zero() {
        return 0;
    };
    if B.is_power_of_two() {
        // `ilog` would also materialize B^log as a heap allocation, only to discard it here
        let bits_per_digit = B.trailing_zeros() as usize;
        (value.bit_len() + bits_per_digit - 1) / bits_per_digit
    } else {
        value.ilog(&UBig::from_word(B)) + 1
    }
}

/// "Left shifting" in given radix, i.e. multiply by a power of radix
#[inline]
pub fn shl_digits<const B: Word>(value: &IBig, exp: usize) -> IBig {
    if exp == 0 {
        return value.clone();
    }

    match B {
        2 => value << exp,
        10 => (value * IBig::from(5).pow(exp)) << exp,
        b if b.is_power_of_two() => value << (exp * b.trailing_zeros() as usize),
        _ => value * base_as_ibig::<B>().pow(exp),
    }
}

#[inline]
pub fn shl_digits_in_place<const B: Word>(value: &mut IBig, exp: usize) {
    if exp != 0 {
        match B {
            2 => *value <<= exp,
            10 => {
                *value *= IBig::from(5).pow(exp);
                *value <<= exp;
            }
            b if b.is_power_of_two() => *value <<= exp * b.trailing_zeros() as usize,
            _ => *value *= base_as_ibig::<B>().pow(exp),
        }
    }
}

// The original Shr operation of IBig follows two's complement representation,
// here we only want to shift the magnitude, so a separate operation is necessary
#[inline]
fn shr_ref(value: &IBig, shift: usize) -> IBig {
    let (sign, words) = value.as_sign_words();
    let n_words = shift / Word::BITS as usize;

    let hi = UBig::from_words(&words[n_words.min(words.len())..]);
    IBig::from_parts(sign, hi >> (shift % Word::BITS as usize))
}

/// "Right shifting" in given radix, i.e. divide by a power of radix
#[inline]
pub fn shr_digits<const B: Word>(value: &IBig, exp: usize) -> IBig {
    if exp == 0 {
        return value.clone();
    }

    match B {
        2 => shr_ref(value, exp),
        10 => shr_ref(value, exp) / IBig::from(5).pow(exp),
        b if b.is_power_of_two() => shr_ref(value, exp * b.trailing_zeros() as usize),
        _ => value / base_as_ibig::<B>().pow(exp),
    }
}

/// Equivalent to value.unsigned_abs().split_bits(n), but returns (hi, lo) and preserving the sign
fn split_bits(value: IBig, n: usize) -> (IBig, IBig) {
    let (sign, mag) = value.into_parts();
    let (lo, hi) = mag.split_bits(n);
    (IBig::from_parts(sign, hi), IBig::from_parts(sign, lo))
}

/// Equivalent to value.unsigned_abs().split_bits(n), but returns (hi, lo) and preserving the sign
fn split_bits_ref(value: &IBig, n: usize) -> (IBig, IBig) {
    debug_assert!(n > 0);
    if value.is_zero() {
        return (IBig::ZERO, IBig::ZERO);
    }

    let (sign, words) = value.as_sign_words();
    let n_words = n / Word::BITS as usize;
    if n_words >= words.len() {
        // shortcut if n is very large
        return (IBig::ZERO, value.clone());
    }

    let mut hi = UBig::from_words(&words[n_words..]);
    hi >>= n % Word::BITS as usize;
    let mut lo = UBig::from_words(&words[..n_words + 1]);
    lo.clear_high_bits(n);

    (IBig::from_parts(sign, hi), IBig::from_parts(sign, lo))
}

/// Same as [split_digits] but take reference of input
#[inline]
pub fn split_digits_ref<const B: Word>(value: &IBig, pos: usize) -> (IBig, IBig) {
    if pos != 0 {
        match B {
            10 => {
                let (q, rem1) = split_bits_ref(value, pos);
                let (q, rem2) = q.div_rem(IBig::from(5).pow(pos));
                let rem = (rem2 << pos) + rem1;
                (q, rem)
            }
            i if i.is_power_of_two() => split_bits_ref(value, pos * i.trailing_zeros() as usize),
            _ => value.div_rem(base_as_ibig::<B>().pow(pos)),
        }
    } else {
        (value.clone(), IBig::ZERO)
    }
}

/// Split the integer at given digit position. Return the high part and low part,
/// and the sign is applied to both parts.
///
/// For example in base 10:
/// * `split_digits(123, 1)` returns `(12, 3)`
/// * `split_digits(-123, 2)` returns `(-1, -23)`
#[inline]
pub fn split_digits<const B: Word>(value: IBig, pos: usize) -> (IBig, IBig) {
    if pos != 0 {
        match B {
            10 => {
                let (q, rem1) = split_bits(value, pos);
                let (q, rem2) = q.div_rem(IBig::from(5).pow(pos));
                let rem = (rem2 << pos) + rem1;
                (q, rem)
            }
            i if i.is_power_of_two() => split_bits(value, pos * i.trailing_zeros() as usize),
            _ => value.div_rem(base_as_ibig::<B>().pow(pos)),
        }
    } else {
        (value, IBig::ZERO)
    }
}

/// If n is a power of base, then return the exponent,
/// otherwise return 0.
///
/// This is a const function replacement of `IBig::ilog`.
pub const fn ilog_exact(n: Word, base: Word) -> u32 {
    if n < base {
        return 0;
    }

    let mut pow = base;
    let mut exp = 1;
    while pow < n {
        pow *= base;
        exp += 1;
    }

    if pow == n {
        exp
    } else {
        0
    }
}

/// `ceil(x) as usize` without `f64::ceil`, which is `std`-only on this crate's
/// MSRV (the `f64` inherent methods only landed in `core` in Rust 1.85). Valid for
/// non-negative `x` within `usize` range.
#[inline]
pub(crate) fn ceil_usize(x: f32) -> usize {
    let i = x as usize;
    if x > i as f32 {
        i + 1
    } else {
        i
    }
}

/// Factor `b` as `newb^a * r`, returning `(a, r)` where `gcd(r, newb) = 1`.
/// Returns `(0, b)` if `newb` does not divide `b`.
pub const fn factor_base(b: Word, newb: Word) -> (usize, Word) {
    let mut a = 0;
    let mut r = b;
    while r % newb == 0 {
        r /= newb;
        a += 1;
    }
    (a, r)
}

// ============================================================================
// log₂ helpers — the fixed-point brackets and bounds behind the base-aware
// magnitude guards (mag.rs's radius scaling, the ×u argument reduction)
// ============================================================================

/// `log₂ BASE` as a `(⌈·, ⌊·⌋, frac_bits)` fixed-point bracket — the same shape as the
/// hand-written base-10 constants it replaces, computed for any base by the classic squaring
/// walk on the normalized significand (one `u128` square per fraction bit). `B` is a const
/// generic, so the whole walk folds into a constant at compile time; the generic-base radius
/// rules cost no more than the base-10 ones. The fraction width shrinks as the integer part
/// grows so the ratio itself stays below `2^63` — that bound is what keeps `|a|·ratio`
/// inside `i128` in `div_scaled` (mag.rs) for every `a` in the `isize` range.
pub const fn log2_base<const B: Word>() -> (u64, u32) {
    // normalize to mant = BASE·2^lz ∈ [2^63, 2^64): ⌊log₂ BASE⌋ = 63 − lz, then walk the
    // fraction bits — square, and each time the value crosses 2 the bit is 1. The walk is
    // pinned to explicit widths (`Word` is `u32` on some targets, whose raw `leading_zeros`
    // would shift the normalization out of `[2^63, 2^64)` and return a wrong log for every
    // non-power-of-two base).
    let int_bits: u32 = Word::BITS - 1 - B.leading_zeros(); // ⌊log₂ BASE⌋ (width-independent)
    let frac_bits: u32 = 62 - int_bits;
    let mut m: u128 = (B as u128) << (63 - int_bits); // ∈ [2^63, 2^64)
    let mut frac: u64 = 0;
    let mut i = frac_bits;
    loop {
        if i == 0 {
            break;
        }
        i -= 1;
        m = (m * m) >> 63;
        if m >> 64 != 0 {
            frac |= 1 << i;
            m >>= 1;
        }
    }
    ((int_bits as u64) << frac_bits | frac, frac_bits)
}

/// A conservative f32 bound of `log₂(u·B^s)` from below: the factors' log2 lower bounds,
/// shaved by the f32 arithmetic's own rounding slack (the `s as f32` conversion, the product
/// and the sum each round at ~2⁻²⁴ relative; 3·10⁻⁷ > 2⁻²² covers them with margin) plus a
/// fixed allowance for the `u` term. Whatever this bound certifies as "below", the true
/// `log₂(u·B^s)` is below too — erring towards the general (always-exact) path, never towards
/// a wrong fast-path verdict.
pub fn log2_u_bs_lb(u_lb: f32, b_lb: f32, s: usize) -> f32 {
    let s_lb = s as f32 * b_lb;
    let raw = u_lb + s_lb;
    raw - (raw * 3e-7 + 1e-4)
}

#[cfg(test)]
mod tests {
    use super::*;
    use dashu_base::{Abs, UnsignedAbs};

    #[test]
    fn test_shr_ref() {
        let a = IBig::from(0x1234567890abcdefu64).pow(12); // 723 bits
        assert_eq!(shr_ref(&a, 10), (&a).abs() >> 10);
        assert_eq!(shr_ref(&a, 100), (&a).abs() >> 100);
        assert_eq!(shr_ref(&a, 1000), (&a).abs() >> 1000);

        let a = IBig::from(-0x1234567890abcdefi64).pow(7); // 422 bits
        assert_eq!(-shr_ref(&a, 10), (&a).abs() >> 10);
        assert_eq!(-shr_ref(&a, 100), (&a).abs() >> 100);
        assert_eq!(-shr_ref(&a, 1000), (&a).abs() >> 1000);
    }

    #[test]
    fn test_split_bits_ref() {
        let a = IBig::from(0x1234567890abcdefu64).pow(12);
        let (hi, lo) = split_bits_ref(&a, 100);
        let (rlo, rhi) = (&a).unsigned_abs().split_bits(100);
        assert_eq!(lo, rlo.into());
        assert_eq!(hi, rhi.into());

        let (hi, lo) = split_bits_ref(&a, 192);
        let (rlo, rhi) = (&a).unsigned_abs().split_bits(192);
        assert_eq!(lo, rlo.into());
        assert_eq!(hi, rhi.into());

        let a = IBig::from(-0x1234567890abcdefi64).pow(7);
        let (hi, lo) = split_bits_ref(&a, 100);
        let (rlo, rhi) = (&a).unsigned_abs().split_bits(100);
        assert_eq!(-lo, rlo.into());
        assert_eq!(-hi, rhi.into());

        let (hi, lo) = split_bits_ref(&a, 192);
        let (rlo, rhi) = (&a).unsigned_abs().split_bits(192);
        assert_eq!(-lo, rlo.into());
        assert_eq!(-hi, rhi.into());
    }
}