clock_bigint/

div.rs

//! Division and modulo operations for BigInt.
//!
//! Implements constant-time binary long division algorithm.

use crate::error::{BigIntError, Result};
use crate::limbs::{limb_add_carry, limb_mul_wide, limb_sub_borrow};
use crate::types::BigIntCore;

#[cfg(feature = "alloc")]
use crate::types::BigInt;
#[cfg(feature = "alloc")]
use alloc::vec;

/// Divide two BigInt values, returning quotient and remainder.
///
/// Computes: a = q × b + r where 0 ≤ |r| < |b|
///
/// Uses binary long division algorithm with constant-time operations.
#[cfg(feature = "alloc")]
pub fn div_rem(a: &BigInt, b: &BigInt) -> Result<(BigInt, BigInt)> {
    // Check for division by zero
    if b.is_zero() {
        return Err(BigIntError::DivisionByZero);
    }

    // Handle zero dividend
    if a.is_zero() {
        let max_limbs = a.max_limbs().max(b.max_limbs());
        let zero = BigInt::new(max_limbs);
        return Ok((zero.clone(), zero));
    }

    // Handle signs: division of signed numbers
    let result_sign = a.sign() != b.sign();
    let a_abs = if a.sign() {
        crate::sub::sub(&BigInt::new(a.max_limbs()), a)?
    } else {
        a.clone()
    };
    let b_abs = if b.sign() {
        crate::sub::sub(&BigInt::new(b.max_limbs()), b)?
    } else {
        b.clone()
    };

    // Perform unsigned division
    let (mut quotient, mut remainder) = div_rem_unsigned(&a_abs, &b_abs)?;

    // Apply correct signs
    if result_sign {
        quotient.set_sign(true);
    }
    if a.sign() {
        remainder.set_sign(true);
    }

    quotient.canonicalize()?;
    remainder.canonicalize()?;

    Ok((quotient, remainder))
}

/// Unsigned division: divides two positive BigInt values using Knuth Algorithm D.
fn div_rem_unsigned(a: &BigInt, b: &BigInt) -> Result<(BigInt, BigInt)> {
    let a_limbs = a.limbs();
    let b_limbs = b.limbs();

    // If |a| < |b|, quotient is 0, remainder is a
    if a_limbs.len() < b_limbs.len()
        || (a_limbs.len() == b_limbs.len()
            && a_limbs[a_limbs.len() - 1] < b_limbs[a_limbs.len() - 1])
    {
        let max_limbs = a.max_limbs().max(b.max_limbs());
        let zero = BigInt::new(max_limbs);
        return Ok((zero, a.clone()));
    }

    // For single limb cases, use built-in division
    if a_limbs.len() == 1 && b_limbs.len() == 1 {
        let max_limbs = a.max_limbs().max(b.max_limbs());
        let mut quotient = BigInt::new(max_limbs);
        let mut remainder = BigInt::new(max_limbs);

        let a_val = a_limbs[0];
        let b_val = b_limbs[0];
        let q = a_val / b_val;
        let r = a_val % b_val;

        quotient.limbs_mut()[0] = q;
        remainder.limbs_mut()[0] = r;

        quotient.canonicalize()?;
        remainder.canonicalize()?;
        return Ok((quotient, remainder));
    }

    // Knuth Algorithm D implementation
    let max_limbs = a.max_limbs().max(b.max_limbs());

    // Step D1: Normalize
    // Find d such that d * v_{n-1} >= 2^{63}
    let v_n1 = b_limbs[b_limbs.len() - 1];
    let d_shift = 63 - (63 - v_n1.leading_zeros());
    let d = if d_shift < 64 {
        1u64 << (63 - v_n1.leading_zeros())
    } else {
        1
    };

    // Normalize u and v by multiplying by d
    let mut u_normalized = vec![0u64; a_limbs.len() + 1]; // Extra space for normalization
    let mut v_normalized = vec![0u64; b_limbs.len()];

    // Multiply a by d
    let mut carry = 0u64;
    for i in 0..a_limbs.len() {
        let prod = limb_mul_wide(a_limbs[i], d);
        let (sum, carry1) = limb_add_carry(prod.lo, carry, 0);
        u_normalized[i] = sum;
        carry = prod.hi + carry1;
    }
    if carry > 0 {
        u_normalized[a_limbs.len()] = carry;
    }

    // Multiply b by d
    carry = 0;
    for i in 0..b_limbs.len() {
        let prod = limb_mul_wide(b_limbs[i], d);
        let (sum, carry1) = limb_add_carry(prod.lo, carry, 0);
        v_normalized[i] = sum;
        carry = prod.hi + carry1;
    }

    // Ensure v_normalized has leading bit set (for proper quotient estimation)
    if v_normalized.len() > 0 && (v_normalized[v_normalized.len() - 1] & (1u64 << 63)) == 0 {
        // Shift everything left by 1
        let mut carry = 0u64;
        for i in 0..u_normalized.len() {
            let new_carry = (u_normalized[i] >> 63) & 1;
            u_normalized[i] = (u_normalized[i] << 1) | carry;
            carry = new_carry;
        }
        carry = 0;
        for i in 0..v_normalized.len() {
            let new_carry = (v_normalized[i] >> 63) & 1;
            v_normalized[i] = (v_normalized[i] << 1) | carry;
            carry = new_carry;
        }
    }

    let m = u_normalized.len() - v_normalized.len();
    let n = v_normalized.len();

    // Create quotient
    let mut quotient = BigInt::new(max_limbs);
    let q_limbs = quotient.limbs_mut();

    // Step D2: Initialize j
    for j in (0..=m).rev() {
        // Step D3: Calculate q̂
        // q̂ = (u[j+n] * 2^64 + u[j+n-1]) / v[n-1]
        let u_jn = if j + n < u_normalized.len() {
            u_normalized[j + n]
        } else {
            0
        };
        let u_jn1 = if j + n - 1 < u_normalized.len() {
            u_normalized[j + n - 1]
        } else {
            0
        };
        let v_nm1 = v_normalized[n - 1];

        // Estimate quotient digit
        let mut q_hat = if v_nm1 == 0 {
            u64::MAX
        } else {
            // (u[j+n] * 2^64 + u[j+n-1]) / v[n-1]
            let num = (u_jn as u128) << 64 | (u_jn1 as u128);
            (num / (v_nm1 as u128)) as u64
        };

        // Adjust q̂ if necessary (Knuth's correction)
        if q_hat > 0 {
            let mut prod_carry = 0u64;
            for i in 0..n {
                let prod = limb_mul_wide(q_hat, v_normalized[i]);
                let (sum, carry1) = limb_add_carry(prod.lo, prod_carry, 0);
                prod_carry = prod.hi + carry1;
                if i + j < u_normalized.len() && sum > u_normalized[i + j] {
                    q_hat -= 1;
                    break;
                } else if i + j < u_normalized.len() && sum < u_normalized[i + j] {
                    break;
                }
            }
        }

        // Step D4: Multiply and subtract
        let mut borrow = 0u64;
        for i in 0..n {
            let idx = j + i;
            if idx >= u_normalized.len() {
                continue;
            }

            let prod = limb_mul_wide(q_hat, v_normalized[i]);
            let (diff, borrow_out) = limb_sub_borrow(u_normalized[idx], prod.lo, borrow);
            u_normalized[idx] = diff;
            borrow = borrow_out + prod.hi;
        }

        // Handle borrow from next limb
        let mut idx = j + n;
        while borrow > 0 && idx < u_normalized.len() {
            let (diff, borrow_out) = limb_sub_borrow(u_normalized[idx], 0, borrow);
            u_normalized[idx] = diff;
            borrow = borrow_out;
            idx += 1;
        }

        // Step D5: Test remainder
        if borrow > 0 {
            // Remainder is negative, add back v
            q_hat -= 1;
            let mut carry = 0u64;
            for i in 0..n {
                let idx = j + i;
                if idx >= u_normalized.len() {
                    continue;
                }
                let (sum, carry_out) = limb_add_carry(u_normalized[idx], v_normalized[i], carry);
                u_normalized[idx] = sum;
                carry = carry_out;
            }
        }

        // Step D6: Store quotient digit
        let q_idx = j;
        let limb_idx = q_idx / 64;
        let bit_idx = q_idx % 64;
        if limb_idx < q_limbs.len() {
            q_limbs[limb_idx] |= q_hat << bit_idx;
        }
    }

    // Step D8: Unnormalize remainder
    let mut remainder = BigInt::new(max_limbs);
    let r_limbs = remainder.limbs_mut();

    // Copy remainder and divide by d
    for i in 0..b_limbs.len() {
        if i < u_normalized.len() {
            r_limbs[i] = u_normalized[i];
        }
    }

    // Divide remainder by d (unnormalize)
    if d > 1 {
        let mut rem = 0u64;
        for i in (0..r_limbs.len()).rev() {
            let num = (rem as u128) << 64 | (r_limbs[i] as u128);
            r_limbs[i] = (num / (d as u128)) as u64;
            rem = (num % (d as u128)) as u64;
        }
    }

    quotient.canonicalize()?;
    remainder.canonicalize()?;

    Ok((quotient, remainder))
}

/// Divide two BigInt values, returning only the quotient.
#[cfg(feature = "alloc")]
pub fn div(a: &BigInt, b: &BigInt) -> Result<BigInt> {
    let (q, _) = div_rem(a, b)?;
    Ok(q)
}

/// Compute the remainder of dividing two BigInt values.
#[cfg(feature = "alloc")]
pub fn rem(a: &BigInt, b: &BigInt) -> Result<BigInt> {
    let (_, remainder) = div_rem(a, b)?;
    Ok(remainder)
}

/// In-place division for dynamic BigInt.
#[cfg(feature = "alloc")]
pub fn div_assign(a: &mut BigInt, b: &BigInt) -> Result<()> {
    let result = div(a, b)?;
    *a = result;
    Ok(())
}

/// In-place remainder for dynamic BigInt.
#[cfg(feature = "alloc")]
pub fn rem_assign(a: &mut BigInt, b: &BigInt) -> Result<()> {
    let result = rem(a, b)?;
    *a = result;
    Ok(())
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_div_by_zero() {
        let a = BigInt::from_u64(10, 10);
        let b = BigInt::from_u64(0, 10);
        assert!(div(&a, &b).is_err());
    }

    #[test]
    fn test_div_zero_dividend() {
        let a = BigInt::from_u64(0, 10);
        let b = BigInt::from_u64(10, 10);
        let (q, r) = div_rem(&a, &b).unwrap();
        assert!(q.is_zero());
        assert!(r.is_zero());
    }

    #[test]
    fn test_div_simple() {
        let a = BigInt::from_u64(100, 10);
        let b = BigInt::from_u64(10, 10);
        let (q, r) = div_rem(&a, &b).unwrap();
        assert_eq!(q.limbs()[0], 10);
        assert_eq!(r.limbs()[0], 0);
    }

    #[test]
    fn test_div_with_remainder() {
        let a = BigInt::from_u64(123, 10);
        let b = BigInt::from_u64(10, 10);
        let (q, r) = div_rem(&a, &b).unwrap();
        assert_eq!(q.limbs()[0], 12);
        assert_eq!(r.limbs()[0], 3);
    }
}