miden-crypto/src/dsa/rpo_falcon512/math/mod.rs
src/dsa/rpo_falcon512/math/mod.rsRaw
1//! Contains different structs and methods related to the Falcon DSA.
2//!
3//! It uses and acknowledges the work in:
4//!
5//! 1. The [reference](https://falcon-sign.info/impl/README.txt.html) implementation by Thomas
6//! Pornin.
7//! 2. The [Rust](https://github.com/aszepieniec/falcon-rust) implementation by Alan Szepieniec.
8use alloc::{string::String, vec::Vec};
9use core::ops::MulAssign;
10
11#[cfg(not(feature = "std"))]
12use num::Float;
13use num::{BigInt, FromPrimitive, One, Zero};
14use num_complex::Complex64;
15use rand::Rng;
16
17use super::MODULUS;
18
19mod fft;
20pub use fft::{CyclotomicFourier, FastFft};
21
22mod field;
23pub use field::FalconFelt;
24
25mod ffsampling;
26pub use ffsampling::{ffldl, ffsampling, gram, normalize_tree, LdlTree};
27
28mod samplerz;
29use self::samplerz::sampler_z;
30
31mod polynomial;
32pub use polynomial::Polynomial;
33
34pub trait Inverse: Copy + Zero + MulAssign + One {
35 /// Gets the inverse of a, or zero if it is zero.
36 fn inverse_or_zero(self) -> Self;
37
38 /// Gets the inverses of a batch of elements, and skip over any that are zero.
39 fn batch_inverse_or_zero(batch: &[Self]) -> Vec<Self> {
40 let mut acc = Self::one();
41 let mut rp: Vec<Self> = Vec::with_capacity(batch.len());
42 for batch_item in batch {
43 if !batch_item.is_zero() {
44 rp.push(acc);
45 acc = *batch_item * acc;
46 } else {
47 rp.push(Self::zero());
48 }
49 }
50 let mut inv = Self::inverse_or_zero(acc);
51 for i in (0..batch.len()).rev() {
52 if !batch[i].is_zero() {
53 rp[i] *= inv;
54 inv *= batch[i];
55 }
56 }
57 rp
58 }
59}
60
61impl Inverse for Complex64 {
62 fn inverse_or_zero(self) -> Self {
63 let modulus = self.re * self.re + self.im * self.im;
64 Complex64::new(self.re / modulus, -self.im / modulus)
65 }
66 fn batch_inverse_or_zero(batch: &[Self]) -> Vec<Self> {
67 batch.iter().map(|&c| Complex64::new(1.0, 0.0) / c).collect()
68 }
69}
70
71impl Inverse for f64 {
72 fn inverse_or_zero(self) -> Self {
73 1.0 / self
74 }
75 fn batch_inverse_or_zero(batch: &[Self]) -> Vec<Self> {
76 batch.iter().map(|&c| 1.0 / c).collect()
77 }
78}
79
80/// Samples 4 small polynomials f, g, F, G such that f * G - g * F = q mod (X^n + 1).
81/// Algorithm 5 (NTRUgen) of the documentation [1, p.34].
82///
83/// [1]: https://falcon-sign.info/falcon.pdf
84pub(crate) fn ntru_gen<R: Rng>(n: usize, rng: &mut R) -> [Polynomial<i16>; 4] {
85 loop {
86 let f = gen_poly(n, rng);
87 let g = gen_poly(n, rng);
88 let f_ntt = f.map(|&i| FalconFelt::new(i)).fft();
89 if f_ntt.coefficients.iter().any(|e| e.is_zero()) {
90 continue;
91 }
92 let gamma = gram_schmidt_norm_squared(&f, &g);
93 if gamma > 1.3689f64 * (MODULUS as f64) {
94 continue;
95 }
96
97 if let Some((capital_f, capital_g)) =
98 ntru_solve(&f.map(|&i| i.into()), &g.map(|&i| i.into()))
99 {
100 return [
101 g,
102 -f,
103 capital_g.map(|i| i.try_into().unwrap()),
104 -capital_f.map(|i| i.try_into().unwrap()),
105 ];
106 }
107 }
108}
109
110/// Solves the NTRU equation. Given f, g in ZZ[X], find F, G in ZZ[X] such that:
111///
112/// f G - g F = q mod (X^n + 1)
113///
114/// Algorithm 6 of the specification [1, p.35].
115///
116/// [1]: https://falcon-sign.info/falcon.pdf
117fn ntru_solve(
118 f: &Polynomial<BigInt>,
119 g: &Polynomial<BigInt>,
120) -> Option<(Polynomial<BigInt>, Polynomial<BigInt>)> {
121 let n = f.coefficients.len();
122 if n == 1 {
123 let (gcd, u, v) = xgcd(&f.coefficients[0], &g.coefficients[0]);
124 if gcd != BigInt::one() {
125 return None;
126 }
127 return Some((
128 (Polynomial::new(vec![-v * BigInt::from_u32(MODULUS as u32).unwrap()])),
129 Polynomial::new(vec![u * BigInt::from_u32(MODULUS as u32).unwrap()]),
130 ));
131 }
132
133 let f_prime = f.field_norm();
134 let g_prime = g.field_norm();
135
136 let (capital_f_prime, capital_g_prime) = ntru_solve(&f_prime, &g_prime)?;
137 let capital_f_prime_xsq = capital_f_prime.lift_next_cyclotomic();
138 let capital_g_prime_xsq = capital_g_prime.lift_next_cyclotomic();
139
140 let f_minx = f.galois_adjoint();
141 let g_minx = g.galois_adjoint();
142
143 let mut capital_f = (capital_f_prime_xsq.karatsuba(&g_minx)).reduce_by_cyclotomic(n);
144 let mut capital_g = (capital_g_prime_xsq.karatsuba(&f_minx)).reduce_by_cyclotomic(n);
145
146 match babai_reduce(f, g, &mut capital_f, &mut capital_g) {
147 Ok(_) => Some((capital_f, capital_g)),
148 Err(_e) => {
149 #[cfg(test)]
150 {
151 panic!("{}", _e);
152 }
153 #[cfg(not(test))]
154 {
155 None
156 }
157 },
158 }
159}
160
161/// Generates a polynomial of degree at most n-1 whose coefficients are distributed according
162/// to a discrete Gaussian with mu = 0 and sigma = 1.17 * sqrt(Q / (2n)).
163fn gen_poly<R: Rng>(n: usize, rng: &mut R) -> Polynomial<i16> {
164 let mu = 0.0;
165 let sigma_star = 1.43300980528773;
166 Polynomial {
167 coefficients: (0..4096)
168 .map(|_| sampler_z(mu, sigma_star, sigma_star - 0.001, rng))
169 .collect::<Vec<i16>>()
170 .chunks(4096 / n)
171 .map(|ch| ch.iter().sum())
172 .collect(),
173 }
174}
175
176/// Computes the Gram-Schmidt norm of B = [[g, -f], [G, -F]] from f and g.
177/// Corresponds to line 9 in algorithm 5 of the spec [1, p.34]
178///
179/// [1]: https://falcon-sign.info/falcon.pdf
180fn gram_schmidt_norm_squared(f: &Polynomial<i16>, g: &Polynomial<i16>) -> f64 {
181 let n = f.coefficients.len();
182 let norm_f_squared = f.l2_norm_squared();
183 let norm_g_squared = g.l2_norm_squared();
184 let gamma1 = norm_f_squared + norm_g_squared;
185
186 let f_fft = f.map(|i| Complex64::new(*i as f64, 0.0)).fft();
187 let g_fft = g.map(|i| Complex64::new(*i as f64, 0.0)).fft();
188 let f_adj_fft = f_fft.map(|c| c.conj());
189 let g_adj_fft = g_fft.map(|c| c.conj());
190 let ffgg_fft = f_fft.hadamard_mul(&f_adj_fft) + g_fft.hadamard_mul(&g_adj_fft);
191 let ffgg_fft_inverse = ffgg_fft.hadamard_inv();
192 let qf_over_ffgg_fft = f_adj_fft.map(|c| c * (MODULUS as f64)).hadamard_mul(&ffgg_fft_inverse);
193 let qg_over_ffgg_fft = g_adj_fft.map(|c| c * (MODULUS as f64)).hadamard_mul(&ffgg_fft_inverse);
194 let norm_f_over_ffgg_squared =
195 qf_over_ffgg_fft.coefficients.iter().map(|c| (c * c.conj()).re).sum::<f64>() / (n as f64);
196 let norm_g_over_ffgg_squared =
197 qg_over_ffgg_fft.coefficients.iter().map(|c| (c * c.conj()).re).sum::<f64>() / (n as f64);
198
199 let gamma2 = norm_f_over_ffgg_squared + norm_g_over_ffgg_squared;
200
201 f64::max(gamma1, gamma2)
202}
203
204/// Reduces the vector (F,G) relative to (f,g). This method follows the python implementation [1].
205/// Note that this algorithm can end up in an infinite loop. (It's one of the things the author
206/// would like to fix.) When this happens, control returns an error (hence the return type) and
207/// generates another keypair with fresh randomness.
208///
209/// Algorithm 7 in the spec [2, p.35]
210///
211/// [1]: https://github.com/tprest/falcon.py
212///
213/// [2]: https://falcon-sign.info/falcon.pdf
214fn babai_reduce(
215 f: &Polynomial<BigInt>,
216 g: &Polynomial<BigInt>,
217 capital_f: &mut Polynomial<BigInt>,
218 capital_g: &mut Polynomial<BigInt>,
219) -> Result<(), String> {
220 let bitsize = |bi: &BigInt| (bi.bits() + 7) & (u64::MAX ^ 7);
221 let n = f.coefficients.len();
222 let size = [
223 f.map(bitsize).fold(0, |a, &b| u64::max(a, b)),
224 g.map(bitsize).fold(0, |a, &b| u64::max(a, b)),
225 53,
226 ]
227 .into_iter()
228 .max()
229 .unwrap();
230 let shift = (size as i64) - 53;
231 let f_adjusted = f
232 .map(|bi| Complex64::new(i64::try_from(bi >> shift).unwrap() as f64, 0.0))
233 .fft();
234 let g_adjusted = g
235 .map(|bi| Complex64::new(i64::try_from(bi >> shift).unwrap() as f64, 0.0))
236 .fft();
237
238 let f_star_adjusted = f_adjusted.map(|c| c.conj());
239 let g_star_adjusted = g_adjusted.map(|c| c.conj());
240 let denominator_fft =
241 f_adjusted.hadamard_mul(&f_star_adjusted) + g_adjusted.hadamard_mul(&g_star_adjusted);
242
243 let mut counter = 0;
244 loop {
245 let capital_size = [
246 capital_f.map(bitsize).fold(0, |a, &b| u64::max(a, b)),
247 capital_g.map(bitsize).fold(0, |a, &b| u64::max(a, b)),
248 53,
249 ]
250 .into_iter()
251 .max()
252 .unwrap();
253
254 if capital_size < size {
255 break;
256 }
257 let capital_shift = (capital_size as i64) - 53;
258 let capital_f_adjusted = capital_f
259 .map(|bi| Complex64::new(i64::try_from(bi >> capital_shift).unwrap() as f64, 0.0))
260 .fft();
261 let capital_g_adjusted = capital_g
262 .map(|bi| Complex64::new(i64::try_from(bi >> capital_shift).unwrap() as f64, 0.0))
263 .fft();
264
265 let numerator = capital_f_adjusted.hadamard_mul(&f_star_adjusted)
266 + capital_g_adjusted.hadamard_mul(&g_star_adjusted);
267 let quotient = numerator.hadamard_div(&denominator_fft).ifft();
268
269 let k = quotient.map(|f| Into::<BigInt>::into(f.re.round() as i64));
270
271 if k.is_zero() {
272 break;
273 }
274 let kf = (k.clone().karatsuba(f))
275 .reduce_by_cyclotomic(n)
276 .map(|bi| bi << (capital_size - size));
277 let kg = (k.clone().karatsuba(g))
278 .reduce_by_cyclotomic(n)
279 .map(|bi| bi << (capital_size - size));
280 *capital_f -= kf;
281 *capital_g -= kg;
282
283 counter += 1;
284 if counter > 1000 {
285 // If we get here, that means that (with high likelihood) we are in an
286 // infinite loop. We know it happens from time to time -- seldomly, but it
287 // does. It would be nice to fix that! But in order to fix it we need to be
288 // able to reproduce it, and for that we need test vectors. So print them
289 // and hope that one day they circle back to the implementor.
290 return Err(format!("Encountered infinite loop in babai_reduce of falcon-rust.\n\\
291 Please help the developer(s) fix it! You can do this by sending them the inputs to the function that caused the behavior:\n\\
292 f: {:?}\n\\
293 g: {:?}\n\\
294 capital_f: {:?}\n\\
295 capital_g: {:?}\n", f.coefficients, g.coefficients, capital_f.coefficients, capital_g.coefficients));
296 }
297 }
298 Ok(())
299}
300
301/// Extended Euclidean algorithm for computing the greatest common divisor (g) and
302/// Bézout coefficients (u, v) for the relation
303///
304/// $$ u a + v b = g . $$
305///
306/// Implementation adapted from Wikipedia [1].
307///
308/// [1]: https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Pseudocode
309fn xgcd(a: &BigInt, b: &BigInt) -> (BigInt, BigInt, BigInt) {
310 let (mut old_r, mut r) = (a.clone(), b.clone());
311 let (mut old_s, mut s) = (BigInt::one(), BigInt::zero());
312 let (mut old_t, mut t) = (BigInt::zero(), BigInt::one());
313
314 while r != BigInt::zero() {
315 let quotient = old_r.clone() / r.clone();
316 (old_r, r) = (r.clone(), old_r.clone() - quotient.clone() * r);
317 (old_s, s) = (s.clone(), old_s.clone() - quotient.clone() * s);
318 (old_t, t) = (t.clone(), old_t.clone() - quotient * t);
319 }
320
321 (old_r, old_s, old_t)
322}