1use crate::modulus::prime::Prime;
  2use crate::modulus::montgomery::{Montgomery, MontgomeryPrecomp};
  3use crate::modulus::shoup::{ShoupPrecomp};
  4use crate::modulus::ONCE;
  5use primality_test::is_prime;
  6use prime_factorization::Factorization;
  7
  8impl Prime<u64>{
  9
 10    /// Returns a new instance of Prime<u64>.
 11    /// Panics if q_base is not a prime > 2 and
 12    /// if q_base^q_power would overflow u64.
 13    pub fn new(q_base: u64, q_power: usize) -> Self{
 14		assert!(is_prime(q_base) && q_base > 2);
 15        Self::new_unchecked(q_base, q_power)
 16	}
 17
 18    /// Returns a new instance of Prime<u64>.
 19    /// Does not check if q_base is a prime > 2.
 20    /// Panics if q_base^q_power would overflow u64.
 21	pub fn new_unchecked(q_base: u64, q_power: usize) -> Self {
 22
 23        let mut q = q_base;
 24        for _i in 1..q_power{
 25            q *= q_base
 26        }
 27
 28        assert!(q.next_power_of_two().ilog2() <= 61);
 29
 30        let mut phi = q_base-1;
 31        for _i in 1..q_power{
 32            phi *= q_base
 33        }
 34
 35        let mut prime: Prime<u64> = Self {
 36            q:q,
 37            q_base:q_base,
 38            q_power:q_power,
 39            factors: Vec::new(),
 40            montgomery:MontgomeryPrecomp::new(q),
 41            shoup:ShoupPrecomp::new(q),
 42            phi:phi,
 43        };
 44
 45        prime.check_factors();
 46
 47        prime
 48
 49    }
 50
 51    pub fn q(&self) -> u64{
 52        self.q
 53    }
 54
 55    pub fn q_base(&self) -> u64{
 56        self.q_base
 57    }
 58
 59    pub fn q_power(&self) -> usize{
 60        self.q_power
 61    }
 62
 63    /// Returns x^exponen mod q.
 64    #[inline(always)]
 65    pub fn pow(&self, x: u64, exponent: u64) -> u64{
 66        let mut y_mont: Montgomery<u64> = self.montgomery.one();
 67        let mut x_mont: Montgomery<u64> = self.montgomery.prepare::<ONCE>(x);
 68        let mut i: u64 = exponent;
 69        while i > 0{
 70            if i & 1 == 1{
 71                self.montgomery.mul_internal_assign::<ONCE>(x_mont, &mut y_mont);
 72            }
 73
 74            self.montgomery.mul_internal_assign::<ONCE>(x_mont, &mut x_mont);
 75
 76            i >>= 1;
 77
 78        }
 79        
 80        self.montgomery.unprepare::<ONCE>(y_mont)
 81    }
 82
 83    /// Returns x^-1 mod q.
 84    /// User must ensure that x is not divisible by q_base.
 85    #[inline(always)]
 86    pub fn inv(&self, x: u64) -> u64{
 87        self.pow(x, self.phi-1)
 88    }
 89}
 90
 91impl Prime<u64>{
 92    /// Returns the smallest nth primitive root of q_base.
 93    pub fn primitive_root(&self) -> u64{
 94
 95        let mut candidate: u64 = 1u64;
 96        let mut not_found: bool = true;
 97
 98        while not_found{
 99
100            candidate += 1;
101
102            for &factor in &self.factors{
103
104                if  Pow(candidate, (self.q_base-1)/factor, self.q_base) == 1{
105                    not_found = true;
106                    break
107                }
108                not_found = false;
109            }
110        }
111
112        if not_found{
113            panic!("failed to find a primitive root for q_base={}", self.q_base)
114        }
115
116        candidate
117    }
118
119    /// Returns an nth primitive root of q = q_base^q_power in Montgomery.
120    pub fn primitive_nth_root(&self, nth_root:u64) -> u64{
121
122        assert!(self.q & (nth_root-1) == 1, "invalid prime: q = {} % nth_root = {} = {} != 1", self.q, nth_root, self.q & (nth_root-1));
123
124        let psi: u64 = self.primitive_root();
125
126        // nth primitive root mod q_base: psi_nth^(prime.q_base-1)/nth_root mod q_base
127        let psi_nth_q_base: u64 = Pow(psi, (self.q_base-1)/nth_root, self.q_base);
128
129        // lifts nth primitive root mod q_base to q = q_base^q_power
130        let psi_nth_q: u64 = self.hensel_lift(psi_nth_q_base, nth_root);
131
132        assert!(self.pow(psi_nth_q, nth_root) == 1, "invalid nth primitive root: psi^nth_root != 1 mod q");
133        assert!(self.pow(psi_nth_q, nth_root>>1) == self.q-1, "invalid nth primitive root: psi^(nth_root/2) != -1 mod q");
134
135        psi_nth_q
136    }
137
138    /// Checks if the field self.factor is populated.
139    /// If not, factorize q_base-1 and populates self.factor.
140    /// If yes, checks that it contains the unique factors of q_base-1.
141    pub fn check_factors(&mut self){
142
143        if self.factors.len() == 0{
144
145            let factors = Factorization::run(self.q_base-1).prime_factor_repr();
146            let mut distincts_factors: Vec<u64> = Vec::with_capacity(factors.len());
147            for factor in factors.iter(){
148                distincts_factors.push(factor.0)
149            }
150            self.factors = distincts_factors
151            
152        }else{
153            let mut q_base: u64 = self.q_base;
154
155            for &factor in &self.factors{
156                if !is_prime(factor){
157                    panic!("invalid factor list: factor {} is not prime", factor)
158                }
159    
160                while q_base%factor != 0{
161                    q_base /= factor
162                }
163            }
164    
165            if q_base != 1{
166                panic!("invalid factor list: does not fully divide q_base: q_base % (all factors) = {}", q_base)
167            }
168        }        
169    }
170
171    /// Returns (psi + a * q_base)^{nth_root} = 1 mod q = q_base^q_power given psi^{nth_root} = 1 mod q_base.
172    /// Panics if psi^{nth_root} != 1 mod q_base.
173    fn hensel_lift(&self, psi: u64, nth_root: u64) -> u64{
174        assert!(Pow(psi, nth_root, self.q_base)==1, "invalid argument psi: psi^nth_root = {} != 1", Pow(psi, nth_root, self.q_base));
175
176        let mut psi_mont: Montgomery<u64> = self.montgomery.prepare::<ONCE>(psi);
177        let nth_root_mont: Montgomery<u64> = self.montgomery.prepare::<ONCE>(nth_root);
178
179        for _i in 1..self.q_power{
180
181            let psi_pow: Montgomery<u64> = self.montgomery.pow(psi_mont, nth_root-1);
182
183            let num: Montgomery<u64> = Montgomery(self.montgomery.one().value() + self.q - self.montgomery.mul_internal::<ONCE>(psi_pow, psi_mont).value());
184
185            let mut den: Montgomery<u64> = self.montgomery.mul_internal::<ONCE>(nth_root_mont, psi_pow);
186
187            den = self.montgomery.pow(den, self.phi-1);
188
189            psi_mont = self.montgomery.add_internal(psi_mont, self.montgomery.mul_internal::<ONCE>(num, den));
190        }
191        
192        self.montgomery.unprepare::<ONCE>(psi_mont)
193    }
194}
195
196/// Returns x^exponent mod q.
197/// This function internally instantiate a new MontgomeryPrecomp<u64>
198/// To be used when called only a few times and if there 
199/// is no Prime instantiated with q.
200pub fn Pow(x:u64, exponent:u64, q:u64) -> u64{
201    let montgomery: MontgomeryPrecomp<u64> = MontgomeryPrecomp::<u64>::new(q);
202    let mut y_mont: Montgomery<u64> = montgomery.one();
203    let mut x_mont: Montgomery<u64> = montgomery.prepare::<ONCE>(x);
204    let mut i: u64 = exponent;
205    while i > 0{
206        if i & 1 == 1{
207            montgomery.mul_internal_assign::<ONCE>(x_mont, &mut y_mont);
208        }
209
210        montgomery.mul_internal_assign::<ONCE>(x_mont, &mut x_mont);
211
212        i >>= 1;
213    }
214    
215    montgomery.unprepare::<ONCE>(y_mont)
216}