Montgomery modular multiplication

monty.rs

// Cut from bigint.rs
pub struct MontyReducer<'a> {
    p: &'a BigUint,
    n: Vec<u32>,
    n0inv: u64
}
fn inv_mod_u32(num: u32) -> u64 {
    assert!( num % 2 != 0 );
    let mut x = 0;
    let mut y = 1;
    let mut u = 1;
    let mut v = 0;

    let mut a: i64 = num as i64;
    let mut b: i64 = (u32::max_value() as i64)+1;
    let mx = b;

    while a != 0 {
        let q = b / a;
        let r = b % a;
        let m = x - u*q;
        let n = y - v*q;
        b = a; a = r;
        x = u; y = v;
        u = m; v = n;
    }

    assert!(b == 1);
    if x < 0 {
        y = y-(num as i64);
        (x + mx) as u64
    } else {
        x as u64
    }
}
impl<'a> MontyReducer<'a> {
    pub fn new(p: &'a BigUint) -> Self {
        let n : Vec<u32> = p.data.clone();
        let n0inv = inv_mod_u32(n[0]);
        MontyReducer { p: p, n: n, n0inv: n0inv }
    }
}
pub fn monty_redc(a: BigUint, mr: &MontyReducer) -> BigUint {
    let mut t = a.data;
    let n = &mr.n;
    let n_size = n.len();
    let old_size = t.len();

    //Give it enough work space
    t.reserve(2*n_size+2-old_size);
    t.extend(repeat(0).take(2*n_size+2-old_size));

    let mask = u32::max_value() as u64;
    let npri = (mask-mr.n0inv)+1;

    for i in 0..n_size {
        let mut c = 0;
        let m = ((t[i] as u64)*npri)&mask;
        for j in 0..n_size {
            let x = (t[i+j] as u64)+m*(n[j] as u64)+c;
            t[i+j] = (x & mask) as u32;
            c = x>>32;
        }
        for j in n_size..2*n_size-i+2 {
            let x = (t[i+j] as u64)+c;
            t[i+j] = (x & mask) as u32;
            c = x>>32;
        }
    }
    let k : Vec<u32> = t.iter().skip(n_size).cloned().collect();
    let ret = BigUint::new(k);
    if &ret < mr.p {
        ret
    } else {
        &ret-mr.p
    }
}
pub fn pow_mod_monty(a: &BigUint, exp: &BigUint, mr: &MontyReducer) -> BigUint{
    let mut r : BigUint = One::one();
    while &r < mr.p {
        r = r << 32;
    }
    let mut e = exp.clone();
    let mut apri = a*&r%mr.p;
    let mut ans = &r%mr.p;
    let zero = Zero::zero();
    while e > zero {
        if e.is_odd() {
            let tmpans = ans*&apri;
            ans = monty_redc(tmpans, mr);
            //assert!(&ans < p)
        }
        apri = pow(apri, 2);
        apri = monty_redc(apri, mr);
        //assert!(&apri < p);
        e = e >> 1;

    }

    monty_redc(ans, mr)
}