miller-rabin.rs

#![feature(step_by)]
use std::io::{BufReader,BufRead};
use std::fs::File;
use std::cmp::Ordering;
use std::collections::{BinaryHeap, HashSet};
use std::thread;
use std::sync::mpsc::channel;

/*
 * Author: Benjamin Brittain
 * Date: 11/23/2015
 *
 * Finds the largest error probability of the Miller-Rabin primality
 * testing algorithm for odd integers between 90001 and 100000
 *
 * Most recent output of miller-rabin.rs
 *
 * time ./target/release/miller-rabin
 * 97921: 12.408% of bases are strong liars | prime? false
 * 90751: 12.397% of bases are strong liars | prime? false
 * 96049: 9.914% of bases are strong liars  | prime? false
 * 96139: 8.257% of bases are strong liars  | prime? false
 * 94657: 5.496% of bases are strong liars  | prime? false
 * 96727: 4.965% of bases are strong liars  | prime? false
 * 91001: 3.571% of bases are strong liars  | prime? false
 * 92509: 3.288% of bases are strong liars  | prime? false
 * 98671: 3.083% of bases are strong liars  | prime? false
 * ./target/release/miller-rabin  116.41s user 0.31s system 792% cpu 14.721 total
 */

fn main() {
    // Load a list of prime numbers into a vector
    // for checking primality
    let file = File::open("10000.txt").unwrap();
    let mut prime_set = HashSet::new();
    for line in BufReader::new(file).lines() {
        let l = line.unwrap();
        let vs: Vec<&str> = l.split_whitespace().collect();
        for x in vs {
            prime_set.insert(x.parse::<u64>().unwrap());
        }
    }
    let mut int_heap = BinaryHeap::new();
    // thread communication channel
    let (tx, rx) = channel();
    // check every number between 90001 and 100000 in increments of two
    for potential_prime in (90001..100000).step_by(2) {
        let tx = tx.clone();
        // spawn a new thread for each prime, speeds code up by a factor of 6
        thread::spawn(move|| {
            let mut mr_prime = 0;
            // check every number between 0 and the number being tested for primality
            // to see if it is false witness
            for base in 0..potential_prime {
                if miller_rabin(potential_prime, base) {
                    mr_prime += 1
                }
            }
            tx.send((potential_prime, mr_prime)).unwrap();
        });
    }

    // get the response from each thread, calculate the error and add it to the sorted heap
    for _ in 0 .. ((100000 - 90001) / 2) {
        let (potential_prime, mr_prime) = rx.recv().unwrap();
        // calculate the number of incorrect bases, then push them onto a binary heap
        // sorted by the error rate
        let error;
        if prime_set.contains(&potential_prime) {
            error = 1.0 - (mr_prime as f32 / (potential_prime - 2) as f32);
            assert!((error as u16) == 0u16 , "error = {}, must be 0", error);
            // floating points are the worst, but error for primes must be 0.
        } else {
            error = mr_prime as f32 / (potential_prime - 2) as f32;
        }
        int_heap.push(IntWrapper { number: potential_prime,
            error: error,
            prime: prime_set.contains(&potential_prime)});
    }

    // Display the top error rate integers by popping them off the heap
    for _ in 1..10 {
        match int_heap.pop() {
            Some(IntWrapper { number, error, prime }) =>{
                println!("{}: {:.3}% of bases are strong liars | prime? {}", number, error * 100f32, prime);
            }
            None => println!("No more numbers were tested"),
        }
    }
}

// change n to the form 2^s·d
fn decompose(n: u64) -> (u64, u64) {
    let mut s = 0;
    let mut m = n;
    while (m & 1) == 0 {
        s += 1;
        m >>= 1;
    }
    (s, n / (1 << s))
}

// modular exponentiation implementation
fn mod_exp(mut a: u64, mut x: u64, modulus: u64) -> u64 {
    let mut remainder = 1;
    while x != 0 {
        if (x & 1) == 1 {
            remainder = a * remainder % modulus;
        }
        x >>= 1;
        a = a * a % modulus;
    }
    remainder
}

fn miller_rabin(n: u64, witness: u64) -> bool {
    if n == 2 || n == 3 {
        return true;
    } else if n <= 1 || (n & 1) == 0 {
        return false;
    }
    let (exponent, remainder) = decompose(n - 1);
    let mut x = mod_exp(witness, remainder, n);
    // If this is true, the number is either a strong liar or a witness
    if x == 1 || x == n - 1 {
        return true;
    }
    for _ in 1..exponent {
        x = mod_exp(x, 2, n);
        if x == 1 {
            // the number is composite for this base
            return false;
        } else if x == n - 1 {
            // the number is probably prime for this base
            return true;
        }
    }
    // the number is composite for this base
    return false;
}

// Language specific stuff to make the Integer Wrapper type sortable in the binary heap
#[derive(Copy, Clone, PartialEq)]
struct IntWrapper {
    number: u64,
    error: f32,
    prime: bool
}
impl Eq for IntWrapper {}
impl Ord for IntWrapper {
    fn cmp(&self, other: &IntWrapper) -> Ordering {
        // sort on error field
        self.error.partial_cmp(&other.error).unwrap()
    }
}
impl PartialOrd for IntWrapper {
    fn partial_cmp(&self, other: &IntWrapper) -> Option<Ordering> {
        Some(self.cmp(other))
    }
}