Skip to main content

falcon_rust/
polynomial.rs

1use falcon_profiler::profiling;
2use num::{One, Zero};
3use shake::{ExtendableOutput, Shake256, Update, XofReader};
4use std::default::Default;
5use std::fmt::{Debug, Display};
6use std::ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub, SubAssign};
7
8use itertools::Itertools;
9
10use crate::falcon_field::Felt;
11use crate::fixed_point::FixedPoint64;
12use crate::inverse::Inverse;
13
14/// Marked pub for benchmarking; not considered part of the public API.
15#[doc(hidden)]
16#[derive(Debug, Clone)]
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<F> Polynomial<F>
31where
32    F: Clone + Neg<Output = F>,
33{
34    /// Compute the Hermitian adjoint of the polynomial f in the
35    /// cyclotomic ring $\mathbb{Q}[X] / \langle X^n + 1 \rangle$ where $n \geq \deg(f)+1$.
36    /// In this structure, the Hermitian adjoint is given by
37    ///
38    /// $$ f*(X) = f_{[0]} + \sum_{i=1}^{n-1} f_{[i]} * X({n-i}) . $$
39    #[allow(dead_code)]
40    pub fn hermitian_adjoint(&self) -> Polynomial<F> {
41        let coefficients = [
42            vec![self.coefficients[0].clone()],
43            self.coefficients
44                .iter()
45                .skip(1)
46                .cloned()
47                .map(|c| -c)
48                .rev()
49                .collect_vec(),
50        ]
51        .concat();
52        Polynomial { coefficients }
53    }
54}
55
56#[profiling]
57fn vector_karatsuba<
58    F: Zero + AddAssign + Mul<Output = F> + Sub<Output = F> + Div<Output = F> + Clone,
59>(
60    left: &[F],
61    right: &[F],
62) -> Vec<F> {
63    let n = left.len();
64    if n <= 8 {
65        let mut product = vec![F::zero(); left.len() + right.len() - 1];
66        for (i, l) in left.iter().enumerate() {
67            for (j, r) in right.iter().enumerate() {
68                product[i + j] += l.clone() * r.clone();
69            }
70        }
71        return product;
72    }
73    let n_over_2 = n / 2;
74    let mut product = vec![F::zero(); 2 * n - 1];
75    let left_lo = &left[0..n_over_2];
76    let right_lo = &right[0..n_over_2];
77    let left_hi = &left[n_over_2..];
78    let right_hi = &right[n_over_2..];
79    let left_sum = left_lo
80        .iter()
81        .zip(left_hi)
82        .map(|(a, b)| a.clone() + b.clone())
83        .collect_vec();
84    let right_sum = right_lo
85        .iter()
86        .zip(right_hi)
87        .map(|(a, b)| a.clone() + b.clone())
88        .collect_vec();
89
90    let prod_lo = vector_karatsuba(left_lo, right_lo);
91    let prod_hi = vector_karatsuba(left_hi, right_hi);
92    let prod_mid = vector_karatsuba(&left_sum, &right_sum)
93        .iter()
94        .zip(prod_lo.iter().zip(prod_hi.iter()))
95        .map(|(s, (l, h))| s.clone() - (l.clone() + h.clone()))
96        .collect_vec();
97
98    for (i, l) in prod_lo.into_iter().enumerate() {
99        product[i] = l;
100    }
101    for (i, m) in prod_mid.into_iter().enumerate() {
102        product[i + n_over_2] += m;
103    }
104    for (i, h) in prod_hi.into_iter().enumerate() {
105        product[i + n] += h
106    }
107    product
108}
109
110#[allow(private_bounds)] // The module is marked `pub` only for benchmarking.
111impl<
112        F: Mul<Output = F> + Sub<Output = F> + AddAssign + Zero + Div<Output = F> + Inverse + Clone,
113    > Polynomial<F>
114{
115    #[profiling]
116    pub fn hadamard_mul(&self, other: &Self) -> Self {
117        Polynomial::new(
118            self.coefficients
119                .iter()
120                .zip(other.coefficients.iter())
121                .map(|(a, b)| *a * *b)
122                .collect_vec(),
123        )
124    }
125    #[profiling]
126    pub fn hadamard_div(&self, other: &Self) -> Self {
127        let other_coefficients_inverse = F::batch_inverse_or_zero(&other.coefficients);
128        Polynomial::new(
129            self.coefficients
130                .iter()
131                .zip(other_coefficients_inverse.iter())
132                .map(|(a, b)| *a * *b)
133                .collect_vec(),
134        )
135    }
136    #[profiling]
137    pub fn hadamard_inv(&self) -> Self {
138        let coefficients_inverse = F::batch_inverse_or_zero(&self.coefficients);
139        Polynomial::new(coefficients_inverse)
140    }
141}
142
143impl<F: Mul<Output = F> + Sub<Output = F> + AddAssign + Zero + Div<Output = F> + Clone>
144    Polynomial<F>
145{
146    /// Multiply two polynomials using Karatsuba's divide-and-conquer algorithm.
147    #[profiling]
148    pub fn karatsuba(&self, other: &Self) -> Self {
149        Polynomial::new(vector_karatsuba(&self.coefficients, &other.coefficients))
150    }
151}
152
153impl<F: Zero + PartialEq + Clone> Polynomial<F> {
154    pub fn degree(&self) -> Option<usize> {
155        if self.coefficients.is_empty() {
156            return None;
157        }
158        let mut max_index = self.coefficients.len() - 1;
159        while self.coefficients[max_index] == F::zero() {
160            if let Some(new_index) = max_index.checked_sub(1) {
161                max_index = new_index;
162            } else {
163                return None;
164            }
165        }
166        Some(max_index)
167    }
168    pub fn lc(&self) -> F {
169        match self.degree() {
170            Some(non_negative_degree) => self.coefficients[non_negative_degree].clone(),
171            None => F::zero(),
172        }
173    }
174}
175
176impl<F: Zero + Clone> Polynomial<F> {
177    pub fn shift(&self, shamt: usize) -> Self {
178        Self {
179            coefficients: [vec![F::zero(); shamt], self.coefficients.clone()].concat(),
180        }
181    }
182
183    pub fn constant(f: F) -> Self {
184        Self {
185            coefficients: vec![f],
186        }
187    }
188
189    pub fn map<G: Clone, C: FnMut(&F) -> G>(&self, closure: C) -> Polynomial<G> {
190        Polynomial::<G>::new(self.coefficients.iter().map(closure).collect_vec())
191    }
192
193    pub fn fold<G, C: FnMut(G, &F) -> G + Clone>(&self, mut initial_value: G, closure: C) -> G {
194        for c in self.coefficients.iter() {
195            initial_value = (closure.clone())(initial_value, c);
196        }
197        initial_value
198    }
199}
200
201/// The following implementations are specific to cyclotomic polynomial rings,
202/// i.e., F[ X ] / <X^n + 1>, and are used extensively in Falcon.
203impl<
204        F: One + Zero + Clone + Neg<Output = F> + MulAssign + AddAssign + Sub<Output = F> + PartialEq,
205    > Polynomial<F>
206{
207    /// Reduce the polynomial by X^n + 1.
208    #[profiling]
209    pub fn reduce_by_cyclotomic(&self, n: usize) -> Self {
210        let mut coefficients = vec![F::zero(); n];
211        let mut sign = -F::one();
212        for (i, c) in self.coefficients.iter().cloned().enumerate() {
213            if i % n == 0 {
214                sign *= -F::one();
215            }
216            coefficients[i % n] += sign.clone() * c;
217        }
218        Polynomial::new(coefficients)
219    }
220}
221
222/// The following implementations are specific to cyclotomic polynomial rings,
223/// i.e., F[ X ] / <X^n + 1>, and are used extensively in Falcon.
224impl<
225        F: One
226            + Zero
227            + Clone
228            + Neg<Output = F>
229            + MulAssign
230            + AddAssign
231            + Div<Output = F>
232            + Sub<Output = F>
233            + PartialEq,
234    > Polynomial<F>
235{
236    /// Compute the multiplicative inverse of the polynomial in the ring
237    /// F[ X ] / <X^n + 1>
238    ///
239    /// This function assumes that F is a field; otherwise the gcd will never end.
240    #[allow(dead_code)]
241    pub(crate) fn cyclotomic_ring_inverse(&self, n: usize) -> Self {
242        let mut cyclotomic_coefficients = vec![F::zero(); n + 1];
243        cyclotomic_coefficients[0] = F::one();
244        cyclotomic_coefficients[n] = F::one();
245        let (_, a, _) = Polynomial::xgcd(self, &Polynomial::new(cyclotomic_coefficients));
246        a
247    }
248
249    /// Compute the field norm of the polynomial as an element of the cyclotomic
250    /// ring  F[ X ] / <X^n + 1 > relative to one of half the size, i.e.,
251    ///  F[ X ] / <X^(n/2) + 1> .
252    ///
253    /// Corresponds to formula 3.25 in the spec [1, p.30].
254    ///
255    /// [1]: https://falcon-sign.info/falcon.pdf
256    ///
257    #[profiling]
258    pub fn field_norm(&self) -> Self {
259        let n = self.coefficients.len();
260        let mut f0_coefficients = vec![F::zero(); n / 2];
261        let mut f1_coefficients = vec![F::zero(); n / 2];
262        for i in 0..n / 2 {
263            f0_coefficients[i] = self.coefficients[2 * i].clone();
264            f1_coefficients[i] = self.coefficients[2 * i + 1].clone();
265        }
266        let f0 = Polynomial::new(f0_coefficients);
267        let f1 = Polynomial::new(f1_coefficients);
268        let f0_squared = (f0.clone() * f0).reduce_by_cyclotomic(n / 2);
269        let f1_squared = (f1.clone() * f1).reduce_by_cyclotomic(n / 2);
270        let x = Polynomial::new(vec![F::zero(), F::one()]);
271        f0_squared - (x * f1_squared).reduce_by_cyclotomic(n / 2)
272    }
273
274    /// Lift an element from a cyclotomic polynomial ring to one of double the
275    /// size. Do this by interleaving zeros.
276    ///
277    #[profiling]
278    pub fn lift_next_cyclotomic(&self) -> Self {
279        let n = self.coefficients.len();
280        let mut coefficients = vec![F::zero(); n * 2];
281        for i in 0..n {
282            coefficients[2 * i] = self.coefficients[i].clone();
283        }
284        Self::new(coefficients)
285    }
286
287    /// Compute the galois adjoint of the polynomial in the cyclotomic ring
288    /// F[ X ] / < X^n + 1 > , which corresponds to flipping the sign of all
289    /// odd coefficients.
290    ///
291    #[profiling]
292    pub fn galois_adjoint(&self) -> Self {
293        Self::new(
294            self.coefficients
295                .iter()
296                .enumerate()
297                .map(|(i, c)| {
298                    if i % 2 == 0 {
299                        c.clone()
300                    } else {
301                        c.clone().neg()
302                    }
303                })
304                .collect_vec(),
305        )
306    }
307}
308
309impl<
310        F: One
311            + Zero
312            + Clone
313            + Neg<Output = F>
314            + MulAssign
315            + AddAssign
316            + Div<Output = F>
317            + Sub<Output = F>
318            + PartialEq,
319    > Polynomial<F>
320{
321    /// Extended Euclidean algorithm for polynomials. Uses the EEA to compute
322    /// the greatest common divisor g and Bezout coefficients u, v such that
323    ///
324    /// $$  u a + v b = 1 . $$
325    ///
326    /// Implementation adapted from Wikipedia [1].
327    ///
328    /// [1]: https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Pseudocode
329    ///
330    #[profiling]
331    pub(crate) fn xgcd(a: &Self, b: &Self) -> (Self, Self, Self) {
332        if a.is_zero() || b.is_zero() {
333            return (Self::zero(), Self::zero(), Self::zero());
334        }
335        let (mut old_r, mut r) = (a.clone(), b.clone());
336        let (mut old_s, mut s) = (Self::one(), Self::zero());
337        let (mut old_t, mut t) = (Self::zero(), Self::one());
338
339        while !r.is_zero() {
340            let quotient = old_r.clone() / r.clone();
341            (old_r, r) = (r.clone(), old_r.clone() - quotient.clone() * r.clone());
342            (old_s, s) = (s.clone(), old_s.clone() - quotient.clone() * s.clone());
343            (old_t, t) = (t.clone(), old_t.clone() - quotient.clone() * t.clone());
344        }
345
346        (old_r, old_s, old_t)
347    }
348}
349
350#[allow(private_bounds)]
351impl<F: Clone + Into<FixedPoint64>> Polynomial<F> {
352    #[allow(dead_code)]
353    #[profiling]
354    pub(crate) fn l2_norm(&self) -> FixedPoint64 {
355        self.l2_norm_squared().sqrt()
356    }
357    #[profiling]
358    pub(crate) fn l2_norm_squared(&self) -> FixedPoint64 {
359        self.coefficients
360            .iter()
361            .map(|i| Into::<FixedPoint64>::into(i.clone()))
362            .map(|i| i * i)
363            .sum::<FixedPoint64>()
364    }
365}
366
367impl<F> PartialEq for Polynomial<F>
368where
369    F: Zero + PartialEq + Clone + AddAssign,
370{
371    fn eq(&self, other: &Self) -> bool {
372        if self.is_zero() && other.is_zero() {
373            true
374        } else if self.is_zero() || other.is_zero() {
375            false
376        } else {
377            let self_degree = self.degree().unwrap();
378            let other_degree = other.degree().unwrap();
379            self.coefficients[0..=self_degree] == other.coefficients[0..=other_degree]
380        }
381    }
382}
383
384impl<F> Eq for Polynomial<F> where F: Zero + PartialEq + Clone + AddAssign {}
385
386impl<F> Add for &Polynomial<F>
387where
388    F: Add<Output = F> + AddAssign + Clone,
389{
390    type Output = Polynomial<F>;
391
392    fn add(self, rhs: Self) -> Self::Output {
393        let coefficients = if self.coefficients.len() >= rhs.coefficients.len() {
394            let mut coefficients = self.coefficients.clone();
395            for (i, c) in rhs.coefficients.iter().enumerate() {
396                coefficients[i] += c.clone();
397            }
398            coefficients
399        } else {
400            let mut coefficients = rhs.coefficients.clone();
401            for (i, c) in self.coefficients.iter().enumerate() {
402                coefficients[i] += c.clone();
403            }
404            coefficients
405        };
406        Self::Output { coefficients }
407    }
408}
409impl<F> Add for Polynomial<F>
410where
411    F: Add<Output = F> + AddAssign + Clone,
412{
413    type Output = Polynomial<F>;
414    fn add(self, rhs: Self) -> Self::Output {
415        let coefficients = if self.coefficients.len() >= rhs.coefficients.len() {
416            let mut coefficients = self.coefficients.clone();
417            for (i, c) in rhs.coefficients.into_iter().enumerate() {
418                coefficients[i] += c;
419            }
420            coefficients
421        } else {
422            let mut coefficients = rhs.coefficients.clone();
423            for (i, c) in self.coefficients.into_iter().enumerate() {
424                coefficients[i] += c;
425            }
426            coefficients
427        };
428        Self::Output { coefficients }
429    }
430}
431
432impl<F> AddAssign for Polynomial<F>
433where
434    F: Add<Output = F> + AddAssign + Clone,
435{
436    fn add_assign(&mut self, rhs: Self) {
437        if self.coefficients.len() >= rhs.coefficients.len() {
438            for (i, c) in rhs.coefficients.into_iter().enumerate() {
439                self.coefficients[i] += c;
440            }
441        } else {
442            let mut coefficients = rhs.coefficients.clone();
443            for (i, c) in self.coefficients.iter().enumerate() {
444                coefficients[i] += c.clone();
445            }
446            self.coefficients = coefficients;
447        }
448    }
449}
450
451impl<F> Sub for &Polynomial<F>
452where
453    F: Sub<Output = F> + Clone + Neg<Output = F> + Add<Output = F> + AddAssign,
454{
455    type Output = Polynomial<F>;
456
457    fn sub(self, rhs: Self) -> Self::Output {
458        self + &(-rhs)
459    }
460}
461impl<F> Sub for Polynomial<F>
462where
463    F: Sub<Output = F> + Clone + Neg<Output = F> + Add<Output = F> + AddAssign,
464{
465    type Output = Polynomial<F>;
466
467    fn sub(self, rhs: Self) -> Self::Output {
468        self + (-rhs)
469    }
470}
471
472impl<F> SubAssign for Polynomial<F>
473where
474    F: Add<Output = F> + Neg<Output = F> + AddAssign + Clone + Sub<Output = F>,
475{
476    fn sub_assign(&mut self, rhs: Self) {
477        self.coefficients = self.clone().sub(rhs).coefficients;
478    }
479}
480
481impl<F: Neg<Output = F> + Clone> Neg for &Polynomial<F> {
482    type Output = Polynomial<F>;
483
484    fn neg(self) -> Self::Output {
485        Self::Output {
486            coefficients: self.coefficients.iter().cloned().map(|a| -a).collect(),
487        }
488    }
489}
490impl<F: Neg<Output = F> + Clone> Neg for Polynomial<F> {
491    type Output = Self;
492
493    fn neg(self) -> Self::Output {
494        Self::Output {
495            coefficients: self.coefficients.iter().cloned().map(|a| -a).collect(),
496        }
497    }
498}
499
500impl<F> Mul for &Polynomial<F>
501where
502    F: Add + AddAssign + Mul<Output = F> + Sub<Output = F> + Zero + PartialEq + Clone,
503{
504    type Output = Polynomial<F>;
505
506    fn mul(self, other: Self) -> Self::Output {
507        if self.is_zero() || other.is_zero() {
508            return Polynomial::<F>::zero();
509        }
510        let mut coefficients =
511            vec![F::zero(); self.coefficients.len() + other.coefficients.len() - 1];
512        for i in 0..self.coefficients.len() {
513            for j in 0..other.coefficients.len() {
514                coefficients[i + j] += self.coefficients[i].clone() * other.coefficients[j].clone();
515            }
516        }
517        Polynomial { coefficients }
518    }
519}
520
521impl<F> Mul for Polynomial<F>
522where
523    F: Add + AddAssign + Mul<Output = F> + Zero + PartialEq + Clone,
524{
525    type Output = Self;
526
527    fn mul(self, other: Self) -> Self::Output {
528        if self.is_zero() || other.is_zero() {
529            return Self::zero();
530        }
531        let mut coefficients =
532            vec![F::zero(); self.coefficients.len() + other.coefficients.len() - 1];
533        for i in 0..self.coefficients.len() {
534            for j in 0..other.coefficients.len() {
535                coefficients[i + j] += self.coefficients[i].clone() * other.coefficients[j].clone();
536            }
537        }
538        Self { coefficients }
539    }
540}
541
542impl<F: Add + Mul<Output = F> + Zero + Clone> Mul<F> for &Polynomial<F> {
543    type Output = Polynomial<F>;
544
545    fn mul(self, other: F) -> Self::Output {
546        Polynomial {
547            coefficients: self
548                .coefficients
549                .iter()
550                .cloned()
551                .map(|i| i * other.clone())
552                .collect_vec(),
553        }
554    }
555}
556impl<F: Add + Mul<Output = F> + Zero + Clone> Mul<F> for Polynomial<F> {
557    type Output = Polynomial<F>;
558
559    fn mul(self, other: F) -> Self::Output {
560        Polynomial {
561            coefficients: self
562                .coefficients
563                .iter()
564                .cloned()
565                .map(|i| i * other.clone())
566                .collect_vec(),
567        }
568    }
569}
570impl<F> One for Polynomial<F>
571where
572    F: Clone + One + PartialEq + Zero + AddAssign,
573{
574    fn one() -> Self {
575        Self {
576            coefficients: vec![F::one()],
577        }
578    }
579}
580
581impl<F> Zero for Polynomial<F>
582where
583    F: Zero + PartialEq + Clone + AddAssign,
584{
585    fn zero() -> Self {
586        Self {
587            coefficients: vec![],
588        }
589    }
590
591    fn is_zero(&self) -> bool {
592        self.degree().is_none()
593    }
594}
595
596impl<F> Div<Polynomial<F>> for Polynomial<F>
597where
598    F: Zero
599        + One
600        + PartialEq
601        + AddAssign
602        + Clone
603        + Mul<Output = F>
604        + MulAssign
605        + Div<Output = F>
606        + Neg<Output = F>
607        + Sub<Output = F>,
608{
609    type Output = Polynomial<F>;
610
611    fn div(self, denominator: Self) -> Self::Output {
612        if denominator.is_zero() {
613            panic!();
614        }
615        if self.is_zero() {
616            Self::zero();
617        }
618        let mut remainder = self.clone();
619        let mut quotient = Polynomial::<F>::zero();
620        while remainder.degree().unwrap() >= denominator.degree().unwrap() {
621            let shift = remainder.degree().unwrap() - denominator.degree().unwrap();
622            let quotient_coefficient = remainder.lc() / denominator.lc();
623            let monomial = Self::constant(quotient_coefficient).shift(shift);
624            quotient += monomial.clone();
625            remainder -= monomial * denominator.clone();
626            if remainder.is_zero() {
627                break;
628            }
629        }
630        quotient
631    }
632}
633
634/// Hash a string to a random polynomial in ZZ[ X ] mod <Phi(X), q>.
635/// Algorithm 3, "HashToPoint" in the spec (page 31).
636#[profiling]
637pub(crate) fn hash_to_point(string: &[u8], n: usize) -> Polynomial<Felt> {
638    const K: u32 = (1u32 << 16) / (Felt::Q as u32);
639
640    let mut hasher = Shake256::default();
641    hasher.update(string);
642    let mut reader = hasher.finalize_xof();
643
644    let mut coefficients: Vec<Felt> = Vec::with_capacity(n);
645    while coefficients.len() != n {
646        let mut randomness = [0u8; 2];
647        reader.read(&mut randomness);
648        // Arabic endianness but so be it
649        let t = ((randomness[0] as u32) << 8) | (randomness[1] as u32);
650        if t < K * (Felt::Q as u32) {
651            coefficients.push(Felt::new((t % (Felt::Q as u32)) as u16));
652        }
653    }
654
655    Polynomial { coefficients }
656}
657
658impl<T: Display> Display for Polynomial<T> {
659    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
660        write!(f, "[{}]", self.coefficients.iter().join(", "))
661    }
662}
663
664#[cfg(test)]
665mod test {
666    use crate::falcon_field::Felt;
667    use crate::polynomial::hash_to_point;
668    use crate::polynomial::Polynomial;
669    use itertools::Itertools;
670    use rand::rng;
671    use rand::RngExt;
672
673    #[test]
674    fn test_hash_to_point_sanity() {
675        // KAT obtained from python script, run in the same directory
676        // as a copy of the official python implementation.
677        // ```
678        // from falcon import SecretKey
679        // sk = SecretKey(512)
680        // sk.hash_to_point(b"", b"")
681        // ```
682        let hash_of_empty = Polynomial {
683            coefficients: [
684                5816, 7463, 2984, 11537, 9019, 4074, 5180, 11040, 4044, 8937, 694, 7042, 9481,
685                10084, 3795, 5677, 5977, 1241, 6332, 2817, 413, 1971, 755, 7241, 6041, 9347, 4136,
686                11948, 9140, 1210, 5150, 1630, 4015, 2390, 2346, 2025, 4272, 10978, 7171, 8764,
687                11920, 888, 12160, 11275, 7043, 10323, 1181, 1873, 5876, 12213, 627, 11319, 5675,
688                8207, 6210, 385, 333, 4581, 1359, 10859, 3346, 3970, 8720, 3640, 8157, 1080, 2794,
689                5769, 11618, 6780, 1734, 6484, 1575, 9433, 10353, 2004, 5921, 5013, 4753, 9865,
690                10931, 6621, 1417, 9804, 12027, 9437, 10657, 3260, 9541, 4967, 12124, 6827, 333,
691                6404, 3498, 6920, 3979, 14, 440, 1293, 8011, 7567, 3899, 3252, 4023, 10727, 11938,
692                957, 2412, 9552, 10409, 8063, 9131, 9835, 10305, 3124, 6303, 12241, 6354, 2540,
693                10113, 10777, 6803, 4879, 10952, 10503, 1728, 5067, 3339, 7045, 11333, 5469, 11062,
694                11666, 5235, 2314, 3345, 2224, 2274, 8060, 4304, 6716, 11595, 1541, 996, 6983, 36,
695                449, 7401, 4987, 9177, 810, 1908, 8650, 7646, 6893, 4919, 1971, 4930, 11763, 201,
696                12223, 9234, 4081, 6199, 12047, 7646, 9639, 6814, 6739, 5279, 4012, 2101, 10707,
697                4241, 12146, 3779, 3999, 3176, 1699, 10294, 5168, 5590, 457, 9709, 6450, 442, 8884,
698                6755, 10995, 10923, 3935, 8499, 3508, 12088, 1115, 11336, 1379, 7557, 4954, 7639,
699                2514, 8672, 6686, 98, 5676, 8800, 5712, 4724, 7724, 3202, 12128, 10940, 10177,
700                9421, 11013, 7372, 8546, 441, 6261, 8779, 2453, 12082, 7922, 5307, 6920, 7726, 823,
701                10561, 1251, 10358, 8383, 10905, 8145, 1733, 1718, 3105, 10756, 6798, 10209, 7976,
702                11148, 9353, 4746, 1089, 11444, 6571, 409, 8381, 10325, 7649, 10042, 5587, 3625,
703                10182, 10494, 228, 4687, 5949, 7995, 12092, 3312, 5339, 5920, 8145, 6796, 1992,
704                3205, 2761, 12199, 11342, 9695, 390, 252, 989, 1385, 12148, 8324, 10694, 3690,
705                3440, 8888, 12238, 9018, 3354, 5859, 6298, 8098, 4388, 3788, 3045, 11095, 2372,
706                10036, 9233, 168, 8500, 3604, 2494, 9854, 5679, 2182, 3350, 7798, 8310, 3544, 948,
707                7646, 7235, 2650, 6008, 4610, 2159, 6884, 10545, 688, 4115, 10312, 4408, 4951,
708                2891, 9791, 1377, 5909, 11147, 11139, 4969, 5158, 350, 1067, 4242, 10820, 1818,
709                6473, 105, 2919, 10892, 7116, 850, 11409, 2652, 6392, 2540, 6892, 8372, 11975,
710                4994, 2621, 2763, 11837, 6132, 11293, 9138, 8769, 10964, 9826, 601, 7007, 9078,
711                10340, 9410, 8746, 10835, 9053, 11010, 5308, 8851, 1976, 11016, 599, 8348, 9876,
712                7100, 1333, 4550, 1735, 4598, 9970, 525, 8320, 1609, 9213, 4178, 484, 10814, 1760,
713                9667, 8369, 2286, 10384, 12139, 24, 1178, 5682, 7074, 3676, 3661, 3322, 1831, 5562,
714                734, 8059, 8750, 6951, 4760, 10395, 1019, 9404, 2923, 6715, 123, 10157, 4892, 7667,
715                1677, 4175, 3455, 12123, 10730, 2000, 8212, 2665, 7088, 8741, 10936, 3172, 225,
716                3867, 5140, 2310, 6453, 2898, 3637, 4580, 113, 5991, 3532, 3363, 11457, 11601,
717                7280, 6792, 11872, 8127, 2192, 10761, 9019, 8197, 8965, 6061, 10799, 988, 10522,
718                1281, 1965, 2716, 9841, 7496, 8456, 5192, 825, 3727, 4664, 7906, 8521, 5901, 10200,
719                5167, 9451, 10825, 12011, 2272, 8698, 8174, 11973, 5155, 6890, 9999, 4391, 12044,
720                1620, 8310, 111, 4481, 9650, 2077, 7691, 7531, 1956, 494, 3297, 1623, 3266, 7018,
721                2031, 6317, 4657, 5206, 2581, 11227, 10508, 4567, 8892, 1363, 6790, 6180, 1588,
722                9776, 11998, 10689, 492, 331,
723            ]
724            .map(Felt::new)
725            .to_vec(),
726        };
727        assert_eq!(hash_of_empty, hash_to_point(&[], 512));
728    }
729
730    #[test]
731    fn test_div() {
732        let mut rng = rng();
733        let n = rng.random_range(1..100);
734        let m = rng.random_range(1..100);
735        let expected_coefficients = (0..n).map(|_| rng.random_range(-5..5)).collect_vec();
736        let cofactor_coefficients = (0..m).map(|_| rng.random_range(-5..5)).collect_vec();
737        let cofactor_polynomial = Polynomial::new(cofactor_coefficients);
738        let product = Polynomial::new(expected_coefficients.clone()) * cofactor_polynomial.clone();
739        let quotient = product / cofactor_polynomial;
740        assert_eq!(Polynomial::new(expected_coefficients), quotient);
741    }
742
743    #[test]
744    fn test_karatsuba() {
745        let mut rng = rng();
746        let n = 64;
747        let coefficients_lhs = (0..n).map(|_| rng.random_range(-5..5) as f64).collect_vec();
748        let lhs = Polynomial::new(coefficients_lhs);
749        let coefficients_rhs = (0..n).map(|_| rng.random_range(-5..5) as f64).collect_vec();
750        let rhs = Polynomial::new(coefficients_rhs);
751
752        let schoolbook = lhs.clone() * rhs.clone();
753        let karatsuba = lhs.karatsuba(&rhs);
754        let difference = schoolbook - karatsuba;
755        assert!(f64::from(difference.l2_norm()) < f64::EPSILON * 100.0);
756    }
757}