poulpy/src/modulus/impl_u64/montgomery.rs
src/modulus/impl_u64/montgomery.rsRaw
1
2use crate::modulus::ReduceOnce;
3use crate::modulus::montgomery::{MontgomeryPrecomp, Montgomery};
4use crate::modulus::barrett::{BarrettPrecomp};
5extern crate test;
6
7/// MontgomeryPrecomp is a set of methods implemented for MontgomeryPrecomp<u64>
8/// enabling Montgomery arithmetic over u64 values.
9impl MontgomeryPrecomp<u64>{
10
11 /// Returns an new instance of MontgomeryPrecomp<u64>.
12 /// This method will fail if gcd(q, 2^64) != 1.
13 #[inline(always)]
14 pub fn new(q: u64) -> MontgomeryPrecomp<u64>{
15 assert!(q & 1 != 0, "Invalid argument: gcd(q={}, radix=2^64) != 1", q);
16 let mut q_inv: u64 = 1;
17 let mut q_pow = q;
18 for _i in 0..63{
19 q_inv = q_inv.wrapping_mul(q_pow);
20 q_pow = q_pow.wrapping_mul(q_pow);
21 }
22 let mut precomp = Self{
23 q: q,
24 barrett: BarrettPrecomp::new(q),
25 q_inv: q_inv,
26 one: Montgomery(0),
27 minus_one: Montgomery(0),
28 };
29
30 precomp.one = precomp.prepare(1);
31 precomp.minus_one = Montgomery(q-precomp.one.value());
32
33 precomp
34 }
35
36 /// Returns 2^64 mod q as a Montgomery<u64>.
37 #[inline(always)]
38 pub fn one(&self) -> Montgomery<u64>{
39 self.one
40 }
41
42 /// Returns (q-1) * 2^64 mod q as a Montgomery<u64>.
43 #[inline(always)]
44 pub fn minus_one(&self) -> Montgomery<u64>{
45 self.minus_one
46 }
47
48 /// Returns lhs * 2^64 mod q as a Montgomery<u64>.
49 #[inline(always)]
50 pub fn prepare(&self, lhs: u64) -> Montgomery<u64>{
51 let mut rhs = Montgomery(0);
52 self.prepare_assign(lhs, &mut rhs);
53 rhs
54 }
55
56 /// Assigns lhs * 2^64 mod q to rhs.
57 #[inline(always)]
58 pub fn prepare_assign(&self, lhs: u64, rhs: &mut Montgomery<u64>){
59 self.prepare_lazy_assign(lhs, rhs);
60 rhs.value_mut().reduce_once_assign(self.q);
61 }
62
63 /// Returns lhs * 2^64 mod q in range [0, 2q-1] as a Montgomery<u64>.
64 #[inline(always)]
65 pub fn prepare_lazy(&self, lhs: u64) -> Montgomery<u64>{
66 let mut rhs = Montgomery(0);
67 self.prepare_lazy_assign(lhs, &mut rhs);
68 rhs
69 }
70
71 /// Assigns lhs * 2^64 mod q in range [0, 2q-1] to rhs.
72 #[inline(always)]
73 pub fn prepare_lazy_assign(&self, lhs: u64, rhs: &mut Montgomery<u64>){
74 let (_, mhi) = lhs.widening_mul(*self.barrett.value_lo());
75 *rhs = Montgomery((lhs.wrapping_mul(*self.barrett.value_hi()).wrapping_add(mhi)).wrapping_mul(self.q).wrapping_neg());
76 }
77
78 /// Returns lhs * (2^64)^-1 mod q as a u64.
79 #[inline(always)]
80 pub fn unprepare(&self, lhs: Montgomery<u64>) -> u64{
81 let mut rhs = 0u64;
82 self.unprepare_assign(lhs, &mut rhs);
83 rhs
84 }
85
86 /// Assigns lhs * (2^64)^-1 mod q to rhs.
87 #[inline(always)]
88 pub fn unprepare_assign(&self, lhs: Montgomery<u64>, rhs: &mut u64){
89 self.unprepare_lazy_assign(lhs, rhs);
90 rhs.reduce_once_assign(self.q);
91
92 }
93
94 /// Returns lhs * (2^64)^-1 mod q in range [0, 2q-1].
95 #[inline(always)]
96 pub fn unprepare_lazy(&self, lhs: Montgomery<u64>) -> u64{
97 let mut rhs = 0u64;
98 self.unprepare_lazy_assign(lhs, &mut rhs);
99 rhs
100 }
101
102 /// Assigns lhs * (2^64)^-1 mod q in range [0, 2q-1] to rhs.
103 #[inline(always)]
104 pub fn unprepare_lazy_assign(&self, lhs: Montgomery<u64>, rhs: &mut u64){
105 let (_, r) = self.q.widening_mul(lhs.value().wrapping_mul(self.q_inv));
106 *rhs = self.q - r
107 }
108
109 /// Returns lhs * rhs * (2^{64})^-1 mod q.
110 #[inline(always)]
111 pub fn mul_external(&self, lhs: Montgomery<u64>, rhs: u64) -> u64{
112 let mut r = self.mul_external_lazy(lhs, rhs);
113 r.reduce_once_assign(self.q);
114 r
115 }
116
117 /// Assigns lhs * rhs * (2^{64})^-1 mod q to rhs.
118 #[inline(always)]
119 pub fn mul_external_assign(&self, lhs: Montgomery<u64>, rhs: &mut u64){
120 self.mul_external_lazy_assign(lhs, rhs);
121 rhs.reduce_once_assign(self.q);
122 }
123
124 /// Returns lhs * rhs * (2^{64})^-1 mod q in range [0, 2q-1].
125 #[inline(always)]
126 pub fn mul_external_lazy(&self, lhs: Montgomery<u64>, rhs: u64) -> u64{
127 let mut result: u64 = rhs;
128 self.mul_external_lazy_assign(lhs, &mut result);
129 result
130 }
131
132 /// Assigns lhs * rhs * (2^{64})^-1 mod q in range [0, 2q-1] to rhs.
133 #[inline(always)]
134 pub fn mul_external_lazy_assign(&self, lhs: Montgomery<u64>, rhs: &mut u64){
135 let (mlo, mhi) = lhs.value().widening_mul(*rhs);
136 let (_, hhi) = self.q.widening_mul(mlo.wrapping_mul(self.q_inv));
137 *rhs = mhi.wrapping_add(self.q - hhi)
138 }
139
140 /// Returns lhs * rhs * (2^{64})^-1 mod q in range [0, 2q-1].
141 #[inline(always)]
142 pub fn mul_internal(&self, lhs: Montgomery<u64>, rhs: Montgomery<u64>) -> Montgomery<u64>{
143 Montgomery(self.mul_external(lhs, *rhs.value()))
144 }
145
146 /// Assigns lhs * rhs * (2^{64})^-1 mod q to rhs.
147 #[inline(always)]
148 pub fn mul_internal_assign(&self, lhs: Montgomery<u64>, rhs: &mut Montgomery<u64>){
149 self.mul_external_assign(lhs, rhs.value_mut());
150 }
151
152 /// Returns lhs * rhs * (2^{64})^-1 mod q in range [0, 2q-1].
153 #[inline(always)]
154 pub fn mul_internal_lazy(&self, lhs: Montgomery<u64>, rhs: Montgomery<u64>) -> Montgomery<u64>{
155 Montgomery(self.mul_external_lazy(lhs, *rhs.value()))
156 }
157
158 /// Assigns lhs * rhs * (2^{64})^-1 mod q in range [0, 2q-1] to rhs.
159 #[inline(always)]
160 pub fn mul_internal_lazy_assign(&self, lhs: Montgomery<u64>, rhs: &mut Montgomery<u64>){
161 self.mul_external_lazy_assign(lhs, rhs.value_mut());
162 }
163
164 #[inline(always)]
165 pub fn add_internal(&self, lhs: Montgomery<u64>, rhs: Montgomery<u64>) -> Montgomery<u64>{
166 Montgomery(self.barrett.reduce(rhs.value() + lhs.value()))
167 }
168
169 /// Assigns lhs + rhs to rhs.
170 #[inline(always)]
171 pub fn add_internal_lazy_assign(&self, lhs: Montgomery<u64>, rhs: &mut Montgomery<u64>){
172 *rhs.value_mut() += lhs.value()
173 }
174
175 /// Assigns lhs + rhs - q if (lhs + rhs) >= q to rhs.
176 #[inline(always)]
177 pub fn add_internal_reduce_once_assign(&self, lhs: Montgomery<u64>, rhs: &mut Montgomery<u64>){
178 self.add_internal_lazy_assign(lhs, rhs);
179 rhs.value_mut().reduce_once_assign(self.q);
180 }
181
182 #[inline(always)]
183 pub fn reduce(&self, lhs: u64) -> u64{
184 self.barrett.reduce(lhs)
185 }
186
187 /// Returns lhs mod q in range [0, 2q-1].
188 #[inline(always)]
189 pub fn reduce_lazy(&self, lhs: u64) -> u64{
190 self.barrett.reduce_lazy(lhs)
191 }
192
193 #[inline(always)]
194 pub fn reduce_assign(&self, lhs: &mut u64){
195 self.barrett.reduce_assign(lhs)
196 }
197
198 /// Returns lhs mod q in range [0, 2q-1].
199 #[inline(always)]
200 pub fn reduce_lazy_assign(&self, lhs: &mut u64){
201 self.barrett.reduce_lazy_assign(lhs)
202 }
203
204 /// Returns (x^exponent) * 2^64 mod q.
205 #[inline(always)]
206 pub fn pow(&self, x: Montgomery<u64>, exponent:u64) -> Montgomery<u64>{
207 let mut y: Montgomery<u64> = self.one();
208 let mut x_mut: Montgomery<u64> = x;
209 let mut i: u64 = exponent;
210 while i > 0{
211 if i & 1 == 1{
212 self.mul_internal_assign(x_mut, &mut y);
213 }
214 self.mul_internal_assign(x_mut, &mut x_mut);
215 i >>= 1;
216 }
217
218 y.value_mut().reduce_once_assign(self.q);
219 y
220 }
221}
222
223/// Returns x^exponent mod q.
224/// This function internally instantiate a new MontgomeryPrecomp<u64>
225/// To be used when called only a few times and if there
226/// is no Prime instantiated with q.
227fn pow(x:u64, exponent:u64, q:u64) -> u64{
228 let montgomery: MontgomeryPrecomp<u64> = MontgomeryPrecomp::<u64>::new(q);
229 let mut y_mont: Montgomery<u64> = montgomery.one();
230 let mut x_mont: Montgomery<u64> = montgomery.prepare(x);
231 while exponent > 0{
232 if exponent & 1 == 1{
233 montgomery.mul_internal_assign(x_mont, &mut y_mont);
234 }
235
236 montgomery.mul_internal_assign(x_mont, &mut x_mont);
237 }
238
239 montgomery.unprepare(y_mont)
240}
241
242#[cfg(test)]
243mod tests {
244 use crate::modulus::montgomery;
245 use super::*;
246 use test::Bencher;
247
248 #[test]
249 fn test_mul_external() {
250 let q: u64 = 0x1fffffffffe00001;
251 let m_precomp = montgomery::MontgomeryPrecomp::new(q);
252 let x: u64 = 0x5f876e514845cc8b;
253 let y: u64 = 0xad726f98f24a761a;
254 let y_mont = m_precomp.prepare(y);
255 assert!(m_precomp.mul_external(y_mont, x) == (x as u128 * y as u128 % q as u128) as u64);
256 }
257
258 #[bench]
259 fn bench_mul_external(b: &mut Bencher){
260 let q: u64 = 0x1fffffffffe00001;
261 let m_precomp = montgomery::MontgomeryPrecomp::new(q);
262 let mut x: u64 = 0x5f876e514845cc8b;
263 let y: u64 = 0xad726f98f24a761a;
264 let y_mont = m_precomp.prepare(y);
265 b.iter(|| m_precomp.mul_external_assign(y_mont, &mut x));
266 }
267}