Upstream GBP, divisor, circuit abstraction, and EC gadgets from FCMP++

crypto/evrf/divisors/src/lib.rs

@@ -0,0 +1,247 @@
+#![cfg_attr(docsrs, feature(doc_auto_cfg))]
+#![doc = include_str!("../README.md")]
+#![deny(missing_docs)]
+#![allow(non_snake_case)]
+
+use group::{
+  ff::{Field, PrimeField},
+  Group,
+};
+
+mod poly;
+pub use poly::*;
+
+#[cfg(test)]
+mod tests;
+
+/// A curve usable with this library.
+pub trait DivisorCurve: Group
+where
+  Self::Scalar: PrimeField,
+{
+  /// An element of the field this curve is defined over.
+  type FieldElement: PrimeField;
+
+  /// The A in the curve equation y^2 = x^3 + A x + B.
+  fn a() -> Self::FieldElement;
+  /// The B in the curve equation y^2 = x^3 + A x + B.
+  fn b() -> Self::FieldElement;
+
+  /// y^2 - x^3 - A x - B
+  ///
+  /// Section 2 of the security proofs define this modulus.
+  ///
+  /// This MUST NOT be overriden.
+  fn divisor_modulus() -> Poly<Self::FieldElement> {
+    Poly {
+      // 0 y**1, 1 y*2
+      y_coefficients: vec![Self::FieldElement::ZERO, Self::FieldElement::ONE],
+      yx_coefficients: vec![],
+      x_coefficients: vec![
+        // - A x
+        -Self::a(),
+        // 0 x^2
+        Self::FieldElement::ZERO,
+        // - x^3
+        -Self::FieldElement::ONE,
+      ],
+      // - B
+      zero_coefficient: -Self::b(),
+    }
+  }
+
+  /// Convert a point to its x and y coordinates.
+  ///
+  /// Returns None if passed the point at infinity.
+  fn to_xy(point: Self) -> Option<(Self::FieldElement, Self::FieldElement)>;
+}
+
+/// Calculate the slope and intercept between two points.
+///
+/// This function panics when `a @ infinity`, `b @ infinity`, `a == b`, or when `a == -b`.
+pub(crate) fn slope_intercept<C: DivisorCurve>(a: C, b: C) -> (C::FieldElement, C::FieldElement) {
+  let (ax, ay) = C::to_xy(a).unwrap();
+  debug_assert_eq!(C::divisor_modulus().eval(ax, ay), C::FieldElement::ZERO);
+  let (bx, by) = C::to_xy(b).unwrap();
+  debug_assert_eq!(C::divisor_modulus().eval(bx, by), C::FieldElement::ZERO);
+  let slope = (by - ay) *
+    Option::<C::FieldElement>::from((bx - ax).invert())
+      .expect("trying to get slope/intercept of points sharing an x coordinate");
+  let intercept = by - (slope * bx);
+  debug_assert!(bool::from((ay - (slope * ax) - intercept).is_zero()));
+  debug_assert!(bool::from((by - (slope * bx) - intercept).is_zero()));
+  (slope, intercept)
+}
+
+// The line interpolating two points.
+fn line<C: DivisorCurve>(a: C, mut b: C) -> Poly<C::FieldElement> {
+  // If they're both the point at infinity, we simply set the line to one
+  if bool::from(a.is_identity() & b.is_identity()) {
+    return Poly {
+      y_coefficients: vec![],
+      yx_coefficients: vec![],
+      x_coefficients: vec![],
+      zero_coefficient: C::FieldElement::ONE,
+    };
+  }
+
+  // If either point is the point at infinity, or these are additive inverses, the line is
+  // `1 * x - x`. The first `x` is a term in the polynomial, the `x` is the `x` coordinate of these
+  // points (of which there is one, as the second point is either at infinity or has a matching `x`
+  // coordinate).
+  if bool::from(a.is_identity() | b.is_identity()) || (a == -b) {
+    let (x, _) = C::to_xy(if !bool::from(a.is_identity()) { a } else { b }).unwrap();
+    return Poly {
+      y_coefficients: vec![],
+      yx_coefficients: vec![],
+      x_coefficients: vec![C::FieldElement::ONE],
+      zero_coefficient: -x,
+    };
+  }
+
+  // If the points are equal, we use the line interpolating the sum of these points with the point
+  // at infinity
+  if a == b {
+    b = -a.double();
+  }
+
+  let (slope, intercept) = slope_intercept::<C>(a, b);
+
+  // Section 4 of the proofs explicitly state the line `L = y - lambda * x - mu`
+  // y - (slope * x) - intercept
+  Poly {
+    y_coefficients: vec![C::FieldElement::ONE],
+    yx_coefficients: vec![],
+    x_coefficients: vec![-slope],
+    zero_coefficient: -intercept,
+  }
+}
+
+/// Create a divisor interpolating the following points.
+///
+/// Returns None if:
+///   - No points were passed in
+///   - The points don't sum to the point at infinity
+///   - A passed in point was the point at infinity
+#[allow(clippy::new_ret_no_self)]
+pub fn new_divisor<C: DivisorCurve>(points: &[C]) -> Option<Poly<C::FieldElement>> {
+  // A single point is either the point at infinity, or this doesn't sum to the point at infinity
+  // Both cause us to return None
+  if points.len() < 2 {
+    None?;
+  }
+  if points.iter().sum::<C>() != C::identity() {
+    None?;
+  }
+
+  // Create the initial set of divisors
+  let mut divs = vec![];
+  let mut iter = points.iter().copied();
+  while let Some(a) = iter.next() {
+    if a == C::identity() {
+      None?;
+    }
+
+    let b = iter.next();
+    if b == Some(C::identity()) {
+      None?;
+    }
+
+    // Draw the line between those points
+    divs.push((a + b.unwrap_or(C::identity()), line::<C>(a, b.unwrap_or(-a))));
+  }
+
+  let modulus = C::divisor_modulus();
+
+  // Pair them off until only one remains
+  while divs.len() > 1 {
+    let mut next_divs = vec![];
+    // If there's an odd amount of divisors, carry the odd one out to the next iteration
+    if (divs.len() % 2) == 1 {
+      next_divs.push(divs.pop().unwrap());
+    }
+
+    while let Some((a, a_div)) = divs.pop() {
+      let (b, b_div) = divs.pop().unwrap();
+
+      // Merge the two divisors
+      let numerator = a_div.mul_mod(b_div, &modulus).mul_mod(line::<C>(a, b), &modulus);
+      let denominator = line::<C>(a, -a).mul_mod(line::<C>(b, -b), &modulus);
+      let (q, r) = numerator.div_rem(&denominator);
+      assert_eq!(r, Poly::zero());
+
+      next_divs.push((a + b, q));
+    }
+
+    divs = next_divs;
+  }
+
+  // Return the unified divisor
+  Some(divs.remove(0).1)
+}
+
+#[cfg(any(test, feature = "ed25519"))]
+mod ed25519 {
+  use group::{
+    ff::{Field, PrimeField},
+    Group, GroupEncoding,
+  };
+  use dalek_ff_group::{FieldElement, EdwardsPoint};
+
+  impl crate::DivisorCurve for EdwardsPoint {
+    type FieldElement = FieldElement;
+
+    // Wei25519 a/b
+    // https://www.ietf.org/archive/id/draft-ietf-lwig-curve-representations-02.pdf E.3
+    fn a() -> Self::FieldElement {
+      let mut be_bytes =
+        hex::decode("2aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa984914a144").unwrap();
+      be_bytes.reverse();
+      let le_bytes = be_bytes;
+      Self::FieldElement::from_repr(le_bytes.try_into().unwrap()).unwrap()
+    }
+    fn b() -> Self::FieldElement {
+      let mut be_bytes =
+        hex::decode("7b425ed097b425ed097b425ed097b425ed097b425ed097b4260b5e9c7710c864").unwrap();
+      be_bytes.reverse();
+      let le_bytes = be_bytes;
+
+      Self::FieldElement::from_repr(le_bytes.try_into().unwrap()).unwrap()
+    }
+
+    // https://www.ietf.org/archive/id/draft-ietf-lwig-curve-representations-02.pdf E.2
+    fn to_xy(point: Self) -> Option<(Self::FieldElement, Self::FieldElement)> {
+      if bool::from(point.is_identity()) {
+        None?;
+      }
+
+      // Extract the y coordinate from the compressed point
+      let mut edwards_y = point.to_bytes();
+      let x_is_odd = edwards_y[31] >> 7;
+      edwards_y[31] &= (1 << 7) - 1;
+      let edwards_y = Self::FieldElement::from_repr(edwards_y).unwrap();
+
+      // Recover the x coordinate
+      let edwards_y_sq = edwards_y * edwards_y;
+      let D = -Self::FieldElement::from(121665u64) *
+        Self::FieldElement::from(121666u64).invert().unwrap();
+      let mut edwards_x = ((edwards_y_sq - Self::FieldElement::ONE) *
+        ((D * edwards_y_sq) + Self::FieldElement::ONE).invert().unwrap())
+      .sqrt()
+      .unwrap();
+      if u8::from(bool::from(edwards_x.is_odd())) != x_is_odd {
+        edwards_x = -edwards_x;
+      }
+
+      // Calculate the x and y coordinates for Wei25519
+      let edwards_y_plus_one = Self::FieldElement::ONE + edwards_y;
+      let one_minus_edwards_y = Self::FieldElement::ONE - edwards_y;
+      let wei_x = (edwards_y_plus_one * one_minus_edwards_y.invert().unwrap()) +
+        (Self::FieldElement::from(486662u64) * Self::FieldElement::from(3u64).invert().unwrap());
+      let c =
+        (-(Self::FieldElement::from(486662u64) + Self::FieldElement::from(2u64))).sqrt().unwrap();
+      let wei_y = c * edwards_y_plus_one * (one_minus_edwards_y * edwards_x).invert().unwrap();
+      Some((wei_x, wei_y))
+    }
+  }
+}

crypto/evrf/divisors/src/poly.rs

@@ -0,0 +1,430 @@
+use core::ops::{Add, Neg, Sub, Mul, Rem};
+
+use zeroize::Zeroize;
+
+use group::ff::PrimeField;
+
+/// A structure representing a Polynomial with x**i, y**i, and y**i * x**j terms.
+#[derive(Clone, PartialEq, Eq, Debug, Zeroize)]
+pub struct Poly<F: PrimeField + From<u64>> {
+  /// c[i] * y ** (i + 1)
+  pub y_coefficients: Vec<F>,
+  /// c[i][j] * y ** (i + 1) x ** (j + 1)
+  pub yx_coefficients: Vec<Vec<F>>,
+  /// c[i] * x ** (i + 1)
+  pub x_coefficients: Vec<F>,
+  /// Coefficient for x ** 0, y ** 0, and x ** 0 y ** 0 (the coefficient for 1)
+  pub zero_coefficient: F,
+}
+
+impl<F: PrimeField + From<u64>> Poly<F> {
+  /// A polynomial for zero.
+  pub fn zero() -> Self {
+    Poly {
+      y_coefficients: vec![],
+      yx_coefficients: vec![],
+      x_coefficients: vec![],
+      zero_coefficient: F::ZERO,
+    }
+  }
+
+  /// The amount of terms in the polynomial.
+  #[allow(clippy::len_without_is_empty)]
+  #[must_use]
+  pub fn len(&self) -> usize {
+    self.y_coefficients.len() +
+      self.yx_coefficients.iter().map(Vec::len).sum::<usize>() +
+      self.x_coefficients.len() +
+      usize::from(u8::from(self.zero_coefficient != F::ZERO))
+  }
+
+  // Remove high-order zero terms, allowing the length of the vectors to equal the amount of terms.
+  pub(crate) fn tidy(&mut self) {
+    let tidy = |vec: &mut Vec<F>| {
+      while vec.last() == Some(&F::ZERO) {
+        vec.pop();
+      }
+    };
+
+    tidy(&mut self.y_coefficients);
+    for vec in self.yx_coefficients.iter_mut() {
+      tidy(vec);
+    }
+    while self.yx_coefficients.last() == Some(&vec![]) {
+      self.yx_coefficients.pop();
+    }
+    tidy(&mut self.x_coefficients);
+  }
+}
+
+impl<F: PrimeField + From<u64>> Add<&Self> for Poly<F> {
+  type Output = Self;
+
+  fn add(mut self, other: &Self) -> Self {
+    // Expand to be the neeeded size
+    while self.y_coefficients.len() < other.y_coefficients.len() {
+      self.y_coefficients.push(F::ZERO);
+    }
+    while self.yx_coefficients.len() < other.yx_coefficients.len() {
+      self.yx_coefficients.push(vec![]);
+    }
+    for i in 0 .. other.yx_coefficients.len() {
+      while self.yx_coefficients[i].len() < other.yx_coefficients[i].len() {
+        self.yx_coefficients[i].push(F::ZERO);
+      }
+    }
+    while self.x_coefficients.len() < other.x_coefficients.len() {
+      self.x_coefficients.push(F::ZERO);
+    }
+
+    // Perform the addition
+    for (i, coeff) in other.y_coefficients.iter().enumerate() {
+      self.y_coefficients[i] += coeff;
+    }
+    for (i, coeffs) in other.yx_coefficients.iter().enumerate() {
+      for (j, coeff) in coeffs.iter().enumerate() {
+        self.yx_coefficients[i][j] += coeff;
+      }
+    }
+    for (i, coeff) in other.x_coefficients.iter().enumerate() {
+      self.x_coefficients[i] += coeff;
+    }
+    self.zero_coefficient += other.zero_coefficient;
+
+    self.tidy();
+    self
+  }
+}
+
+impl<F: PrimeField + From<u64>> Neg for Poly<F> {
+  type Output = Self;
+
+  fn neg(mut self) -> Self {
+    for y_coeff in self.y_coefficients.iter_mut() {
+      *y_coeff = -*y_coeff;
+    }
+    for yx_coeffs in self.yx_coefficients.iter_mut() {
+      for yx_coeff in yx_coeffs.iter_mut() {
+        *yx_coeff = -*yx_coeff;
+      }
+    }
+    for x_coeff in self.x_coefficients.iter_mut() {
+      *x_coeff = -*x_coeff;
+    }
+    self.zero_coefficient = -self.zero_coefficient;
+
+    self
+  }
+}
+
+impl<F: PrimeField + From<u64>> Sub for Poly<F> {
+  type Output = Self;
+
+  fn sub(self, other: Self) -> Self {
+    self + &-other
+  }
+}
+
+impl<F: PrimeField + From<u64>> Mul<F> for Poly<F> {
+  type Output = Self;
+
+  fn mul(mut self, scalar: F) -> Self {
+    if scalar == F::ZERO {
+      return Poly::zero();
+    }
+
+    for y_coeff in self.y_coefficients.iter_mut() {
+      *y_coeff *= scalar;
+    }
+    for coeffs in self.yx_coefficients.iter_mut() {
+      for coeff in coeffs.iter_mut() {
+        *coeff *= scalar;
+      }
+    }
+    for x_coeff in self.x_coefficients.iter_mut() {
+      *x_coeff *= scalar;
+    }
+    self.zero_coefficient *= scalar;
+    self
+  }
+}
+
+impl<F: PrimeField + From<u64>> Poly<F> {
+  #[must_use]
+  fn shift_by_x(mut self, power_of_x: usize) -> Self {
+    if power_of_x == 0 {
+      return self;
+    }
+
+    // Shift up every x coefficient
+    for _ in 0 .. power_of_x {
+      self.x_coefficients.insert(0, F::ZERO);
+      for yx_coeffs in &mut self.yx_coefficients {
+        yx_coeffs.insert(0, F::ZERO);
+      }
+    }
+
+    // Move the zero coefficient
+    self.x_coefficients[power_of_x - 1] = self.zero_coefficient;
+    self.zero_coefficient = F::ZERO;
+
+    // Move the y coefficients
+    // Start by creating yx coefficients with the necessary powers of x
+    let mut yx_coefficients_to_push = vec![];
+    while yx_coefficients_to_push.len() < power_of_x {
+      yx_coefficients_to_push.push(F::ZERO);
+    }
+    // Now, ensure the yx coefficients has the slots for the y coefficients we're moving
+    while self.yx_coefficients.len() < self.y_coefficients.len() {
+      self.yx_coefficients.push(yx_coefficients_to_push.clone());
+    }
+    // Perform the move
+    for (i, y_coeff) in self.y_coefficients.drain(..).enumerate() {
+      self.yx_coefficients[i][power_of_x - 1] = y_coeff;
+    }
+
+    self
+  }
+
+  #[must_use]
+  fn shift_by_y(mut self, power_of_y: usize) -> Self {
+    if power_of_y == 0 {
+      return self;
+    }
+
+    // Shift up every y coefficient
+    for _ in 0 .. power_of_y {
+      self.y_coefficients.insert(0, F::ZERO);
+      self.yx_coefficients.insert(0, vec![]);
+    }
+
+    // Move the zero coefficient
+    self.y_coefficients[power_of_y - 1] = self.zero_coefficient;
+    self.zero_coefficient = F::ZERO;
+
+    // Move the x coefficients
+    self.yx_coefficients[power_of_y - 1] = self.x_coefficients;
+    self.x_coefficients = vec![];
+
+    self
+  }
+}
+
+impl<F: PrimeField + From<u64>> Mul for Poly<F> {
+  type Output = Self;
+
+  fn mul(self, other: Self) -> Self {
+    let mut res = self.clone() * other.zero_coefficient;
+
+    for (i, y_coeff) in other.y_coefficients.iter().enumerate() {
+      let scaled = self.clone() * *y_coeff;
+      res = res + &scaled.shift_by_y(i + 1);
+    }
+
+    for (y_i, yx_coeffs) in other.yx_coefficients.iter().enumerate() {
+      for (x_i, yx_coeff) in yx_coeffs.iter().enumerate() {
+        let scaled = self.clone() * *yx_coeff;
+        res = res + &scaled.shift_by_y(y_i + 1).shift_by_x(x_i + 1);
+      }
+    }
+
+    for (i, x_coeff) in other.x_coefficients.iter().enumerate() {
+      let scaled = self.clone() * *x_coeff;
+      res = res + &scaled.shift_by_x(i + 1);
+    }
+
+    res.tidy();
+    res
+  }
+}
+
+impl<F: PrimeField + From<u64>> Poly<F> {
+  /// Perform multiplication mod `modulus`.
+  #[must_use]
+  pub fn mul_mod(self, other: Self, modulus: &Self) -> Self {
+    ((self % modulus) * (other % modulus)) % modulus
+  }
+
+  /// Perform division, returning the result and remainder.
+  ///
+  /// Panics upon division by zero, with undefined behavior if a non-tidy divisor is used.
+  #[must_use]
+  pub fn div_rem(self, divisor: &Self) -> (Self, Self) {
+    // The leading y coefficient and associated x coefficient.
+    let leading_y = |poly: &Self| -> (_, _) {
+      if poly.y_coefficients.len() > poly.yx_coefficients.len() {
+        (poly.y_coefficients.len(), 0)
+      } else if !poly.yx_coefficients.is_empty() {
+        (poly.yx_coefficients.len(), poly.yx_coefficients.last().unwrap().len())
+      } else {
+        (0, poly.x_coefficients.len())
+      }
+    };
+
+    let (div_y, div_x) = leading_y(divisor);
+    // If this divisor is actually a scalar, don't perform long division
+    if (div_y == 0) && (div_x == 0) {
+      return (self * divisor.zero_coefficient.invert().unwrap(), Poly::zero());
+    }
+
+    // Remove leading terms until the value is less than the divisor
+    let mut quotient: Poly<F> = Poly::zero();
+    let mut remainder = self.clone();
+    loop {
+      // If there's nothing left to divide, return
+      if remainder == Poly::zero() {
+        break;
+      }
+
+      let (rem_y, rem_x) = leading_y(&remainder);
+      if (rem_y < div_y) || (rem_x < div_x) {
+        break;
+      }
+
+      let get = |poly: &Poly<F>, y_pow: usize, x_pow: usize| -> F {
+        if (y_pow == 0) && (x_pow == 0) {
+          poly.zero_coefficient
+        } else if x_pow == 0 {
+          poly.y_coefficients[y_pow - 1]
+        } else if y_pow == 0 {
+          poly.x_coefficients[x_pow - 1]
+        } else {
+          poly.yx_coefficients[y_pow - 1][x_pow - 1]
+        }
+      };
+      let coeff_numerator = get(&remainder, rem_y, rem_x);
+      let coeff_denominator = get(divisor, div_y, div_x);
+
+      // We want coeff_denominator scaled by x to equal coeff_numerator
+      // x * d = n
+      // n / d = x
+      let mut quotient_term = Poly::zero();
+      // Because this is the coefficient for the leading term of a tidied polynomial, it must be
+      // non-zero
+      quotient_term.zero_coefficient = coeff_numerator * coeff_denominator.invert().unwrap();
+
+      // Add the necessary yx powers
+      let delta_y = rem_y - div_y;
+      let delta_x = rem_x - div_x;
+      let quotient_term = quotient_term.shift_by_y(delta_y).shift_by_x(delta_x);
+
+      let to_remove = quotient_term.clone() * divisor.clone();
+      debug_assert_eq!(get(&to_remove, rem_y, rem_x), coeff_numerator);
+
+      remainder = remainder - to_remove;
+      quotient = quotient + &quotient_term;
+    }
+    debug_assert_eq!((quotient.clone() * divisor.clone()) + &remainder, self);
+
+    (quotient, remainder)
+  }
+}
+
+impl<F: PrimeField + From<u64>> Rem<&Self> for Poly<F> {
+  type Output = Self;
+
+  fn rem(self, modulus: &Self) -> Self {
+    self.div_rem(modulus).1
+  }
+}
+
+impl<F: PrimeField + From<u64>> Poly<F> {
+  /// Evaluate this polynomial with the specified x/y values.
+  ///
+  /// Panics on polynomials with terms whose powers exceed 2**64.
+  #[must_use]
+  pub fn eval(&self, x: F, y: F) -> F {
+    let mut res = self.zero_coefficient;
+    for (pow, coeff) in
+      self.y_coefficients.iter().enumerate().map(|(i, v)| (u64::try_from(i + 1).unwrap(), v))
+    {
+      res += y.pow([pow]) * coeff;
+    }
+    for (y_pow, coeffs) in
+      self.yx_coefficients.iter().enumerate().map(|(i, v)| (u64::try_from(i + 1).unwrap(), v))
+    {
+      let y_pow = y.pow([y_pow]);
+      for (x_pow, coeff) in
+        coeffs.iter().enumerate().map(|(i, v)| (u64::try_from(i + 1).unwrap(), v))
+      {
+        res += y_pow * x.pow([x_pow]) * coeff;
+      }
+    }
+    for (pow, coeff) in
+      self.x_coefficients.iter().enumerate().map(|(i, v)| (u64::try_from(i + 1).unwrap(), v))
+    {
+      res += x.pow([pow]) * coeff;
+    }
+    res
+  }
+
+  /// Differentiate a polynomial, reduced by a modulus with a leading y term y**2 x**0, by x and y.
+  ///
+  /// This function panics if a y**2 term is present within the polynomial.
+  #[must_use]
+  pub fn differentiate(&self) -> (Poly<F>, Poly<F>) {
+    assert!(self.y_coefficients.len() <= 1);
+    assert!(self.yx_coefficients.len() <= 1);
+
+    // Differentation by x practically involves:
+    // - Dropping everything without an x component
+    // - Shifting everything down a power of x
+    // - Multiplying the new coefficient by the power it prior was used with
+    let diff_x = {
+      let mut diff_x = Poly {
+        y_coefficients: vec![],
+        yx_coefficients: vec![],
+        x_coefficients: vec![],
+        zero_coefficient: F::ZERO,
+      };
+      if !self.x_coefficients.is_empty() {
+        let mut x_coeffs = self.x_coefficients.clone();
+        diff_x.zero_coefficient = x_coeffs.remove(0);
+        diff_x.x_coefficients = x_coeffs;
+
+        let mut prior_x_power = F::from(2);
+        for x_coeff in &mut diff_x.x_coefficients {
+          *x_coeff *= prior_x_power;
+          prior_x_power += F::ONE;
+        }
+      }
+
+      if !self.yx_coefficients.is_empty() {
+        let mut yx_coeffs = self.yx_coefficients[0].clone();
+        diff_x.y_coefficients = vec![yx_coeffs.remove(0)];
+        diff_x.yx_coefficients = vec![yx_coeffs];
+
+        let mut prior_x_power = F::from(2);
+        for yx_coeff in &mut diff_x.yx_coefficients[0] {
+          *yx_coeff *= prior_x_power;
+          prior_x_power += F::ONE;
+        }
+      }
+
+      diff_x.tidy();
+      diff_x
+    };
+
+    // Differentation by y is trivial
+    // It's the y coefficient as the zero coefficient, and the yx coefficients as the x
+    // coefficients
+    // This is thanks to any y term over y^2 being reduced out
+    let diff_y = Poly {
+      y_coefficients: vec![],
+      yx_coefficients: vec![],
+      x_coefficients: self.yx_coefficients.first().cloned().unwrap_or(vec![]),
+      zero_coefficient: self.y_coefficients.first().cloned().unwrap_or(F::ZERO),
+    };
+
+    (diff_x, diff_y)
+  }
+
+  /// Normalize the x coefficient to 1.
+  ///
+  /// Panics if there is no x coefficient to normalize or if it cannot be normalized to 1.
+  #[must_use]
+  pub fn normalize_x_coefficient(self) -> Self {
+    let scalar = self.x_coefficients[0].invert().unwrap();
+    self * scalar
+  }
+}

crypto/evrf/divisors/src/tests/mod.rs

@@ -0,0 +1,247 @@
+use rand_core::OsRng;
+
+use group::{ff::Field, Group, Curve};
+use dalek_ff_group::EdwardsPoint;
+use pasta_curves::{
+  arithmetic::{Coordinates, CurveAffine},
+  Ep, Fp,
+};
+
+use crate::{DivisorCurve, Poly, new_divisor};
+
+impl DivisorCurve for Ep {
+  type FieldElement = Fp;
+
+  fn a() -> Self::FieldElement {
+    Self::FieldElement::ZERO
+  }
+  fn b() -> Self::FieldElement {
+    Self::FieldElement::from(5u64)
+  }
+
+  fn to_xy(point: Self) -> Option<(Self::FieldElement, Self::FieldElement)> {
+    Option::<Coordinates<_>>::from(point.to_affine().coordinates())
+      .map(|coords| (*coords.x(), *coords.y()))
+  }
+}
+
+// Equation 4 in the security proofs
+fn check_divisor<C: DivisorCurve>(points: Vec<C>) {
+  // Create the divisor
+  let divisor = new_divisor::<C>(&points).unwrap();
+  let eval = |c| {
+    let (x, y) = C::to_xy(c).unwrap();
+    divisor.eval(x, y)
+  };
+
+  // Decide challgenges
+  let c0 = C::random(&mut OsRng);
+  let c1 = C::random(&mut OsRng);
+  let c2 = -(c0 + c1);
+  let (slope, intercept) = crate::slope_intercept::<C>(c0, c1);
+
+  let mut rhs = <C as DivisorCurve>::FieldElement::ONE;
+  for point in points {
+    let (x, y) = C::to_xy(point).unwrap();
+    rhs *= intercept - (y - (slope * x));
+  }
+  assert_eq!(eval(c0) * eval(c1) * eval(c2), rhs);
+}
+
+fn test_divisor<C: DivisorCurve>() {
+  for i in 1 ..= 255 {
+    println!("Test iteration {i}");
+
+    // Select points
+    let mut points = vec![];
+    for _ in 0 .. i {
+      points.push(C::random(&mut OsRng));
+    }
+    points.push(-points.iter().sum::<C>());
+    println!("Points {}", points.len());
+
+    // Perform the original check
+    check_divisor(points.clone());
+
+    // Create the divisor
+    let divisor = new_divisor::<C>(&points).unwrap();
+
+    // For a divisor interpolating 256 points, as one does when interpreting a 255-bit discrete log
+    // with the result of its scalar multiplication against a fixed generator, the lengths of the
+    // yx/x coefficients shouldn't supersede the following bounds
+    assert!((divisor.yx_coefficients.first().unwrap_or(&vec![]).len()) <= 126);
+    assert!((divisor.x_coefficients.len() - 1) <= 127);
+    assert!(
+      (1 + divisor.yx_coefficients.first().unwrap_or(&vec![]).len() +
+        (divisor.x_coefficients.len() - 1) +
+        1) <=
+        255
+    );
+
+    // Decide challgenges
+    let c0 = C::random(&mut OsRng);
+    let c1 = C::random(&mut OsRng);
+    let c2 = -(c0 + c1);
+    let (slope, intercept) = crate::slope_intercept::<C>(c0, c1);
+
+    // Perform the Logarithmic derivative check
+    {
+      let dx_over_dz = {
+        let dx = Poly {
+          y_coefficients: vec![],
+          yx_coefficients: vec![],
+          x_coefficients: vec![C::FieldElement::ZERO, C::FieldElement::from(3)],
+          zero_coefficient: C::a(),
+        };
+
+        let dy = Poly {
+          y_coefficients: vec![C::FieldElement::from(2)],
+          yx_coefficients: vec![],
+          x_coefficients: vec![],
+          zero_coefficient: C::FieldElement::ZERO,
+        };
+
+        let dz = (dy.clone() * -slope) + &dx;
+
+        // We want dx/dz, and dz/dx is equal to dy/dx - slope
+        // Sagemath claims this, dy / dz, is the proper inverse
+        (dy, dz)
+      };
+
+      {
+        let sanity_eval = |c| {
+          let (x, y) = C::to_xy(c).unwrap();
+          dx_over_dz.0.eval(x, y) * dx_over_dz.1.eval(x, y).invert().unwrap()
+        };
+        let sanity = sanity_eval(c0) + sanity_eval(c1) + sanity_eval(c2);
+        // This verifies the dx/dz polynomial is correct
+        assert_eq!(sanity, C::FieldElement::ZERO);
+      }
+
+      // Logarithmic derivative check
+      let test = |divisor: Poly<_>| {
+        let (dx, dy) = divisor.differentiate();
+
+        let lhs = |c| {
+          let (x, y) = C::to_xy(c).unwrap();
+
+          let n_0 = (C::FieldElement::from(3) * (x * x)) + C::a();
+          let d_0 = (C::FieldElement::from(2) * y).invert().unwrap();
+          let p_0_n_0 = n_0 * d_0;
+
+          let n_1 = dy.eval(x, y);
+          let first = p_0_n_0 * n_1;
+
+          let second = dx.eval(x, y);
+
+          let d_1 = divisor.eval(x, y);
+
+          let fraction_1_n = first + second;
+          let fraction_1_d = d_1;
+
+          let fraction_2_n = dx_over_dz.0.eval(x, y);
+          let fraction_2_d = dx_over_dz.1.eval(x, y);
+
+          fraction_1_n * fraction_2_n * (fraction_1_d * fraction_2_d).invert().unwrap()
+        };
+        let lhs = lhs(c0) + lhs(c1) + lhs(c2);
+
+        let mut rhs = C::FieldElement::ZERO;
+        for point in &points {
+          let (x, y) = <C as DivisorCurve>::to_xy(*point).unwrap();
+          rhs += (intercept - (y - (slope * x))).invert().unwrap();
+        }
+
+        assert_eq!(lhs, rhs);
+      };
+      // Test the divisor and the divisor with a normalized x coefficient
+      test(divisor.clone());
+      test(divisor.normalize_x_coefficient());
+    }
+  }
+}
+
+fn test_same_point<C: DivisorCurve>() {
+  let mut points = vec![C::random(&mut OsRng)];
+  points.push(points[0]);
+  points.push(-points.iter().sum::<C>());
+  check_divisor(points);
+}
+
+fn test_subset_sum_to_infinity<C: DivisorCurve>() {
+  // Internally, a binary tree algorithm is used
+  // This executes the first pass to end up with [0, 0] for further reductions
+  {
+    let mut points = vec![C::random(&mut OsRng)];
+    points.push(-points[0]);
+
+    let next = C::random(&mut OsRng);
+    points.push(next);
+    points.push(-next);
+    check_divisor(points);
+  }
+
+  // This executes the first pass to end up with [0, X, -X, 0]
+  {
+    let mut points = vec![C::random(&mut OsRng)];
+    points.push(-points[0]);
+
+    let x_1 = C::random(&mut OsRng);
+    let x_2 = C::random(&mut OsRng);
+    points.push(x_1);
+    points.push(x_2);
+
+    points.push(-x_1);
+    points.push(-x_2);
+
+    let next = C::random(&mut OsRng);
+    points.push(next);
+    points.push(-next);
+    check_divisor(points);
+  }
+}
+
+#[test]
+fn test_divisor_pallas() {
+  test_divisor::<Ep>();
+  test_same_point::<Ep>();
+  test_subset_sum_to_infinity::<Ep>();
+}
+
+#[test]
+fn test_divisor_ed25519() {
+  // Since we're implementing Wei25519 ourselves, check the isomorphism works as expected
+  {
+    let incomplete_add = |p1, p2| {
+      let (x1, y1) = EdwardsPoint::to_xy(p1).unwrap();
+      let (x2, y2) = EdwardsPoint::to_xy(p2).unwrap();
+
+      // mmadd-1998-cmo
+      let u = y2 - y1;
+      let uu = u * u;
+      let v = x2 - x1;
+      let vv = v * v;
+      let vvv = v * vv;
+      let R = vv * x1;
+      let A = uu - vvv - R.double();
+      let x3 = v * A;
+      let y3 = (u * (R - A)) - (vvv * y1);
+      let z3 = vvv;
+
+      // Normalize from XYZ to XY
+      let x3 = x3 * z3.invert().unwrap();
+      let y3 = y3 * z3.invert().unwrap();
+
+      // Edwards addition -> Wei25519 coordinates should be equivalent to Wei25519 addition
+      assert_eq!(EdwardsPoint::to_xy(p1 + p2).unwrap(), (x3, y3));
+    };
+
+    for _ in 0 .. 256 {
+      incomplete_add(EdwardsPoint::random(&mut OsRng), EdwardsPoint::random(&mut OsRng));
+    }
+  }
+
+  test_divisor::<EdwardsPoint>();
+  test_same_point::<EdwardsPoint>();
+  test_subset_sum_to_infinity::<EdwardsPoint>();
+}

crypto/evrf/divisors/src/tests/poly.rs

@@ -0,0 +1,129 @@
+use group::ff::Field;
+use pasta_curves::Ep;
+
+use crate::{DivisorCurve, Poly};
+
+type F = <Ep as DivisorCurve>::FieldElement;
+
+#[test]
+fn test_poly() {
+  let zero = F::ZERO;
+  let one = F::ONE;
+
+  {
+    let mut poly = Poly::zero();
+    poly.y_coefficients = vec![zero, one];
+
+    let mut modulus = Poly::zero();
+    modulus.y_coefficients = vec![one];
+    assert_eq!(poly % &modulus, Poly::zero());
+  }
+
+  {
+    let mut poly = Poly::zero();
+    poly.y_coefficients = vec![zero, one];
+
+    let mut squared = Poly::zero();
+    squared.y_coefficients = vec![zero, zero, zero, one];
+    assert_eq!(poly.clone() * poly.clone(), squared);
+  }
+
+  {
+    let mut a = Poly::zero();
+    a.zero_coefficient = F::from(2u64);
+
+    let mut b = Poly::zero();
+    b.zero_coefficient = F::from(3u64);
+
+    let mut res = Poly::zero();
+    res.zero_coefficient = F::from(6u64);
+    assert_eq!(a.clone() * b.clone(), res);
+
+    b.y_coefficients = vec![F::from(4u64)];
+    res.y_coefficients = vec![F::from(8u64)];
+    assert_eq!(a.clone() * b.clone(), res);
+    assert_eq!(b.clone() * a.clone(), res);
+
+    a.x_coefficients = vec![F::from(5u64)];
+    res.x_coefficients = vec![F::from(15u64)];
+    res.yx_coefficients = vec![vec![F::from(20u64)]];
+    assert_eq!(a.clone() * b.clone(), res);
+    assert_eq!(b * a.clone(), res);
+
+    // res is now 20xy + 8*y + 15*x + 6
+    // res ** 2 =
+    //   400*x^2*y^2 + 320*x*y^2 + 64*y^2 + 600*x^2*y + 480*x*y + 96*y + 225*x^2 + 180*x + 36
+
+    let mut squared = Poly::zero();
+    squared.y_coefficients = vec![F::from(96u64), F::from(64u64)];
+    squared.yx_coefficients =
+      vec![vec![F::from(480u64), F::from(600u64)], vec![F::from(320u64), F::from(400u64)]];
+    squared.x_coefficients = vec![F::from(180u64), F::from(225u64)];
+    squared.zero_coefficient = F::from(36u64);
+    assert_eq!(res.clone() * res, squared);
+  }
+}
+
+#[test]
+fn test_differentation() {
+  let random = || F::random(&mut OsRng);
+
+  let input = Poly {
+    y_coefficients: vec![random()],
+    yx_coefficients: vec![vec![random()]],
+    x_coefficients: vec![random(), random(), random()],
+    zero_coefficient: random(),
+  };
+  let (diff_x, diff_y) = input.differentiate();
+  assert_eq!(
+    diff_x,
+    Poly {
+      y_coefficients: vec![input.yx_coefficients[0][0]],
+      yx_coefficients: vec![],
+      x_coefficients: vec![
+        F::from(2) * input.x_coefficients[1],
+        F::from(3) * input.x_coefficients[2]
+      ],
+      zero_coefficient: input.x_coefficients[0],
+    }
+  );
+  assert_eq!(
+    diff_y,
+    Poly {
+      y_coefficients: vec![],
+      yx_coefficients: vec![],
+      x_coefficients: vec![input.yx_coefficients[0][0]],
+      zero_coefficient: input.y_coefficients[0],
+    }
+  );
+
+  let input = Poly {
+    y_coefficients: vec![random()],
+    yx_coefficients: vec![vec![random(), random()]],
+    x_coefficients: vec![random(), random(), random(), random()],
+    zero_coefficient: random(),
+  };
+  let (diff_x, diff_y) = input.differentiate();
+  assert_eq!(
+    diff_x,
+    Poly {
+      y_coefficients: vec![input.yx_coefficients[0][0]],
+      yx_coefficients: vec![vec![F::from(2) * input.yx_coefficients[0][1]]],
+      x_coefficients: vec![
+        F::from(2) * input.x_coefficients[1],
+        F::from(3) * input.x_coefficients[2],
+        F::from(4) * input.x_coefficients[3],
+      ],
+      zero_coefficient: input.x_coefficients[0],
+    }
+  );
+  assert_eq!(
+    diff_y,
+    Poly {
+      y_coefficients: vec![],
+      yx_coefficients: vec![],
+      x_coefficients: vec![input.yx_coefficients[0][0], input.yx_coefficients[0][1]],
+      zero_coefficient: input.y_coefficients[0],
+    }
+  );
+}