Division By constants

main.rs

// https://stackoverflow.com/questions/28003334/modulus-optimization-in-c
pub struct MU128 {
    pub magic: u128,
    pub do_add: bool,
    pub shift: i32,
}

fn magic_u128(d: u128) -> MU128 {
    debug_assert!(d > 1);

    const MAX_UNSIGNED: u128 = u128::max_value();
    const MAX_SIGNED: u128 = i128::max_value() as u128;
    const MIN_SIGNED: u128 = i128::min_value() as u128;

    let mut do_add: bool = false;
    let mut p: i32 = 127;
    let nc: u128 = MAX_UNSIGNED - u128::wrapping_neg(d) % d;
    let mut q1: u128 = MIN_SIGNED / nc;
    let mut r1: u128 = MIN_SIGNED - q1 * nc;
    let mut q2: u128 = MAX_SIGNED / d;
    let mut r2: u128 = MAX_SIGNED - q2 * d;
    loop {
        p = p + 1;
        if r1 >= nc - r1 {
            q1 = u128::wrapping_add(u128::wrapping_mul(2, q1), 1);
            r1 = u128::wrapping_sub(u128::wrapping_mul(2, r1), nc);
        } else {
            q1 = u128::wrapping_mul(2, q1);
            r1 = 2 * r1;
        }
        if r2 + 1 >= d - r2 {
            if q2 >= MAX_SIGNED {
                do_add = true;
            }
            q2 = 2 * q2 + 1;
            r2 = u128::wrapping_sub(u128::wrapping_add(u128::wrapping_mul(2, r2), 1), d);
        } else {
            if q2 >= MIN_SIGNED {
                do_add = true;
            }
            q2 = u128::wrapping_mul(2, q2);
            r2 = 2 * r2 + 1;
        }
        let delta: u128 = d - 1 - r2;
        if !(p < 256 && (q1 < delta || (q1 == delta && r1 == 0))) {
            break;
        }
    }

    MU128 {
        magic: q2 + 1,
        do_add,
        shift: p - 128,
    }
}

// https://stackoverflow.com/questions/28868367/getting-the-high-part-of-64-bit-integer-multiplication
fn u128_mulhi(a: u128, b: u128) -> u128 {
    let a_lo = a as u64 as u128;
    let a_hi = (a >> 64);
    let b_lo = b as u64 as u128;
    let b_hi = (b >> 64);

    let ab_hi = a_hi * b_hi;
    let ab_mid = a_hi * b_lo;
    let ab_lo = a_lo * b_lo;

    let m = a_lo * b_hi + (ab_lo >> 64);

    let m_lo = m as u64;

    ab_hi + (m >> 64) + (ab_mid + m_lo as u128 >> 64)
}

fn main() {
    let result = magic_u128(10000000000000000); // 10^16
    println!("magic: {}, add: {}, shift: {}", result.magic, result.do_add, result.shift);
    let mul_result = u128_mulhi(411849172299223580994, result.magic) >> result.shift;
    println!("411849172299223580994/10^16={}", mul_result);
}