Skip to content

Instantly share code, notes, and snippets.

@benbrittain
Created November 23, 2015 15:41
Show Gist options
  • Save benbrittain/efe514e5752d4a14f4a4 to your computer and use it in GitHub Desktop.
Save benbrittain/efe514e5752d4a14f4a4 to your computer and use it in GitHub Desktop.
#![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))
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment