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}