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

miden_crypto/dsa/falcon512_rpo/math/
polynomial.rs
