kzg 0.6.0

use core::borrow::Borrow;
use core::cmp::{Eq, PartialEq};
use core::ops::{Add, AddAssign, Mul, Sub, SubAssign};
use pairing::group::ff::PrimeField;
use std::iter::Iterator;

#[cfg(feature = "serde_support")]
use std::fmt;
#[cfg(feature = "serde_support")]
use serde::{Serialize, Deserialize, Serializer, Deserializer, ser::SerializeStruct, de::Unexpected, de::Visitor, de::Error as SerdeError, de::SeqAccess, de::MapAccess, self};
#[cfg(feature = "serde_support")]
use pairing::group::ff::PrimeFieldBits;


use crate::ft::EvaluationDomain;

const FFT_MUL_THRESHOLD: usize = 128;

#[derive(Clone, Debug)]
pub struct Polynomial<S: PrimeField> {
    pub degree: usize,
    pub coeffs: Vec<S>,
}

impl<S: PrimeField> PartialEq<Polynomial<S>> for Polynomial<S> {
    fn eq(&self, other: &Self) -> bool {
        if self.degree() != other.degree() {
            false
        } else {
            self.coeffs
                .iter()
                .zip(other.coeffs.iter())
                .all(|(l, r)| l == r)
        }
    }
}

impl<S: PrimeField> Eq for Polynomial<S> {}

impl<S: PrimeField> Polynomial<S> {
    pub fn is_zero(&self) -> bool {
        self.degree() == 0 && self.coeffs[0] == S::zero()
    }

    pub fn new_zero() -> Polynomial<S> {
        Polynomial {
            degree: 0,
            coeffs: vec![S::zero()],
        }
    }

    pub fn from_scalar(scalar: S) -> Polynomial<S> {
        Polynomial {
            degree: 0,
            coeffs: vec![scalar]
        }
    }

    pub fn new_monic_of_degree(degree: usize) -> Polynomial<S> {
        Polynomial {
            degree,
            coeffs: vec![S::one(); degree + 1]
        }
    }

    pub fn new_single_term(degree: usize) -> Polynomial<S> {
        let mut coeffs = vec![S::zero(); degree + 1];
        coeffs[degree] = S::one();
        Polynomial {
            degree,
            coeffs
        }
    }

    pub fn new_zero_with_size(cap: usize) -> Polynomial<S> {
        Polynomial {
            degree: 0,
            coeffs: vec![S::zero(); cap],
        }
    }

    pub fn new(coeffs: Vec<S>) -> Polynomial<S> {
        // figure out what the initial degree is
        let degree = Self::compute_degree(&coeffs, coeffs.len() - 1);
        Polynomial { degree, coeffs }
    }

    /// note: use this carefully, as setting the degree incorrect can lead to the degree being inconsistent
    pub fn new_from_coeffs(coeffs: Vec<S>, degree: usize) -> Polynomial<S> {
        Polynomial { degree, coeffs }
    }

    pub fn compute_degree(coeffs: &Vec<S>, upper_bound: usize) -> usize {
        let mut i = upper_bound;
        loop {
            if i == 0 {
                break 0;
            } else if coeffs[i] != S::zero() {
                break i;
            }

            i -= 1;
        }
    }

    pub fn truncate(&mut self, degree: usize) {
        self.degree = degree;
        self.coeffs.truncate(degree + 1);
    }

    pub fn reverse(&mut self) {
        self.coeffs.truncate(self.num_coeffs());
        self.coeffs.reverse();
    }

    pub fn shrink_degree(&mut self) {
        let degree = Self::compute_degree(&self.coeffs, self.degree);
        self.degree = degree;
    }

    pub fn fixup_degree(&mut self) {
        let degree = Self::compute_degree(&self.coeffs, self.coeffs.len() - 1);
        self.degree = degree;
    }

    pub fn lead(&self) -> S {
        self.coeffs[self.degree]
    }

    pub fn constant(&self) -> S {
        self.coeffs[0]
    }

    pub fn num_coeffs(&self) -> usize {
        self.degree + 1
    }

    pub fn degree(&self) -> usize {
        self.degree
    }

    pub fn iter_coeffs(&self) -> impl Iterator<Item = &S> {
        self.coeffs.iter().take(self.num_coeffs())
    }

    pub fn eval(&self, x: S) -> S {
        let mut res = self.coeffs[self.degree()];

        for i in (0..self.degree()).rev() {
            res *= x;
            res += self.coeffs[i];
        }

        res
    }

    pub fn fft_mul(&self, other: &Polynomial<S>) -> Polynomial<S> {
        let n = self.num_coeffs();
        let k = other.num_coeffs();
        let mut lhs = self.coeffs.clone();
        let mut rhs = other.coeffs.clone();
        lhs.resize(n + k, S::zero());
        rhs.resize(n + k, S::zero());

        let mut lhs = EvaluationDomain::from_coeffs(lhs).unwrap();
        let mut rhs = EvaluationDomain::from_coeffs(rhs).unwrap();

        lhs.fft();
        rhs.fft();
        lhs.mul_assign(&rhs);
        lhs.ifft();
        lhs.into()
    }

    pub fn best_mul(&self, other: &Polynomial<S>) -> Polynomial<S> {
        if self.degree() < FFT_MUL_THRESHOLD || other.degree() < FFT_MUL_THRESHOLD {
            self.clone() * other.clone()
        } else {
            self.fft_mul(&other)
        }
    }

    pub fn long_division(&self, divisor: &Self) -> (Polynomial<S>, Option<Polynomial<S>>) {
        if self.is_zero() {
            (Self::new_zero(), None)
        } else if divisor.is_zero() {
            panic!("divisor must not be zero!")
        } else if self.degree < divisor.degree() {
            (Self::new_zero(), Some(self.clone()))
        } else {
            let mut remainder = self.clone();
            let mut quotient = Polynomial::new_from_coeffs(
                vec![S::zero(); self.degree() - divisor.degree() + 1],
                self.degree() - divisor.degree(),
            );

            // inverse guaranteed to succeed because divisor isn't 0.
            let lead_inverse = divisor.lead().invert().unwrap();
            while !remainder.is_zero() && remainder.degree() >= divisor.degree() {
                let factor = remainder.lead() * lead_inverse;
                let i = remainder.degree() - divisor.degree();
                quotient.coeffs[i] = factor;

                for (j, &coeff) in divisor.iter_coeffs().enumerate() {
                    remainder.coeffs[i + j] -= coeff * factor;
                }

                remainder.shrink_degree();
            }

            if remainder.is_zero() {
                (quotient, None)
            } else {
                (quotient, Some(remainder))
            }
        }
    }

    pub fn fft_div(&self, divisor: &Self) -> (Polynomial<S>, Option<Polynomial<S>>) {
        let m = self.degree();
        let n = divisor.degree();

        let mut a_rev = self.clone();
        let mut b_rev = divisor.clone();

        a_rev.reverse();
        b_rev.reverse();

        let inv = b_rev.invert(m - n);
        
        let q_rev = a_rev.best_mul(&inv);
        let mut q = q_rev.clone();
        q.truncate(m - n);
        q.reverse();

        let r = self - &divisor.best_mul(&q);
        if r.is_zero() {
            (q, None)
        } else {
            (q, Some(r))
        }
    }

    /// computes the first degree + 1 terms of the formal series 1/f(x)
    /// panics if coeffs[0] == 0 or lead_coeff == 0
    pub fn invert(&self, degree: usize) -> Polynomial<S> {
        if degree == 0 {
            Polynomial::new_from_coeffs(vec![self.coeffs[0].invert().unwrap()], 0)
        } else {
            let c = self.invert(degree / 2);
            let mut res = c.best_mul(&Polynomial::from_scalar(2.into()).sub(&c.best_mul(self)));
            res.truncate(degree);
            res
        }
    }

    pub fn multi_eval(&self, xs: &[S]) -> Vec<S> {
        assert!(xs.len() > self.degree());
        let tree = SubProductTree::new_from_points(xs);
        tree.eval(xs.as_ref(), self)
    }

    /// Performs lagrange interpolation on a pre-computed sub-product tree of xs.
    /// `tree` must be the same as the result when calling SubProductTree::new_from_points(xs)
    pub fn lagrange_interpolation_with_tree(xs: &[S], ys: &[S], tree: &SubProductTree<S>) -> Polynomial<S> {
        assert_eq!(xs.len(), ys.len());

        if xs.len() == 1 {
            let coeffs = vec![ys[0] - xs[0], S::one()];
            return Polynomial::new_from_coeffs(coeffs, 1);
        }

        let mut m_prime = tree.product.clone();
        for i in 1..m_prime.num_coeffs() {
            m_prime.coeffs[i] *= S::from(i as u64);
        }
        m_prime.coeffs.remove(0);
        m_prime.degree -= 1;


        let cs: Vec<S> = m_prime.multi_eval(xs).iter().enumerate().map(|(i, c)| ys[i] * c.invert().unwrap()).collect();

        tree.linear_mod_combination(cs.as_slice())
    }

    pub fn lagrange_interpolation(xs: &[S], ys: &[S]) -> Polynomial<S> {
        assert_eq!(xs.len(), ys.len());

        if xs.len() == 1 {
            let coeffs = vec![ys[0] - xs[0], S::one()];
            return Polynomial::new_from_coeffs(coeffs, 1);
        }

        // let xs = pad_to_power_of_two(xs);
        // let ys = pad_to_power_of_two(ys);
        let tree = SubProductTree::new_from_points(xs);

        let mut m_prime = tree.product.clone();
        for i in 1..m_prime.num_coeffs() {
            m_prime.coeffs[i] *= S::from(i as u64);
        }
        m_prime.coeffs.remove(0);
        m_prime.degree -= 1;


        let cs: Vec<S> = m_prime.multi_eval(xs).iter().enumerate().map(|(i, c)| ys[i] * c.invert().unwrap()).collect();

        tree.linear_mod_combination(cs.as_slice())
    }

    pub fn scalar_multiplication(mut self, rhs: S) -> Polynomial<S> {
        for i in 0..self.num_coeffs() {
            self.coeffs[i] *= rhs;
        }
        self
    }
}

pub struct SubProductTree<S: PrimeField> {
    pub product: Polynomial<S>,
    pub left: Option<Box<SubProductTree<S>>>,
    pub right: Option<Box<SubProductTree<S>>>
}

impl<S: PrimeField> SubProductTree<S> {
    pub fn new_from_points(xs: &[S]) -> SubProductTree<S> {
        match xs.len() {
            1 => SubProductTree {
                product: Polynomial::new_from_coeffs(vec![-xs[0], S::one()], 1),
                left: None,
                right: None
            },
            n => {
                let left = SubProductTree::new_from_points(&xs[..n / 2]);
                let right = SubProductTree::new_from_points(&xs[n / 2..]);
                SubProductTree {
                    product: left.product.best_mul(&right.product),
                    left: Some(Box::new(left)),
                    right: Some(Box::new(right))
                }
            }
        }
    }

    pub fn eval(&self, xs: &[S], f: &Polynomial<S>) -> Vec<S> {
        let n = xs.len();

        if n == 1 {
            let y = f.eval(xs[0]);
            vec![y]
        } else {

            let left = self.left.as_ref().unwrap();
            let right = self.right.as_ref().unwrap();

            let (_, r0) = f.long_division(&left.product);
            let (_, r1) = f.long_division(&right.product);

            let mut l0 = left.eval(&xs[..n/2], &r0.unwrap());
            let l1 = right.eval(&xs[n/2..], &r1.unwrap());

            l0.extend(l1);
            l0
        }
    }

    pub fn linear_mod_combination(&self, cs: &[S]) -> Polynomial<S> {
        let n = cs.len();

        if n == 1 {
            Polynomial::new_from_coeffs(vec![cs[0]], 0)
        } else {
            let left = self.left.as_ref().unwrap();
            let right = self.right.as_ref().unwrap();

            let l = left.linear_mod_combination(&cs[..n/2]);
            let r = right.linear_mod_combination(&cs[n/2..]);
            
            right.product.best_mul(&l) + left.product.best_mul(&r)
        }
    }
}

fn op_tree_inner<T, F, O>(left: usize, size: usize, get_elem: &F, op: &O) -> T
where
    F: Fn(usize) -> T,
    O: Fn(T, T) -> T,
{
    assert!(size > 0);
    if size == 1 {
        get_elem(left)
    } else if size == 2 {
        op(get_elem(left), get_elem(left + 1))
    } else {
        let mid = left + (size / 2);
        op(
            op_tree_inner(left, size / 2, get_elem, op),
            op_tree_inner(mid, size - (size / 2), get_elem, op),
        )
    }
}

pub fn op_tree<T, F, O>(size: usize, get_elem: &F, op: &O) -> T
where
    F: Fn(usize) -> T,
    O: Fn(T, T) -> T,
{
    op_tree_inner(0, size, get_elem, op)
}

impl<'a, S: PrimeField> Add for &'a Polynomial<S> {
    type Output = Polynomial<S>;

    fn add(self, rhs: Self) -> Self::Output {
        let (mut res, shorter) = if rhs.degree() > self.degree {
            (rhs.clone(), self)
        } else {
            (self.clone(), rhs)
        };

        for i in 0..shorter.degree() {
            res.coeffs[i] += shorter.coeffs[i];
        }

        res
    }
}

impl<S: PrimeField> Add for Polynomial<S> {
    type Output = Polynomial<S>;

    fn add(self, rhs: Self) -> Self::Output {
        let (mut res, shorter) = if rhs.degree() > self.degree() {
            (rhs, self)
        } else {
            (self, rhs)
        };

        for i in 0..shorter.num_coeffs() {
            res.coeffs[i] += shorter.coeffs[i];
        }

        res
    }
}

impl<S: PrimeField, R: Borrow<Polynomial<S>>> AddAssign<R> for Polynomial<S> {
    fn add_assign(&mut self, rhs: R) {
        let rhs = rhs.borrow();
        for i in 0..rhs.num_coeffs() {
            self.coeffs[i] += rhs.coeffs[i];
        }

        if self.degree() < rhs.degree() {
            self.degree = rhs.degree();
        }
    }
}

impl<'a, S: PrimeField> Sub for &'a Polynomial<S> {
    type Output = Polynomial<S>;

    fn sub(self, rhs: Self) -> Self::Output {
        let mut res = self.clone();
        if rhs.num_coeffs() > self.num_coeffs() {
            res.coeffs.resize(rhs.num_coeffs(), S::zero());
            res.degree = rhs.degree();
        }

        for i in 0..rhs.num_coeffs() {
            res.coeffs[i] -= rhs.coeffs[i];
        }

        res.shrink_degree();
        res
    }
}

impl<S: PrimeField, R: Borrow<Polynomial<S>>> SubAssign<R> for Polynomial<S> {
    fn sub_assign(&mut self, rhs: R) {
        let rhs = rhs.borrow();
        for i in 0..rhs.num_coeffs() {
            self.coeffs[i] -= rhs.coeffs[i];
        }

        self.fixup_degree()
    }
}

impl<S: PrimeField> Mul<Polynomial<S>> for Polynomial<S> {
    type Output = Polynomial<S>;

    fn mul(self, rhs: Self) -> Self::Output {
        let mut res = Polynomial::new_zero_with_size(self.degree() + rhs.degree() + 1);
        for i in 0..self.num_coeffs() {
            for j in 0..rhs.num_coeffs() {
                res.coeffs[i + j] += self.coeffs[i] * rhs.coeffs[j];
            }
        }

        res.degree = self.degree() + rhs.degree();
        res
    }
}

#[cfg(all(feature = "serde_support", any(feature = "b12_381")))]
#[derive(Debug, Clone)]
pub struct SerializablePolynomial<S: PrimeFieldBits> {
    degree: usize,
    coeffs: Vec<SerializablePrimeField<S>>,
}

#[cfg(all(feature = "serde_support", any(feature = "b12_381")))]
#[derive(Debug, Clone)]
pub struct SerializablePrimeField<S: PrimeField>(S);

#[cfg(all(feature = "serde_support", any(feature = "b12_381")))]
impl<S: PrimeField> From<S> for SerializablePrimeField<S> {
    fn from(s: S) -> SerializablePrimeField<S> {
        SerializablePrimeField(s)
    }
}

#[cfg(all(feature = "serde_support", any(feature = "b12_381")))]
impl<S: PrimeField> SerializablePrimeField<S> {
    fn into_inner(self) -> S {
        self.0
    }
}

#[cfg(all(feature = "serde_support", any(feature = "b12_381")))]
impl<F: PrimeField> Serialize for SerializablePrimeField<F> {
    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
    where
        S: Serializer
    {
        serializer.serialize_bytes(self.0.to_repr().as_ref())
    }
}


#[cfg(all(feature = "serde_support", feature = "b12_381"))]
use bls12_381::Scalar;

#[cfg(all(feature = "serde_support", feature = "b12_381"))]
impl<'de> Deserialize<'de> for SerializablePrimeField<Scalar> {
    fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>
    where
        D: Deserializer<'de>
    {
        struct PrimeFieldVisitor;

        impl<'de> Visitor<'de> for PrimeFieldVisitor {
            type Value = Scalar;

            fn expecting(&self, formatter: &mut fmt::Formatter) -> fmt::Result {
                formatter.write_str("canonical byte representation of a prime field element")
            }

            fn visit_bytes<E>(self, v: &[u8]) -> Result<Self::Value, E>
            where
                E: SerdeError
            {
                if v.len() != 32 {
                    return Err(SerdeError::invalid_value(Unexpected::Bytes(v), &self))
                }

                let mut encoding = [0; 32];
                encoding.as_mut().copy_from_slice(&v[0..32]);
                let s = Scalar::from_bytes(&encoding);
                if s.is_none().into() {
                    Err(SerdeError::invalid_value(Unexpected::Bytes(&encoding), &self))
                } else {
                    Ok(s.unwrap())
                }
            }
        }

        let inner = deserializer.deserialize_bytes(PrimeFieldVisitor)?;
        Ok(inner.into())
    }
}

#[cfg(all(feature = "serde_support", feature = "b12_381"))]
impl From<Polynomial<Scalar>> for SerializablePolynomial<Scalar> {
    fn from(inner: Polynomial<Scalar>) -> Self {
        SerializablePolynomial {
            degree: inner.degree,
            // TODO: is there a way to do this without actually iterating?
            // TODO: since 'x' is a newtype will this whole loop get compiled away?
            coeffs: inner.coeffs.into_iter().map(|x| x.into()).collect()
        }
    }
}

#[cfg(all(feature = "serde_support", feature = "b12_381"))]
impl From<SerializablePolynomial<Scalar>> for Polynomial<Scalar> {
    fn from(inner: SerializablePolynomial<Scalar>) -> Self {
        Polynomial {
            degree: inner.degree,
            // TODO: is there a way to do this without actually iterating?
            // TODO: since 'x' is a newtype will this whole loop get compiled away?
            coeffs: inner.coeffs.into_iter().map(|x| x.into_inner()).collect()
        }
    }
}

#[cfg(all(feature = "serde_support", feature = "b12_381"))]
impl SerializablePolynomial<Scalar> {
    pub fn from_inner_ref(inner: &Polynomial<Scalar>) -> Self {
        SerializablePolynomial {
            degree: inner.degree,
            // TODO: is there a way to do this without actually iterating?
            // TODO: since 'x' is a newtype will this whole loop get compiled away?
            coeffs: inner.coeffs.iter().cloned().map(|x| x.into()).collect()
        }
    }

    pub fn to_inner_ref(&self) -> Polynomial<Scalar> {
        Polynomial {
            degree: self.degree,
            coeffs: self.coeffs.iter().cloned().map(|x| x.into_inner()).collect()
        }
    }
}

#[cfg(all(feature = "serde_support", feature = "b12_381"))]
impl Serialize for SerializablePolynomial<Scalar> {
    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
    where
        S: Serializer,
    {
        let mut state = serializer.serialize_struct("SerializablePolynomial<bls12_381::Scalar>", 2)?;
        state.serialize_field("degree", &self.degree)?;
        state.serialize_field("coeffs", &self.coeffs)?;
        state.end()
    }
}

#[cfg(all(feature = "serde_support", feature = "b12_381"))]
impl<'de> Deserialize<'de> for SerializablePolynomial<Scalar> {
    fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>
    where
        D: Deserializer<'de>
    {
        #[derive(Deserialize)]
        #[serde(field_identifier, rename_all = "lowercase")]
        enum Field { Degree, Coeffs }

        struct SerializablePolynomialVisitor;

        impl<'de> Visitor<'de> for SerializablePolynomialVisitor {
            type Value = SerializablePolynomial<Scalar>;


            fn expecting(&self, formatter: &mut fmt::Formatter) -> fmt::Result {
                formatter.write_str("SerializablePolynomial<bls12_381::Scalar>")
            }

            fn visit_seq<V>(self, mut seq: V) -> Result<SerializablePolynomial<Scalar>, V::Error>
            where
                V: SeqAccess<'de>,
            {
                let degree = seq.next_element()?
                    .ok_or_else(|| SerdeError::invalid_length(0, &self))?;
                let coeffs = seq.next_element()?
                    .ok_or_else(|| SerdeError::invalid_length(1, &self))?;
                Ok(SerializablePolynomial { degree, coeffs })
            }

            fn visit_map<V>(self, mut map: V) -> Result<SerializablePolynomial<Scalar>, V::Error>
            where
                V: MapAccess<'de>,
            {
                let mut degree = None;
                let mut coeffs = None;
                while let Some(key) = map.next_key()? {
                    match key {
                        Field::Degree => {
                            if degree.is_some() {
                                return Err(SerdeError::duplicate_field("degree"));
                            }
                            degree = Some(map.next_value()?);
                        }
                        Field::Coeffs => {
                            if coeffs.is_some() {
                                return Err(SerdeError::duplicate_field("coeffs"));
                            }
                            coeffs = Some(map.next_value()?);
                        }
                    }
                }
                let degree = degree.ok_or_else(|| SerdeError::missing_field("secs"))?;
                let coeffs = coeffs.ok_or_else(|| SerdeError::missing_field("nanos"))?;
                Ok(SerializablePolynomial { degree, coeffs })
            }

        }

        const FIELDS: &'static [&'static str] = &["degree", "coeffs"];
        deserializer.deserialize_struct("SerializablePolynomial<bls12_381::Scalar>", FIELDS, SerializablePolynomialVisitor)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use bls12_381::{Bls12, Scalar};

    #[test]
    fn test_long_division() {
        // test cases taken from https://tutorial.math.lamar.edu/Solutions/Alg/DividingPolynomials

        // 3x^4 - 5x^2 + 3 / x + 2 = 3x^3 - 6x^2 + 7x - 14 r 31
        let x = Polynomial::new(vec![
            3.into(),
            Scalar::zero(),
            -Scalar::from(5),
            Scalar::zero(),
            3.into(),
        ]);
        let y = Polynomial::new(vec![
            2.into(),
            Scalar::one(),
            Scalar::zero(),
            Scalar::zero(),
            Scalar::zero(),
        ]);

        let (q, r) = x.long_division(&y);
        assert!(r.is_some());
        assert_eq!(
            r.unwrap(),
            Polynomial::new(vec![
                31.into(),
                Scalar::zero(),
                Scalar::zero(),
                Scalar::zero(),
                Scalar::zero()
            ])
        );
        assert_eq!(
            q,
            Polynomial::new(vec![
                -Scalar::from(14),
                7.into(),
                -Scalar::from(6),
                3.into(),
                Scalar::zero(),
            ])
        );

        // x^3 + 2x^2 - 3x + 4 / x - 7 = x^2 + 9x + 60 r 424
        let x = Polynomial::new(vec![4.into(), -Scalar::from(3), 2.into(), Scalar::one()]);
        let y = Polynomial::new(vec![
            -Scalar::from(7),
            Scalar::one(),
            Scalar::zero(),
            Scalar::zero(),
        ]);

        let (q, r) = x.long_division(&y);
        assert!(r.is_some());
        assert_eq!(
            r.unwrap(),
            Polynomial::new(vec![
                424.into(),
                Scalar::zero(),
                Scalar::zero(),
                Scalar::zero(),
            ])
        );
        assert_eq!(
            q,
            Polynomial::new(vec![60.into(), 9.into(), Scalar::one(), Scalar::zero(),])
        );

        // x^3 + 6x^2 + 13x + 10 / x + 2 = x^2 + 4x + 5 r 0
        let x = Polynomial::new(vec![10.into(), 13.into(), 6.into(), Scalar::one()]);
        let y = Polynomial::new(vec![
            Scalar::from(2),
            Scalar::one(),
            Scalar::zero(),
            Scalar::zero(),
        ]);

        let (q, r) = x.long_division(&y);
        assert!(r.is_none());
        assert_eq!(
            q,
            Polynomial::new(vec![5.into(), 4.into(), Scalar::one(), Scalar::zero(),])
        );
    }

    #[test]
    fn test_fft_division() {
        // test cases taken from https://tutorial.math.lamar.edu/Solutions/Alg/DividingPolynomials

        // 3x^4 - 5x^2 + 3 / x + 2 = 3x^3 - 6x^2 + 7x - 14 r 31
        let x = Polynomial::new(vec![
            3.into(),
            Scalar::zero(),
            -Scalar::from(5),
            Scalar::zero(),
            3.into(),
        ]);
        let y = Polynomial::new(vec![
            2.into(),
            Scalar::one(),
            Scalar::zero(),
            Scalar::zero(),
            Scalar::zero(),
        ]);

        let (q, r) = x.fft_div(&y);
        assert!(r.is_some());
        assert_eq!(
            r.unwrap(),
            Polynomial::new(vec![
                31.into(),
                Scalar::zero(),
                Scalar::zero(),
                Scalar::zero(),
                Scalar::zero()
            ])
        );
        assert_eq!(
            q,
            Polynomial::new(vec![
                -Scalar::from(14),
                7.into(),
                -Scalar::from(6),
                3.into(),
                Scalar::zero(),
            ])
        );

        // x^3 + 2x^2 - 3x + 4 / x - 7 = x^2 + 9x + 60 r 424
        let x = Polynomial::new(vec![4.into(), -Scalar::from(3), 2.into(), Scalar::one()]);
        let y = Polynomial::new(vec![
            -Scalar::from(7),
            Scalar::one(),
            Scalar::zero(),
            Scalar::zero(),
        ]);

        let (q, r) = x.fft_div(&y);
        assert!(r.is_some());
        assert_eq!(
            r.unwrap(),
            Polynomial::new(vec![
                424.into(),
                Scalar::zero(),
                Scalar::zero(),
                Scalar::zero(),
            ])
        );
        assert_eq!(
            q,
            Polynomial::new(vec![60.into(), 9.into(), Scalar::one(), Scalar::zero(),])
        );

        // x^3 + 6x^2 + 13x + 10 / x + 2 = x^2 + 4x + 5 r 0
        let x = Polynomial::new(vec![10.into(), 13.into(), 6.into(), Scalar::one()]);
        let y = Polynomial::new(vec![
            Scalar::from(2),
            Scalar::one(),
            Scalar::zero(),
            Scalar::zero(),
        ]);

        let (q, r) = x.fft_div(&y);
        assert!(r.is_none());
        assert_eq!(
            q,
            Polynomial::new(vec![5.into(), 4.into(), Scalar::one(), Scalar::zero(),])
        );
    }

    #[test]
    fn test_eval_basic() {
        // y(x) = x^5 + 4x^3 + 7x^2 + 34
        let polynomial = Polynomial::new(vec![
            34.into(),
            Scalar::zero(),
            7.into(),
            4.into(),
            Scalar::zero(),
            Scalar::one(),
        ]);

        // y(0) = 34
        assert_eq!(polynomial.eval(Scalar::zero()), 34.into());
        // y(1) = 46
        assert_eq!(polynomial.eval(Scalar::one()), 46.into());
        // y(5) = 3834
        assert_eq!(polynomial.eval(5.into()), 3834.into());
    }

    fn verify_tree(tree: &SubProductTree<Scalar>) {
        if tree.left.is_some() && tree.right.is_some() {
            assert!(
                tree.product == tree.left.as_ref().unwrap().product.best_mul(&tree.right.as_ref().unwrap().product)
            );
        }
    }

    #[test]
    fn test_new_subproduct_tree() {
        let xs = [Scalar::from(2), Scalar::from(5), Scalar::from(7), Scalar::from(90), Scalar::from(111), Scalar::from(31), Scalar::from(29)];

        let tree = SubProductTree::new_from_points(&xs);
        verify_tree(&tree);

        let xs = [Scalar::from(2), Scalar::from(5), Scalar::from(7), Scalar::from(90), Scalar::from(111)];
        let tree = SubProductTree::new_from_points(&xs);
        verify_tree(&tree);
    }

    #[test]
    fn test_fast_multi_eval() {
        let polynomial: Polynomial<Scalar> = Polynomial::new(
            vec![2, 5, 7, 90, 111]
                .into_iter()
                .map(|x| x.into())
                .collect(),
        );

        let xs: Vec<Scalar> = vec![1, 2, 3, 4, 5, 6, 7, 8].into_iter().map(|x| x.into()).collect();
        
        let mut fast = polynomial.multi_eval(xs.as_slice());
        fast.truncate(xs.len());
        let mut slow = Vec::new();

        for i in 0..xs.len() {
           let slow_y = polynomial.eval(xs[i]);
           slow.push(slow_y);
        }

        let slow: Vec<Scalar> =  xs.iter().map(|x| polynomial.eval(*x)).collect();
        assert!(fast == slow);
    }

    #[test]
    fn test_interpolation() {
        let xs: Vec<Scalar> = vec![2].into_iter().map(|x| x.into()).collect();
        let ys: Vec<Scalar> = vec![8].into_iter().map(|x| x.into()).collect();

        let interpolation = Polynomial::lagrange_interpolation(xs.as_slice(), ys.as_slice());

        for (&x, &y) in xs.iter().zip(ys.iter()) {
            assert_eq!(interpolation.eval(x), y);
        }

        let xs: Vec<Scalar> = vec![2, 5, 7, 90, 111, 31, 29]
            .into_iter()
            .map(|x| x.into())
            .collect();
        let ys: Vec<Scalar> = vec![8, 1, 43, 2, 87, 122, 13]
            .into_iter()
            .map(|x| x.into())
            .collect();
        let interpolation = Polynomial::lagrange_interpolation(xs.as_slice(), ys.as_slice());

        for (&x, &y) in xs.iter().zip(ys.iter()) {
            assert_eq!(interpolation.eval(x), y);
        }
    }

    #[cfg(feature = "serde_support")]
    use bincode::{serialize, deserialize};

    #[cfg(all(feature = "serde_support", feature = "b12_381"))]
    #[test]
    fn test_polynomial_serialization() {
        let f = Polynomial::new(vec![
            3.into(),
            -Scalar::from(9)
            -Scalar::from(5),
            120.into(),
            4.into(),
        ]);
        
        let serable = SerializablePolynomial::from_inner_ref(&f);
        let ser = serialize(&serable).unwrap();
        let de: SerializablePolynomial<Scalar> = deserialize(ser.as_slice()).unwrap();
        let f_de: Polynomial<Scalar> = de.into();
        assert_eq!(f, f_de);
    }
}