diff --git a/shamirsecretsharing-rs/Cargo.toml b/shamirsecretsharing-rs/Cargo.toml
index 410275d..248218c 100644
--- a/shamirsecretsharing-rs/Cargo.toml
+++ b/shamirsecretsharing-rs/Cargo.toml
@@ -9,3 +9,12 @@ rand = "0.6.5"
 num = "0.2.0"
 num-bigint = {version = "0.2.2", features = ["rand"]}
 num-traits = "0.2.8"
+
+[dev-dependencies]
+criterion = "0.2"
+
+# Criterion benchmarks
+[[bench]]
+path = "./benchmarks/benches.rs"
+name = "benches"
+harness = false
\ No newline at end of file
diff --git a/shamirsecretsharing-rs/benchmarks/benches.rs b/shamirsecretsharing-rs/benchmarks/benches.rs
new file mode 100644
index 0000000..dd68381
--- /dev/null
+++ b/shamirsecretsharing-rs/benchmarks/benches.rs
@@ -0,0 +1,38 @@
+#[macro_use]
+extern crate criterion;
+extern crate shamirsecretsharing_rs;
+extern crate num_bigint;
+
+use criterion::{Criterion, Benchmark};
+use shamirsecretsharing_rs::*;
+use num_bigint::BigInt;
+
+use std::str::FromStr;
+
+
+mod mod_inv_benches {
+    use super::*;
+
+    pub fn bench_modular_inv(c: &mut Criterion) {
+
+        let modul1 = BigInt::from_str("7237005577332262213973186563042994240857116359379907606001950938285454250989").unwrap();
+        let d1 = BigInt::from_str("182687704666362864775460604089535377456991567872").unwrap();
+
+        let modul2 = BigInt::from_str("7237005577332262213973186563042994240857116359379907606001950938285454250989").unwrap();
+        let d2 = BigInt::from_str("182687704666362864775460604089535377456991567872").unwrap();
+
+        c.bench(
+            "Modular Inverse",
+            Benchmark::new("Kalinski Modular inverse", move |b| b.iter(|| kalinski_inv(&d1, &modul1)))
+        );
+
+        c.bench(
+            "Modular Inverse",
+            Benchmark::new("Standard Mod Inv", move |b| b.iter(|| mod_inverse(d2.clone(), modul2.clone())))
+        );
+    }
+}
+
+criterion_group!(benches,
+                mod_inv_benches::bench_modular_inv);
+criterion_main!(benches);
\ No newline at end of file
diff --git a/shamirsecretsharing-rs/src/lib.rs b/shamirsecretsharing-rs/src/lib.rs
index a8a56c7..f361ac4 100644
--- a/shamirsecretsharing-rs/src/lib.rs
+++ b/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};
@@ -72,7 +75,7 @@ fn unpack_shares(s: Vec<[BigInt;2]>) -> (Vec<BigInt>, Vec<BigInt>) {
     (shares, is)
 }
 
-fn mod_inverse(a: BigInt, module: BigInt) -> BigInt {
+pub fn mod_inverse(a: BigInt, module: BigInt) -> BigInt {
     // TODO search biguint impl of mod_inv
     let mut mn = (module.clone(), a);
     let mut xy: (BigInt, BigInt) = (Zero::zero(), One::one());
@@ -89,6 +92,94 @@ 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 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 << 1;
+                },
+                // u isn't even but v is even
+                (false, true, _, _) => {
+
+                    v = v >> 1;
+                    r = &r << 1;
+                },
+                // u and v aren't even and u > v
+                (false, false, true, _) => {
+
+                    u = &u - &v;
+                    u = u >> 1;
+                    r = &r + &s;
+                    s = &s << 1;
+                },
+                // u and v aren't even and v > u
+                (false, false, false, true) => {
+
+                    v = &v - &u;
+                    v = v >> 1;
+                    s = &r + &s;
+                    r = &r << 1;
+                },
+                (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 {
+            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();
@@ -118,7 +209,7 @@ pub fn lagrange_interpolation(p: &BigInt, shares_packed: Vec<[BigInt;2]>) -> Big
     }
     let modinv_mul: BigInt;
     if res_d != Zero::zero() {
-        let modinv = mod_inverse(res_d, p.clone());
+        let modinv = kalinski_inv(&res_d, &p);
         modinv_mul = res_n * modinv;
     } else {
         modinv_mul = res_n;
@@ -130,7 +221,9 @@ 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,29 @@ 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);
+
+        // Modulo: 2^252 + 27742317777372353535851937790883648493
+        // Tested: 182687704666362864775460604089535377456991567872
+        // Expected for: inverse_mod(a, l) computed on SageMath:
+        // `7155219595916845557842258654134856828180378438239419449390401977965479867845`.
+        let modul3 = BigInt::from_str("7237005577332262213973186563042994240857116359379907606001950938285454250989").unwrap();
+        let d = BigInt::from_str("182687704666362864775460604089535377456991567872").unwrap();
+        let res4 = kalinski_inv(&d, &modul3); 
+        let expected4 = BigInt::from_str("7155219595916845557842258654134856828180378438239419449390401977965479867845").unwrap();
+        assert_eq!(expected4, res4);
+    }
 }