diff --git a/src/lib.rs b/src/lib.rs index 38f5b39..e1bc79f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,5 +1,14 @@ use reikna::totient::totient; +fn gcd(mut a: i64, mut b: i64) -> i64 { + while b != 0 { + let temp = b; + b = a % b; + a = temp; + } + a.abs() +} + // Modular arithmetic functions using i64 fn mod_add(a: i64, b: i64, p: i64) -> i64 { (a + b) % p @@ -24,6 +33,7 @@ pub fn mod_exp(mut base: i64, mut exp: i64, p: i64) -> i64 { //compute the modular inverse of a modulo p using Fermat's little theorem, p not necessarily prime fn mod_inv(a: i64, p: i64) -> i64 { + assert!(gcd(a, p) == 1, "{} and {} are not coprime", a, p); mod_exp(a, totient(p as u64) as i64 - 1, p) // order of mult. group is Euler's totient function } diff --git a/src/test.rs b/src/test.rs index c0b0d4a..3553a79 100644 --- a/src/test.rs +++ b/src/test.rs @@ -50,7 +50,7 @@ mod tests { #[test] fn test_polymul_ntt_prime_power_modulus() { - let modulus: i64 = (17 as i64).pow(4); // modulus p^k or 2*p^k + let modulus: i64 = (17 as i64).pow(4); // modulus p^k let root: i64 = 3; // Primitive root of unity let n: usize = 8; // Length of the NTT (must be a power of 2) let omega = omega(root, modulus, n); // n-th root of unity @@ -70,4 +70,5 @@ mod tests { // Ensure both methods produce the same result assert_eq!(c_std, c_fast, "The results of polymul and polymul_ntt do not match"); } + }