diff --git a/examples/main.rs b/examples/main.rs new file mode 100644 index 0000000..634257e --- /dev/null +++ b/examples/main.rs @@ -0,0 +1,22 @@ +extern crate math; +use math::modulus::prime::Prime; +use math::dft::ntt::Table; + +fn main() { + // Example usage of `Prime` + let q_base: u64 = 65537; // Example prime base + let q_power: u64 = 2; // Example power + let mut prime_instance: Prime = Prime::::new(q_base, q_power); + + // Display the fields of `Prime` to verify + println!("Prime instance created:"); + println!("q: {}", prime_instance.q()); + println!("q_base: {}", prime_instance.q_base()); + println!("q_power: {}", prime_instance.q_power()); + + let n: u64 = 1024; + let nth_root: u64 = n<<1; + + let ntt_table: Table<'_, u64> = Table::::new(&mut prime_instance, n, nth_root); + +} \ No newline at end of file diff --git a/src/dft.rs b/src/dft.rs index 27eaa65..221c7d6 100644 --- a/src/dft.rs +++ b/src/dft.rs @@ -1 +1 @@ -pub(crate) mod primitive_root; \ No newline at end of file +pub mod ntt; \ No newline at end of file diff --git a/src/dft/ntt.rs b/src/dft/ntt.rs new file mode 100644 index 0000000..ae35c27 --- /dev/null +++ b/src/dft/ntt.rs @@ -0,0 +1,49 @@ +use crate::modulus::montgomery::Montgomery; +use crate::modulus::prime::Prime; + +pub struct Table<'a, O>{ + prime:&'a Prime, + pub psi_forward_rev:Vec>, + psi_backward_rev: Vec>, + n_inv: Montgomery, +} + +impl<'a> Table<'a, u64> { + pub fn new(prime: &'a mut Prime, n: u64, nth_root: u64)->Self{ + + assert!(n&(n-1) == 0, "invalid argument: n = {} is not a power of two", n); + assert!(n&(n-1) == 0, "invalid argument: nth_root = {} is not a power of two", nth_root); + assert!(n < nth_root, "invalid argument: n = {} cannot be greater or equal to nth_root = {}", n, nth_root); + + let psi: u64 = prime.primitive_nth_root(nth_root); + + let psi_mont: Montgomery = prime.montgomery.prepare(psi); + let psi_inv_mont: Montgomery = prime.montgomery.pow(psi_mont, prime.phi-1); + + let mut psi_forward_rev: Vec> = vec![Montgomery(0); (nth_root >> 1) as usize]; + let mut psi_backward_rev: Vec> = vec![Montgomery(0); (nth_root >> 1) as usize]; + + psi_forward_rev[0] = prime.montgomery.one(); + psi_backward_rev[0] = prime.montgomery.one(); + + let log_nth_root_half: usize = (usize::MAX - ((nth_root>>1 as usize)-1).leading_zeros() as usize) as usize; + + for i in 1..(nth_root>>1) as usize{ + + let i_rev_prev: usize = (i-1).reverse_bits() >> (usize::MAX - log_nth_root_half) as usize; + let i_rev_next: usize = i.reverse_bits() >> (usize::MAX - log_nth_root_half) as usize; + + psi_forward_rev[i_rev_next] = prime.montgomery.mul_internal(psi_forward_rev[i_rev_prev], psi_mont); + psi_backward_rev[i_rev_next] = prime.montgomery.mul_internal(psi_backward_rev[i_rev_prev], psi_inv_mont); + } + + let n_inv: Montgomery = prime.montgomery.pow(prime.montgomery.prepare(nth_root>>1), prime.phi-1); + + Self{ + prime: prime, + psi_forward_rev: psi_forward_rev, + psi_backward_rev: psi_backward_rev, + n_inv: n_inv, + } + } +} \ No newline at end of file diff --git a/src/dft/primitive_root.rs b/src/dft/primitive_root.rs deleted file mode 100644 index e69de29..0000000 diff --git a/src/lib.rs b/src/lib.rs index afb5ba9..ef496bd 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,4 +1,5 @@ #![feature(bigint_helper_methods)] #![feature(test)] -pub mod modulus; \ No newline at end of file +pub mod modulus; +pub mod dft; \ No newline at end of file diff --git a/src/modulus.rs b/src/modulus.rs index 581f1d4..7592558 100644 --- a/src/modulus.rs +++ b/src/modulus.rs @@ -1,6 +1,6 @@ -pub(crate) mod prime; -pub(crate) mod montgomery; -pub(crate) mod barrett; +pub mod prime; +pub mod barrett; +pub mod montgomery; trait ReduceOnce{ fn reduce_once_assign(&mut self, q: O); @@ -8,12 +8,13 @@ trait ReduceOnce{ } impl ReduceOnce for u64{ + #[inline(always)] fn reduce_once_assign(&mut self, q: u64){ if *self >= q{ *self -= q } } - + #[inline(always)] fn reduce_once(&self, q:u64) -> u64{ if *self >= q { *self - q @@ -21,4 +22,4 @@ impl ReduceOnce for u64{ *self } } -} \ No newline at end of file +} diff --git a/src/modulus/barrett.rs b/src/modulus/barrett.rs index 9a2b541..2a625f3 100644 --- a/src/modulus/barrett.rs +++ b/src/modulus/barrett.rs @@ -1,19 +1,25 @@ +use crate::modulus::ReduceOnce; + use num_bigint::BigUint; use num_traits::cast::ToPrimitive; #[derive(Clone, Copy, Debug, PartialEq, Eq)] -pub struct BarrettPrecomp(O, O); +pub struct BarrettPrecomp{ + q: O, + lo:O, + hi:O, +} impl BarrettPrecomp{ #[inline(always)] pub fn value_hi(&self) -> &O{ - &self.1 + &self.hi } #[inline(always)] pub fn value_lo(&self) -> &O{ - &self.0 + &self.lo } } @@ -23,6 +29,21 @@ impl BarrettPrecomp{ big_r = big_r / BigUint::from(q); let lo = (&big_r & BigUint::from(u64::MAX)).to_u64().unwrap(); let hi = (big_r >> 64u64).to_u64().unwrap(); - Self(lo, hi) + Self{q, lo, hi} + } + + /// Returns lhs mod q. + #[inline(always)] + pub fn reduce(&self, lhs: u64) -> u64{ + let mut r: u64 = self.reduce_lazy(lhs); + r.reduce_once_assign(self.q); + r + } + + /// Returns lhs mod q in range [0, 2q-1]. + #[inline(always)] + pub fn reduce_lazy(&self, lhs: u64) -> u64{ + let (_, mhi) = lhs.widening_mul(self.hi); + lhs - mhi.wrapping_mul(self.q) } } \ No newline at end of file diff --git a/src/modulus/prime_generation.rs b/src/modulus/generation.rs similarity index 95% rename from src/modulus/prime_generation.rs rename to src/modulus/generation.rs index a4d43af..246dfe9 100644 --- a/src/modulus/prime_generation.rs +++ b/src/modulus/generation.rs @@ -3,16 +3,8 @@ use crate::modulus::prime; use prime::Prime; use primality_test::is_prime; -pub struct NTTFriendlyPrimesGenerator{ - size: f64, - next_prime: u64, - prev_prime: u64, - nth_root: u64, - check_next_prime: bool, - check_prev_prime: bool, -} - -impl NTTFriendlyPrimesGenerator { +impl NTTFriendlyPrimesGenerator { + pub fn new(bit_size: u64, nth_root: u64) -> Self{ let mut check_next_prime: bool = true; let mut check_prev_prime: bool = true; diff --git a/src/modulus/montgomery.rs b/src/modulus/montgomery.rs index c42d477..f9366c4 100644 --- a/src/modulus/montgomery.rs +++ b/src/modulus/montgomery.rs @@ -6,7 +6,7 @@ extern crate test; /// Montgomery is a generic struct storing /// an element in the Montgomery domain. #[derive(Clone, Copy, Debug, PartialEq, Eq)] -pub struct Montgomery(O); +pub struct Montgomery(pub O); /// Implements helper methods on the struct Montgomery. impl Montgomery{ @@ -31,8 +31,10 @@ impl Montgomery{ #[derive(Clone, Copy, Debug, PartialEq, Eq)] pub struct MontgomeryPrecomp{ q: O, - q_barrett: BarrettPrecomp, + barrett: BarrettPrecomp, q_inv: O, + one: Montgomery, + minus_one: Montgomery, } /// MontgomeryPrecomp is a set of methods implemented for MontgomeryPrecomp @@ -42,7 +44,7 @@ impl MontgomeryPrecomp{ /// Returns an new instance of MontgomeryPrecomp. /// This method will fail if gcd(q, 2^64) != 1. #[inline(always)] - fn new(q: u64) -> MontgomeryPrecomp{ + pub fn new(q: u64) -> MontgomeryPrecomp{ assert!(q & 1 != 0, "Invalid argument: gcd(q={}, radix=2^64) != 1", q); let mut q_inv: u64 = 1; let mut q_pow = q; @@ -50,42 +52,96 @@ impl MontgomeryPrecomp{ q_inv = q_inv.wrapping_mul(q_pow); q_pow = q_pow.wrapping_mul(q_pow); } - Self{ q: q, q_barrett: BarrettPrecomp::new(q), q_inv: q_inv} + let mut precomp = Self{ + q: q, + barrett: BarrettPrecomp::new(q), + q_inv: q_inv, + one: Montgomery(0), + minus_one: Montgomery(0), + }; + + precomp.one = precomp.prepare(1); + precomp.minus_one = Montgomery(q-precomp.one.value()); + + precomp } - /// Returns (lhs<<64)%q as a Montgomery. + /// Returns 2^64 mod q as a Montgomery. #[inline(always)] - fn prepare(&self, lhs: u64) -> Montgomery{ + pub fn one(&self) -> Montgomery{ + self.one + } + + /// Returns (q-1) * 2^64 mod q as a Montgomery. + #[inline(always)] + pub fn minus_one(&self) -> Montgomery{ + self.minus_one + } + + /// Returns lhs * 2^64 mod q as a Montgomery. + #[inline(always)] + pub fn prepare(&self, lhs: u64) -> Montgomery{ let mut rhs = Montgomery(0); self.prepare_assign(lhs, &mut rhs); rhs } - /// Assigns (lhs<<64)%q to rhs. + /// Assigns lhs * 2^64 mod q to rhs. #[inline(always)] - fn prepare_assign(&self, lhs: u64, rhs: &mut Montgomery){ + pub fn prepare_assign(&self, lhs: u64, rhs: &mut Montgomery){ self.prepare_lazy_assign(lhs, rhs); rhs.value_mut().reduce_once_assign(self.q); } - /// Returns (lhs<<64)%q in range [0, 2q-1] as a Montgomery. + /// Returns lhs * 2^64 mod q in range [0, 2q-1] as a Montgomery. #[inline(always)] - fn prepare_lazy(&self, lhs: u64) -> Montgomery{ + pub fn prepare_lazy(&self, lhs: u64) -> Montgomery{ let mut rhs = Montgomery(0); self.prepare_lazy_assign(lhs, &mut rhs); rhs } - /// Assigns (lhs<<64)%q in range [0, 2q-1] to rhs. + /// Assigns lhs * 2^64 mod q in range [0, 2q-1] to rhs. #[inline(always)] - fn prepare_lazy_assign(&self, lhs: u64, rhs: &mut Montgomery){ - let (_, mhi) = lhs.widening_mul(*self.q_barrett.value_lo()); - *rhs = Montgomery((lhs.wrapping_mul(*self.q_barrett.value_hi()).wrapping_add(mhi)).wrapping_mul(self.q).wrapping_neg()); + pub fn prepare_lazy_assign(&self, lhs: u64, rhs: &mut Montgomery){ + let (_, mhi) = lhs.widening_mul(*self.barrett.value_lo()); + *rhs = Montgomery((lhs.wrapping_mul(*self.barrett.value_hi()).wrapping_add(mhi)).wrapping_mul(self.q).wrapping_neg()); + } + + /// Returns lhs * (2^64)^-1 mod q as a u64. + #[inline(always)] + pub fn unprepare(&self, lhs: Montgomery) -> u64{ + let mut rhs = 0u64; + self.unprepare_assign(lhs, &mut rhs); + rhs + } + + /// Assigns lhs * (2^64)^-1 mod q to rhs. + #[inline(always)] + pub fn unprepare_assign(&self, lhs: Montgomery, rhs: &mut u64){ + self.unprepare_lazy_assign(lhs, rhs); + rhs.reduce_once_assign(self.q); + + } + + /// Returns lhs * (2^64)^-1 mod q in range [0, 2q-1]. + #[inline(always)] + pub fn unprepare_lazy(&self, lhs: Montgomery) -> u64{ + let mut rhs = 0u64; + self.unprepare_lazy_assign(lhs, &mut rhs); + rhs + } + + /// Assigns lhs * (2^64)^-1 mod q in range [0, 2q-1] to rhs. + #[inline(always)] + pub fn unprepare_lazy_assign(&self, lhs: Montgomery, rhs: &mut u64){ + let (_, r) = self.q.widening_mul(lhs.value().wrapping_mul(self.q_inv)); + *rhs = self.q - r } /// Returns lhs * rhs * (2^{64})^-1 mod q. #[inline(always)] - fn mul_external(&self, lhs: Montgomery, rhs: u64) -> u64{ + pub fn mul_external(&self, lhs: Montgomery, rhs: u64) -> u64{ let mut r = self.mul_external_lazy(lhs, rhs); r.reduce_once_assign(self.q); r @@ -93,14 +149,14 @@ impl MontgomeryPrecomp{ /// Assigns lhs * rhs * (2^{64})^-1 mod q to rhs. #[inline(always)] - fn mul_external_assign(&self, lhs: Montgomery, rhs: &mut u64){ + pub fn mul_external_assign(&self, lhs: Montgomery, rhs: &mut u64){ self.mul_external_lazy_assign(lhs, rhs); rhs.reduce_once_assign(self.q); } /// Returns lhs * rhs * (2^{64})^-1 mod q in range [0, 2q-1]. #[inline(always)] - fn mul_external_lazy(&self, lhs: Montgomery, rhs: u64) -> u64{ + pub fn mul_external_lazy(&self, lhs: Montgomery, rhs: u64) -> u64{ let mut result = rhs; self.mul_external_lazy_assign(lhs, &mut result); result @@ -108,7 +164,7 @@ impl MontgomeryPrecomp{ /// Assigns lhs * rhs * (2^{64})^-1 mod q in range [0, 2q-1] to rhs. #[inline(always)] - fn mul_external_lazy_assign(&self, lhs: Montgomery, rhs: &mut u64){ + pub fn mul_external_lazy_assign(&self, lhs: Montgomery, rhs: &mut u64){ let (mlo, mhi) = lhs.value().widening_mul(*rhs); let (_, hhi) = self.q.widening_mul(mlo.wrapping_mul(self.q_inv)); *rhs = mhi.wrapping_add(self.q - hhi) @@ -116,29 +172,83 @@ impl MontgomeryPrecomp{ /// Returns lhs * rhs * (2^{64})^-1 mod q in range [0, 2q-1]. #[inline(always)] - fn mul_internal(&self, lhs: Montgomery, rhs: Montgomery) -> Montgomery{ + pub fn mul_internal(&self, lhs: Montgomery, rhs: Montgomery) -> Montgomery{ Montgomery(self.mul_external(lhs, *rhs.value())) } /// Assigns lhs * rhs * (2^{64})^-1 mod q to rhs. #[inline(always)] - fn mul_internal_assign(&self, lhs: Montgomery, rhs: &mut Montgomery){ + pub fn mul_internal_assign(&self, lhs: Montgomery, rhs: &mut Montgomery){ self.mul_external_assign(lhs, rhs.value_mut()); } /// Returns lhs * rhs * (2^{64})^-1 mod q in range [0, 2q-1]. #[inline(always)] - fn mul_internal_lazy(&self, lhs: Montgomery, rhs: Montgomery) -> Montgomery{ + pub fn mul_internal_lazy(&self, lhs: Montgomery, rhs: Montgomery) -> Montgomery{ Montgomery(self.mul_external_lazy(lhs, *rhs.value())) } /// Assigns lhs * rhs * (2^{64})^-1 mod q in range [0, 2q-1] to rhs. #[inline(always)] - fn mul_internal_lazy_assign(&self, lhs: Montgomery, rhs: &mut Montgomery){ + pub fn mul_internal_lazy_assign(&self, lhs: Montgomery, rhs: &mut Montgomery){ self.mul_external_lazy_assign(lhs, rhs.value_mut()); } + + #[inline(always)] + pub fn add_internal(&self, lhs: Montgomery, rhs: Montgomery) -> Montgomery{ + Montgomery(self.barrett.reduce(rhs.value() + lhs.value())) + } + + /// Assigns lhs + rhs to rhs. + #[inline(always)] + pub fn add_internal_lazy_assign(&self, lhs: Montgomery, rhs: &mut Montgomery){ + *rhs.value_mut() += lhs.value() + } + + /// Assigns lhs + rhs - q if (lhs + rhs) >= q to rhs. + #[inline(always)] + pub fn add_internal_reduce_once_assign(&self, lhs: Montgomery, rhs: &mut Montgomery){ + self.add_internal_lazy_assign(lhs, rhs); + rhs.value_mut().reduce_once_assign(self.q); + } + + /// Returns (x^exponent) * 2^64 mod q. + #[inline(always)] + pub fn pow(&self, x: Montgomery, exponent:u64) -> Montgomery{ + let mut y: Montgomery = self.one(); + let mut x_mut: Montgomery = x; + let mut i: u64 = exponent; + while i > 0{ + if i & 1 == 1{ + self.mul_internal_assign(x_mut, &mut y); + } + self.mul_internal_assign(x_mut, &mut x_mut); + i >>= 1; + } + + y.value_mut().reduce_once_assign(self.q); + y + } } +/// Returns x^exponent mod q. +/// This function internally instantiate a new MontgomeryPrecomp +/// To be used when called only a few times and if there +/// is no Prime instantiated with q. +fn pow(x:u64, exponent:u64, q:u64) -> u64{ + let montgomery: MontgomeryPrecomp = MontgomeryPrecomp::::new(q); + let mut y_mont: Montgomery = montgomery.one(); + let mut x_mont: Montgomery = montgomery.prepare(x); + while exponent > 0{ + if exponent & 1 == 1{ + montgomery.mul_internal_assign(x_mont, &mut y_mont); + } + + montgomery.mul_internal_assign(x_mont, &mut x_mont); + } + + montgomery.unprepare(y_mont) +} #[cfg(test)] mod tests { diff --git a/src/modulus/prime.rs b/src/modulus/prime.rs index dad4ded..cd3f211 100644 --- a/src/modulus/prime.rs +++ b/src/modulus/prime.rs @@ -1,78 +1,230 @@ +use crate::modulus::montgomery::{MontgomeryPrecomp, Montgomery}; +use crate::dft::ntt::Table; use primality_test::is_prime; use prime_factorization::Factorization; pub struct Prime { - q: O, /// q_base^q_powers - q_base: O, - q_powers: O, - factors: Vec, /// distinct factors of q-1 - nth_root: O, + pub q: O, /// q_base^q_powers + pub q_base: O, + pub q_power: O, + pub factors: Vec, /// distinct factors of q-1 + pub montgomery: MontgomeryPrecomp, + pub phi: O, +} + +pub struct NTTFriendlyPrimesGenerator{ + size: f64, + next_prime: O, + prev_prime: O, + check_next_prime: bool, + check_prev_prime: bool, } impl Prime{ + + /// Returns a new instance of Prime. + /// Panics if q_base is not a prime > 2 and + /// if q_base^q_power would overflow u64. pub fn new(q_base: u64, q_power: u64) -> Self{ - assert!(is_prime(q) && q > 2); - assert!() - - q_exp - for i in 0..q_power{ - - } - - Self::new_unchecked(q) + assert!(is_prime(q_base) && q_base > 2); + Self::new_unchecked(q_base, q_power) } - pub fn new_unchecked(q: u64) -> Self { + /// Returns a new instance of Prime. + /// Does not check if q_base is a prime > 2. + /// Panics if q_base^q_power would overflow u64. + pub fn new_unchecked(q_base: u64, q_power: u64) -> Self { + + let mut q = q_base; + for _i in 1..q_power{ + q *= q_base + } + assert!(q.next_power_of_two().ilog2() <= 61); + + let mut phi = q_base-1; + for _i in 1..q_power{ + phi *= q_base + } + Self { - q, + q:q, + q_base:q_base, + q_power:q_power, + factors: Vec::new(), + montgomery:MontgomeryPrecomp::new(q), + phi:phi, } } - /// Returns returns Phi(BaseModulus^BaseModulusPower) - pub fn phi() -> u64 { - + pub fn q(&self) -> u64{ + self.q } - /// Returns the smallest primitive root. The unique factors - /// can be given as argument to avoid factorization of q-1. - pub fn primitive_root(&self) -> u64{ - if self.factors.len() != 0{ - self.check_factors(); - }else{ - let factors = Factorization::run(q).prime_factor_repr(); + pub fn q_base(&self) -> u64{ + self.q_base + } + + pub fn q_power(&self) -> u64{ + self.q_power + } + + /// Returns x^exponen mod q. + #[inline(always)] + pub fn pow(&self, x: u64, exponent: u64) -> u64{ + let mut y_mont: Montgomery = self.montgomery.one(); + let mut x_mont: Montgomery = self.montgomery.prepare(x); + let mut i: u64 = exponent; + while i > 0{ + if i & 1 == 1{ + self.montgomery.mul_internal_assign(x_mont, &mut y_mont); + } + + self.montgomery.mul_internal_assign(x_mont, &mut x_mont); + + i >>= 1; + + + } + + self.montgomery.unprepare(y_mont) + } + + /// Returns x^-1 mod q. + /// User must ensure that x is not divisible by q_base. + #[inline(always)] + pub fn inv(&self, x: u64) -> u64{ + self.pow(x, self.phi) + } +} + +/// Returns x^exponent mod q. +/// This function internally instantiate a new MontgomeryPrecomp +/// To be used when called only a few times and if there +/// is no Prime instantiated with q. +pub fn Pow(x:u64, exponent:u64, q:u64) -> u64{ + let montgomery: MontgomeryPrecomp = MontgomeryPrecomp::::new(q); + let mut y_mont: Montgomery = montgomery.one(); + let mut x_mont: Montgomery = montgomery.prepare(x); + let mut i: u64 = exponent; + while i > 0{ + if i & 1 == 1{ + montgomery.mul_internal_assign(x_mont, &mut y_mont); + } + + montgomery.mul_internal_assign(x_mont, &mut x_mont); + + i >>= 1; + } + + montgomery.unprepare(y_mont) +} + +impl Prime{ + /// Returns the smallest nth primitive root of q_base. + pub fn primitive_root(&mut self) -> u64{ + + self.check_factors(); + + let mut candidate: u64 = 1u64; + let mut not_found: bool = true; + + while not_found{ + + candidate += 1; + + for &factor in &self.factors{ + + if Pow(candidate, (self.q_base-1)/factor, self.q_base) == 1{ + not_found = true; + break + } + not_found = false; + } + } + + if not_found{ + panic!("failed to find a primitive root for q_base={}", self.q_base) + } + + candidate + } + + /// Returns an nth primitive root of q = q_base^q_power in Montgomery. + pub fn primitive_nth_root(&mut self, nth_root:u64) -> u64{ + + assert!(self.q & (nth_root-1) == 1, "invalid prime: q = {} % nth_root = {} = {} != 1", self.q, nth_root, self.q & (nth_root-1)); + + let psi: u64 = self.primitive_root(); + + // nth primitive root mod q_base: psi_nth^(prime.q_base-1)/nth_root mod q_base + let psi_nth_q_base: u64 = Pow(psi, (self.q_base-1)/nth_root, self.q_base); + + // lifts nth primitive root mod q_base to q = q_base^q_power + let psi_nth_q: u64 = self.hensel_lift(psi_nth_q_base, nth_root); + + assert!(self.pow(psi_nth_q, nth_root) == 1, "invalid nth primitive root: psi^nth_root != 1 mod q"); + assert!(self.pow(psi_nth_q, nth_root>>1) == self.q-1, "invalid nth primitive root: psi^(nth_root/2) != -1 mod q"); + + psi_nth_q + } + + /// Checks if the field self.factor is populated. + /// If not, factorize q_base-1 and populates self.factor. + /// If yes, checks that it contains the unique factors of q_base-1. + pub fn check_factors(&mut self){ + + if self.factors.len() == 0{ + + let factors = Factorization::run(self.q_base-1).prime_factor_repr(); let mut distincts_factors: Vec = Vec::with_capacity(factors.len()); for factor in factors.iter(){ distincts_factors.push(factor.0) } self.factors = distincts_factors - } + + }else{ + let mut q_base: u64 = self.q_base; - let log_nth_root = 64 - self.q.leading_zeros() as usize; - - 0 - } - - pub fn check_factors(&self){ - - if self.factors.len() == 0{ - panic!("invalid factor list: empty") - } - - let mut q = self.q; - - for &factor in &self.factors{ - if !is_prime(factor){ - panic!("invalid factor list: factor {} is not prime", factor) + for &factor in &self.factors{ + if !is_prime(factor){ + panic!("invalid factor list: factor {} is not prime", factor) + } + + while q_base%factor != 0{ + q_base /= factor + } } - - while q%factor != 0{ - q /= factor + + if q_base != 1{ + panic!("invalid factor list: does not fully divide q_base: q_base % (all factors) = {}", q_base) } } - if q != 1{ - panic!("invalid factor list: does not fully divide q: q % (alll factors) = {}", q) - } + } -} \ No newline at end of file + + /// Returns (psi + a * q_base)^{nth_root} = 1 mod q = q_base^q_power given psi^{nth_root} = 1 mod q_base. + /// Panics if psi^{nth_root} != 1 mod q_base. + fn hensel_lift(&self, psi: u64, nth_root: u64) -> u64{ + assert!(Pow(psi, nth_root, self.q_base)==1, "invalid argument psi: psi^nth_root = {} != 1", Pow(psi, nth_root, self.q_base)); + + let mut psi_mont: Montgomery = self.montgomery.prepare(psi); + let nth_root_mont: Montgomery = self.montgomery.prepare(nth_root); + + for _i in 1..self.q_power{ + + let psi_pow: Montgomery = self.montgomery.pow(psi_mont, nth_root-1); + + let num: Montgomery = Montgomery(self.montgomery.one().value() + self.q - self.montgomery.mul_internal(psi_pow, psi_mont).value()); + + let mut den: Montgomery = self.montgomery.mul_internal(nth_root_mont, psi_pow); + + den = self.montgomery.pow(den, self.phi-1); + + psi_mont = self.montgomery.add_internal(psi_mont, self.montgomery.mul_internal(num, den)); + } + + self.montgomery.unprepare(psi_mont) + } +}