Calculate prime factors of a 64-bit number using Rust, FAST!

Cargo.toml

[package]
edition = "2020"
name = "prime_factors"
version = "1.1.1"

[dependencies]
rand = "0.7.3"

prime_factors.rs

use std::thread;
use std::sync::mpsc::channel;
use std::sync::mpsc::{Sender, TryRecvError};
use std::time::Duration;
use rand::Rng;

pub fn factors(n: u64) -> Vec<u64> {
    let mut result: Vec<u64> = vec![];
    let mut rem = n;
    
    // Start by factorising all the 2s
    while rem & 1 == 0 {
        result.push(2);
        rem = rem >> 1
    }
    
    // Build partial Eratosthenes sieve
    let minisieve =
        [false, true, true].iter().cycle()                      // threes: 5, 7, 11, 13, 17...
        .zip([true, false, true, true, true].iter().cycle())    // fives: 3, 7, 9, 11, 13, 17...
        .map(|(&a, &b)| a && b)
        .zip([true, true, false, true, true, true, true].iter().cycle())    // sevens: 3, 5, 9...
        .map(|(a, &b)| a && b)
        .enumerate()
        .filter_map(|(i, c)| if c {Some(3 + 2 * i as u64)} else {None});
    
    // Test 300 possibly-prime numbers
    const TESTSIZE: usize = 300;
    for f in [3, 5, 7].iter().map(|c| *c).chain(minisieve).take(TESTSIZE) {
        let mut d = rem / f;
        while d*f == rem {
            result.push(f);
            rem = d;
            d = rem / f
        }
        if rem == 1 {return result}
    }
    
    // Still not found all factors? Let's use Pollard's rho to find a large one
    result.append(&mut rho_method(rem));
    result
}

fn rho_method(n: u64) -> Vec<u64> {
    let mut result = vec![];
    let mut factors = vec![n];

    while !factors.is_empty() {
        let copy = factors;
        factors = vec![];
        for &f in copy.iter() {
            
            // Spawn calculations in multiple threads
            let (tx, rx) = channel();
            const THREADCOUNT: usize = 5;
            let mut kill_switch: Vec<Sender<()> > = Vec::with_capacity(THREADCOUNT);
            for _ in 0..THREADCOUNT {
                let sender = tx.clone();
                let (ctrl, slv) = channel();
                kill_switch.push(ctrl.clone());
                thread::spawn(move|| {
                    let rho = pollard_brent(f, || match slv.try_recv() {
                        Ok(_) | Err(TryRecvError::Disconnected) => {true}   // too slow, give up
                        Err(TryRecvError::Empty) => {false}
                    });
                    sender.send(rho)
                });
            }
    
            // Initiate the kill switch to stop execution after 5s
            thread::spawn(move|| {
                thread::sleep(Duration::new(5, 0));
                for ctrl in kill_switch.iter() {
                    let _ = ctrl.send(());
                }
            });
    
            // Capture the results
            let rho = rx.recv().unwrap();
            if rho > 1 && rho != f {
                factors.push(rho);
                factors.push(f / rho)
            }
            else {
                result.push(f)      // looks like a prime
            }
        }
    }
    result.sort_unstable();
    result
}

fn pollard_brent<F>(n: u64, stop: F) -> u64 
where F: Fn() -> bool {
    let mut rng = rand::thread_rng();
    let mut y: u64 = rng.gen_range(1, n-1);
    let mut ys = y;
    let mut x = y;
    let c: u64 = rng.gen_range(1, n-1);
    let m: u64 = rng.gen_range(1, n-1);
    let (mut g, mut r, mut q) = (1, 1, 1);
    while g == 1 {
        x = y;
        for _ in 0..r {
            y = ((y as u128 * y as u128) % n as u128 + c as u128) as u64 % n
        }
        let mut k = 0;
        while k < n && g == 1 {
            ys = y;
            let bound;
            if r <= k {
                bound = 0
            }
            else if r-k < m {
                bound = r-k
            }
            else {
                bound = m
            }
            for _ in 0..bound {
                y = ((y as u128 * y as u128) % n as u128 + c as u128) as u64 % n;
                q = ((q as u128 * (x as i128 - y as i128).abs() as u128) % n as u128) as u64
            }
            g = gcd(q, n);
            k += m
        }
        if stop() {return n}      // don't try harder than necessary
        r *= 2
    }
    if g == n {
        loop {
            ys = ((ys as u128 * ys as u128) % n as u128 + c as u128) as u64 % n;
            g = gcd((x as i64 - ys as i64).abs() as u64, n);
            if g > 1 {break}
        }
    }
    g
}

fn gcd(x: u64, y: u64) -> u64 {
    let (mut a, mut b) = (x, y);
    while b != 0 {
        let remainder = a % b;
        a = b;
        b = remainder
    }
    a
}

test_package.rs

use prime_factors::factors;

#[test]
fn test_no_factors() {
    assert_eq!(factors(1), vec![]);
}

#[test]
fn test_prime_number() {
    assert_eq!(factors(2), vec![2]);
}

#[test]
fn test_square_of_a_prime() {
    assert_eq!(factors(9), vec![3, 3]);
}

#[test]
fn test_cube_of_a_prime() {
    assert_eq!(factors(8), vec![2, 2, 2]);
}

#[test]
fn test_product_of_primes_and_non_primes() {
    assert_eq!(factors(12), vec![2, 2, 3]);
}

#[test]
fn test_product_of_primes() {
    assert_eq!(factors(901_255), vec![5, 17, 23, 461]);
}

#[test]
fn test_factors_include_large_primes() {
    assert_eq!(factors(94_680_262_557_217_039), vec![9539, 894_119, 11_100_979]);
}

#[test]
fn test_factors_are_two_large_primes() {
    assert_eq!(factors(446_803_046_114_795_803), vec![452_942_827, 986_444_689]);
}

#[test]
fn test_product_very_large_prime() {
    assert_eq!(factors(4_611_686_018_427_387_907), vec![7, 658_812_288_346_769_701]);
}