arnaucube
arnaucube/poulpyPublic Self-hosted repository
poulpy/math/src/modulus/impl_u64/montgomery.rs
math/src/modulus/impl_u64/montgomery.rsRaw
  1
  2use crate::modulus::ReduceOnce;
  3use crate::modulus::montgomery::{MontgomeryPrecomp, Montgomery};
  4use crate::modulus::barrett::BarrettPrecomp;
  5use crate::modulus::{REDUCEMOD, NONE, ONCE, TWICE, FOURTIMES, BARRETT, BARRETTLAZY};
  6extern crate test;
  7
  8/// MontgomeryPrecomp is a set of methods implemented for MontgomeryPrecomp<u64>
  9/// enabling Montgomery arithmetic over u64 values.
 10impl MontgomeryPrecomp<u64>{
 11
 12    /// Returns an new instance of MontgomeryPrecomp<u64>.
 13    /// This method will fail if gcd(q, 2^64) != 1.
 14    #[inline(always)]
 15    pub fn new(q: u64) -> MontgomeryPrecomp<u64>{
 16        assert!(q & 1 != 0, "Invalid argument: gcd(q={}, radix=2^64) != 1", q);
 17        let mut q_inv: u64 = 1;
 18        let mut q_pow = q;
 19        for _i in 0..63{
 20            q_inv = q_inv.wrapping_mul(q_pow);
 21            q_pow = q_pow.wrapping_mul(q_pow);
 22        }
 23        let mut precomp = Self{ 
 24            q: q,
 25            two_q: q<<1,
 26            four_q: q<<2,
 27            barrett: BarrettPrecomp::new(q), 
 28            q_inv: q_inv,
 29            one: 0,
 30            minus_one:0,
 31        };
 32
 33        precomp.one = precomp.prepare::<ONCE>(1);
 34        precomp.minus_one = q-precomp.one;
 35
 36        precomp
 37    }
 38
 39    /// Returns 2^64 mod q as a Montgomery<u64>.
 40    #[inline(always)]
 41    pub fn one(&self) -> Montgomery<u64>{
 42        self.one
 43    }
 44
 45    /// Returns (q-1) * 2^64 mod q as a Montgomery<u64>.
 46    #[inline(always)]
 47    pub fn minus_one(&self) -> Montgomery<u64>{
 48        self.minus_one
 49    }
 50
 51    /// Applies a modular reduction on x based on REDUCE:
 52    /// - LAZY: no modular reduction.
 53    /// - ONCE: subtracts q if x >= q.
 54    /// - FULL: maps x to x mod q using Barrett reduction.
 55    #[inline(always)]
 56    pub fn reduce<const REDUCE:REDUCEMOD>(&self, x: u64) -> u64{
 57        let mut r: u64 = x;
 58        self.reduce_assign::<REDUCE>(&mut r);
 59        r
 60    }
 61
 62    /// Applies a modular reduction on x based on REDUCE:
 63    /// - LAZY: no modular reduction.
 64    /// - ONCE: subtracts q if x >= q.
 65    /// - FULL: maps x to x mod q using Barrett reduction.
 66    #[inline(always)]
 67    pub fn reduce_assign<const REDUCE:REDUCEMOD>(&self, x: &mut u64){
 68        match REDUCE {
 69            NONE =>{},
 70            ONCE =>{x.reduce_once_assign(self.q)},
 71            TWICE=>{x.reduce_once_assign(self.two_q)},
 72            FOURTIMES =>{x.reduce_once_assign(self.four_q)},
 73            BARRETT =>{self.barrett.reduce_assign(x)},
 74            BARRETTLAZY =>{self.barrett.reduce_lazy_assign(x)},
 75            _ => unreachable!("invalid REDUCE argument")
 76        }
 77    }
 78
 79    /// Returns lhs * 2^64 mod q as a Montgomery<u64>.
 80    #[inline(always)]
 81    pub fn prepare<const REDUCE:REDUCEMOD>(&self, lhs: u64) -> Montgomery<u64>{
 82        let mut rhs: u64 = 0;
 83        self.prepare_assign::<REDUCE>(lhs, &mut rhs);
 84        rhs
 85    }
 86
 87    /// Assigns lhs * 2^64 mod q to rhs.
 88    #[inline(always)]
 89    pub fn prepare_assign<const REDUCE:REDUCEMOD>(&self, lhs: u64, rhs: &mut Montgomery<u64>){
 90        let (_, mhi) = lhs.widening_mul(*self.barrett.value_lo());
 91        *rhs = (lhs.wrapping_mul(*self.barrett.value_hi()).wrapping_add(mhi)).wrapping_mul(self.q).wrapping_neg();
 92        self.reduce_assign::<REDUCE>(rhs); 
 93    }
 94
 95    /// Returns lhs * (2^64)^-1 mod q as a u64.
 96    #[inline(always)]
 97    pub fn unprepare<const REDUCE:REDUCEMOD>(&self, lhs: Montgomery<u64>) -> u64{
 98        let mut rhs = 0u64;
 99        self.unprepare_assign::<REDUCE>(lhs, &mut rhs);
100        rhs
101    }
102
103    /// Assigns lhs * (2^64)^-1 mod q to rhs.
104    #[inline(always)]
105    pub fn unprepare_assign<const REDUCE:REDUCEMOD>(&self, lhs: Montgomery<u64>, rhs: &mut u64){
106        let (_, r) = self.q.widening_mul(lhs.wrapping_mul(self.q_inv));
107        *rhs = self.reduce::<REDUCE>(self.q.wrapping_sub(r));        
108    }
109
110    /// Returns lhs * rhs * (2^{64})^-1 mod q.
111    #[inline(always)]
112    pub fn mul_external<const REDUCE:REDUCEMOD>(&self, lhs: Montgomery<u64>, rhs: u64) -> u64{
113        let mut r: u64 = rhs;
114        self.mul_external_assign::<REDUCE>(lhs, &mut r);
115        r
116    }
117
118    /// Assigns lhs * rhs * (2^{64})^-1 mod q to rhs.
119    #[inline(always)]
120    pub fn mul_external_assign<const REDUCE:REDUCEMOD>(&self, lhs: Montgomery<u64>, rhs: &mut u64){
121        let (mlo, mhi) = lhs.widening_mul(*rhs);
122        let (_, hhi) = self.q.widening_mul(mlo.wrapping_mul(self.q_inv));
123        *rhs = self.reduce::<REDUCE>(mhi.wrapping_sub(hhi).wrapping_add(self.q));
124    }
125
126    /// Returns lhs * rhs * (2^{64})^-1 mod q in range [0, 2q-1].
127    #[inline(always)]
128    pub fn mul_internal<const REDUCE:REDUCEMOD>(&self, lhs: Montgomery<u64>, rhs: Montgomery<u64>) -> Montgomery<u64>{
129        self.mul_external::<REDUCE>(lhs, rhs)
130    }
131
132    /// Assigns lhs * rhs * (2^{64})^-1 mod q to rhs.
133    #[inline(always)]
134    pub fn mul_internal_assign<const REDUCE:REDUCEMOD>(&self, lhs: Montgomery<u64>, rhs: &mut Montgomery<u64>){
135        self.mul_external_assign::<REDUCE>(lhs, rhs);
136    }
137
138    #[inline(always)]
139    pub fn add_internal(&self, lhs: Montgomery<u64>, rhs: Montgomery<u64>) -> Montgomery<u64>{
140        self.barrett.reduce(rhs + lhs)
141    }
142
143    /// Assigns lhs + rhs to rhs.
144    #[inline(always)]
145    pub fn add_internal_lazy_assign(&self, lhs: Montgomery<u64>, rhs: &mut Montgomery<u64>){
146        *rhs += lhs
147    }
148
149    /// Assigns lhs + rhs - q if (lhs + rhs) >= q to rhs.
150    #[inline(always)]
151    pub fn add_internal_reduce_once_assign<const LAZY:bool>(&self, lhs: Montgomery<u64>, rhs: &mut Montgomery<u64>){
152        self.add_internal_lazy_assign(lhs, rhs);
153        rhs.reduce_once_assign(self.q);
154    }
155
156    /// Returns lhs mod q in range [0, 2q-1].
157    #[inline(always)]
158    pub fn reduce_lazy_assign(&self, lhs: &mut u64){
159        self.barrett.reduce_lazy_assign(lhs)
160    }
161
162    /// Returns (x^exponent) * 2^64 mod q.
163    #[inline(always)]
164    pub fn pow(&self, x: Montgomery<u64>, exponent:u64) -> Montgomery<u64>{
165        let mut y: Montgomery<u64> = self.one();
166        let mut x_mut: Montgomery<u64> = x;
167        let mut i: u64 = exponent;
168        while i > 0{
169            if i & 1 == 1{
170                self.mul_internal_assign::<ONCE>(x_mut, &mut y);
171            }
172            self.mul_internal_assign::<ONCE>(x_mut, &mut x_mut);
173            i >>= 1;
174        }
175
176        y.reduce_once_assign(self.q);
177        y
178    }
179}
180
181#[cfg(test)]
182mod tests {
183    use crate::modulus::montgomery;
184    use super::*;
185    use test::Bencher;
186
187    #[test]
188    fn test_mul_external() {
189        let q: u64 = 0x1fffffffffe00001;
190    	let m_precomp = montgomery::MontgomeryPrecomp::new(q);
191        let x: u64 = 0x5f876e514845cc8b;
192        let y: u64 = 0xad726f98f24a761a;
193        let y_mont = m_precomp.prepare::<ONCE>(y);
194    	assert!(m_precomp.mul_external::<ONCE>(y_mont, x) == (x as u128 * y as u128 % q as u128) as u64);
195    }
196
197    #[bench]
198    fn bench_mul_external(b: &mut Bencher){
199        let q: u64 = 0x1fffffffffe00001;
200    	let m_precomp = montgomery::MontgomeryPrecomp::new(q);
201        let mut x: u64 = 0x5f876e514845cc8b;
202        let y: u64 = 0xad726f98f24a761a;
203        let y_mont = m_precomp.prepare::<ONCE>(y);
204        b.iter(|| m_precomp.mul_external_assign::<ONCE>(y_mont, &mut x));
205    }
206}