1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123

use std::u64;

extern crate rand;
use rand::distributions::{Distribution, Uniform};
 
extern crate modular_arithmetic;
use modular_arithmetic::{mut_even_power, split_odd_even, jacobi_symbol};

mod pseudoprime;
use pseudoprime::{miller_rabin_witness, lucas_test};



//TODO: ensure # rounds < n; infinite loop otherwise
//true => composite | false => maybe prime
pub fn miller_rabin(n: u64, rounds: u8) -> bool {
	//manually check divisibility properties for small primes
	//improves efficiency
	match trial_division_test(&n) {
		Some(b) => return b,
		None => ()
	}
	
	let mut tested_numbers = Vec::with_capacity(rounds as usize);
	//uniform sample space
	let mut rng = rand::thread_rng();
	let sample_space = Uniform::from(2..n-1);

	//split n-1 into odd component and even power in 2 steps
	let mut odd_component = n-1;
	let power2 = mut_even_power(&mut odd_component); //mutate n-1 into odd component, return even power

	//loop until #rounds different numbers have been tried
	while tested_numbers.len() < rounds as usize {

		let test_num = sample_space.sample(&mut rng);
		//excludes numbers already tried
		if !tested_numbers.contains(&test_num){
			//run a single miller rabin test 
			let test_result = miller_rabin_witness(n, odd_component, power2, test_num);
			if test_result {
				return true
			} else {
				tested_numbers.push(test_num);
			}
		}
	}
	return false
}





//returns true if n is definitely composite and false if n is likely to be a prime
pub fn baillie_psw(n: u64) -> bool {
	//perform trial division
	match trial_division_test(&n) {
		Some(b) => return b,
		None => ()
	}

	//check if n is a square manually
	let possible_square = (n as f64).sqrt() as u64;
	if possible_square * possible_square == n {
		return true
	}

	if (possible_square + 1) * (possible_square + 1) == n {
		return true
	}

	//base 2 strong pseudoprime test
	let (odd, even_pow) = split_odd_even(n-1);

	if miller_rabin_witness(n, odd, even_pow, 2) {
		return true
	}

	//find D with jacobi symbol (D/n) = -1
	//guaranteed to terminate
	let mut d:i64 = 5;
	let mut d_is_neg = false;

	loop {
		if jacobi_symbol(d, n) == -1 {
			break()
		}

		if d_is_neg == true {
			d *= -1;
			d += 2;
			d_is_neg = !d_is_neg;
		} else {
			d += 2;
			d += -1;
			d_is_neg = !d_is_neg;
		}
	}

	let p = 1;
	let q = (1-d)/4;



	//lucas probable prime test with D, P Q and n
	return lucas_test(n, odd, even_pow as u32, p, q)
}


//checks if n is divisible by one of the primes in the collection
fn trial_division_test(n: &u64) -> Option<bool> {
	let primes = vec![2,3,5,7,11];
	let primes_iter = primes.iter();
	for prime in primes_iter {
		if n == prime {
			return Some(false)
		} else if n % prime == 0 {
			return Some(true)
		}
	}
	return None
}