machine_factor/

mfactor.rs

use machine_prime::PRIME_TABLE;
/// Deterministic Random Bit Generator (xorshift)
pub const fn drbg(mut x: u64) -> u64 {
    x ^= x.wrapping_shr(12);
    x ^= x.wrapping_shl(25);
    x ^= x.wrapping_shr(27);
    x.wrapping_mul(0x2545F4914F6CDD1D)
}

/// GCD(A,B) where B is odd
#[no_mangle]
pub const extern "C" fn partial_gcd(mut a: u64, mut b: u64) -> u64 {
    if a == 0 {
        return b;
    }

    let self_two_factor = a.trailing_zeros();

    a >>= self_two_factor;

    loop {
        if b > a {
            let interim = b;
            b = a;
            a = interim;
        }
        a -= b;

        if a == 0 {
            return b;
        }
        a >>= a.trailing_zeros();
    }
}

const fn poly_eval(x: u64, subtrahend: u64, inv: u64, n: u64) -> u64 {
    machine_prime::mont_sub(machine_prime::mont_prod(x, x, inv, n), subtrahend, n)
}


// This function only returns prime factors. Ideally you would also return the composite factors
// However in practice it is actually rare to return a composite factor, it is just 
// as likely as not finding any factor at all.
const fn pollard_brent(base: u64, inv: u64, subtrahend: u64, n: u64) -> Option<u64> {
    let m = 128;
    let mut r = 1;
    let mut q = 1;
    let mut g = 1;
    let mut ys = 1;
    let mut y = base;
    let mut x = y;
    let mut cycle = 0;
    //let mut i = 0;

    while cycle < 17 {
        cycle += 1;
        x = y;

        let mut yloop = 0;

        while yloop < r {
            yloop += 1;
            y = poly_eval(y, subtrahend, inv, n);
        }

        let mut k = 0;

        loop {
            let mut i = 0;

            while i < m * cycle {
                if i >= r - k {
                    break;
                }

                y = poly_eval(y, subtrahend, inv, n);
                q = machine_prime::mont_prod(q, x.abs_diff(y), inv, n);
                i += 1;
            } // end loop

            ys = y;
            g = partial_gcd(q, n);
            k += m;
            if k >= r || g != 1 {
                break;
            }
        }

        r <<= 1;
        if g != 1 {
            break;
        }
    }

    if g == n {
        while g == 1 {
            ys = poly_eval(ys, subtrahend, inv,n);
            g = partial_gcd(x.abs_diff(ys), n);
        }
    }
    
    if g != 1 && g != n && machine_prime::is_prime_wc(g) {
        return Some(g);
    }
    None
}


/// Returns some prime factor of an 64-bit integer
///
/// This function uses the Pollard-rho algorithm and mostly exists for FFI
#[no_mangle]
pub const extern "C" fn get_factor(n: u64) -> u64 {
    let inv = machine_prime::mul_inv2(n);
    let one = machine_prime::one_mont(n);
    let base = machine_prime::two_mont(one, n);
    
    match pollard_brent(base, inv, one, n) {
        Some(factor) => return factor,
        None => {
            // if x^2 -1 failed try x^2+1
            // No particular reason except to reuse some values
            let coef = n.wrapping_sub(one);
            match pollard_brent(base, inv, coef, n) {
                Some(factor) =>return factor,
                None => {
                    // Loop that has a roughly 0.5 probability of factoring each iteration
                    let mut param = drbg(n);
                    loop {
                        let rand_base = param % (n - 3) + 3; // We don't bother mapping to montgomery ring
                        match pollard_brent(rand_base, inv, one, n) {
                            Some(factor) => return factor,
                            None => param = drbg(param),
                        }
                    }
                }
            }
        }
    }
}





/// Factorisation of an integer
///
/// Representation of the factors in the form p^k p2^k
#[repr(C)]
pub struct Factorization {
    pub factors: [u64; 15],
    pub powers: [u8; 15],
    pub len: usize,
}

impl Factorization {
    const fn new() -> Self {
        Factorization {
            factors: [0u64; 15],
            powers: [0u8; 15],
            len: 0,
        }
    }
}

const fn revert(n: u64) -> u64{
   let n16 = n as u16;
   let mut est : u16 = 3u16.wrapping_mul(n as u16) ^ 2;
   est = 2u16.wrapping_sub(est.wrapping_mul(n16)).wrapping_mul(est);
   2u16.wrapping_sub(est.wrapping_mul(n16)).wrapping_mul(est) as u64
}

// 
/// Complete factorization of N
#[no_mangle]
pub const extern "C" fn factorize(mut n: u64) -> Factorization {
    let mut t = Factorization::new();

    let mut idx = 0usize;

    if n == 0 {
        return t;
    }
    if n == 1 {
        t.factors[0] = 1;
        t.powers[1] = 1;
        t.len = 1;
        return t;
    }

    let twofactor = n.trailing_zeros();

    if twofactor != 0 {
        t.factors[0] = 2u64;
        t.powers[0] = twofactor as u8;
        n >>= twofactor;
        idx += 1;
    }

    let mut primeidx = 0usize;
    
    
    'trial : while primeidx < 256 && n > 1{
    
     let pinv = PRIME_TABLE[primeidx];
     let bound = PRIME_TABLE[primeidx+1];
     let mut tmp = n;
     let mut prod = n.wrapping_mul(pinv);
        
     if prod <= bound{
       t.factors[idx]=revert(pinv);
       let mut power = 1u8;
       
       loop{
         tmp = prod;
        prod = prod.wrapping_mul(pinv);
         if prod > bound{
           break;
         }
         power+=1;
         if prod == 1{
            break 'trial;
         }
       }
      
       t.powers[idx]=power;
       idx+=1;
     }
     primeidx+=2;
     n=tmp;
    }

    if n == 1 {
        t.len = idx;
        return t;
    }

    if machine_prime::is_prime_wc(n) {
        t.factors[idx] = n;
        t.powers[idx] = 1;
        idx += 1;
        t.len = idx;
        return t;
    }
    while n != 1 {
        let k = get_factor(n);
        t.factors[idx] = k;
        let mut count = 0u8;
        while n % k == 0 {
            count += 1;
            n /= k;
        }
        t.powers[idx] = count;
        idx += 1;
        if n == 1 {
            t.len = idx;
            return t;
        }
        if machine_prime::is_prime_wc(n) {
            t.factors[idx] = n;
            t.powers[idx] = 1;
            idx += 1;
            t.len = idx;
            return t;
        }
    }
    t
}