First whack at prime numbers and such.

src/cryptonum/mod.rs

@@ -8,53 +8,4 @@ mod unsigned;
 pub use self::signed::SCN;
 pub use self::unsigned::UCN;
 
-use std::ops::Neg;
-
-pub fn modinv(e: &UCN, phi: &UCN) -> UCN {
-    let (_, mut x, _) = extended_euclidean(&e, &phi);
-    let int_phi = SCN::from(phi.clone());
-    while x.is_negative() {
-        x = x + &int_phi;
-    }
-    x.into()
-}
-
-fn extended_euclidean(a: &UCN, b: &UCN) -> (SCN, SCN, SCN) {
-    let posinta = SCN::from(a.clone());
-    let posintb = SCN::from(b.clone());
-    let (d, x, y) = egcd(posinta, posintb);
-
-    if d.is_negative() {
-        (d.neg(), x.neg(), y.neg())
-    } else {
-        (d, x, y)
-    }
-}
-
-fn egcd(a: SCN, b: SCN) -> (SCN, SCN, SCN) {
-    let mut s     = SCN::zero();
-    let mut old_s = SCN::from(1 as u8);
-    let mut t     = SCN::from(1 as u8);
-    let mut old_t = SCN::zero();
-    let mut r     = b;
-    let mut old_r = a;
-
-    while !r.is_zero() {
-        let quotient = old_r.clone() / r.clone();
-
-        let prov_r = r.clone();
-        let prov_s = s.clone();
-        let prov_t = t.clone();
-
-        r = old_r - (r * &quotient);
-        s = old_s - (s * &quotient);
-        t = old_t - (t * &quotient);
-
-        old_r = prov_r;
-        old_s = prov_s;
-        old_t = prov_t;
-    }
-
-    (old_r, old_s, old_t)
-}
 

src/cryptonum/signed.rs

@@ -37,6 +37,33 @@ impl SCN {
             self.negative = false;
         }
     }
+
+    pub fn egcd(self, b: SCN) -> (SCN, SCN, SCN) {
+        let mut s     = SCN::zero();
+        let mut old_s = SCN::from(1 as u8);
+        let mut t     = SCN::from(1 as u8);
+        let mut old_t = SCN::zero();
+        let mut r     = b;
+        let mut old_r = self;
+
+        while !r.is_zero() {
+            let quotient = old_r.clone() / r.clone();
+
+            let prov_r = r.clone();
+            let prov_s = s.clone();
+            let prov_t = t.clone();
+
+            r = old_r - (r * &quotient);
+            s = old_s - (s * &quotient);
+            t = old_t - (t * &quotient);
+
+            old_r = prov_r;
+            old_s = prov_s;
+            old_t = prov_t;
+        }
+
+        (old_r, old_s, old_t)
+    }
 }
 
 impl fmt::UpperHex for SCN {
@@ -280,8 +307,6 @@ mod test {
         }
         fn division_identity(x: SCN) -> bool {
             let result = &x / &one();
-            println!("\nx: {:?}", x);
-            println!("result: {:?}", result);
             result == x
         }
 

src/cryptonum/unsigned.rs

@@ -1,4 +1,6 @@
+use cryptonum::signed::SCN;
 use num::{BigUint,ToPrimitive,Zero};
+use rand::Rng;
 use std::fmt;
 use std::fmt::Write;
 use std::cmp::Ordering;
@@ -10,15 +12,70 @@ pub struct UCN {
     pub(crate) contents: Vec<u64>
 }
 
+static SMALL_PRIMES: [u32; 310] = [
+      2,     3,     5,     7,    11,    13,    17,    19,    23,    29,
+     31,    37,    41,    43,    47,    53,    59,    61,    67,    71,
+     73,    79,    83,    89,    97,   101,   103,   107,   109,   113,
+    127,   131,   137,   139,   149,   151,   157,   163,   167,   173,
+    179,   181,   191,   193,   197,   199,   211,   223,   227,   229,
+    233,   239,   241,   251,   257,   263,   269,   271,   277,   281,
+    283,   293,   307,   311,   313,   317,   331,   337,   347,   349,
+    353,   359,   367,   373,   379,   383,   389,   397,   401,   409,
+    419,   421,   431,   433,   439,   443,   449,   457,   461,   463,
+    467,   479,   487,   491,   499,   503,   509,   521,   523,   541,
+    547,   557,   563,   569,   571,   577,   587,   593,   599,   601,
+    607,   613,   617,   619,   631,   641,   643,   647,   653,   659,
+    661,   673,   677,   683,   691,   701,   709,   719,   727,   733,
+    739,   743,   751,   757,   761,   769,   773,   787,   797,   809,
+    811,   821,   823,   827,   829,   839,   853,   857,   859,   863,
+    877,   881,   883,   887,   907,   911,   919,   929,   937,   941,
+    947,   953,   967,   971,   977,   983,   991,   997,  1009,  1013,
+   1019,  1021,  1031,  1033,  1039,  1049,  1051,  1061,  1063,  1069,
+   1087,  1091,  1093,  1097,  1103,  1109,  1117,  1123,  1129,  1151,
+   1153,  1163,  1171,  1181,  1187,  1193,  1201,  1213,  1217,  1223,
+   1229,  1231,  1237,  1249,  1259,  1277,  1279,  1283,  1289,  1291,
+   1297,  1301,  1303,  1307,  1319,  1321,  1327,  1361,  1367,  1373,
+   1381,  1399,  1409,  1423,  1427,  1429,  1433,  1439,  1447,  1451,
+   1453,  1459,  1471,  1481,  1483,  1487,  1489,  1493,  1499,  1511,
+   1523,  1531,  1543,  1549,  1553,  1559,  1567,  1571,  1579,  1583,
+   1597,  1601,  1607,  1609,  1613,  1619,  1621,  1627,  1637,  1657,
+   1663,  1667,  1669,  1693,  1697,  1699,  1709,  1721,  1723,  1733,
+   1741,  1747,  1753,  1759,  1777,  1783,  1787,  1789,  1801,  1811,
+   1823,  1831,  1847,  1861,  1867,  1871,  1873,  1877,  1879,  1889,
+   1901,  1907,  1913,  1931,  1933,  1949,  1951,  1973,  1979,  1987,
+   1993,  1997,  1999,  2003,  2011,  2017,  2027,  2029,  2039,  2053];
+
+
+
 impl UCN {
     pub fn zero() -> UCN {
         UCN{ contents: vec![] }
     }
 
+    pub fn bits(&self) -> usize {
+        self.contents.len() * 64
+    }
+
     pub fn is_zero(&self) -> bool {
         self.contents.len() == 0
     }
 
+    pub fn is_odd(&self) -> bool {
+        if self.contents.len() == 0 {
+            false
+        } else {
+            (self.contents[0] & 1) == 1
+        }
+    }
+
+    pub fn is_even(&self) -> bool {
+        if self.contents.len() == 0 {
+            false
+        } else {
+            (self.contents[0] & 1) == 0
+        }
+    }
+
     fn clean(&mut self) {
         loop {
             match self.contents.pop() {
@@ -65,6 +122,145 @@ impl UCN {
         res.clean();
         res
     }
+
+    pub fn modinv(&self, phi: &UCN) -> UCN {
+        let (_, mut x, _) = self.extended_euclidean(&phi);
+        let int_phi = SCN::from(phi.clone());
+        while x.is_negative() {
+            x = x + &int_phi;
+        }
+        x.into()
+    }
+
+    fn extended_euclidean(&self, b: &UCN) -> (SCN, SCN, SCN) {
+        let posinta = SCN::from(self.clone());
+        let posintb = SCN::from(b.clone());
+        let (d, x, y) = posinta.egcd(posintb);
+
+        if d.is_negative() {
+            (d.neg(), x.neg(), y.neg())
+        } else {
+            (d, x, y)
+        }
+    }
+
+    pub fn generate_prime<F,G>(rng: &mut G,
+                               bitlen: usize,
+                               iterations: usize,
+                               ok_value: F)
+        -> UCN
+      where
+        G: Rng,
+        F: Fn(&UCN) -> bool
+    {
+        let one = UCN::from(1 as u8);
+        let topbit = &one << (bitlen - 1);
+
+        loop {
+            let base = random_number(rng, bitlen);
+            let candidate = base | &topbit | &one;
+
+            if !ok_value(&candidate) {
+                continue;
+            }
+
+            if candidate.probably_prime(rng, iterations) {
+                return candidate
+            }
+        }
+    }
+
+    pub fn probably_prime<G: Rng>(&self, g: &mut G, iters: usize) -> bool {
+        for tester in SMALL_PRIMES.iter() {
+            if (self % UCN::from(*tester)).is_zero() {
+                return false;
+            }
+        }
+        miller_rabin(g, &self, iters)
+    }
+
+    pub fn modexp(&self, e: &UCN, m: &UCN) -> UCN {
+        let mut b = self.clone() % m;
+        let mut eprime = e.clone();
+        let mut result = UCN::from(1 as u8);
+
+        loop {
+            if eprime.is_zero() {
+                return result;
+            }
+
+            if eprime.is_odd() {
+                result = (result * &b) % m;
+            }
+
+            b = (&b * &b) % m;
+            eprime >>= 1;
+        }
+    }
+}
+
+fn miller_rabin<G: Rng>(g: &mut G, n: &UCN, iters: usize) -> bool {
+    let one = UCN::from(1 as u8);
+    let two = UCN::from(2 as u8);
+    let nm1 = n - &one;
+    // Quoth Wikipedia:
+    // write n - 1 as 2^r*d with d odd by factoring powers of 2 from n - 1
+    let mut d = nm1.clone();
+    let mut r = 0;
+    while d.is_even() {
+        d >>= 1;
+        r += 1;
+        assert!(r < n.bits());
+    }
+    // WitnessLoop: repeat k times
+    'WitnessLoop: for _k in 0..iters {
+        // pick a random integer a in the range [2, n - 2]
+        let a = random_in_range(g, &two, &nm1);
+        // x <- a^d mod n
+        let mut x = a.modexp(&d, &n);
+        // if x = 1 or x = n - 1 then
+        if (&x == &one) || (&x == &nm1) {
+            // continue WitnessLoop
+            continue 'WitnessLoop;
+        }
+        // repeat r - 1 times:
+        for _i in 0..r {
+            // x <- x^2 mod n
+            x = x.modexp(&two, &n);
+            // if x = 1 then
+            if &x == &one {
+                // return composite
+                return false;
+            }
+            // if x = n - 1 then
+            if &x == &nm1 {
+                // continue WitnessLoop
+                continue 'WitnessLoop;
+            }
+        }
+        // return composite
+        return false;
+    }
+    // return probably prime
+    true
+}
+
+fn random_in_range<G: Rng>(rng: &mut G, min: &UCN, max: &UCN) -> UCN {
+    let bitlen = ((max.bits() + 31) / 32) * 32;
+    loop {
+        let candidate = random_number(rng, bitlen);
+
+        if (&candidate >= min) && (&candidate < max) {
+            return candidate;
+        }
+    }
+}
+
+fn random_number<G: Rng>(rng: &mut G, bitlen: usize) -> UCN {
+    assert!(bitlen % 32 == 0);
+    let wordlen = bitlen / 32;
+    let components = rng.gen_iter().take(wordlen).collect();
+    UCN{ contents: components }
 }
 
 impl fmt::UpperHex for UCN {
@@ -1026,6 +1222,30 @@ mod test {
         gold_test("tests/mod_tests.txt", |x,y| x % y);
     }
 
+    #[test]
+    fn modexp_tests() {
+        let mut file = File::open("tests/modpow_tests.txt").unwrap();
+        let mut contents = String::new();
+        file.read_to_string(&mut contents).unwrap();
+        let mut iter = contents.lines();
+
+        while let Some(xstr) = iter.next() {
+            let ystr = iter.next().unwrap();
+            let zstr = iter.next().unwrap();
+            let rstr = iter.next().unwrap();
+
+            assert!(xstr.starts_with("x: "));
+            assert!(ystr.starts_with("y: "));
+            assert!(zstr.starts_with("z: "));
+            assert!(rstr.starts_with("r: "));
+            let x = UCN::from_str(&xstr[3..]);
+            let y = UCN::from_str(&ystr[3..]);
+            let z = UCN::from_str(&zstr[3..]);
+            let r = UCN::from_str(&rstr[3..]);
+            assert_eq!(x.modexp(&y, &z), r);
+        }
+    }
+
     quickcheck! {
         fn and_over_or_distribution(a: UCN, b: UCN, c: UCN) -> bool {
             (&a & (&b | &c)) == ((&a & &b) | (&a & &c))