lambdaworks_math/polynomial/
mod.rs

1use super::field::element::FieldElement;
2use crate::field::traits::{IsField, IsPrimeField, IsSubFieldOf};
3use alloc::string::{String, ToString};
4use alloc::{borrow::ToOwned, format, vec, vec::Vec};
5use core::{fmt::Display, ops};
6pub mod dense_multilinear_poly;
7mod error;
8pub mod sparse_multilinear_poly;
9
10/// Represents the polynomial c_0 + c_1 * X + c_2 * X^2 + ... + c_n * X^n
11/// as a vector of coefficients `[c_0, c_1, ... , c_n]`
12#[derive(Debug, Clone, PartialEq, Eq)]
13pub struct Polynomial<FE> {
14    pub coefficients: Vec<FE>,
15}
16
17impl<F: IsField> Polynomial<FieldElement<F>> {
18    /// Creates a new polynomial with the given coefficients
19    pub fn new(coefficients: &[FieldElement<F>]) -> Self {
20        // Removes trailing zero coefficients at the end
21        let mut unpadded_coefficients = coefficients
22            .iter()
23            .rev()
24            .skip_while(|x| **x == FieldElement::zero())
25            .cloned()
26            .collect::<Vec<FieldElement<F>>>();
27        unpadded_coefficients.reverse();
28        Polynomial {
29            coefficients: unpadded_coefficients,
30        }
31    }
32
33    pub fn new_monomial(coefficient: FieldElement<F>, degree: usize) -> Self {
34        let mut coefficients = vec![FieldElement::zero(); degree];
35        coefficients.push(coefficient);
36        Self::new(&coefficients)
37    }
38
39    pub fn zero() -> Self {
40        Self::new(&[])
41    }
42
43    /// Returns a polynomial that interpolates the points with x coordinates and y coordinates given by
44    /// `xs` and `ys`.
45    /// `xs` and `ys` must be the same length, and `xs` values should be unique. If not, panics.
46    pub fn interpolate(
47        xs: &[FieldElement<F>],
48        ys: &[FieldElement<F>],
49    ) -> Result<Self, InterpolateError> {
50        // TODO: try to use the type system to avoid this assert
51        if xs.len() != ys.len() {
52            return Err(InterpolateError::UnequalLengths(xs.len(), ys.len()));
53        }
54        if xs.is_empty() {
55            return Ok(Polynomial::new(&[]));
56        }
57
58        let mut denominators = Vec::with_capacity(xs.len() * (xs.len() - 1) / 2);
59        let mut indexes = Vec::with_capacity(xs.len());
60
61        let mut idx = 0;
62
63        for (i, xi) in xs.iter().enumerate().skip(1) {
64            indexes.push(idx);
65            for xj in xs.iter().take(i) {
66                if xi == xj {
67                    return Err(InterpolateError::NonUniqueXs);
68                }
69                denominators.push(xi - xj);
70                idx += 1;
71            }
72        }
73
74        FieldElement::inplace_batch_inverse(&mut denominators).unwrap();
75
76        let mut result = Polynomial::zero();
77
78        for (i, y) in ys.iter().enumerate() {
79            let mut y_term = Polynomial::new(&[y.clone()]);
80            for (j, x) in xs.iter().enumerate() {
81                if i == j {
82                    continue;
83                }
84                let denominator = if i > j {
85                    denominators[indexes[i - 1] + j].clone()
86                } else {
87                    -&denominators[indexes[j - 1] + i]
88                };
89                let denominator_poly = Polynomial::new(&[denominator]);
90                let numerator = Polynomial::new(&[-x, FieldElement::one()]);
91                y_term = y_term.mul_with_ref(&(numerator * denominator_poly));
92            }
93            result = result + y_term;
94        }
95        Ok(result)
96    }
97
98    pub fn evaluate<E>(&self, x: &FieldElement<E>) -> FieldElement<E>
99    where
100        E: IsField,
101        F: IsSubFieldOf<E>,
102    {
103        self.coefficients
104            .iter()
105            .rev()
106            .fold(FieldElement::zero(), |acc, coeff| {
107                coeff + acc * x.to_owned()
108            })
109    }
110
111    pub fn evaluate_slice(&self, input: &[FieldElement<F>]) -> Vec<FieldElement<F>> {
112        input.iter().map(|x| self.evaluate(x)).collect()
113    }
114
115    pub fn degree(&self) -> usize {
116        if self.coefficients.is_empty() {
117            0
118        } else {
119            self.coefficients.len() - 1
120        }
121    }
122
123    pub fn leading_coefficient(&self) -> FieldElement<F> {
124        if let Some(coefficient) = self.coefficients.last() {
125            coefficient.clone()
126        } else {
127            FieldElement::zero()
128        }
129    }
130
131    /// Returns coefficients of the polynomial as an array
132    /// \[c_0, c_1, c_2, ..., c_n\]
133    /// that represents the polynomial
134    /// c_0 + c_1 * X + c_2 * X^2 + ... + c_n * X^n
135    pub fn coefficients(&self) -> &[FieldElement<F>] {
136        &self.coefficients
137    }
138
139    pub fn coeff_len(&self) -> usize {
140        self.coefficients().len()
141    }
142
143    /// Returns the derivative of the polynomial with respect to x.
144    pub fn differentiate(&self) -> Self {
145        let degree = self.degree();
146        if degree == 0 {
147            return Polynomial::zero();
148        }
149        let mut derivative = Vec::with_capacity(degree);
150        for (i, coeff) in self.coefficients().iter().enumerate().skip(1) {
151            derivative.push(FieldElement::<F>::from(i as u64) * coeff);
152        }
153        Polynomial::new(&derivative)
154    }
155
156    /// Computes quotient with `x - b` in place.
157    pub fn ruffini_division_inplace(&mut self, b: &FieldElement<F>) {
158        let mut c = FieldElement::zero();
159        for coeff in self.coefficients.iter_mut().rev() {
160            *coeff = &*coeff + b * &c;
161            core::mem::swap(coeff, &mut c);
162        }
163        self.coefficients.pop();
164    }
165
166    pub fn ruffini_division<L>(&self, b: &FieldElement<L>) -> Polynomial<FieldElement<L>>
167    where
168        L: IsField,
169        F: IsSubFieldOf<L>,
170    {
171        if let Some(c) = self.coefficients.last() {
172            let mut c = c.clone().to_extension();
173            let mut coefficients = Vec::with_capacity(self.degree());
174            for coeff in self.coefficients.iter().rev().skip(1) {
175                coefficients.push(c.clone());
176                c = coeff + c * b;
177            }
178            coefficients = coefficients.into_iter().rev().collect();
179            Polynomial::new(&coefficients)
180        } else {
181            Polynomial::zero()
182        }
183    }
184
185    /// Computes quotient and remainder of polynomial division.
186    ///
187    /// Output: (quotient, remainder)
188    pub fn long_division_with_remainder(self, dividend: &Self) -> (Self, Self) {
189        if dividend.degree() > self.degree() {
190            (Polynomial::zero(), self)
191        } else {
192            let mut n = self;
193            let mut q: Vec<FieldElement<F>> = vec![FieldElement::zero(); n.degree() + 1];
194            let denominator = dividend.leading_coefficient().inv().unwrap();
195            while n != Polynomial::zero() && n.degree() >= dividend.degree() {
196                let new_coefficient = n.leading_coefficient() * &denominator;
197                q[n.degree() - dividend.degree()] = new_coefficient.clone();
198                let d = dividend.mul_with_ref(&Polynomial::new_monomial(
199                    new_coefficient,
200                    n.degree() - dividend.degree(),
201                ));
202                n = n - d;
203            }
204            (Polynomial::new(&q), n)
205        }
206    }
207
208    /// Extended Euclidean Algorithm for polynomials.
209    ///
210    /// This method computes the extended greatest common divisor (GCD) of two polynomials `self` and `y`.
211    /// It returns a tuple of three elements: `(a, b, g)` such that `a * self + b * y = g`, where `g` is the
212    /// greatest common divisor of `self` and `y`.
213    pub fn xgcd(&self, y: &Self) -> (Self, Self, Self) {
214        let one = Polynomial::new(&[FieldElement::one()]);
215        let zero = Polynomial::zero();
216        let (mut old_r, mut r) = (self.clone(), y.clone());
217        let (mut old_s, mut s) = (one.clone(), zero.clone());
218        let (mut old_t, mut t) = (zero.clone(), one.clone());
219
220        while r != Polynomial::zero() {
221            let quotient = old_r.clone().div_with_ref(&r);
222            old_r = old_r - &quotient * &r;
223            core::mem::swap(&mut old_r, &mut r);
224            old_s = old_s - &quotient * &s;
225            core::mem::swap(&mut old_s, &mut s);
226            old_t = old_t - &quotient * &t;
227            core::mem::swap(&mut old_t, &mut t);
228        }
229
230        let lcinv = old_r.leading_coefficient().inv().unwrap();
231        (
232            old_s.scale_coeffs(&lcinv),
233            old_t.scale_coeffs(&lcinv),
234            old_r.scale_coeffs(&lcinv),
235        )
236    }
237
238    pub fn div_with_ref(self, dividend: &Self) -> Self {
239        let (quotient, _remainder) = self.long_division_with_remainder(dividend);
240        quotient
241    }
242
243    pub fn mul_with_ref(&self, factor: &Self) -> Self {
244        let degree = self.degree() + factor.degree();
245        let mut coefficients = vec![FieldElement::zero(); degree + 1];
246
247        if self.coefficients.is_empty() || factor.coefficients.is_empty() {
248            Polynomial::new(&[FieldElement::zero()])
249        } else {
250            for i in 0..=factor.degree() {
251                for j in 0..=self.degree() {
252                    coefficients[i + j] += &factor.coefficients[i] * &self.coefficients[j];
253                }
254            }
255            Polynomial::new(&coefficients)
256        }
257    }
258
259    pub fn scale<S: IsSubFieldOf<F>>(&self, factor: &FieldElement<S>) -> Self {
260        let scaled_coefficients = self
261            .coefficients
262            .iter()
263            .zip(core::iter::successors(Some(FieldElement::one()), |x| {
264                Some(x * factor)
265            }))
266            .map(|(coeff, power)| power * coeff)
267            .collect();
268        Self {
269            coefficients: scaled_coefficients,
270        }
271    }
272
273    pub fn scale_coeffs(&self, factor: &FieldElement<F>) -> Self {
274        let scaled_coefficients = self
275            .coefficients
276            .iter()
277            .map(|coeff| factor * coeff)
278            .collect();
279        Self {
280            coefficients: scaled_coefficients,
281        }
282    }
283
284    /// Returns a vector of polynomials [p₀, p₁, ..., p_{d-1}], where d is `number_of_parts`, such that `self` equals
285    /// p₀(Xᵈ) + Xp₁(Xᵈ) + ... + X^(d-1)p_{d-1}(Xᵈ).
286    ///
287    /// Example: if d = 2 and `self` is 3 X^3 + X^2 + 2X + 1, then `poly.break_in_parts(2)`
288    /// returns a vector with two polynomials `(p₀, p₁)`, where p₀ = X + 1 and p₁ = 3X + 2.
289    pub fn break_in_parts(&self, number_of_parts: usize) -> Vec<Self> {
290        let coef = self.coefficients();
291        let mut parts: Vec<Self> = Vec::with_capacity(number_of_parts);
292        for i in 0..number_of_parts {
293            let coeffs: Vec<_> = coef
294                .iter()
295                .skip(i)
296                .step_by(number_of_parts)
297                .cloned()
298                .collect();
299            parts.push(Polynomial::new(&coeffs));
300        }
301        parts
302    }
303
304    pub fn to_extension<L: IsField>(self) -> Polynomial<FieldElement<L>>
305    where
306        F: IsSubFieldOf<L>,
307    {
308        Polynomial {
309            coefficients: self
310                .coefficients
311                .into_iter()
312                .map(|x| x.to_extension::<L>())
313                .collect(),
314        }
315    }
316}
317
318impl<F: IsPrimeField> Polynomial<FieldElement<F>> {
319    // Print the polynomial as a string ready to be used in SageMath, or just for pretty printing.
320    pub fn print_as_sage_poly(&self, var_name: Option<char>) -> String {
321        let var_name = var_name.unwrap_or('x');
322        if self.coefficients.is_empty()
323            || self.coefficients.len() == 1 && self.coefficients[0] == FieldElement::zero()
324        {
325            return String::new();
326        }
327
328        let mut string = String::new();
329        let zero = FieldElement::<F>::zero();
330
331        for (i, coeff) in self.coefficients.iter().rev().enumerate() {
332            if *coeff == zero {
333                continue;
334            }
335
336            let coeff_str = coeff.representative().to_string();
337
338            if i == self.coefficients.len() - 1 {
339                string.push_str(&coeff_str);
340            } else if i == self.coefficients.len() - 2 {
341                string.push_str(&format!("{}*{} + ", coeff_str, var_name));
342            } else {
343                string.push_str(&format!(
344                    "{}*{}^{} + ",
345                    coeff_str,
346                    var_name,
347                    self.coefficients.len() - 1 - i
348                ));
349            }
350        }
351
352        string
353    }
354}
355
356pub fn pad_with_zero_coefficients_to_length<F: IsField>(
357    pa: &mut Polynomial<FieldElement<F>>,
358    n: usize,
359) {
360    pa.coefficients.resize(n, FieldElement::zero());
361}
362
363/// Pads polynomial representations with minimum number of zeros to match lengths.
364pub fn pad_with_zero_coefficients<L: IsField, F: IsSubFieldOf<L>>(
365    pa: &Polynomial<FieldElement<F>>,
366    pb: &Polynomial<FieldElement<L>>,
367) -> (Polynomial<FieldElement<F>>, Polynomial<FieldElement<L>>) {
368    let mut pa = pa.clone();
369    let mut pb = pb.clone();
370
371    if pa.coefficients.len() > pb.coefficients.len() {
372        pad_with_zero_coefficients_to_length(&mut pb, pa.coefficients.len());
373    } else {
374        pad_with_zero_coefficients_to_length(&mut pa, pb.coefficients.len());
375    }
376    (pa, pb)
377}
378
379pub fn compose<F>(
380    poly_1: &Polynomial<FieldElement<F>>,
381    poly_2: &Polynomial<FieldElement<F>>,
382) -> Polynomial<FieldElement<F>>
383where
384    F: IsField,
385{
386    let max_degree: u64 = (poly_1.degree() * poly_2.degree()) as u64;
387
388    let mut interpolation_points = vec![];
389    for i in 0_u64..max_degree + 1 {
390        interpolation_points.push(FieldElement::<F>::from(i));
391    }
392
393    let values: Vec<_> = interpolation_points
394        .iter()
395        .map(|value| {
396            let intermediate_value = poly_2.evaluate(value);
397            poly_1.evaluate(&intermediate_value)
398        })
399        .collect();
400
401    Polynomial::interpolate(interpolation_points.as_slice(), values.as_slice())
402        .expect("xs and ys have equal length and xs are unique")
403}
404
405impl<F, L> ops::Add<&Polynomial<FieldElement<L>>> for &Polynomial<FieldElement<F>>
406where
407    L: IsField,
408    F: IsSubFieldOf<L>,
409{
410    type Output = Polynomial<FieldElement<L>>;
411
412    fn add(self, a_polynomial: &Polynomial<FieldElement<L>>) -> Self::Output {
413        let (pa, pb) = pad_with_zero_coefficients(self, a_polynomial);
414        let iter_coeff_pa = pa.coefficients.iter();
415        let iter_coeff_pb = pb.coefficients.iter();
416        let new_coefficients = iter_coeff_pa.zip(iter_coeff_pb).map(|(x, y)| x + y);
417        let new_coefficients_vec = new_coefficients.collect::<Vec<FieldElement<L>>>();
418        Polynomial::new(&new_coefficients_vec)
419    }
420}
421
422impl<F, L> ops::Add<Polynomial<FieldElement<L>>> for Polynomial<FieldElement<F>>
423where
424    L: IsField,
425    F: IsSubFieldOf<L>,
426{
427    type Output = Polynomial<FieldElement<L>>;
428
429    fn add(self, a_polynomial: Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
430        &self + &a_polynomial
431    }
432}
433
434impl<F, L> ops::Add<&Polynomial<FieldElement<L>>> for Polynomial<FieldElement<F>>
435where
436    L: IsField,
437    F: IsSubFieldOf<L>,
438{
439    type Output = Polynomial<FieldElement<L>>;
440
441    fn add(self, a_polynomial: &Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
442        &self + a_polynomial
443    }
444}
445
446impl<F, L> ops::Add<Polynomial<FieldElement<L>>> for &Polynomial<FieldElement<F>>
447where
448    L: IsField,
449    F: IsSubFieldOf<L>,
450{
451    type Output = Polynomial<FieldElement<L>>;
452
453    fn add(self, a_polynomial: Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
454        self + &a_polynomial
455    }
456}
457
458impl<F: IsField> ops::Neg for &Polynomial<FieldElement<F>> {
459    type Output = Polynomial<FieldElement<F>>;
460
461    fn neg(self) -> Polynomial<FieldElement<F>> {
462        let neg = self
463            .coefficients
464            .iter()
465            .map(|x| -x)
466            .collect::<Vec<FieldElement<F>>>();
467        Polynomial::new(&neg)
468    }
469}
470
471impl<F: IsField> ops::Neg for Polynomial<FieldElement<F>> {
472    type Output = Polynomial<FieldElement<F>>;
473
474    fn neg(self) -> Polynomial<FieldElement<F>> {
475        -&self
476    }
477}
478
479impl<F, L> ops::Sub<&Polynomial<FieldElement<L>>> for &Polynomial<FieldElement<F>>
480where
481    L: IsField,
482    F: IsSubFieldOf<L>,
483{
484    type Output = Polynomial<FieldElement<L>>;
485
486    fn sub(self, substrahend: &Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
487        self + (-substrahend)
488    }
489}
490
491impl<F, L> ops::Sub<Polynomial<FieldElement<L>>> for Polynomial<FieldElement<F>>
492where
493    L: IsField,
494    F: IsSubFieldOf<L>,
495{
496    type Output = Polynomial<FieldElement<L>>;
497
498    fn sub(self, substrahend: Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
499        &self - &substrahend
500    }
501}
502
503impl<F, L> ops::Sub<&Polynomial<FieldElement<L>>> for Polynomial<FieldElement<F>>
504where
505    L: IsField,
506    F: IsSubFieldOf<L>,
507{
508    type Output = Polynomial<FieldElement<L>>;
509
510    fn sub(self, substrahend: &Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
511        &self - substrahend
512    }
513}
514
515impl<F, L> ops::Sub<Polynomial<FieldElement<L>>> for &Polynomial<FieldElement<F>>
516where
517    L: IsField,
518    F: IsSubFieldOf<L>,
519{
520    type Output = Polynomial<FieldElement<L>>;
521
522    fn sub(self, substrahend: Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
523        self - &substrahend
524    }
525}
526
527impl<F> ops::Div<Polynomial<FieldElement<F>>> for Polynomial<FieldElement<F>>
528where
529    F: IsField,
530{
531    type Output = Polynomial<FieldElement<F>>;
532
533    fn div(self, dividend: Polynomial<FieldElement<F>>) -> Polynomial<FieldElement<F>> {
534        self.div_with_ref(&dividend)
535    }
536}
537
538impl<F: IsField> ops::Mul<&Polynomial<FieldElement<F>>> for &Polynomial<FieldElement<F>> {
539    type Output = Polynomial<FieldElement<F>>;
540    fn mul(self, factor: &Polynomial<FieldElement<F>>) -> Polynomial<FieldElement<F>> {
541        self.mul_with_ref(factor)
542    }
543}
544
545impl<F: IsField> ops::Mul<Polynomial<FieldElement<F>>> for Polynomial<FieldElement<F>> {
546    type Output = Polynomial<FieldElement<F>>;
547    fn mul(self, factor: Polynomial<FieldElement<F>>) -> Polynomial<FieldElement<F>> {
548        &self * &factor
549    }
550}
551
552impl<F: IsField> ops::Mul<Polynomial<FieldElement<F>>> for &Polynomial<FieldElement<F>> {
553    type Output = Polynomial<FieldElement<F>>;
554    fn mul(self, factor: Polynomial<FieldElement<F>>) -> Polynomial<FieldElement<F>> {
555        self * &factor
556    }
557}
558
559impl<F: IsField> ops::Mul<&Polynomial<FieldElement<F>>> for Polynomial<FieldElement<F>> {
560    type Output = Polynomial<FieldElement<F>>;
561    fn mul(self, factor: &Polynomial<FieldElement<F>>) -> Polynomial<FieldElement<F>> {
562        &self * factor
563    }
564}
565
566/* Operations between Polynomials and field elements */
567/* Multiplication field element at left */
568impl<F, L> ops::Mul<FieldElement<F>> for Polynomial<FieldElement<L>>
569where
570    L: IsField,
571    F: IsSubFieldOf<L>,
572{
573    type Output = Polynomial<FieldElement<L>>;
574
575    fn mul(self, multiplicand: FieldElement<F>) -> Polynomial<FieldElement<L>> {
576        let new_coefficients = self
577            .coefficients
578            .iter()
579            .map(|value| &multiplicand * value)
580            .collect();
581        Polynomial {
582            coefficients: new_coefficients,
583        }
584    }
585}
586
587impl<F, L> ops::Mul<&FieldElement<F>> for &Polynomial<FieldElement<L>>
588where
589    L: IsField,
590    F: IsSubFieldOf<L>,
591{
592    type Output = Polynomial<FieldElement<L>>;
593
594    fn mul(self, multiplicand: &FieldElement<F>) -> Polynomial<FieldElement<L>> {
595        self.clone() * multiplicand.clone()
596    }
597}
598
599impl<F, L> ops::Mul<FieldElement<F>> for &Polynomial<FieldElement<L>>
600where
601    L: IsField,
602    F: IsSubFieldOf<L>,
603{
604    type Output = Polynomial<FieldElement<L>>;
605
606    fn mul(self, multiplicand: FieldElement<F>) -> Polynomial<FieldElement<L>> {
607        self * &multiplicand
608    }
609}
610
611impl<F, L> ops::Mul<&FieldElement<F>> for Polynomial<FieldElement<L>>
612where
613    L: IsField,
614    F: IsSubFieldOf<L>,
615{
616    type Output = Polynomial<FieldElement<L>>;
617
618    fn mul(self, multiplicand: &FieldElement<F>) -> Polynomial<FieldElement<L>> {
619        &self * multiplicand
620    }
621}
622
623/* Multiplication field element at right */
624impl<F, L> ops::Mul<&Polynomial<FieldElement<L>>> for &FieldElement<F>
625where
626    L: IsField,
627    F: IsSubFieldOf<L>,
628{
629    type Output = Polynomial<FieldElement<L>>;
630
631    fn mul(self, multiplicand: &Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
632        multiplicand * self
633    }
634}
635
636impl<F, L> ops::Mul<Polynomial<FieldElement<L>>> for &FieldElement<F>
637where
638    L: IsField,
639    F: IsSubFieldOf<L>,
640{
641    type Output = Polynomial<FieldElement<L>>;
642
643    fn mul(self, multiplicand: Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
644        &multiplicand * self
645    }
646}
647
648impl<F, L> ops::Mul<&Polynomial<FieldElement<L>>> for FieldElement<F>
649where
650    L: IsField,
651    F: IsSubFieldOf<L>,
652{
653    type Output = Polynomial<FieldElement<L>>;
654
655    fn mul(self, multiplicand: &Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
656        multiplicand * self
657    }
658}
659
660impl<F, L> ops::Mul<Polynomial<FieldElement<L>>> for FieldElement<F>
661where
662    L: IsField,
663    F: IsSubFieldOf<L>,
664{
665    type Output = Polynomial<FieldElement<L>>;
666
667    fn mul(self, multiplicand: Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
668        &multiplicand * &self
669    }
670}
671
672/* Addition field element at left */
673impl<F, L> ops::Add<&FieldElement<F>> for &Polynomial<FieldElement<L>>
674where
675    L: IsField,
676    F: IsSubFieldOf<L>,
677{
678    type Output = Polynomial<FieldElement<L>>;
679
680    fn add(self, other: &FieldElement<F>) -> Polynomial<FieldElement<L>> {
681        Polynomial::new_monomial(other.clone(), 0) + self
682    }
683}
684
685impl<F, L> ops::Add<FieldElement<F>> for Polynomial<FieldElement<L>>
686where
687    L: IsField,
688    F: IsSubFieldOf<L>,
689{
690    type Output = Polynomial<FieldElement<L>>;
691
692    fn add(self, other: FieldElement<F>) -> Polynomial<FieldElement<L>> {
693        &self + &other
694    }
695}
696
697impl<F, L> ops::Add<FieldElement<F>> for &Polynomial<FieldElement<L>>
698where
699    L: IsField,
700    F: IsSubFieldOf<L>,
701{
702    type Output = Polynomial<FieldElement<L>>;
703
704    fn add(self, other: FieldElement<F>) -> Polynomial<FieldElement<L>> {
705        self + &other
706    }
707}
708
709impl<F, L> ops::Add<&FieldElement<F>> for Polynomial<FieldElement<L>>
710where
711    L: IsField,
712    F: IsSubFieldOf<L>,
713{
714    type Output = Polynomial<FieldElement<L>>;
715
716    fn add(self, other: &FieldElement<F>) -> Polynomial<FieldElement<L>> {
717        &self + other
718    }
719}
720
721/* Addition field element at right */
722impl<F, L> ops::Add<&Polynomial<FieldElement<L>>> for &FieldElement<F>
723where
724    L: IsField,
725    F: IsSubFieldOf<L>,
726{
727    type Output = Polynomial<FieldElement<L>>;
728
729    fn add(self, other: &Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
730        Polynomial::new_monomial(self.clone(), 0) + other
731    }
732}
733
734impl<F, L> ops::Add<Polynomial<FieldElement<L>>> for FieldElement<F>
735where
736    L: IsField,
737    F: IsSubFieldOf<L>,
738{
739    type Output = Polynomial<FieldElement<L>>;
740
741    fn add(self, other: Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
742        &self + &other
743    }
744}
745
746impl<F, L> ops::Add<Polynomial<FieldElement<L>>> for &FieldElement<F>
747where
748    L: IsField,
749    F: IsSubFieldOf<L>,
750{
751    type Output = Polynomial<FieldElement<L>>;
752
753    fn add(self, other: Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
754        self + &other
755    }
756}
757
758impl<F, L> ops::Add<&Polynomial<FieldElement<L>>> for FieldElement<F>
759where
760    L: IsField,
761    F: IsSubFieldOf<L>,
762{
763    type Output = Polynomial<FieldElement<L>>;
764
765    fn add(self, other: &Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
766        &self + other
767    }
768}
769
770/* Substraction field element at left */
771impl<F, L> ops::Sub<&FieldElement<F>> for &Polynomial<FieldElement<L>>
772where
773    L: IsField,
774    F: IsSubFieldOf<L>,
775{
776    type Output = Polynomial<FieldElement<L>>;
777
778    fn sub(self, other: &FieldElement<F>) -> Polynomial<FieldElement<L>> {
779        -Polynomial::new_monomial(other.clone(), 0) + self
780    }
781}
782
783impl<F, L> ops::Sub<FieldElement<F>> for Polynomial<FieldElement<L>>
784where
785    L: IsField,
786    F: IsSubFieldOf<L>,
787{
788    type Output = Polynomial<FieldElement<L>>;
789
790    fn sub(self, other: FieldElement<F>) -> Polynomial<FieldElement<L>> {
791        &self - &other
792    }
793}
794
795impl<F, L> ops::Sub<FieldElement<F>> for &Polynomial<FieldElement<L>>
796where
797    L: IsField,
798    F: IsSubFieldOf<L>,
799{
800    type Output = Polynomial<FieldElement<L>>;
801
802    fn sub(self, other: FieldElement<F>) -> Polynomial<FieldElement<L>> {
803        self - &other
804    }
805}
806
807impl<F, L> ops::Sub<&FieldElement<F>> for Polynomial<FieldElement<L>>
808where
809    L: IsField,
810    F: IsSubFieldOf<L>,
811{
812    type Output = Polynomial<FieldElement<L>>;
813
814    fn sub(self, other: &FieldElement<F>) -> Polynomial<FieldElement<L>> {
815        &self - other
816    }
817}
818
819/* Substraction field element at right */
820impl<F, L> ops::Sub<&Polynomial<FieldElement<L>>> for &FieldElement<F>
821where
822    L: IsField,
823    F: IsSubFieldOf<L>,
824{
825    type Output = Polynomial<FieldElement<L>>;
826
827    fn sub(self, other: &Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
828        Polynomial::new_monomial(self.clone(), 0) - other
829    }
830}
831
832impl<F, L> ops::Sub<Polynomial<FieldElement<L>>> for FieldElement<F>
833where
834    L: IsField,
835    F: IsSubFieldOf<L>,
836{
837    type Output = Polynomial<FieldElement<L>>;
838
839    fn sub(self, other: Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
840        &self - &other
841    }
842}
843
844impl<F, L> ops::Sub<Polynomial<FieldElement<L>>> for &FieldElement<F>
845where
846    L: IsField,
847    F: IsSubFieldOf<L>,
848{
849    type Output = Polynomial<FieldElement<L>>;
850
851    fn sub(self, other: Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
852        self - &other
853    }
854}
855
856impl<F, L> ops::Sub<&Polynomial<FieldElement<L>>> for FieldElement<F>
857where
858    L: IsField,
859    F: IsSubFieldOf<L>,
860{
861    type Output = Polynomial<FieldElement<L>>;
862
863    fn sub(self, other: &Polynomial<FieldElement<L>>) -> Polynomial<FieldElement<L>> {
864        &self - other
865    }
866}
867
868#[derive(Debug)]
869pub enum InterpolateError {
870    UnequalLengths(usize, usize),
871    NonUniqueXs,
872}
873
874impl Display for InterpolateError {
875    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
876        match self {
877            InterpolateError::UnequalLengths(x, y) => {
878                write!(f, "xs and ys must be the same length. Got: {x} != {y}")
879            }
880            InterpolateError::NonUniqueXs => write!(f, "xs values should be unique."),
881        }
882    }
883}
884
885#[cfg(feature = "std")]
886impl std::error::Error for InterpolateError {}
887
888#[cfg(test)]
889mod tests {
890    use crate::field::fields::u64_prime_field::U64PrimeField;
891
892    // Some of these tests work when the finite field has order greater than 2.
893    use super::*;
894    const ORDER: u64 = 23;
895    type F = U64PrimeField<ORDER>;
896    type FE = FieldElement<F>;
897
898    fn polynomial_a() -> Polynomial<FE> {
899        Polynomial::new(&[FE::new(1), FE::new(2), FE::new(3)])
900    }
901
902    fn polynomial_minus_a() -> Polynomial<FE> {
903        Polynomial::new(&[FE::new(ORDER - 1), FE::new(ORDER - 2), FE::new(ORDER - 3)])
904    }
905
906    fn polynomial_b() -> Polynomial<FE> {
907        Polynomial::new(&[FE::new(3), FE::new(4), FE::new(5)])
908    }
909
910    fn polynomial_a_plus_b() -> Polynomial<FE> {
911        Polynomial::new(&[FE::new(4), FE::new(6), FE::new(8)])
912    }
913
914    fn polynomial_b_minus_a() -> Polynomial<FE> {
915        Polynomial::new(&[FE::new(2), FE::new(2), FE::new(2)])
916    }
917
918    #[test]
919    fn adding_a_and_b_equals_a_plus_b() {
920        assert_eq!(polynomial_a() + polynomial_b(), polynomial_a_plus_b());
921    }
922
923    #[test]
924    fn adding_a_and_a_plus_b_does_not_equal_b() {
925        assert_ne!(polynomial_a() + polynomial_a_plus_b(), polynomial_b());
926    }
927
928    #[test]
929    fn add_5_to_0_is_5() {
930        let p1 = Polynomial::new(&[FE::new(5)]);
931        let p2 = Polynomial::new(&[FE::new(0)]);
932        assert_eq!(p1 + p2, Polynomial::new(&[FE::new(5)]));
933    }
934
935    #[test]
936    fn add_0_to_5_is_5() {
937        let p1 = Polynomial::new(&[FE::new(5)]);
938        let p2 = Polynomial::new(&[FE::new(0)]);
939        assert_eq!(p2 + p1, Polynomial::new(&[FE::new(5)]));
940    }
941
942    #[test]
943    fn negating_0_returns_0() {
944        let p1 = Polynomial::new(&[FE::new(0)]);
945        assert_eq!(-p1, Polynomial::new(&[FE::new(0)]));
946    }
947
948    #[test]
949    fn negating_a_is_equal_to_minus_a() {
950        assert_eq!(-polynomial_a(), polynomial_minus_a());
951    }
952
953    #[test]
954    fn negating_a_is_not_equal_to_a() {
955        assert_ne!(-polynomial_a(), polynomial_a());
956    }
957
958    #[test]
959    fn substracting_5_5_gives_0() {
960        let p1 = Polynomial::new(&[FE::new(5)]);
961        let p2 = Polynomial::new(&[FE::new(5)]);
962        let p3 = Polynomial::new(&[FE::new(0)]);
963        assert_eq!(p1 - p2, p3);
964    }
965
966    #[test]
967    fn substracting_b_and_a_equals_b_minus_a() {
968        assert_eq!(polynomial_b() - polynomial_a(), polynomial_b_minus_a());
969    }
970
971    #[test]
972    fn constructor_removes_zeros_at_the_end_of_polynomial() {
973        let p1 = Polynomial::new(&[FE::new(3), FE::new(4), FE::new(0)]);
974        assert_eq!(p1.coefficients, &[FE::new(3), FE::new(4)]);
975    }
976
977    #[test]
978    fn pad_with_zero_coefficients_returns_polynomials_with_zeros_until_matching_size() {
979        let p1 = Polynomial::new(&[FE::new(3), FE::new(4)]);
980        let p2 = Polynomial::new(&[FE::new(3)]);
981
982        assert_eq!(p2.coefficients, &[FE::new(3)]);
983        let (pp1, pp2) = pad_with_zero_coefficients(&p1, &p2);
984        assert_eq!(pp1, p1);
985        assert_eq!(pp2.coefficients, &[FE::new(3), FE::new(0)]);
986    }
987
988    #[test]
989    fn multiply_5_and_0_is_0() {
990        let p1 = Polynomial::new(&[FE::new(5)]);
991        let p2 = Polynomial::new(&[FE::new(0)]);
992        assert_eq!(p1 * p2, Polynomial::new(&[FE::new(0)]));
993    }
994
995    #[test]
996    fn multiply_0_and_x_is_0() {
997        let p1 = Polynomial::new(&[FE::new(0)]);
998        let p2 = Polynomial::new(&[FE::new(0), FE::new(1)]);
999        assert_eq!(p1 * p2, Polynomial::new(&[FE::new(0)]));
1000    }
1001
1002    #[test]
1003    fn multiply_2_by_3_is_6() {
1004        let p1 = Polynomial::new(&[FE::new(2)]);
1005        let p2 = Polynomial::new(&[FE::new(3)]);
1006        assert_eq!(p1 * p2, Polynomial::new(&[FE::new(6)]));
1007    }
1008
1009    #[test]
1010    fn multiply_2xx_3x_3_times_x_4() {
1011        let p1 = Polynomial::new(&[FE::new(3), FE::new(3), FE::new(2)]);
1012        let p2 = Polynomial::new(&[FE::new(4), FE::new(1)]);
1013        assert_eq!(
1014            p1 * p2,
1015            Polynomial::new(&[FE::new(12), FE::new(15), FE::new(11), FE::new(2)])
1016        );
1017    }
1018
1019    #[test]
1020    fn multiply_x_4_times_2xx_3x_3() {
1021        let p1 = Polynomial::new(&[FE::new(3), FE::new(3), FE::new(2)]);
1022        let p2 = Polynomial::new(&[FE::new(4), FE::new(1)]);
1023        assert_eq!(
1024            p2 * p1,
1025            Polynomial::new(&[FE::new(12), FE::new(15), FE::new(11), FE::new(2)])
1026        );
1027    }
1028
1029    #[test]
1030    fn division_works() {
1031        let p1 = Polynomial::new(&[FE::new(1), FE::new(3)]);
1032        let p2 = Polynomial::new(&[FE::new(1), FE::new(3)]);
1033        let p3 = p1.mul_with_ref(&p2);
1034        assert_eq!(p3 / p2, p1);
1035    }
1036
1037    #[test]
1038    fn division_by_zero_degree_polynomial_works() {
1039        let four = FE::new(4);
1040        let two = FE::new(2);
1041        let p1 = Polynomial::new(&[four, four]);
1042        let p2 = Polynomial::new(&[two]);
1043        assert_eq!(Polynomial::new(&[two, two]), p1 / p2);
1044    }
1045
1046    #[test]
1047    fn evaluate_constant_polynomial_returns_constant() {
1048        let three = FE::new(3);
1049        let p = Polynomial::new(&[three]);
1050        assert_eq!(p.evaluate(&FE::new(10)), three);
1051    }
1052
1053    #[test]
1054    fn evaluate_slice() {
1055        let three = FE::new(3);
1056        let p = Polynomial::new(&[three]);
1057        let ret = p.evaluate_slice(&[FE::new(10), FE::new(15)]);
1058        assert_eq!(ret, [three, three]);
1059    }
1060
1061    #[test]
1062    fn create_degree_0_new_monomial() {
1063        assert_eq!(
1064            Polynomial::new_monomial(FE::new(3), 0),
1065            Polynomial::new(&[FE::new(3)])
1066        );
1067    }
1068
1069    #[test]
1070    fn zero_poly_evals_0_in_3() {
1071        assert_eq!(
1072            Polynomial::new_monomial(FE::new(0), 0).evaluate(&FE::new(3)),
1073            FE::new(0)
1074        );
1075    }
1076
1077    #[test]
1078    fn evaluate_degree_1_new_monomial() {
1079        let two = FE::new(2);
1080        let four = FE::new(4);
1081        let p = Polynomial::new_monomial(two, 1);
1082        assert_eq!(p.evaluate(&two), four);
1083    }
1084
1085    #[test]
1086    fn evaluate_degree_2_monomyal() {
1087        let two = FE::new(2);
1088        let eight = FE::new(8);
1089        let p = Polynomial::new_monomial(two, 2);
1090        assert_eq!(p.evaluate(&two), eight);
1091    }
1092
1093    #[test]
1094    fn evaluate_3_term_polynomial() {
1095        let p = Polynomial::new(&[FE::new(3), -FE::new(2), FE::new(4)]);
1096        assert_eq!(p.evaluate(&FE::new(2)), FE::new(15));
1097    }
1098
1099    #[test]
1100    fn simple_interpolating_polynomial_by_hand_works() {
1101        let denominator = Polynomial::new(&[FE::new(1) * (FE::new(2) - FE::new(4)).inv().unwrap()]);
1102        let numerator = Polynomial::new(&[-FE::new(4), FE::new(1)]);
1103        let interpolating = numerator * denominator;
1104        assert_eq!(
1105            (FE::new(2) - FE::new(4)) * (FE::new(1) * (FE::new(2) - FE::new(4)).inv().unwrap()),
1106            FE::new(1)
1107        );
1108        assert_eq!(interpolating.evaluate(&FE::new(2)), FE::new(1));
1109        assert_eq!(interpolating.evaluate(&FE::new(4)), FE::new(0));
1110    }
1111
1112    #[test]
1113    fn interpolate_x_2_y_3() {
1114        let p = Polynomial::interpolate(&[FE::new(2)], &[FE::new(3)]).unwrap();
1115        assert_eq!(FE::new(3), p.evaluate(&FE::new(2)));
1116    }
1117
1118    #[test]
1119    fn interpolate_x_0_2_y_3_4() {
1120        let p =
1121            Polynomial::interpolate(&[FE::new(0), FE::new(2)], &[FE::new(3), FE::new(4)]).unwrap();
1122        assert_eq!(FE::new(3), p.evaluate(&FE::new(0)));
1123        assert_eq!(FE::new(4), p.evaluate(&FE::new(2)));
1124    }
1125
1126    #[test]
1127    fn interpolate_x_2_5_7_y_10_19_43() {
1128        let p = Polynomial::interpolate(
1129            &[FE::new(2), FE::new(5), FE::new(7)],
1130            &[FE::new(10), FE::new(19), FE::new(43)],
1131        )
1132        .unwrap();
1133
1134        assert_eq!(FE::new(10), p.evaluate(&FE::new(2)));
1135        assert_eq!(FE::new(19), p.evaluate(&FE::new(5)));
1136        assert_eq!(FE::new(43), p.evaluate(&FE::new(7)));
1137    }
1138
1139    #[test]
1140    fn interpolate_x_0_0_y_1_1() {
1141        let p =
1142            Polynomial::interpolate(&[FE::new(0), FE::new(1)], &[FE::new(0), FE::new(1)]).unwrap();
1143
1144        assert_eq!(FE::new(0), p.evaluate(&FE::new(0)));
1145        assert_eq!(FE::new(1), p.evaluate(&FE::new(1)));
1146    }
1147
1148    #[test]
1149    fn interpolate_x_0_y_0() {
1150        let p = Polynomial::interpolate(&[FE::new(0)], &[FE::new(0)]).unwrap();
1151        assert_eq!(FE::new(0), p.evaluate(&FE::new(0)));
1152    }
1153
1154    #[test]
1155    fn composition_works() {
1156        let p = Polynomial::new(&[FE::new(0), FE::new(2)]);
1157        let q = Polynomial::new(&[FE::new(0), FE::new(0), FE::new(1)]);
1158        assert_eq!(
1159            compose(&p, &q),
1160            Polynomial::new(&[FE::new(0), FE::new(0), FE::new(2)])
1161        );
1162    }
1163
1164    #[test]
1165    fn break_in_parts() {
1166        // p = 3 X^3 + X^2 + 2X + 1
1167        let p = Polynomial::new(&[FE::new(1), FE::new(2), FE::new(1), FE::new(3)]);
1168        let p0_expected = Polynomial::new(&[FE::new(1), FE::new(1)]);
1169        let p1_expected = Polynomial::new(&[FE::new(2), FE::new(3)]);
1170        let parts = p.break_in_parts(2);
1171        assert_eq!(parts.len(), 2);
1172        let p0 = &parts[0];
1173        let p1 = &parts[1];
1174        assert_eq!(p0, &p0_expected);
1175        assert_eq!(p1, &p1_expected);
1176    }
1177
1178    use alloc::format;
1179    use proptest::prelude::*;
1180    proptest! {
1181        #[test]
1182        fn ruffini_inplace_equals_division(p in any::<Vec<u64>>(), b in any::<u64>()) {
1183            let p: Vec<_> = p.into_iter().map(FE::from).collect();
1184            let mut p = Polynomial::new(&p);
1185            let b = FE::from(b);
1186
1187            let p_ref = p.clone();
1188            let m = Polynomial::new_monomial(FE::one(), 1) - b;
1189
1190            p.ruffini_division_inplace(&b);
1191            prop_assert_eq!(p, p_ref / m);
1192        }
1193    }
1194
1195    proptest! {
1196        #[test]
1197        fn ruffini_inplace_equals_ruffini(p in any::<Vec<u64>>(), b in any::<u64>()) {
1198            let p: Vec<_> = p.into_iter().map(FE::from).collect();
1199            let mut p = Polynomial::new(&p);
1200            let b = FE::from(b);
1201            let q = p.ruffini_division(&b);
1202            p.ruffini_division_inplace(&b);
1203            prop_assert_eq!(q, p);
1204        }
1205    }
1206    #[test]
1207    fn test_xgcd() {
1208        // Case 1: Simple polynomials
1209        let p1 = Polynomial::new(&[FE::new(1), FE::new(0), FE::new(1)]); // x^2 + 1
1210        let p2 = Polynomial::new(&[FE::new(1), FE::new(1)]); // x + 1
1211        let (a, b, g) = p1.xgcd(&p2);
1212        // Check that a * p1 + b * p2 = g
1213        let lhs = a.mul_with_ref(&p1) + b.mul_with_ref(&p2);
1214        assert_eq!(a, Polynomial::new(&[FE::new(12)]));
1215        assert_eq!(b, Polynomial::new(&[FE::new(12), FE::new(11)]));
1216        assert_eq!(lhs, g);
1217        assert_eq!(g, Polynomial::new(&[FE::new(1)]));
1218
1219        // x^2-1 :
1220        let p3 = Polynomial::new(&[FE::new(ORDER - 1), FE::new(0), FE::new(1)]);
1221        // x^3-x = x(x^2-1)
1222        let p4 = Polynomial::new(&[FE::new(0), FE::new(ORDER - 1), FE::new(0), FE::new(1)]);
1223        let (a, b, g) = p3.xgcd(&p4);
1224
1225        let lhs = a.mul_with_ref(&p3) + b.mul_with_ref(&p4);
1226        assert_eq!(a, Polynomial::new(&[FE::new(1)]));
1227        assert_eq!(b, Polynomial::zero());
1228        assert_eq!(lhs, g);
1229        assert_eq!(g, p3);
1230    }
1231
1232    #[test]
1233    fn test_differentiate() {
1234        // 3x^2 + 2x + 42
1235        let px = Polynomial::new(&[FE::new(42), FE::new(2), FE::new(3)]);
1236        // 6x + 2
1237        let dpdx = px.differentiate();
1238        assert_eq!(dpdx, Polynomial::new(&[FE::new(2), FE::new(6)]));
1239
1240        // 128
1241        let px = Polynomial::new(&[FE::new(128)]);
1242        // 0
1243        let dpdx = px.differentiate();
1244        assert_eq!(dpdx, Polynomial::new(&[FE::new(0)]));
1245    }
1246
1247    #[test]
1248    fn test_print_as_sage_poly() {
1249        let p = Polynomial::new(&[FE::new(1), FE::new(2), FE::new(3)]);
1250        assert_eq!(p.print_as_sage_poly(None), "3*x^2 + 2*x + 1");
1251    }
1252}
lambdaworks_math/polynomial/mod.rs

lambdaworks_math/polynomial/
mod.rs
