miden-crypto/src/dsa/rpo_falcon512/math/polynomial.rs
src/dsa/rpo_falcon512/math/polynomial.rsRaw
1use super::{field::FalconFelt, Inverse};
2use crate::dsa::rpo_falcon512::{MODULUS, N};
3use crate::Felt;
4use alloc::vec::Vec;
5use core::default::Default;
6use core::fmt::Debug;
7use core::ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub, SubAssign};
8use num::{One, Zero};
9
10#[derive(Debug, Clone, Default)]
11pub struct Polynomial<F> {
12 pub coefficients: Vec<F>,
13}
14
15impl<F> Polynomial<F>
16where
17 F: Clone,
18{
19 pub fn new(coefficients: Vec<F>) -> Self {
20 Self { coefficients }
21 }
22}
23
24impl<
25 F: Mul<Output = F> + Sub<Output = F> + AddAssign + Zero + Div<Output = F> + Clone + Inverse,
26 > Polynomial<F>
27{
28 pub fn hadamard_mul(&self, other: &Self) -> Self {
29 Polynomial::new(
30 self.coefficients
31 .iter()
32 .zip(other.coefficients.iter())
33 .map(|(a, b)| *a * *b)
34 .collect(),
35 )
36 }
37 pub fn hadamard_div(&self, other: &Self) -> Self {
38 let other_coefficients_inverse = F::batch_inverse_or_zero(&other.coefficients);
39 Polynomial::new(
40 self.coefficients
41 .iter()
42 .zip(other_coefficients_inverse.iter())
43 .map(|(a, b)| *a * *b)
44 .collect(),
45 )
46 }
47
48 pub fn hadamard_inv(&self) -> Self {
49 let coefficients_inverse = F::batch_inverse_or_zero(&self.coefficients);
50 Polynomial::new(coefficients_inverse)
51 }
52}
53
54impl<F: Zero + PartialEq + Clone> Polynomial<F> {
55 pub fn degree(&self) -> Option<usize> {
56 if self.coefficients.is_empty() {
57 return None;
58 }
59 let mut max_index = self.coefficients.len() - 1;
60 while self.coefficients[max_index] == F::zero() {
61 if let Some(new_index) = max_index.checked_sub(1) {
62 max_index = new_index;
63 } else {
64 return None;
65 }
66 }
67 Some(max_index)
68 }
69
70 pub fn lc(&self) -> F {
71 match self.degree() {
72 Some(non_negative_degree) => self.coefficients[non_negative_degree].clone(),
73 None => F::zero(),
74 }
75 }
76}
77
78/// The following implementations are specific to cyclotomic polynomial rings,
79/// i.e., F\[ X \] / <X^n + 1>, and are used extensively in Falcon.
80impl<
81 F: One
82 + Zero
83 + Clone
84 + Neg<Output = F>
85 + MulAssign
86 + AddAssign
87 + Div<Output = F>
88 + Sub<Output = F>
89 + PartialEq,
90 > Polynomial<F>
91{
92 /// Reduce the polynomial by X^n + 1.
93 pub fn reduce_by_cyclotomic(&self, n: usize) -> Self {
94 let mut coefficients = vec![F::zero(); n];
95 let mut sign = -F::one();
96 for (i, c) in self.coefficients.iter().cloned().enumerate() {
97 if i % n == 0 {
98 sign *= -F::one();
99 }
100 coefficients[i % n] += sign.clone() * c;
101 }
102 Polynomial::new(coefficients)
103 }
104
105 /// Computes the field norm of the polynomial as an element of the cyclotomic ring
106 /// F\[ X \] / <X^n + 1 > relative to one of half the size, i.e., F\[ X \] / <X^(n/2) + 1> .
107 ///
108 /// Corresponds to formula 3.25 in the spec [1, p.30].
109 ///
110 /// [1]: https://falcon-sign.info/falcon.pdf
111 pub fn field_norm(&self) -> Self {
112 let n = self.coefficients.len();
113 let mut f0_coefficients = vec![F::zero(); n / 2];
114 let mut f1_coefficients = vec![F::zero(); n / 2];
115 for i in 0..n / 2 {
116 f0_coefficients[i] = self.coefficients[2 * i].clone();
117 f1_coefficients[i] = self.coefficients[2 * i + 1].clone();
118 }
119 let f0 = Polynomial::new(f0_coefficients);
120 let f1 = Polynomial::new(f1_coefficients);
121 let f0_squared = (f0.clone() * f0).reduce_by_cyclotomic(n / 2);
122 let f1_squared = (f1.clone() * f1).reduce_by_cyclotomic(n / 2);
123 let x = Polynomial::new(vec![F::zero(), F::one()]);
124 f0_squared - (x * f1_squared).reduce_by_cyclotomic(n / 2)
125 }
126
127 /// Lifts an element from a cyclotomic polynomial ring to one of double the size.
128 pub fn lift_next_cyclotomic(&self) -> Self {
129 let n = self.coefficients.len();
130 let mut coefficients = vec![F::zero(); n * 2];
131 for i in 0..n {
132 coefficients[2 * i] = self.coefficients[i].clone();
133 }
134 Self::new(coefficients)
135 }
136
137 /// Computes the galois adjoint of the polynomial in the cyclotomic ring F\[ X \] / < X^n + 1 > ,
138 /// which corresponds to f(x^2).
139 pub fn galois_adjoint(&self) -> Self {
140 Self::new(
141 self.coefficients
142 .iter()
143 .enumerate()
144 .map(|(i, c)| if i % 2 == 0 { c.clone() } else { c.clone().neg() })
145 .collect(),
146 )
147 }
148}
149
150impl<F: Clone + Into<f64>> Polynomial<F> {
151 pub(crate) fn l2_norm_squared(&self) -> f64 {
152 self.coefficients
153 .iter()
154 .map(|i| Into::<f64>::into(i.clone()))
155 .map(|i| i * i)
156 .sum::<f64>()
157 }
158}
159
160impl<F> PartialEq for Polynomial<F>
161where
162 F: Zero + PartialEq + Clone + AddAssign,
163{
164 fn eq(&self, other: &Self) -> bool {
165 if self.is_zero() && other.is_zero() {
166 true
167 } else if self.is_zero() || other.is_zero() {
168 false
169 } else {
170 let self_degree = self.degree().unwrap();
171 let other_degree = other.degree().unwrap();
172 self.coefficients[0..=self_degree] == other.coefficients[0..=other_degree]
173 }
174 }
175}
176
177impl<F> Eq for Polynomial<F> where F: Zero + PartialEq + Clone + AddAssign {}
178
179impl<F> Add for &Polynomial<F>
180where
181 F: Add<Output = F> + AddAssign + Clone,
182{
183 type Output = Polynomial<F>;
184
185 fn add(self, rhs: Self) -> Self::Output {
186 let coefficients = if self.coefficients.len() >= rhs.coefficients.len() {
187 let mut coefficients = self.coefficients.clone();
188 for (i, c) in rhs.coefficients.iter().enumerate() {
189 coefficients[i] += c.clone();
190 }
191 coefficients
192 } else {
193 let mut coefficients = rhs.coefficients.clone();
194 for (i, c) in self.coefficients.iter().enumerate() {
195 coefficients[i] += c.clone();
196 }
197 coefficients
198 };
199 Self::Output { coefficients }
200 }
201}
202
203impl<F> Add for Polynomial<F>
204where
205 F: Add<Output = F> + AddAssign + Clone,
206{
207 type Output = Polynomial<F>;
208 fn add(self, rhs: Self) -> Self::Output {
209 let coefficients = if self.coefficients.len() >= rhs.coefficients.len() {
210 let mut coefficients = self.coefficients.clone();
211 for (i, c) in rhs.coefficients.into_iter().enumerate() {
212 coefficients[i] += c;
213 }
214 coefficients
215 } else {
216 let mut coefficients = rhs.coefficients.clone();
217 for (i, c) in self.coefficients.into_iter().enumerate() {
218 coefficients[i] += c;
219 }
220 coefficients
221 };
222 Self::Output { coefficients }
223 }
224}
225
226impl<F> AddAssign for Polynomial<F>
227where
228 F: Add<Output = F> + AddAssign + Clone,
229{
230 fn add_assign(&mut self, rhs: Self) {
231 if self.coefficients.len() >= rhs.coefficients.len() {
232 for (i, c) in rhs.coefficients.into_iter().enumerate() {
233 self.coefficients[i] += c;
234 }
235 } else {
236 let mut coefficients = rhs.coefficients.clone();
237 for (i, c) in self.coefficients.iter().enumerate() {
238 coefficients[i] += c.clone();
239 }
240 self.coefficients = coefficients;
241 }
242 }
243}
244
245impl<F> Sub for &Polynomial<F>
246where
247 F: Sub<Output = F> + Clone + Neg<Output = F> + Add<Output = F> + AddAssign,
248{
249 type Output = Polynomial<F>;
250
251 fn sub(self, rhs: Self) -> Self::Output {
252 self + &(-rhs)
253 }
254}
255
256impl<F> Sub for Polynomial<F>
257where
258 F: Sub<Output = F> + Clone + Neg<Output = F> + Add<Output = F> + AddAssign,
259{
260 type Output = Polynomial<F>;
261
262 fn sub(self, rhs: Self) -> Self::Output {
263 self + (-rhs)
264 }
265}
266
267impl<F> SubAssign for Polynomial<F>
268where
269 F: Add<Output = F> + Neg<Output = F> + AddAssign + Clone + Sub<Output = F>,
270{
271 fn sub_assign(&mut self, rhs: Self) {
272 self.coefficients = self.clone().sub(rhs).coefficients;
273 }
274}
275
276impl<F: Neg<Output = F> + Clone> Neg for &Polynomial<F> {
277 type Output = Polynomial<F>;
278
279 fn neg(self) -> Self::Output {
280 Self::Output {
281 coefficients: self.coefficients.iter().cloned().map(|a| -a).collect(),
282 }
283 }
284}
285
286impl<F: Neg<Output = F> + Clone> Neg for Polynomial<F> {
287 type Output = Self;
288
289 fn neg(self) -> Self::Output {
290 Self::Output {
291 coefficients: self.coefficients.iter().cloned().map(|a| -a).collect(),
292 }
293 }
294}
295
296impl<F> Mul for &Polynomial<F>
297where
298 F: Add + AddAssign + Mul<Output = F> + Sub<Output = F> + Zero + PartialEq + Clone,
299{
300 type Output = Polynomial<F>;
301
302 fn mul(self, other: Self) -> Self::Output {
303 if self.is_zero() || other.is_zero() {
304 return Polynomial::<F>::zero();
305 }
306 let mut coefficients =
307 vec![F::zero(); self.coefficients.len() + other.coefficients.len() - 1];
308 for i in 0..self.coefficients.len() {
309 for j in 0..other.coefficients.len() {
310 coefficients[i + j] += self.coefficients[i].clone() * other.coefficients[j].clone();
311 }
312 }
313 Polynomial { coefficients }
314 }
315}
316
317impl<F> Mul for Polynomial<F>
318where
319 F: Add + AddAssign + Mul<Output = F> + Zero + PartialEq + Clone,
320{
321 type Output = Self;
322
323 fn mul(self, other: Self) -> Self::Output {
324 if self.is_zero() || other.is_zero() {
325 return Self::zero();
326 }
327 let mut coefficients =
328 vec![F::zero(); self.coefficients.len() + other.coefficients.len() - 1];
329 for i in 0..self.coefficients.len() {
330 for j in 0..other.coefficients.len() {
331 coefficients[i + j] += self.coefficients[i].clone() * other.coefficients[j].clone();
332 }
333 }
334 Self { coefficients }
335 }
336}
337
338impl<F: Add + Mul<Output = F> + Zero + Clone> Mul<F> for &Polynomial<F> {
339 type Output = Polynomial<F>;
340
341 fn mul(self, other: F) -> Self::Output {
342 Polynomial {
343 coefficients: self.coefficients.iter().cloned().map(|i| i * other.clone()).collect(),
344 }
345 }
346}
347
348impl<F: Add + Mul<Output = F> + Zero + Clone> Mul<F> for Polynomial<F> {
349 type Output = Polynomial<F>;
350
351 fn mul(self, other: F) -> Self::Output {
352 Polynomial {
353 coefficients: self.coefficients.iter().cloned().map(|i| i * other.clone()).collect(),
354 }
355 }
356}
357
358impl<F: Mul<Output = F> + Sub<Output = F> + AddAssign + Zero + Div<Output = F> + Clone>
359 Polynomial<F>
360{
361 /// Multiply two polynomials using Karatsuba's divide-and-conquer algorithm.
362 pub fn karatsuba(&self, other: &Self) -> Self {
363 Polynomial::new(vector_karatsuba(&self.coefficients, &other.coefficients))
364 }
365}
366
367impl<F> One for Polynomial<F>
368where
369 F: Clone + One + PartialEq + Zero + AddAssign,
370{
371 fn one() -> Self {
372 Self { coefficients: vec![F::one()] }
373 }
374}
375
376impl<F> Zero for Polynomial<F>
377where
378 F: Zero + PartialEq + Clone + AddAssign,
379{
380 fn zero() -> Self {
381 Self { coefficients: vec![] }
382 }
383
384 fn is_zero(&self) -> bool {
385 self.degree().is_none()
386 }
387}
388
389impl<F: Zero + Clone> Polynomial<F> {
390 pub fn shift(&self, shamt: usize) -> Self {
391 Self {
392 coefficients: [vec![F::zero(); shamt], self.coefficients.clone()].concat(),
393 }
394 }
395
396 pub fn constant(f: F) -> Self {
397 Self { coefficients: vec![f] }
398 }
399
400 pub fn map<G: Clone, C: FnMut(&F) -> G>(&self, closure: C) -> Polynomial<G> {
401 Polynomial::<G>::new(self.coefficients.iter().map(closure).collect())
402 }
403
404 pub fn fold<G, C: FnMut(G, &F) -> G + Clone>(&self, mut initial_value: G, closure: C) -> G {
405 for c in self.coefficients.iter() {
406 initial_value = (closure.clone())(initial_value, c);
407 }
408 initial_value
409 }
410}
411
412impl<F> Div<Polynomial<F>> for Polynomial<F>
413where
414 F: Zero
415 + One
416 + PartialEq
417 + AddAssign
418 + Clone
419 + Mul<Output = F>
420 + MulAssign
421 + Div<Output = F>
422 + Neg<Output = F>
423 + Sub<Output = F>,
424{
425 type Output = Polynomial<F>;
426
427 fn div(self, denominator: Self) -> Self::Output {
428 if denominator.is_zero() {
429 panic!();
430 }
431 if self.is_zero() {
432 Self::zero();
433 }
434 let mut remainder = self.clone();
435 let mut quotient = Polynomial::<F>::zero();
436 while remainder.degree().unwrap() >= denominator.degree().unwrap() {
437 let shift = remainder.degree().unwrap() - denominator.degree().unwrap();
438 let quotient_coefficient = remainder.lc() / denominator.lc();
439 let monomial = Self::constant(quotient_coefficient).shift(shift);
440 quotient += monomial.clone();
441 remainder -= monomial * denominator.clone();
442 if remainder.is_zero() {
443 break;
444 }
445 }
446 quotient
447 }
448}
449
450fn vector_karatsuba<
451 F: Zero + AddAssign + Mul<Output = F> + Sub<Output = F> + Div<Output = F> + Clone,
452>(
453 left: &[F],
454 right: &[F],
455) -> Vec<F> {
456 let n = left.len();
457 if n <= 8 {
458 let mut product = vec![F::zero(); left.len() + right.len() - 1];
459 for (i, l) in left.iter().enumerate() {
460 for (j, r) in right.iter().enumerate() {
461 product[i + j] += l.clone() * r.clone();
462 }
463 }
464 return product;
465 }
466 let n_over_2 = n / 2;
467 let mut product = vec![F::zero(); 2 * n - 1];
468 let left_lo = &left[0..n_over_2];
469 let right_lo = &right[0..n_over_2];
470 let left_hi = &left[n_over_2..];
471 let right_hi = &right[n_over_2..];
472 let left_sum: Vec<F> =
473 left_lo.iter().zip(left_hi).map(|(a, b)| a.clone() + b.clone()).collect();
474 let right_sum: Vec<F> =
475 right_lo.iter().zip(right_hi).map(|(a, b)| a.clone() + b.clone()).collect();
476
477 let prod_lo = vector_karatsuba(left_lo, right_lo);
478 let prod_hi = vector_karatsuba(left_hi, right_hi);
479 let prod_mid: Vec<F> = vector_karatsuba(&left_sum, &right_sum)
480 .iter()
481 .zip(prod_lo.iter().zip(prod_hi.iter()))
482 .map(|(s, (l, h))| s.clone() - (l.clone() + h.clone()))
483 .collect();
484
485 for (i, l) in prod_lo.into_iter().enumerate() {
486 product[i] = l;
487 }
488 for (i, m) in prod_mid.into_iter().enumerate() {
489 product[i + n_over_2] += m;
490 }
491 for (i, h) in prod_hi.into_iter().enumerate() {
492 product[i + n] += h
493 }
494 product
495}
496
497impl From<Polynomial<FalconFelt>> for Polynomial<Felt> {
498 fn from(item: Polynomial<FalconFelt>) -> Self {
499 let res: Vec<Felt> =
500 item.coefficients.iter().map(|a| Felt::from(a.value() as u16)).collect();
501 Polynomial::new(res)
502 }
503}
504
505impl From<&Polynomial<FalconFelt>> for Polynomial<Felt> {
506 fn from(item: &Polynomial<FalconFelt>) -> Self {
507 let res: Vec<Felt> =
508 item.coefficients.iter().map(|a| Felt::from(a.value() as u16)).collect();
509 Polynomial::new(res)
510 }
511}
512
513impl From<Polynomial<i16>> for Polynomial<FalconFelt> {
514 fn from(item: Polynomial<i16>) -> Self {
515 let res: Vec<FalconFelt> = item.coefficients.iter().map(|&a| FalconFelt::new(a)).collect();
516 Polynomial::new(res)
517 }
518}
519
520impl From<&Polynomial<i16>> for Polynomial<FalconFelt> {
521 fn from(item: &Polynomial<i16>) -> Self {
522 let res: Vec<FalconFelt> = item.coefficients.iter().map(|&a| FalconFelt::new(a)).collect();
523 Polynomial::new(res)
524 }
525}
526
527impl From<Vec<i16>> for Polynomial<FalconFelt> {
528 fn from(item: Vec<i16>) -> Self {
529 let res: Vec<FalconFelt> = item.iter().map(|&a| FalconFelt::new(a)).collect();
530 Polynomial::new(res)
531 }
532}
533
534impl From<&Vec<i16>> for Polynomial<FalconFelt> {
535 fn from(item: &Vec<i16>) -> Self {
536 let res: Vec<FalconFelt> = item.iter().map(|&a| FalconFelt::new(a)).collect();
537 Polynomial::new(res)
538 }
539}
540
541impl Polynomial<FalconFelt> {
542 pub fn norm_squared(&self) -> u64 {
543 self.coefficients
544 .iter()
545 .map(|&i| i.balanced_value() as i64)
546 .map(|i| (i * i) as u64)
547 .sum::<u64>()
548 }
549
550 // PUBLIC ACCESSORS
551 // --------------------------------------------------------------------------------------------
552
553 /// Returns the coefficients of this polynomial as field elements.
554 pub fn to_elements(&self) -> Vec<Felt> {
555 self.coefficients.iter().map(|&a| Felt::from(a.value() as u16)).collect()
556 }
557
558 // POLYNOMIAL OPERATIONS
559 // --------------------------------------------------------------------------------------------
560
561 /// Multiplies two polynomials over Z_p\[x\] without reducing modulo p. Given that the degrees
562 /// of the input polynomials are less than 512 and their coefficients are less than the modulus
563 /// q equal to 12289, the resulting product polynomial is guaranteed to have coefficients less
564 /// than the Miden prime.
565 ///
566 /// Note that this multiplication is not over Z_p\[x\]/(phi).
567 pub fn mul_modulo_p(a: &Self, b: &Self) -> [u64; 1024] {
568 let mut c = [0; 2 * N];
569 for i in 0..N {
570 for j in 0..N {
571 c[i + j] += a.coefficients[i].value() as u64 * b.coefficients[j].value() as u64;
572 }
573 }
574
575 c
576 }
577
578 /// Reduces a polynomial, that is the product of two polynomials over Z_p\[x\], modulo
579 /// the irreducible polynomial phi. This results in an element in Z_p\[x\]/(phi).
580 pub fn reduce_negacyclic(a: &[u64; 1024]) -> Self {
581 let mut c = [FalconFelt::zero(); N];
582 let modulus = MODULUS as u16;
583 for i in 0..N {
584 let ai = a[N + i] % modulus as u64;
585 let neg_ai = (modulus - ai as u16) % modulus;
586
587 let bi = (a[i] % modulus as u64) as u16;
588 c[i] = FalconFelt::new(((neg_ai + bi) % modulus) as i16);
589 }
590
591 Self::new(c.to_vec())
592 }
593}
594
595// TESTS
596// ================================================================================================
597
598#[cfg(test)]
599mod tests {
600 use super::{FalconFelt, Polynomial, N};
601
602 #[test]
603 fn test_negacyclic_reduction() {
604 let coef1: [u8; N] = rand_utils::rand_array();
605 let coef2: [u8; N] = rand_utils::rand_array();
606
607 let poly1 = Polynomial::new(coef1.iter().map(|&a| FalconFelt::new(a as i16)).collect());
608 let poly2 = Polynomial::new(coef2.iter().map(|&a| FalconFelt::new(a as i16)).collect());
609 let prod = poly1.clone() * poly2.clone();
610
611 assert_eq!(
612 prod.reduce_by_cyclotomic(N),
613 Polynomial::reduce_negacyclic(&Polynomial::mul_modulo_p(&poly1, &poly2))
614 );
615 }
616}