[WIP] Implement kalinski modular inverse and tests

Implemented Kalinski Modular inverse algorithm.
See: https://www.researchgate.net/publication/3387259_Improved_Montgomery_modular_inverse_algorithm

Also implemented test for big and small Bignums.

One test is giving results over the Montgomery Domain, so
a little bit of research has to be done.

shamirsecretsharing-rs/src/lib.rs

@@ -3,8 +3,11 @@ extern crate num;
 extern crate num_bigint;
 extern crate num_traits;
 
+use std::str::FromStr;
+
 use num_bigint::RandBigInt;
 use num::pow::pow;
+use num::Integer;
 
 
 use num_bigint::{BigInt, ToBigInt};
@@ -89,6 +92,95 @@ fn mod_inverse(a: BigInt, module: BigInt) -> BigInt {
     xy.0
 }
 
+/// Compute `a^-1 (mod l)` using the the Kalinski implementation
+/// of the Montgomery Modular Inverse algorithm.
+/// B. S. Kaliski Jr. - The  Montgomery  inverse  and  its  applica-tions.
+/// IEEE Transactions on Computers, 44(8):1064–1065, August-1995
+pub fn kalinski_inv(a: &BigInt, modulo: &BigInt) -> BigInt {
+    // This Phase I indeed is the Binary GCD algorithm , a version o Stein's algorithm
+    // which tries to remove the expensive division operation away from the Classical
+    // Euclidean GDC algorithm replacing it for Bit-shifting, subtraction and comparaison.
+    // 
+    // Output = `a^(-1) * 2^k (mod l)` where `k = log2(modulo) == Number of bits`.
+    // 
+    // Stein, J.: Computational problems associated with Racah algebra.J. Comput. Phys.1, 397–405 (1967).
+    let phase1 = |a: &BigInt| -> (BigInt, u64) {
+        assert!(a != &BigInt::zero());
+        let p = modulo;
+        let mut u = modulo.clone();
+        let mut v = a.clone();
+        let mut r = BigInt::zero();
+        let mut s = BigInt::one();
+        let two = BigInt::from(2u64); 
+        let mut k = 0u64;
+
+        while v > BigInt::zero() {
+            match(u.is_even(), v.is_even(), u > v, v >= u) {
+                // u is even
+                (true, _, _, _) => {
+
+                    u = u >> 1;
+                    s = &s * &two;
+                },
+                // u isn't even but v is even
+                (false, true, _, _) => {
+
+                    v = v >> 1;
+                    r = &r * &two;
+                },
+                // u and v aren't even and u > v
+                (false, false, true, _) => {
+
+                    u = &u - &v;
+                    u = u >> 1;
+                    r = &r + &s;
+                    s = &s * &two;
+                },
+                // u and v aren't even and v > u
+                (false, false, false, true) => {
+
+                    v = &v - &u;
+                    v = v >> 1;
+                    s = &r + &s;
+                    r = &r * &two;
+                },
+                (false, false, false, false) => panic!("Unexpected error has ocurred."),
+            }
+            k += 1;
+        }
+        if &r > p {
+            r = &r - p;
+        }
+        ((p - &r), k)
+    };
+
+    // Phase II performs some adjustments to obtain
+    // the Montgomery inverse.
+    // 
+    // We implement it as a clousure to be able to grap the 
+    // kalinski_inv scope to get `modulo` variable.
+    let phase2 = |r: &BigInt, k: &u64| -> BigInt {
+        let mut rr = r.clone();
+        let _p = modulo;
+
+        for _i in 0..(k - modulo.bits() as u64) {
+            match rr.is_even() {
+                true => {
+                    rr = rr >> 1;
+                },
+                false => {
+                    rr = (rr + modulo) >> 1;
+                }
+            }
+        }
+        rr 
+    };
+
+    let (r, z) = phase1(&a.clone());
+
+    phase2(&r, &z)
+}
+
 pub fn lagrange_interpolation(p: &BigInt, shares_packed: Vec<[BigInt;2]>) -> BigInt {
     let mut res_n: BigInt = Zero::zero();
     let mut res_d: BigInt = Zero::zero();
@@ -131,6 +223,7 @@ pub fn lagrange_interpolation(p: &BigInt, shares_packed: Vec<[BigInt;2]>) -> Big
 #[cfg(test)]
 mod tests {
     use super::*;
+    use std::str::FromStr;
 
     #[test]
     fn test_create_and_lagrange_interpolation() {
@@ -151,4 +244,37 @@ mod tests {
         println!("original secret: {:?}", k.to_string());
         assert_eq!(k, r);
     }
+
+    #[test]
+    fn kalinski_modular_inverse() {
+        let modul1 = BigInt::from(127u64);
+
+        let a = BigInt::from(79u64);
+        let res1 = kalinski_inv(&a, &modul1);
+        let expected1 = BigInt::from(82u64);
+        assert_eq!(res1, expected1);
+
+        let b = BigInt::from(50u64);
+        let res2 = kalinski_inv(&b, &modul1);
+        let expected2 = BigInt::from(94u64);
+        assert_eq!(res2, expected2);
+
+        // Big numbers testing.
+        // C = 19.
+        // modul2 = 2^255 - 19.
+        let modul2 = BigInt::from_str("57896044618658097711785492504343953926634992332820282019728792003956564819949").unwrap();
+        let c = BigInt::from_str("19").unwrap();
+        let res3 = kalinski_inv(&c, &modul2);
+        let expected3 = BigInt::from_str("1").unwrap();
+        assert_eq!(res3, expected3);
+
+        /*// D = 182687704666362864775460604089535377456991567872.
+        // modul3 = 2^252 + 27742317777372353535851937790883648493.
+        let modul3 = BigInt::from_str("7237005577332262213973186563042994240857116359379907606001950938285454250989").unwrap();
+        let d = BigInt::from_str("182687704666362864775460604089535377456991567872").unwrap();
+        let res4 = kalinski_inv(&d, &modul3); 
+        println!("RES ON IMPL: {}", res4);
+        let expected4 = BigInt::from_str("7155219595916845557842258654134856828180378438239419449390401977965479867845").unwrap();
+        assert_eq!(expected4, res4);*/
+    }
 }