miden_crypto/dsa/falcon512_poseidon2/math/
polynomial.rs

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

miden_crypto/dsa/falcon512_poseidon2/math/
polynomial.rs
