machine_factor/

wfactor.rs

use crate::mfactor::{drbg, get_factor};
use machine_prime::PRIME_TABLE_128;

/// GCD(A,B) for 128-bit integers, where B is odd
#[no_mangle]
pub const extern "C" fn partial_gcd_128(mut a: u128, mut b: u128) -> u128 {
    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;
        }
        if b.leading_zeros() != a.leading_zeros(){
         a%=b;  
        }
        else{
        a -= b;
        }
        if a == 0 {
            return b;
        }
        a >>= a.trailing_zeros();
    }
}

const fn poly_eval_128(x: u128, subtrahend: u128, inv: u128, n: u128) -> u128 {
    machine_prime::mont_sub_128(machine_prime::mont_sqr_128(x, inv, n), subtrahend, n)
}

const fn pollard_brent_128(base: u128, inv: u128, subtrahend: u128, n: u128) -> Option<u128> {
    let m = 512;
    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;
 
    while cycle < 33 {
        cycle += 1;
        x = y;

        let mut yloop = 0;

        while yloop < r {
            yloop += 1;
            y = poly_eval_128(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_128(y, subtrahend, inv, n);
                q = machine_prime::mont_prod_128(q, x.abs_diff(y), inv, n);
                i += 1;
            } // end loop

            ys = y;
            g = partial_gcd_128(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_128(ys, subtrahend, inv, n);
            g = partial_gcd_128(x.abs_diff(ys), n);
        }
    }
    if g != 1 && g != n && machine_prime::is_prime_wc_128(g) {
        return Some(g);
    }
    None
}

/// Returns some prime factor of an 128-bit integer
#[no_mangle]
pub  const extern "C" fn get_factor_128(n: u128) -> u128 {
    if n < 1u128 << 64 {
        return get_factor(n as u64) as u128;
    }

    let inv = machine_prime::mul_inv2_128(n);
    let one = machine_prime::one_mont_128(n);
    let base = machine_prime::two_mont_128(one, n);

    match pollard_brent_128(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_128(base, inv, coef, n) {
                Some(factor) => return factor,
                None => {
                    // Loop that has a roughly 0.5 probability of factoring each iteration
                    // The probability of being unable to factor a composite is 1/2^period and the period is likely very large
                    // So the risk of error is vanishingly small
                    let mut param = drbg(n as u64);
                    loop {
                        let rand_base = (param as u128) % (n - 3) + 3;
                        match pollard_brent_128(rand_base, inv, one, n) {
                            Some(factor) => return factor,
                            None => param = drbg(param),
                        }
                    }
                }
            }
        }
    }
}

/// Factorization of a 128-bit integer
#[repr(C)]
pub struct Factorization128 {
    pub factors: [u128; 26],
    pub powers: [u8; 26],
    pub len: usize,
}

impl Factorization128 {
    const fn new() -> Self {
        Factorization128 {
            factors: [0u128; 26],
            powers: [0u8; 26],
            len: 0,
        }
    }
}

const fn revert(n: u128) -> u128{
   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 u128
}


/// Complete factorization of 128-bit N
#[no_mangle]
pub const extern "C" fn factorize_128(mut n: u128) -> Factorization128 {
    if n < 1u128 << 64 {
        let mut zero = Factorization128::new();

        let t = crate::mfactor::factorize(n as u64);

        let mut i = 0;

        while i < t.len {
            zero.powers[i] = t.powers[i];
            zero.factors[i] = t.factors[i] as u128;
            i += 1;
        }
        zero.len = t.len;
        return zero;
    }

    let mut t = Factorization128::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] = 2u128;
        t.powers[0] = twofactor as u8;
        n >>= twofactor;
        idx += 1;
    }
    
    let mut primeidx = 0usize;
    
    
    'trial : while primeidx < 256 && n > 1{
    
     let pinv = PRIME_TABLE_128[primeidx];
     let bound = PRIME_TABLE_128[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_128(n) {
        t.factors[idx] = n;
        t.powers[idx] = 1;
        idx += 1;
        t.len = idx;
        return t;
    }
    while n != 1 {
        let k = get_factor_128(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_128(n) {
            t.factors[idx] = n;
            t.powers[idx] = 1;
            idx += 1;
            t.len = idx;
            return t;
        }
    }
    t
}