1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149

// TODO: Naming?
#![allow(clippy::module_name_repetitions)]
#[cfg(feature = "std")]
use rayon::prelude::*;
use std::prelude::v1::*;
use zkp_mmap_vec::MmapVec;
#[cfg(feature = "std")]
use zkp_primefield::fft::{fft_cofactor_permuted_out, permute_index};
use zkp_primefield::FieldElement;

#[derive(PartialEq, Clone)]
pub struct DensePolynomial(MmapVec<FieldElement>);

// We normally don't want to spill thousands of coefficients in the logs.
#[cfg(feature = "std")]
impl std::fmt::Debug for DensePolynomial {
    fn fmt(&self, fmt: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        write!(fmt, "DensePolynomial(degree = {:?})", self.degree())
    }
}

impl DensePolynomial {
    pub fn from_mmap_vec(coefficients: MmapVec<FieldElement>) -> Self {
        assert!(coefficients.len().is_power_of_two());
        Self(coefficients)
    }

    // Coefficients are in order of ascending degree. E.g. &[1, 2] corresponds to
    // the polynomial f(x) = 1 + 2x.
    pub fn new(coefficients: &[FieldElement]) -> Self {
        assert!(coefficients.len().is_power_of_two());
        let mut vec = MmapVec::with_capacity(coefficients.len());
        vec.extend_from_slice(coefficients);
        Self(vec)
    }

    pub fn zeros(size: usize) -> Self {
        assert!(size.is_power_of_two());
        let mut vec = MmapVec::with_capacity(size);
        vec.resize(size, FieldElement::ZERO);
        Self(vec)
    }

    // Note that the length of a polynomial is not its degree, because the leading
    // coefficient of a DensePolynomial can be zero.
    pub fn len(&self) -> usize {
        self.0.len()
    }

    pub fn is_empty(&self) -> bool {
        self.0.is_empty()
    }

    pub fn coefficients(&self) -> &[FieldElement] {
        &self.0
    }

    // TODO: The zero polynomial is assigned a degree of 0, but it is
    // more correctly left undefined or sometimes assigned `-1` or `-∞`.
    pub fn degree(&self) -> usize {
        let mut degree = self.len() - 1;
        while self.0[degree] == FieldElement::ZERO && degree > 0 {
            degree -= 1;
        }
        degree
    }

    pub fn evaluate(&self, x: &FieldElement) -> FieldElement {
        let mut result = FieldElement::ZERO;
        for coefficient in self.0.iter().rev() {
            result *= x;
            result += coefficient;
        }
        result
    }

    #[cfg(feature = "std")]
    pub fn low_degree_extension(&self, blowup: usize) -> MmapVec<FieldElement> {
        // TODO: shift polynomial by FieldElement::GENERATOR outside of this function.
        const SHIFT_FACTOR: FieldElement = FieldElement::GENERATOR;
        let length = self.len() * blowup;
        let generator =
            FieldElement::root(length).expect("No generator for extended_domain_length.");

        // FieldElement is safe to initialize zero (which maps to zero)
        #[allow(unsafe_code)]
        let mut result: MmapVec<FieldElement> = unsafe { MmapVec::zero_initialized(length) };

        // Compute cosets in parallel
        result
            .as_mut_slice()
            .par_chunks_mut(self.len())
            .enumerate()
            .for_each(|(i, slice)| {
                let cofactor = &SHIFT_FACTOR * generator.pow(permute_index(blowup, i));
                fft_cofactor_permuted_out(&cofactor, &self.coefficients(), slice);
            });
        result
    }

    /// Divide out a point and add the scaled result to target.
    ///
    /// target += c * (P(X) - P(z)) / (X - z)
    /// See: https://en.wikipedia.org/wiki/Synthetic_division
    pub fn divide_out_point_into(&self, z: &FieldElement, c: &FieldElement, target: &mut Self) {
        let mut remainder = FieldElement::ZERO;
        for (coefficient, target) in self.0.iter().rev().zip(target.0.iter_mut().rev()) {
            *target += c * &remainder;
            remainder *= z;
            remainder += coefficient;
        }
    }
}

#[cfg(test)]
use quickcheck::{Arbitrary, Gen};
#[cfg(test)]
impl Arbitrary for DensePolynomial {
    fn arbitrary<G: Gen>(g: &mut G) -> Self {
        let mut coefficients = Vec::<FieldElement>::arbitrary(g);
        let length = coefficients.len();
        coefficients.extend_from_slice(&vec![
            FieldElement::ZERO;
            length.next_power_of_two() - length
        ]);
        assert!(coefficients.len().is_power_of_two());
        Self::new(&coefficients)
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    fn dense_polynomial(coefficients: &[isize]) -> DensePolynomial {
        DensePolynomial::new(
            &coefficients
                .iter()
                .map(|c| FieldElement::from(*c))
                .collect::<Vec<_>>(),
        )
    }

    #[test]
    fn example_evaluate() {
        let p = dense_polynomial(&[1, 0, 0, 2]);
        assert_eq!(p.evaluate(&FieldElement::from(2)), FieldElement::from(17));
    }
}