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};
11
12use super::{Inverse, field::FalconFelt};
13use crate::{
14    Felt,
15    dsa::falcon512_poseidon2::{MODULUS, N},
16    utils::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            max_index = max_index.checked_sub(1)?;
77        }
78        Some(max_index)
79    }
80
81    /// Returns the leading coefficient of the polynomial.
82    pub fn lc(&self) -> F {
83        match self.degree() {
84            Some(non_negative_degree) => self.coefficients[non_negative_degree].clone(),
85            None => F::zero(),
86        }
87    }
88}
89
90/// The following implementations are specific to cyclotomic polynomial rings,
91/// i.e., F\[ X \] / <X^n + 1>, and are used extensively in Falcon.
92impl<
93    F: One
94        + Zero
95        + Clone
96        + Neg<Output = F>
97        + MulAssign
98        + AddAssign
99        + Div<Output = F>
100        + Sub<Output = F>
101        + PartialEq,
102> Polynomial<F>
103{
104    /// Reduce the polynomial by X^n + 1.
105    pub fn reduce_by_cyclotomic(&self, n: usize) -> Self {
106        let mut coefficients = vec![F::zero(); n];
107        let mut sign = -F::one();
108        for (i, c) in self.coefficients.iter().cloned().enumerate() {
109            if i.is_multiple_of(n) {
110                sign *= -F::one();
111            }
112            coefficients[i % n] += sign.clone() * c;
113        }
114        Polynomial::new(coefficients)
115    }
116
117    /// Computes the field norm of the polynomial as an element of the cyclotomic ring
118    ///  F\[ X \] / <X^n + 1 > relative to one of half the size, i.e., F\[ X \] / <X^(n/2) + 1> .
119    ///
120    /// Corresponds to formula 3.25 in the spec [1, p.30].
121    ///
122    /// [1]: <https://falcon-sign.info/falcon.pdf>
123    pub fn field_norm(&self) -> Self {
124        let n = self.coefficients.len();
125        let mut f0_coefficients = vec![F::zero(); n / 2];
126        let mut f1_coefficients = vec![F::zero(); n / 2];
127        for i in 0..n / 2 {
128            f0_coefficients[i] = self.coefficients[2 * i].clone();
129            f1_coefficients[i] = self.coefficients[2 * i + 1].clone();
130        }
131        let f0 = Polynomial::new(f0_coefficients);
132        let f1 = Polynomial::new(f1_coefficients);
133        let f0_squared = (f0.clone() * f0).reduce_by_cyclotomic(n / 2);
134        let f1_squared = (f1.clone() * f1).reduce_by_cyclotomic(n / 2);
135        let x = Polynomial::new(vec![F::zero(), F::one()]);
136        f0_squared - (x * f1_squared).reduce_by_cyclotomic(n / 2)
137    }
138
139    /// Lifts an element from a cyclotomic polynomial ring to one of double the size.
140    pub fn lift_next_cyclotomic(&self) -> Self {
141        let n = self.coefficients.len();
142        let mut coefficients = vec![F::zero(); n * 2];
143        for i in 0..n {
144            coefficients[2 * i] = self.coefficients[i].clone();
145        }
146        Self::new(coefficients)
147    }
148
149    /// Computes the galois adjoint of the polynomial in the cyclotomic ring F\[ X \] / < X^n + 1 >
150    /// , which corresponds to f(x^2).
151    pub fn galois_adjoint(&self) -> Self {
152        Self::new(
153            self.coefficients
154                .iter()
155                .enumerate()
156                .map(|(i, c)| {
157                    if i.is_multiple_of(2) {
158                        c.clone()
159                    } else {
160                        c.clone().neg()
161                    }
162                })
163                .collect(),
164        )
165    }
166}
167
168impl<F: Clone + Into<f64>> Polynomial<F> {
169    pub(crate) fn l2_norm_squared(&self) -> f64 {
170        self.coefficients
171            .iter()
172            .map(|i| Into::<f64>::into(i.clone()))
173            .map(|i| i * i)
174            .sum::<f64>()
175    }
176}
177
178impl<F> PartialEq for Polynomial<F>
179where
180    F: Zero + PartialEq + Clone + AddAssign,
181{
182    fn eq(&self, other: &Self) -> bool {
183        if self.is_zero() && other.is_zero() {
184            true
185        } else if self.is_zero() || other.is_zero() {
186            false
187        } else {
188            let self_degree = self.degree().unwrap();
189            let other_degree = other.degree().unwrap();
190            self.coefficients[0..=self_degree] == other.coefficients[0..=other_degree]
191        }
192    }
193}
194
195impl<F> Eq for Polynomial<F> where F: Zero + PartialEq + Clone + AddAssign {}
196
197impl<F> Add for &Polynomial<F>
198where
199    F: Add<Output = F> + AddAssign + Clone,
200{
201    type Output = Polynomial<F>;
202
203    fn add(self, rhs: Self) -> Self::Output {
204        let coefficients = if self.coefficients.len() >= rhs.coefficients.len() {
205            let mut coefficients = self.coefficients.clone();
206            for (i, c) in rhs.coefficients.iter().enumerate() {
207                coefficients[i] += c.clone();
208            }
209            coefficients
210        } else {
211            let mut coefficients = rhs.coefficients.clone();
212            for (i, c) in self.coefficients.iter().enumerate() {
213                coefficients[i] += c.clone();
214            }
215            coefficients
216        };
217        Self::Output { coefficients }
218    }
219}
220
221impl<F> Add for Polynomial<F>
222where
223    F: Add<Output = F> + AddAssign + Clone,
224{
225    type Output = Polynomial<F>;
226    fn add(self, rhs: Self) -> Self::Output {
227        let coefficients = if self.coefficients.len() >= rhs.coefficients.len() {
228            let mut coefficients = self.coefficients;
229            for (i, c) in rhs.coefficients.into_iter().enumerate() {
230                coefficients[i] += c;
231            }
232            coefficients
233        } else {
234            let mut coefficients = rhs.coefficients;
235            for (i, c) in self.coefficients.into_iter().enumerate() {
236                coefficients[i] += c;
237            }
238            coefficients
239        };
240        Self::Output { coefficients }
241    }
242}
243
244impl<F> AddAssign for Polynomial<F>
245where
246    F: Add<Output = F> + AddAssign + Clone,
247{
248    fn add_assign(&mut self, rhs: Self) {
249        if self.coefficients.len() >= rhs.coefficients.len() {
250            for (i, c) in rhs.coefficients.into_iter().enumerate() {
251                self.coefficients[i] += c;
252            }
253        } else {
254            let mut coefficients = rhs.coefficients;
255            for (i, c) in self.coefficients.iter().enumerate() {
256                coefficients[i] += c.clone();
257            }
258            self.coefficients = coefficients;
259        }
260    }
261}
262
263impl<F> Sub for &Polynomial<F>
264where
265    F: Sub<Output = F> + Clone + Neg<Output = F> + Add<Output = F> + AddAssign,
266{
267    type Output = Polynomial<F>;
268
269    fn sub(self, rhs: Self) -> Self::Output {
270        self + &(-rhs)
271    }
272}
273
274impl<F> Sub for Polynomial<F>
275where
276    F: Sub<Output = F> + Clone + Neg<Output = F> + Add<Output = F> + AddAssign,
277{
278    type Output = Polynomial<F>;
279
280    fn sub(self, rhs: Self) -> Self::Output {
281        self + (-rhs)
282    }
283}
284
285impl<F> SubAssign for Polynomial<F>
286where
287    F: Add<Output = F> + Neg<Output = F> + AddAssign + Clone + Sub<Output = F>,
288{
289    fn sub_assign(&mut self, rhs: Self) {
290        self.coefficients = self.clone().sub(rhs).coefficients;
291    }
292}
293
294impl<F: Neg<Output = F> + Clone> Neg for &Polynomial<F> {
295    type Output = Polynomial<F>;
296
297    fn neg(self) -> Self::Output {
298        Self::Output {
299            coefficients: self.coefficients.iter().cloned().map(|a| -a).collect(),
300        }
301    }
302}
303
304impl<F: Neg<Output = F> + Clone> Neg for Polynomial<F> {
305    type Output = Self;
306
307    fn neg(self) -> Self::Output {
308        Self::Output {
309            coefficients: self.coefficients.iter().cloned().map(|a| -a).collect(),
310        }
311    }
312}
313
314impl<F> Mul for &Polynomial<F>
315where
316    F: Add + AddAssign + Mul<Output = F> + Sub<Output = F> + Zero + PartialEq + Clone,
317{
318    type Output = Polynomial<F>;
319
320    fn mul(self, other: Self) -> Self::Output {
321        if self.is_zero() || other.is_zero() {
322            return Polynomial::<F>::zero();
323        }
324        let mut coefficients =
325            vec![F::zero(); self.coefficients.len() + other.coefficients.len() - 1];
326        for i in 0..self.coefficients.len() {
327            for j in 0..other.coefficients.len() {
328                coefficients[i + j] += self.coefficients[i].clone() * other.coefficients[j].clone();
329            }
330        }
331        Polynomial { coefficients }
332    }
333}
334
335impl<F> Mul for Polynomial<F>
336where
337    F: Add + AddAssign + Mul<Output = F> + Zero + PartialEq + Clone,
338{
339    type Output = Self;
340
341    fn mul(self, other: Self) -> Self::Output {
342        if self.is_zero() || other.is_zero() {
343            return Self::zero();
344        }
345        let mut coefficients =
346            vec![F::zero(); self.coefficients.len() + other.coefficients.len() - 1];
347        for i in 0..self.coefficients.len() {
348            for j in 0..other.coefficients.len() {
349                coefficients[i + j] += self.coefficients[i].clone() * other.coefficients[j].clone();
350            }
351        }
352        Self { coefficients }
353    }
354}
355
356impl<F: Add + Mul<Output = F> + Zero + Clone> Mul<F> for &Polynomial<F> {
357    type Output = Polynomial<F>;
358
359    fn mul(self, other: F) -> Self::Output {
360        Polynomial {
361            coefficients: self.coefficients.iter().cloned().map(|i| i * other.clone()).collect(),
362        }
363    }
364}
365
366impl<F: Add + Mul<Output = F> + Zero + Clone> Mul<F> for Polynomial<F> {
367    type Output = Polynomial<F>;
368
369    fn mul(self, other: F) -> Self::Output {
370        Polynomial {
371            coefficients: self.coefficients.iter().cloned().map(|i| i * other.clone()).collect(),
372        }
373    }
374}
375
376impl<F: Mul<Output = F> + Sub<Output = F> + AddAssign + Zero + Div<Output = F> + Clone>
377    Polynomial<F>
378{
379    /// Multiply two polynomials using Karatsuba's divide-and-conquer algorithm.
380    pub fn karatsuba(&self, other: &Self) -> Self {
381        Polynomial::new(vector_karatsuba(&self.coefficients, &other.coefficients))
382    }
383}
384
385impl<F> One for Polynomial<F>
386where
387    F: Clone + One + PartialEq + Zero + AddAssign,
388{
389    fn one() -> Self {
390        Self { coefficients: vec![F::one()] }
391    }
392}
393
394impl<F> Zero for Polynomial<F>
395where
396    F: Zero + PartialEq + Clone + AddAssign,
397{
398    fn zero() -> Self {
399        Self { coefficients: vec![] }
400    }
401
402    fn is_zero(&self) -> bool {
403        self.degree().is_none()
404    }
405}
406
407impl<F: Zero + Clone> Polynomial<F> {
408    /// Shifts the polynomial by the specified amount (adds leading zeros).
409    pub fn shift(&self, shamt: usize) -> Self {
410        Self {
411            coefficients: [vec![F::zero(); shamt], self.coefficients.clone()].concat(),
412        }
413    }
414
415    /// Creates a constant polynomial with a single coefficient.
416    pub fn constant(f: F) -> Self {
417        Self { coefficients: vec![f] }
418    }
419
420    /// Applies a function to each coefficient and returns a new polynomial.
421    pub fn map<G: Clone, C: FnMut(&F) -> G>(&self, closure: C) -> Polynomial<G> {
422        Polynomial::<G>::new(self.coefficients.iter().map(closure).collect())
423    }
424
425    /// Folds the coefficients using the provided function and initial value.
426    pub fn fold<G, C: FnMut(G, &F) -> G + Clone>(&self, mut initial_value: G, closure: C) -> G {
427        for c in self.coefficients.iter() {
428            initial_value = (closure.clone())(initial_value, c);
429        }
430        initial_value
431    }
432}
433
434impl<F> Div<Polynomial<F>> for Polynomial<F>
435where
436    F: Zero
437        + One
438        + PartialEq
439        + AddAssign
440        + Clone
441        + Mul<Output = F>
442        + MulAssign
443        + Div<Output = F>
444        + Neg<Output = F>
445        + Sub<Output = F>,
446{
447    type Output = Polynomial<F>;
448
449    fn div(self, denominator: Self) -> Self::Output {
450        if denominator.is_zero() {
451            panic!();
452        }
453        if self.is_zero() {
454            Self::zero();
455        }
456        let mut remainder = self;
457        let mut quotient = Polynomial::<F>::zero();
458        while remainder.degree().unwrap() >= denominator.degree().unwrap() {
459            let shift = remainder.degree().unwrap() - denominator.degree().unwrap();
460            let quotient_coefficient = remainder.lc() / denominator.lc();
461            let monomial = Self::constant(quotient_coefficient).shift(shift);
462            quotient += monomial.clone();
463            remainder -= monomial * denominator.clone();
464            if remainder.is_zero() {
465                break;
466            }
467        }
468        quotient
469    }
470}
471
472fn vector_karatsuba<
473    F: Zero + AddAssign + Mul<Output = F> + Sub<Output = F> + Div<Output = F> + Clone,
474>(
475    left: &[F],
476    right: &[F],
477) -> Vec<F> {
478    let n = left.len();
479    if n <= 8 {
480        let mut product = vec![F::zero(); left.len() + right.len() - 1];
481        for (i, l) in left.iter().enumerate() {
482            for (j, r) in right.iter().enumerate() {
483                product[i + j] += l.clone() * r.clone();
484            }
485        }
486        return product;
487    }
488    let n_over_2 = n / 2;
489    let mut product = vec![F::zero(); 2 * n - 1];
490    let left_lo = &left[0..n_over_2];
491    let right_lo = &right[0..n_over_2];
492    let left_hi = &left[n_over_2..];
493    let right_hi = &right[n_over_2..];
494    let left_sum: Vec<F> =
495        left_lo.iter().zip(left_hi).map(|(a, b)| a.clone() + b.clone()).collect();
496    let right_sum: Vec<F> =
497        right_lo.iter().zip(right_hi).map(|(a, b)| a.clone() + b.clone()).collect();
498
499    let prod_lo = vector_karatsuba(left_lo, right_lo);
500    let prod_hi = vector_karatsuba(left_hi, right_hi);
501    let prod_mid: Vec<F> = vector_karatsuba(&left_sum, &right_sum)
502        .iter()
503        .zip(prod_lo.iter().zip(prod_hi.iter()))
504        .map(|(s, (l, h))| s.clone() - (l.clone() + h.clone()))
505        .collect();
506
507    for (i, l) in prod_lo.into_iter().enumerate() {
508        product[i] = l;
509    }
510    for (i, m) in prod_mid.into_iter().enumerate() {
511        product[i + n_over_2] += m;
512    }
513    for (i, h) in prod_hi.into_iter().enumerate() {
514        product[i + n] += h
515    }
516    product
517}
518
519impl From<Polynomial<FalconFelt>> for Polynomial<Felt> {
520    fn from(item: Polynomial<FalconFelt>) -> Self {
521        let res: Vec<Felt> =
522            item.coefficients.iter().map(|a| Felt::from_u16(a.value() as u16)).collect();
523        Polynomial::new(res)
524    }
525}
526
527impl From<&Polynomial<FalconFelt>> for Polynomial<Felt> {
528    fn from(item: &Polynomial<FalconFelt>) -> Self {
529        let res: Vec<Felt> =
530            item.coefficients.iter().map(|a| Felt::from_u16(a.value() as u16)).collect();
531        Polynomial::new(res)
532    }
533}
534
535impl From<Polynomial<i16>> for Polynomial<FalconFelt> {
536    fn from(item: Polynomial<i16>) -> Self {
537        let res: Vec<FalconFelt> = item.coefficients.iter().map(|&a| FalconFelt::new(a)).collect();
538        Polynomial::new(res)
539    }
540}
541
542impl From<&Polynomial<i16>> for Polynomial<FalconFelt> {
543    fn from(item: &Polynomial<i16>) -> Self {
544        let res: Vec<FalconFelt> = item.coefficients.iter().map(|&a| FalconFelt::new(a)).collect();
545        Polynomial::new(res)
546    }
547}
548
549impl From<Vec<i16>> for Polynomial<FalconFelt> {
550    fn from(item: Vec<i16>) -> Self {
551        let res: Vec<FalconFelt> = item.iter().map(|&a| FalconFelt::new(a)).collect();
552        Polynomial::new(res)
553    }
554}
555
556impl From<&Vec<i16>> for Polynomial<FalconFelt> {
557    fn from(item: &Vec<i16>) -> Self {
558        let res: Vec<FalconFelt> = item.iter().map(|&a| FalconFelt::new(a)).collect();
559        Polynomial::new(res)
560    }
561}
562
563impl Polynomial<FalconFelt> {
564    /// Computes the squared L2 norm of the polynomial.
565    pub fn norm_squared(&self) -> u64 {
566        self.coefficients
567            .iter()
568            .map(|&i| i.balanced_value() as i64)
569            .map(|i| (i * i) as u64)
570            .sum::<u64>()
571    }
572
573    // PUBLIC ACCESSORS
574    // --------------------------------------------------------------------------------------------
575
576    /// Returns the coefficients of this polynomial as field elements.
577    pub fn to_elements(&self) -> Vec<Felt> {
578        self.coefficients.iter().map(|&a| Felt::from_u16(a.value() as u16)).collect()
579    }
580
581    /// Returns the coefficients of this polynomial as balanced signed values.
582    pub fn to_balanced_values(&self) -> Vec<i16> {
583        self.coefficients.iter().copied().map(FalconFelt::balanced_value).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.clone()
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
