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