Use crate num-bigint in miller rabin

Cargo.toml

@@ -4,4 +4,7 @@ version = "0.1.0"
 edition = "2021"
 
 [dependencies]
+num-bigint = { version = "0.4.4", features = ["rand"] }
+num-integer = "0.1.46"
+num-traits = "0.2.18"
 rand = "0.8.5"

src/primality_test/miller_rabin.rs

@@ -3,6 +3,7 @@ mod utils;
 
 use self::composite_evidence::CompositeEvidence;
 use crate::traits::PrimalityTest;
+use num_bigint::BigUint;
 
 pub struct MillerRabin;
 
@@ -14,14 +15,14 @@ impl PrimalityTest for MillerRabin {
         if p < 2 || p % 2 == 0 {
             return false;
         }
-        miller_rabin(p, 10)
+        miller_rabin(&BigUint::from(p), 10)
     }
 }
 
-fn miller_rabin(p: u128, trials: usize) -> bool {
+fn miller_rabin(p: &BigUint, trials: usize) -> bool {
     let evidence = CompositeEvidence::new(p);
-    let likely_prime = |witness| !evidence.witnessed_by(witness);
-    utils::RandomIntegers::new(2..p - 1)
+    let likely_prime = |witness| !evidence.witnessed_by(&witness);
+    utils::RandomIntegers::new(BigUint::from(2u8)..p - 1u8)
         .take(trials)
         .all(likely_prime)
 }

src/primality_test/miller_rabin/composite_evidence.rs

@@ -1,57 +1,64 @@
 use super::utils;
+use num_bigint::BigUint;
+use num_traits::One;
 
-pub struct CompositeEvidence {
-    n: u128,
+pub struct CompositeEvidence<'a> {
+    n: &'a BigUint,
     n_minus_1: Decomposed,
 }
 
-impl CompositeEvidence {
-    pub fn new(n: u128) -> Self {
-        let n_minus_1 = Decomposed::new(n - 1);
+impl<'a> CompositeEvidence<'a> {
+    pub fn new(n: &'a BigUint) -> Self {
+        let n_minus_1 = Decomposed::new(n - 1u8);
         Self { n, n_minus_1 }
     }
 
-    pub fn witnessed_by(&self, witness: u128) -> bool {
+    pub fn witnessed_by(&self, witness: &BigUint) -> bool {
         match self.raise_to_n_minus_1_mod_n(witness) {
             Ok(result) => fails_fermats_condition(result),
             Err(FoundNonTrivialSqrtOf1) => true,
         }
     }
 
-    fn raise_to_n_minus_1_mod_n(&self, base: u128) -> ExponentiationResult {
-        let odd_factor_in_exp = self.n_minus_1.odd_factor;
-        let mut result = utils::modular_exponentiation(base, odd_factor_in_exp, self.n);
+    fn raise_to_n_minus_1_mod_n(&self, base: &BigUint) -> ExponentiationResult {
+        let odd_factor_in_exp = &self.n_minus_1.odd_factor;
+        let mut result = base.modpow(odd_factor_in_exp, &self.n);
         for _ in 0..self.n_minus_1.exponent_of_2 {
-            if utils::is_nontrivial_sqrt_of_1(result, self.n) {
+            if self.is_nontrivial_sqrt_of_1(&result) {
                 return Err(FoundNonTrivialSqrtOf1);
             }
-            result = utils::modular_exponentiation(result, 2, self.n);
+            result = result.modpow(&BigUint::from(2u8), &self.n);
         }
         Ok(RaisedToNMinus1ModN(result))
     }
+
+    pub fn is_nontrivial_sqrt_of_1(&self, solution: &BigUint) -> bool {
+        let squared = solution.modpow(&BigUint::from(2u8), &self.n);
+        squared == BigUint::one() && solution != &BigUint::one() && solution != &(self.n - 1u8)
+    }
 }
 
 fn fails_fermats_condition(r: RaisedToNMinus1ModN) -> bool {
-    r.0 != 1
+    !r.0.is_one()
 }
 
 type ExponentiationResult = Result<RaisedToNMinus1ModN, FoundNonTrivialSqrtOf1>;
 
-struct RaisedToNMinus1ModN(u128);
+struct RaisedToNMinus1ModN(BigUint);
 
 struct FoundNonTrivialSqrtOf1;
 
 struct Decomposed {
     exponent_of_2: u32,
-    odd_factor: u128,
+    odd_factor: BigUint,
 }
 
 impl Decomposed {
     /// Decomposes `number` into `exponent_of_2` and `odd_factor`,
     /// where `number = 2^exponent_of_2 * odd_factor`.
-    pub fn new(number: u128) -> Self {
-        let exponent_of_2 = utils::highest_power_of_2_divisor(number);
-        let odd_factor = number / 2u128.pow(exponent_of_2);
+    pub fn new(number: BigUint) -> Self {
+        let exponent_of_2 = utils::highest_power_of_2_divisor(&number);
+        let odd_factor = number / BigUint::from(2u8).pow(exponent_of_2);
         Self {
             exponent_of_2,
             odd_factor,

src/primality_test/miller_rabin/utils.rs

@@ -1,52 +1,35 @@
-use rand::Rng;
+use num_bigint::{BigUint, RandBigInt};
+use num_integer::Integer;
 use std::ops::Range;
 
 pub struct RandomIntegers {
-    range: Range<u128>,
+    lo: BigUint,
+    hi: BigUint,
 }
 
 impl RandomIntegers {
-    pub fn new(range: Range<u128>) -> Self {
-        Self { range }
+    pub fn new(range: Range<BigUint>) -> Self {
+        Self {
+            lo: range.start,
+            hi: range.end,
+        }
     }
 }
 
 impl Iterator for RandomIntegers {
-    type Item = u128;
+    type Item = BigUint;
 
     fn next(&mut self) -> Option<Self::Item> {
-        Some(rand::thread_rng().gen_range(self.range.clone()))
-    }
-}
-
-pub fn modular_exponentiation(base: u128, exp: u128, modulus: u128) -> u128 {
-    if modulus == 1 {
-        return 0;
+        Some(rand::thread_rng().gen_biguint_range(&self.lo, &self.hi))
     }
-    let mut base = base % modulus;
-    let mut exp = exp;
-    let mut result = 1;
-    while exp > 0 {
-        if exp % 2 == 1 {
-            result = (result * base) % modulus; // If the exponent is odd, multiply the result by the base
-        }
-        exp >>= 1; // (divide by 2, dropping any remainder)
-        base = (base * base) % modulus;
-    }
-    result
 }
 
-pub fn highest_power_of_2_divisor(base: u128) -> u32 {
+pub fn highest_power_of_2_divisor(base: &BigUint) -> u32 {
     let mut exp = 0;
-    let mut base = base;
-    while base % 2 == 0 {
+    let mut base = base.clone();
+    while base.is_even() {
         exp += 1;
-        base /= 2;
+        base /= 2u8;
     }
     exp
 }
-
-pub fn is_nontrivial_sqrt_of_1(solution: u128, number: u128) -> bool {
-    let squared = (solution * solution) % number;
-    squared == 1 && solution != 1 && solution != number - 1
-}