fat muldiv

muldiv.rs

fn lo128(a: u128) -> u128 {
    a & ((1<<64)-1)
}

/// Returns the most significant 64 bits of a
fn hi128(a: u128) -> u128 {
    a >> 64
}

/// Returns 2^128 - a (two's complement)
fn neg128(a: u128) -> u128 {
    (!a).wrapping_add(1)
}

/// Returns 2^128 / a
fn div128(a: u128) -> Option<u128> {
    Some((neg128(a).checked_div(a)?).wrapping_add(1))
}

/// Returns 2^128 % a
fn mod128(a: u128) -> Option<u128> {
    neg128(a).checked_rem(a)
}

/// Returns a+b (wrapping)
fn add256(a0: u128, a1: u128, b0: u128, b1: u128) -> (u128, u128) {
    let (r0, overflow) = a0.overflowing_add(b0);
    let r1 = a1.wrapping_add(b1).wrapping_add(overflow as u128);
    (r0, r1)
}

/// Returns a*b (in 256 bits)
fn mul128(a: u128, b: u128) -> (u128, u128) {
    // Split a and b into hi and lo 64-bit parts
    // a*b = (a1<<64 + a0)*(b1<<64 + b0)
    // = (a1*b1)<<128 + (a1*b0 + b1*a0)<<64 + a0*b0
    let (a0, a1) = (lo128(a), hi128(a));
    let (b0, b1) = (lo128(b), hi128(b));
    let (x, y) = (a1*b0, b1*a0);
    
    let (r0, r1) = (a0*b0, a1*b1);
    let (r0, r1) = add256(r0, r1, lo128(x)<<64, hi128(x));
    let (r0, r1) = add256(r0, r1, lo128(y)<<64, hi128(y));
    (r0, r1)
}

/// Returns a/b
fn div256_128(mut a0: u128, mut a1: u128, b: u128) -> Option<(u128, u128)> {
    if b == 1 {
        return Some((a0, a1));
    }
    // Calculate a/b
    // = (a0<<128 + a1)/b
    //   let (q, r) = (div128(b), mod128(b));
    // = (a1*(q*b + r)) + a0)/b
    // = (a1*q*b + a1*r + a0)/b
    // = (a1*r + a0)/b + a1*q
    let (q, r) = (div128(b)?, mod128(b)?);
    
    // x = current result
    // a = next number
    // t = temp
    let (mut x0, mut x1) = (0, 0);
    while a1 != 0 {
        // x += a1*q
        let (t0, t1) = mul128(a1, q);
        let (new_x0, new_x1) = add256(x0, x1, t0, t1);
        x0 = new_x0;
        x1 = new_x1;
        // a = a1*r + a0
        let (t0, t1) = mul128(a1, r);
        let (new_a0, new_a1) = add256(t0, t1, a0, 0);
        a0 = new_a0;
        a1 = new_a1;
    }
    
    Some(add256(x0, x1, a0.checked_div(b)?, 0))
}

/// Returns a*b/c (wrapping to 128 bits)
pub fn muldiv128(a: u128, b: u128, c: u128) -> Option<u128> {
    let (t0, t1) = mul128(a, b);
    Some(div256_128(t0, t1, c)?.0)
}