[dft]: working NTT roots generation with prime power

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<u64>`
+    let q_base: u64 = 65537;      // Example prime base
+    let q_power: u64 = 2;     // Example power
+    let mut prime_instance: Prime<u64> = Prime::<u64>::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::<u64>::new(&mut prime_instance, n, nth_root);
+
+}

src/dft.rs

@@ -1 +1 @@
-pub(crate) mod primitive_root;
+pub mod ntt;

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<O>,
+    pub psi_forward_rev:Vec<Montgomery<O>>,
+    psi_backward_rev: Vec<Montgomery<O>>,
+    n_inv: Montgomery<O>,
+}
+
+impl<'a> Table<'a, u64> {
+    pub fn new(prime: &'a mut Prime<u64>, 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<u64> = prime.montgomery.prepare(psi);
+        let psi_inv_mont: Montgomery<u64> = prime.montgomery.pow(psi_mont, prime.phi-1);
+        
+        let mut psi_forward_rev: Vec<Montgomery<u64>> = vec![Montgomery(0); (nth_root >> 1) as usize];
+        let mut psi_backward_rev: Vec<Montgomery<u64>> = 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<u64> = 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,
+        }
+    }
+}

src/lib.rs

@@ -2,3 +2,4 @@
 #![feature(test)]
 
 pub mod modulus;
+pub mod dft;

src/modulus.rs

@@ -1,6 +1,6 @@
-pub(crate) mod prime;
+pub mod prime;
-pub(crate) mod montgomery;
+pub mod barrett;
-pub(crate) mod barrett;
+pub mod montgomery;
 
 trait ReduceOnce<O>{
     fn reduce_once_assign(&mut self, q: O);
@@ -8,12 +8,13 @@ trait ReduceOnce<O>{
 }
 
 impl ReduceOnce<u64> 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

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, O);
+pub struct BarrettPrecomp<O>{
+    q: O,
+    lo:O,
+    hi:O,
+}
 
 impl<O> BarrettPrecomp<O>{
 
     #[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<u64>{
         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)
     }
 }

src/modulus/generation.rs

@@ -3,16 +3,8 @@ use crate::modulus::prime;
 use prime::Prime;
 use primality_test::is_prime;
 
-pub struct NTTFriendlyPrimesGenerator{
+impl NTTFriendlyPrimesGenerator<u64> {
-	size: f64,
-	next_prime: u64,
-	prev_prime: u64,
-	nth_root: u64,
-	check_next_prime: bool,
-	check_prev_prime: bool,
-}
 	
-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;

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>(O);
+pub struct Montgomery<O>(pub O);
 
 /// Implements helper methods on the struct Montgomery<O>.
 impl<O> Montgomery<O>{
@@ -31,8 +31,10 @@ impl<O> Montgomery<O>{
 #[derive(Clone, Copy, Debug, PartialEq, Eq)]
 pub struct MontgomeryPrecomp<O>{
     q: O,
-    q_barrett: BarrettPrecomp<O>,
+    barrett: BarrettPrecomp<O>,
     q_inv: O,
+    one: Montgomery<O>,
+    minus_one: Montgomery<O>,
 }
 
 /// MontgomeryPrecomp is a set of methods implemented for MontgomeryPrecomp<u64>
@@ -42,7 +44,7 @@ impl MontgomeryPrecomp<u64>{
     /// Returns an new instance of MontgomeryPrecomp<u64>.
     /// This method will fail if gcd(q, 2^64) != 1.
     #[inline(always)]
-    fn new(q: u64) -> MontgomeryPrecomp<u64>{
+    pub fn new(q: u64) -> MontgomeryPrecomp<u64>{
         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<u64>{
             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<u64>.
+    /// Returns 2^64 mod q as a Montgomery<u64>.
     #[inline(always)]
-    fn prepare(&self, lhs: u64) -> Montgomery<u64>{
+    pub fn one(&self) -> Montgomery<u64>{
+        self.one
+    }
+
+    /// Returns (q-1) * 2^64 mod q as a Montgomery<u64>.
+    #[inline(always)]
+    pub fn minus_one(&self) -> Montgomery<u64>{
+        self.minus_one
+    }
+
+    /// Returns lhs * 2^64 mod q as a Montgomery<u64>.
+    #[inline(always)]
+    pub fn prepare(&self, lhs: u64) -> Montgomery<u64>{
         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<u64>){
+    pub fn prepare_assign(&self, lhs: u64, rhs: &mut Montgomery<u64>){
         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<u64>.
+    /// Returns lhs * 2^64 mod q in range [0, 2q-1] as a Montgomery<u64>.
     #[inline(always)]
-    fn prepare_lazy(&self, lhs: u64) -> Montgomery<u64>{
+    pub fn prepare_lazy(&self, lhs: u64) -> Montgomery<u64>{
         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<u64>){
+    pub fn prepare_lazy_assign(&self, lhs: u64, rhs: &mut Montgomery<u64>){
-        let (_, mhi) = lhs.widening_mul(*self.q_barrett.value_lo());
+        let (_, mhi) = lhs.widening_mul(*self.barrett.value_lo());
-        *rhs = Montgomery((lhs.wrapping_mul(*self.q_barrett.value_hi()).wrapping_add(mhi)).wrapping_mul(self.q).wrapping_neg());
+        *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>) -> 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<u64>, 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>) -> 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<u64>, 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<u64>, rhs: u64) -> u64{
+    pub fn mul_external(&self, lhs: Montgomery<u64>, 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<u64>{
 
     /// Assigns lhs * rhs * (2^{64})^-1 mod q to rhs.
     #[inline(always)]
-    fn mul_external_assign(&self, lhs: Montgomery<u64>, rhs: &mut u64){
+    pub fn mul_external_assign(&self, lhs: Montgomery<u64>, 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<u64>, rhs: u64) -> u64{
+    pub fn mul_external_lazy(&self, lhs: Montgomery<u64>, rhs: u64) -> u64{
         let mut result = rhs;
         self.mul_external_lazy_assign(lhs, &mut result);
         result
@@ -108,7 +164,7 @@ impl MontgomeryPrecomp<u64>{
 
     /// 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<u64>, rhs: &mut u64){
+    pub fn mul_external_lazy_assign(&self, lhs: Montgomery<u64>, 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<u64>{
 
     /// Returns lhs * rhs * (2^{64})^-1 mod q in range [0, 2q-1].
     #[inline(always)]
-    fn mul_internal(&self, lhs: Montgomery<u64>, rhs: Montgomery<u64>) -> Montgomery<u64>{
+    pub fn mul_internal(&self, lhs: Montgomery<u64>, rhs: Montgomery<u64>) -> Montgomery<u64>{
         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<u64>, rhs: &mut Montgomery<u64>){
+    pub fn mul_internal_assign(&self, lhs: Montgomery<u64>, rhs: &mut Montgomery<u64>){
         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<u64>, rhs: Montgomery<u64>) -> Montgomery<u64>{
+    pub fn mul_internal_lazy(&self, lhs: Montgomery<u64>, rhs: Montgomery<u64>) -> Montgomery<u64>{
         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<u64>, rhs: &mut Montgomery<u64>){
+    pub fn mul_internal_lazy_assign(&self, lhs: Montgomery<u64>, rhs: &mut Montgomery<u64>){
         self.mul_external_lazy_assign(lhs, rhs.value_mut());
     }
+
+    #[inline(always)]
+    pub fn add_internal(&self, lhs: Montgomery<u64>, rhs: Montgomery<u64>) -> Montgomery<u64>{
+        Montgomery(self.barrett.reduce(rhs.value() + lhs.value()))
+    }
+
+    /// Assigns lhs + rhs to rhs.
+    #[inline(always)]
+    pub fn add_internal_lazy_assign(&self, lhs: Montgomery<u64>, rhs: &mut Montgomery<u64>){
+        *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<u64>, rhs: &mut Montgomery<u64>){
+        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<u64>, exponent:u64) -> Montgomery<u64>{
+        let mut y: Montgomery<u64> = self.one();
+        let mut x_mut: Montgomery<u64> = 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<u64>
+/// 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<u64> = MontgomeryPrecomp::<u64>::new(q);
+    let mut y_mont: Montgomery<u64> = montgomery.one();
+    let mut x_mont: Montgomery<u64> = 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 {

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<O> {
-	q: O, /// q_base^q_powers
+	pub q: O, /// q_base^q_powers
-    q_base: O,
+    pub q_base: O,
-    q_powers: O,
+    pub q_power: O,
-    factors: Vec<O>, /// distinct factors of q-1
+    pub factors: Vec<O>, /// distinct factors of q-1
-    nth_root: O,
+    pub montgomery: MontgomeryPrecomp<O>,
+    pub phi: O,
+}
+
+pub struct NTTFriendlyPrimesGenerator<O>{
+	size: f64,
+	next_prime: O,
+	prev_prime: O,
+	check_next_prime: bool,
+	check_prev_prime: bool,
 }
 
 impl Prime<u64>{
+
+    /// Returns a new instance of Prime<u64>.
+    /// 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!(is_prime(q_base) && q_base > 2);
-        assert!()
+        Self::new_unchecked(q_base, q_power)
-
-        q_exp
-        for i in 0..q_power{
-
 	}
 
-        Self::new_unchecked(q)
+    /// Returns a new instance of Prime<u64>.
+    /// 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
         }
 
-	pub fn new_unchecked(q: u64) -> Self {
         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 q(&self) -> u64{
-    pub fn phi() -> u64 {
+        self.q
+    }
+
+    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<u64> = self.montgomery.one();
+        let mut x_mont: Montgomery<u64> = 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;
+
             
         }
         
-    /// Returns the smallest primitive root. The unique factors
+        self.montgomery.unprepare(y_mont)
-    /// can be given as argument to avoid factorization of q-1.
+    }
-    pub fn primitive_root(&self) -> u64{
+
-        if self.factors.len() != 0{
+    /// 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<u64>
+/// 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<u64> = MontgomeryPrecomp::<u64>::new(q);
+    let mut y_mont: Montgomery<u64> = montgomery.one();
+    let mut x_mont: Montgomery<u64> = 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<u64>{
+    /// Returns the smallest nth primitive root of q_base.
+    pub fn primitive_root(&mut self) -> u64{
+
         self.check_factors();
-        }else{
+
-            let factors = Factorization::run(q).prime_factor_repr();
+        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<u64> = Vec::with_capacity(factors.len());
             for factor in factors.iter(){
                 distincts_factors.push(factor.0)
             }
             self.factors = distincts_factors
-        }
             
-        let log_nth_root = 64 - self.q.leading_zeros() as usize;
+        }else{
-
+            let mut q_base: u64 = self.q_base;
-        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)
                 }
     
-            while q%factor != 0{
+                while q_base%factor != 0{
-                q /= factor
+                    q_base /= factor
                 }
             }
     
-        if q != 1{
+            if q_base != 1{
-            panic!("invalid factor list: does not fully divide q: q % (alll factors) = {}", q)
+                panic!("invalid factor list: does not fully divide q_base: q_base % (all factors) = {}", q_base)
             }
         }
+
+        
+    }
+
+    /// 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<u64> = self.montgomery.prepare(psi);
+        let nth_root_mont: Montgomery<u64> = self.montgomery.prepare(nth_root);
+
+        for _i in 1..self.q_power{
+
+            let psi_pow: Montgomery<u64> = self.montgomery.pow(psi_mont, nth_root-1);
+
+            let num: Montgomery<u64> = Montgomery(self.montgomery.one().value() + self.q - self.montgomery.mul_internal(psi_pow, psi_mont).value());
+
+            let mut den: Montgomery<u64> = 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)
+    }
 }