first commit

src/curve/curve_adds.rs

@@ -0,0 +1,158 @@
+use core::ops::Add;
+
+use plonky2::field::ops::Square;
+use plonky2::field::types::Field;
+
+use crate::curve::curve_types::{AffinePoint, Curve, ProjectivePoint};
+
+impl<C: Curve> Add<ProjectivePoint<C>> for ProjectivePoint<C> {
+    type Output = ProjectivePoint<C>;
+
+    fn add(self, rhs: ProjectivePoint<C>) -> Self::Output {
+        let ProjectivePoint {
+            x: x1,
+            y: y1,
+            z: z1,
+        } = self;
+        let ProjectivePoint {
+            x: x2,
+            y: y2,
+            z: z2,
+        } = rhs;
+
+        if z1 == C::BaseField::ZERO {
+            return rhs;
+        }
+        if z2 == C::BaseField::ZERO {
+            return self;
+        }
+
+        let x1z2 = x1 * z2;
+        let y1z2 = y1 * z2;
+        let x2z1 = x2 * z1;
+        let y2z1 = y2 * z1;
+
+        // Check if we're doubling or adding inverses.
+        if x1z2 == x2z1 {
+            if y1z2 == y2z1 {
+                // TODO: inline to avoid redundant muls.
+                return self.double();
+            }
+            if y1z2 == -y2z1 {
+                return ProjectivePoint::ZERO;
+            }
+        }
+
+        // From https://www.hyperelliptic.org/EFD/g1p/data/shortw/projective/addition/add-1998-cmo-2
+        let z1z2 = z1 * z2;
+        let u = y2z1 - y1z2;
+        let uu = u.square();
+        let v = x2z1 - x1z2;
+        let vv = v.square();
+        let vvv = v * vv;
+        let r = vv * x1z2;
+        let a = uu * z1z2 - vvv - r.double();
+        let x3 = v * a;
+        let y3 = u * (r - a) - vvv * y1z2;
+        let z3 = vvv * z1z2;
+        ProjectivePoint::nonzero(x3, y3, z3)
+    }
+}
+
+impl<C: Curve> Add<AffinePoint<C>> for ProjectivePoint<C> {
+    type Output = ProjectivePoint<C>;
+
+    fn add(self, rhs: AffinePoint<C>) -> Self::Output {
+        let ProjectivePoint {
+            x: x1,
+            y: y1,
+            z: z1,
+        } = self;
+        let AffinePoint {
+            x: x2,
+            y: y2,
+            zero: zero2,
+        } = rhs;
+
+        if z1 == C::BaseField::ZERO {
+            return rhs.to_projective();
+        }
+        if zero2 {
+            return self;
+        }
+
+        let x2z1 = x2 * z1;
+        let y2z1 = y2 * z1;
+
+        // Check if we're doubling or adding inverses.
+        if x1 == x2z1 {
+            if y1 == y2z1 {
+                // TODO: inline to avoid redundant muls.
+                return self.double();
+            }
+            if y1 == -y2z1 {
+                return ProjectivePoint::ZERO;
+            }
+        }
+
+        // From https://www.hyperelliptic.org/EFD/g1p/data/shortw/projective/addition/madd-1998-cmo
+        let u = y2z1 - y1;
+        let uu = u.square();
+        let v = x2z1 - x1;
+        let vv = v.square();
+        let vvv = v * vv;
+        let r = vv * x1;
+        let a = uu * z1 - vvv - r.double();
+        let x3 = v * a;
+        let y3 = u * (r - a) - vvv * y1;
+        let z3 = vvv * z1;
+        ProjectivePoint::nonzero(x3, y3, z3)
+    }
+}
+
+impl<C: Curve> Add<AffinePoint<C>> for AffinePoint<C> {
+    type Output = ProjectivePoint<C>;
+
+    fn add(self, rhs: AffinePoint<C>) -> Self::Output {
+        let AffinePoint {
+            x: x1,
+            y: y1,
+            zero: zero1,
+        } = self;
+        let AffinePoint {
+            x: x2,
+            y: y2,
+            zero: zero2,
+        } = rhs;
+
+        if zero1 {
+            return rhs.to_projective();
+        }
+        if zero2 {
+            return self.to_projective();
+        }
+
+        // Check if we're doubling or adding inverses.
+        if x1 == x2 {
+            if y1 == y2 {
+                return self.to_projective().double();
+            }
+            if y1 == -y2 {
+                return ProjectivePoint::ZERO;
+            }
+        }
+
+        // From https://www.hyperelliptic.org/EFD/g1p/data/shortw/projective/addition/mmadd-1998-cmo
+        let u = y2 - y1;
+        let uu = u.square();
+        let v = x2 - x1;
+        let vv = v.square();
+        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;
+        ProjectivePoint::nonzero(x3, y3, z3)
+    }
+}

src/curve/curve_msm.rs

@@ -0,0 +1,265 @@
+use alloc::vec::Vec;
+
+use itertools::Itertools;
+use plonky2::field::types::{Field, PrimeField};
+use plonky2_maybe_rayon::*;
+
+use crate::curve::curve_summation::affine_multisummation_best;
+use crate::curve::curve_types::{AffinePoint, Curve, ProjectivePoint};
+
+/// In Yao's method, we compute an affine summation for each digit. In a parallel setting, it would
+/// be easiest to assign individual summations to threads, but this would be sub-optimal because
+/// multi-summations can be more efficient than repeating individual summations (see
+/// `affine_multisummation_best`). Thus we divide digits into large chunks, and assign chunks of
+/// digits to threads. Note that there is a delicate balance here, as large chunks can result in
+/// uneven distributions of work among threads.
+const DIGITS_PER_CHUNK: usize = 80;
+
+#[derive(Clone, Debug)]
+pub struct MsmPrecomputation<C: Curve> {
+    /// For each generator (in the order they were passed to `msm_precompute`), contains a vector
+    /// of powers, i.e. [(2^w)^i] for i < DIGITS.
+    // TODO: Use compressed coordinates here.
+    powers_per_generator: Vec<Vec<AffinePoint<C>>>,
+
+    /// The window size.
+    w: usize,
+}
+
+pub fn msm_precompute<C: Curve>(
+    generators: &[ProjectivePoint<C>],
+    w: usize,
+) -> MsmPrecomputation<C> {
+    MsmPrecomputation {
+        powers_per_generator: generators
+            .into_par_iter()
+            .map(|&g| precompute_single_generator(g, w))
+            .collect(),
+        w,
+    }
+}
+
+fn precompute_single_generator<C: Curve>(g: ProjectivePoint<C>, w: usize) -> Vec<AffinePoint<C>> {
+    let digits = (C::ScalarField::BITS + w - 1) / w;
+    let mut powers: Vec<ProjectivePoint<C>> = Vec::with_capacity(digits);
+    powers.push(g);
+    for i in 1..digits {
+        let mut power_i_proj = powers[i - 1];
+        for _j in 0..w {
+            power_i_proj = power_i_proj.double();
+        }
+        powers.push(power_i_proj);
+    }
+    ProjectivePoint::batch_to_affine(&powers)
+}
+
+pub fn msm_parallel<C: Curve>(
+    scalars: &[C::ScalarField],
+    generators: &[ProjectivePoint<C>],
+    w: usize,
+) -> ProjectivePoint<C> {
+    let precomputation = msm_precompute(generators, w);
+    msm_execute_parallel(&precomputation, scalars)
+}
+
+pub fn msm_execute<C: Curve>(
+    precomputation: &MsmPrecomputation<C>,
+    scalars: &[C::ScalarField],
+) -> ProjectivePoint<C> {
+    assert_eq!(precomputation.powers_per_generator.len(), scalars.len());
+    let w = precomputation.w;
+    let digits = (C::ScalarField::BITS + w - 1) / w;
+    let base = 1 << w;
+
+    // This is a variant of Yao's method, adapted to the multi-scalar setting. Because we use
+    // extremely large windows, the repeated scans in Yao's method could be more expensive than the
+    // actual group operations. To avoid this, we store a multimap from each possible digit to the
+    // positions in which that digit occurs in the scalars. These positions have the form (i, j),
+    // where i is the index of the generator and j is an index into the digits of the scalar
+    // associated with that generator.
+    let mut digit_occurrences: Vec<Vec<(usize, usize)>> = Vec::with_capacity(digits);
+    for _i in 0..base {
+        digit_occurrences.push(Vec::new());
+    }
+    for (i, scalar) in scalars.iter().enumerate() {
+        let digits = to_digits::<C>(scalar, w);
+        for (j, &digit) in digits.iter().enumerate() {
+            digit_occurrences[digit].push((i, j));
+        }
+    }
+
+    let mut y = ProjectivePoint::ZERO;
+    let mut u = ProjectivePoint::ZERO;
+
+    for digit in (1..base).rev() {
+        for &(i, j) in &digit_occurrences[digit] {
+            u = u + precomputation.powers_per_generator[i][j];
+        }
+        y = y + u;
+    }
+
+    y
+}
+
+pub fn msm_execute_parallel<C: Curve>(
+    precomputation: &MsmPrecomputation<C>,
+    scalars: &[C::ScalarField],
+) -> ProjectivePoint<C> {
+    assert_eq!(precomputation.powers_per_generator.len(), scalars.len());
+    let w = precomputation.w;
+    let digits = (C::ScalarField::BITS + w - 1) / w;
+    let base = 1 << w;
+
+    // This is a variant of Yao's method, adapted to the multi-scalar setting. Because we use
+    // extremely large windows, the repeated scans in Yao's method could be more expensive than the
+    // actual group operations. To avoid this, we store a multimap from each possible digit to the
+    // positions in which that digit occurs in the scalars. These positions have the form (i, j),
+    // where i is the index of the generator and j is an index into the digits of the scalar
+    // associated with that generator.
+    let mut digit_occurrences: Vec<Vec<(usize, usize)>> = Vec::with_capacity(digits);
+    for _i in 0..base {
+        digit_occurrences.push(Vec::new());
+    }
+    for (i, scalar) in scalars.iter().enumerate() {
+        let digits = to_digits::<C>(scalar, w);
+        for (j, &digit) in digits.iter().enumerate() {
+            digit_occurrences[digit].push((i, j));
+        }
+    }
+
+    // For each digit, we add up the powers associated with all occurrences that digit.
+    let digits: Vec<usize> = (0..base).collect();
+    let digit_acc: Vec<ProjectivePoint<C>> = digits
+        .par_chunks(DIGITS_PER_CHUNK)
+        .flat_map(|chunk| {
+            let summations: Vec<Vec<AffinePoint<C>>> = chunk
+                .iter()
+                .map(|&digit| {
+                    digit_occurrences[digit]
+                        .iter()
+                        .map(|&(i, j)| precomputation.powers_per_generator[i][j])
+                        .collect()
+                })
+                .collect();
+            affine_multisummation_best(summations)
+        })
+        .collect();
+    // println!("Computing the per-digit summations (in parallel) took {}s", start.elapsed().as_secs_f64());
+
+    let mut y = ProjectivePoint::ZERO;
+    let mut u = ProjectivePoint::ZERO;
+    for digit in (1..base).rev() {
+        u = u + digit_acc[digit];
+        y = y + u;
+    }
+    // println!("Final summation (sequential) {}s", start.elapsed().as_secs_f64());
+    y
+}
+
+pub(crate) fn to_digits<C: Curve>(x: &C::ScalarField, w: usize) -> Vec<usize> {
+    let scalar_bits = C::ScalarField::BITS;
+    let num_digits = (scalar_bits + w - 1) / w;
+
+    // Convert x to a bool array.
+    let x_canonical: Vec<_> = x
+        .to_canonical_biguint()
+        .to_u64_digits()
+        .iter()
+        .cloned()
+        .pad_using(scalar_bits / 64, |_| 0)
+        .collect();
+    let mut x_bits = Vec::with_capacity(scalar_bits);
+    for i in 0..scalar_bits {
+        x_bits.push((x_canonical[i / 64] >> (i as u64 % 64) & 1) != 0);
+    }
+
+    let mut digits = Vec::with_capacity(num_digits);
+    for i in 0..num_digits {
+        let mut digit = 0;
+        for j in ((i * w)..((i + 1) * w).min(scalar_bits)).rev() {
+            digit <<= 1;
+            digit |= x_bits[j] as usize;
+        }
+        digits.push(digit);
+    }
+    digits
+}
+
+#[cfg(test)]
+mod tests {
+    use alloc::vec;
+
+    use num::BigUint;
+    use plonky2::field::secp256k1_scalar::Secp256K1Scalar;
+
+    use super::*;
+    use crate::curve::secp256k1::Secp256K1;
+
+    #[test]
+    fn test_to_digits() {
+        let x_canonical = [
+            0b10101010101010101010101010101010,
+            0b10101010101010101010101010101010,
+            0b11001100110011001100110011001100,
+            0b11001100110011001100110011001100,
+            0b11110000111100001111000011110000,
+            0b11110000111100001111000011110000,
+            0b00001111111111111111111111111111,
+            0b11111111111111111111111111111111,
+        ];
+        let x = Secp256K1Scalar::from_noncanonical_biguint(BigUint::from_slice(&x_canonical));
+        assert_eq!(x.to_canonical_biguint().to_u32_digits(), x_canonical);
+        assert_eq!(
+            to_digits::<Secp256K1>(&x, 17),
+            vec![
+                0b01010101010101010,
+                0b10101010101010101,
+                0b01010101010101010,
+                0b11001010101010101,
+                0b01100110011001100,
+                0b00110011001100110,
+                0b10011001100110011,
+                0b11110000110011001,
+                0b01111000011110000,
+                0b00111100001111000,
+                0b00011110000111100,
+                0b11111111111111110,
+                0b01111111111111111,
+                0b11111111111111000,
+                0b11111111111111111,
+                0b1,
+            ]
+        );
+    }
+
+    #[test]
+    fn test_msm() {
+        let w = 5;
+
+        let generator_1 = Secp256K1::GENERATOR_PROJECTIVE;
+        let generator_2 = generator_1 + generator_1;
+        let generator_3 = generator_1 + generator_2;
+
+        let scalar_1 = Secp256K1Scalar::from_noncanonical_biguint(BigUint::from_slice(&[
+            11111111, 22222222, 33333333, 44444444,
+        ]));
+        let scalar_2 = Secp256K1Scalar::from_noncanonical_biguint(BigUint::from_slice(&[
+            22222222, 22222222, 33333333, 44444444,
+        ]));
+        let scalar_3 = Secp256K1Scalar::from_noncanonical_biguint(BigUint::from_slice(&[
+            33333333, 22222222, 33333333, 44444444,
+        ]));
+
+        let generators = vec![generator_1, generator_2, generator_3];
+        let scalars = vec![scalar_1, scalar_2, scalar_3];
+
+        let precomputation = msm_precompute(&generators, w);
+        let result_msm = msm_execute(&precomputation, &scalars);
+
+        let result_naive = Secp256K1::convert(scalar_1) * generator_1
+            + Secp256K1::convert(scalar_2) * generator_2
+            + Secp256K1::convert(scalar_3) * generator_3;
+
+        assert_eq!(result_msm, result_naive);
+    }
+}

src/curve/curve_multiplication.rs

@@ -0,0 +1,100 @@
+use alloc::vec::Vec;
+use core::ops::Mul;
+
+use plonky2::field::types::{Field, PrimeField};
+
+use crate::curve::curve_types::{Curve, CurveScalar, ProjectivePoint};
+
+const WINDOW_BITS: usize = 4;
+const BASE: usize = 1 << WINDOW_BITS;
+
+fn digits_per_scalar<C: Curve>() -> usize {
+    (C::ScalarField::BITS + WINDOW_BITS - 1) / WINDOW_BITS
+}
+
+/// Precomputed state used for scalar x ProjectivePoint multiplications,
+/// specific to a particular generator.
+#[derive(Clone)]
+pub struct MultiplicationPrecomputation<C: Curve> {
+    /// [(2^w)^i] g for each i < digits_per_scalar.
+    powers: Vec<ProjectivePoint<C>>,
+}
+
+impl<C: Curve> ProjectivePoint<C> {
+    pub fn mul_precompute(&self) -> MultiplicationPrecomputation<C> {
+        let num_digits = digits_per_scalar::<C>();
+        let mut powers = Vec::with_capacity(num_digits);
+        powers.push(*self);
+        for i in 1..num_digits {
+            let mut power_i = powers[i - 1];
+            for _j in 0..WINDOW_BITS {
+                power_i = power_i.double();
+            }
+            powers.push(power_i);
+        }
+
+        MultiplicationPrecomputation { powers }
+    }
+
+    #[must_use]
+    pub fn mul_with_precomputation(
+        &self,
+        scalar: C::ScalarField,
+        precomputation: MultiplicationPrecomputation<C>,
+    ) -> Self {
+        // Yao's method; see https://koclab.cs.ucsb.edu/teaching/ecc/eccPapers/Doche-ch09.pdf
+        let precomputed_powers = precomputation.powers;
+
+        let digits = to_digits::<C>(&scalar);
+
+        let mut y = ProjectivePoint::ZERO;
+        let mut u = ProjectivePoint::ZERO;
+        let mut all_summands = Vec::new();
+        for j in (1..BASE).rev() {
+            let mut u_summands = Vec::new();
+            for (i, &digit) in digits.iter().enumerate() {
+                if digit == j as u64 {
+                    u_summands.push(precomputed_powers[i]);
+                }
+            }
+            all_summands.push(u_summands);
+        }
+
+        let all_sums: Vec<ProjectivePoint<C>> = all_summands
+            .iter()
+            .cloned()
+            .map(|vec| vec.iter().fold(ProjectivePoint::ZERO, |a, &b| a + b))
+            .collect();
+        for i in 0..all_sums.len() {
+            u = u + all_sums[i];
+            y = y + u;
+        }
+        y
+    }
+}
+
+impl<C: Curve> Mul<ProjectivePoint<C>> for CurveScalar<C> {
+    type Output = ProjectivePoint<C>;
+
+    fn mul(self, rhs: ProjectivePoint<C>) -> Self::Output {
+        let precomputation = rhs.mul_precompute();
+        rhs.mul_with_precomputation(self.0, precomputation)
+    }
+}
+
+#[allow(clippy::assertions_on_constants)]
+fn to_digits<C: Curve>(x: &C::ScalarField) -> Vec<u64> {
+    debug_assert!(
+        64 % WINDOW_BITS == 0,
+        "For simplicity, only power-of-two window sizes are handled for now"
+    );
+    let digits_per_u64 = 64 / WINDOW_BITS;
+    let mut digits = Vec::with_capacity(digits_per_scalar::<C>());
+    for limb in x.to_canonical_biguint().to_u64_digits() {
+        for j in 0..digits_per_u64 {
+            digits.push((limb >> (j * WINDOW_BITS) as u64) % BASE as u64);
+        }
+    }
+
+    digits
+}

src/curve/curve_summation.rs

@@ -0,0 +1,238 @@
+use alloc::vec;
+use alloc::vec::Vec;
+use core::iter::Sum;
+
+use plonky2::field::ops::Square;
+use plonky2::field::types::Field;
+
+use crate::curve::curve_types::{AffinePoint, Curve, ProjectivePoint};
+
+impl<C: Curve> Sum<AffinePoint<C>> for ProjectivePoint<C> {
+    fn sum<I: Iterator<Item = AffinePoint<C>>>(iter: I) -> ProjectivePoint<C> {
+        let points: Vec<_> = iter.collect();
+        affine_summation_best(points)
+    }
+}
+
+impl<C: Curve> Sum for ProjectivePoint<C> {
+    fn sum<I: Iterator<Item = ProjectivePoint<C>>>(iter: I) -> ProjectivePoint<C> {
+        iter.fold(ProjectivePoint::ZERO, |acc, x| acc + x)
+    }
+}
+
+pub fn affine_summation_best<C: Curve>(summation: Vec<AffinePoint<C>>) -> ProjectivePoint<C> {
+    let result = affine_multisummation_best(vec![summation]);
+    debug_assert_eq!(result.len(), 1);
+    result[0]
+}
+
+pub fn affine_multisummation_best<C: Curve>(
+    summations: Vec<Vec<AffinePoint<C>>>,
+) -> Vec<ProjectivePoint<C>> {
+    let pairwise_sums: usize = summations.iter().map(|summation| summation.len() / 2).sum();
+
+    // This threshold is chosen based on data from the summation benchmarks.
+    if pairwise_sums < 70 {
+        affine_multisummation_pairwise(summations)
+    } else {
+        affine_multisummation_batch_inversion(summations)
+    }
+}
+
+/// Adds each pair of points using an affine + affine = projective formula, then adds up the
+/// intermediate sums using a projective formula.
+pub fn affine_multisummation_pairwise<C: Curve>(
+    summations: Vec<Vec<AffinePoint<C>>>,
+) -> Vec<ProjectivePoint<C>> {
+    summations
+        .into_iter()
+        .map(affine_summation_pairwise)
+        .collect()
+}
+
+/// Adds each pair of points using an affine + affine = projective formula, then adds up the
+/// intermediate sums using a projective formula.
+pub fn affine_summation_pairwise<C: Curve>(points: Vec<AffinePoint<C>>) -> ProjectivePoint<C> {
+    let mut reduced_points: Vec<ProjectivePoint<C>> = Vec::new();
+    for chunk in points.chunks(2) {
+        match chunk.len() {
+            1 => reduced_points.push(chunk[0].to_projective()),
+            2 => reduced_points.push(chunk[0] + chunk[1]),
+            _ => panic!(),
+        }
+    }
+    // TODO: Avoid copying (deref)
+    reduced_points
+        .iter()
+        .fold(ProjectivePoint::ZERO, |sum, x| sum + *x)
+}
+
+/// Computes several summations of affine points by applying an affine group law, except that the
+/// divisions are batched via Montgomery's trick.
+pub fn affine_summation_batch_inversion<C: Curve>(
+    summation: Vec<AffinePoint<C>>,
+) -> ProjectivePoint<C> {
+    let result = affine_multisummation_batch_inversion(vec![summation]);
+    debug_assert_eq!(result.len(), 1);
+    result[0]
+}
+
+/// Computes several summations of affine points by applying an affine group law, except that the
+/// divisions are batched via Montgomery's trick.
+pub fn affine_multisummation_batch_inversion<C: Curve>(
+    summations: Vec<Vec<AffinePoint<C>>>,
+) -> Vec<ProjectivePoint<C>> {
+    let mut elements_to_invert = Vec::new();
+
+    // For each pair of points, (x1, y1) and (x2, y2), that we're going to add later, we want to
+    // invert either y (if the points are equal) or x1 - x2 (otherwise). We will use these later.
+    for summation in &summations {
+        let n = summation.len();
+        // The special case for n=0 is to avoid underflow.
+        let range_end = if n == 0 { 0 } else { n - 1 };
+
+        for i in (0..range_end).step_by(2) {
+            let p1 = summation[i];
+            let p2 = summation[i + 1];
+            let AffinePoint {
+                x: x1,
+                y: y1,
+                zero: zero1,
+            } = p1;
+            let AffinePoint {
+                x: x2,
+                y: _y2,
+                zero: zero2,
+            } = p2;
+
+            if zero1 || zero2 || p1 == -p2 {
+                // These are trivial cases where we won't need any inverse.
+            } else if p1 == p2 {
+                elements_to_invert.push(y1.double());
+            } else {
+                elements_to_invert.push(x1 - x2);
+            }
+        }
+    }
+
+    let inverses: Vec<C::BaseField> =
+        C::BaseField::batch_multiplicative_inverse(&elements_to_invert);
+
+    let mut all_reduced_points = Vec::with_capacity(summations.len());
+    let mut inverse_index = 0;
+    for summation in summations {
+        let n = summation.len();
+        let mut reduced_points = Vec::with_capacity((n + 1) / 2);
+
+        // The special case for n=0 is to avoid underflow.
+        let range_end = if n == 0 { 0 } else { n - 1 };
+
+        for i in (0..range_end).step_by(2) {
+            let p1 = summation[i];
+            let p2 = summation[i + 1];
+            let AffinePoint {
+                x: x1,
+                y: y1,
+                zero: zero1,
+            } = p1;
+            let AffinePoint {
+                x: x2,
+                y: y2,
+                zero: zero2,
+            } = p2;
+
+            let sum = if zero1 {
+                p2
+            } else if zero2 {
+                p1
+            } else if p1 == -p2 {
+                AffinePoint::ZERO
+            } else {
+                // It's a non-trivial case where we need one of the inverses we computed earlier.
+                let inverse = inverses[inverse_index];
+                inverse_index += 1;
+
+                if p1 == p2 {
+                    // This is the doubling case.
+                    let mut numerator = x1.square().triple();
+                    if C::A.is_nonzero() {
+                        numerator += C::A;
+                    }
+                    let quotient = numerator * inverse;
+                    let x3 = quotient.square() - x1.double();
+                    let y3 = quotient * (x1 - x3) - y1;
+                    AffinePoint::nonzero(x3, y3)
+                } else {
+                    // This is the general case. We use the incomplete addition formulas 4.3 and 4.4.
+                    let quotient = (y1 - y2) * inverse;
+                    let x3 = quotient.square() - x1 - x2;
+                    let y3 = quotient * (x1 - x3) - y1;
+                    AffinePoint::nonzero(x3, y3)
+                }
+            };
+            reduced_points.push(sum);
+        }
+
+        // If n is odd, the last point was not part of a pair.
+        if n % 2 == 1 {
+            reduced_points.push(summation[n - 1]);
+        }
+
+        all_reduced_points.push(reduced_points);
+    }
+
+    // We should have consumed all of the inverses from the batch computation.
+    debug_assert_eq!(inverse_index, inverses.len());
+
+    // Recurse with our smaller set of points.
+    affine_multisummation_best(all_reduced_points)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::curve::secp256k1::Secp256K1;
+
+    #[test]
+    fn test_pairwise_affine_summation() {
+        let g_affine = Secp256K1::GENERATOR_AFFINE;
+        let g2_affine = (g_affine + g_affine).to_affine();
+        let g3_affine = (g_affine + g_affine + g_affine).to_affine();
+        let g2_proj = g2_affine.to_projective();
+        let g3_proj = g3_affine.to_projective();
+        assert_eq!(
+            affine_summation_pairwise::<Secp256K1>(vec![g_affine, g_affine]),
+            g2_proj
+        );
+        assert_eq!(
+            affine_summation_pairwise::<Secp256K1>(vec![g_affine, g2_affine]),
+            g3_proj
+        );
+        assert_eq!(
+            affine_summation_pairwise::<Secp256K1>(vec![g_affine, g_affine, g_affine]),
+            g3_proj
+        );
+        assert_eq!(
+            affine_summation_pairwise::<Secp256K1>(vec![]),
+            ProjectivePoint::ZERO
+        );
+    }
+
+    #[test]
+    fn test_pairwise_affine_summation_batch_inversion() {
+        let g = Secp256K1::GENERATOR_AFFINE;
+        let g_proj = g.to_projective();
+        assert_eq!(
+            affine_summation_batch_inversion::<Secp256K1>(vec![g, g]),
+            g_proj + g_proj
+        );
+        assert_eq!(
+            affine_summation_batch_inversion::<Secp256K1>(vec![g, g, g]),
+            g_proj + g_proj + g_proj
+        );
+        assert_eq!(
+            affine_summation_batch_inversion::<Secp256K1>(vec![]),
+            ProjectivePoint::ZERO
+        );
+    }
+}

src/curve/curve_types.rs

@@ -0,0 +1,286 @@
+use alloc::vec::Vec;
+use core::fmt::Debug;
+use core::hash::{Hash, Hasher};
+use core::ops::Neg;
+
+use plonky2::field::ops::Square;
+use plonky2::field::types::{Field, PrimeField};
+use serde::{Deserialize, Serialize};
+
+// To avoid implementation conflicts from associated types,
+// see https://github.com/rust-lang/rust/issues/20400
+pub struct CurveScalar<C: Curve>(pub <C as Curve>::ScalarField);
+
+/// A short Weierstrass curve.
+pub trait Curve: 'static + Sync + Sized + Copy + Debug {
+    type BaseField: PrimeField;
+    type ScalarField: PrimeField;
+
+    const A: Self::BaseField;
+    const B: Self::BaseField;
+
+    const GENERATOR_AFFINE: AffinePoint<Self>;
+
+    const GENERATOR_PROJECTIVE: ProjectivePoint<Self> = ProjectivePoint {
+        x: Self::GENERATOR_AFFINE.x,
+        y: Self::GENERATOR_AFFINE.y,
+        z: Self::BaseField::ONE,
+    };
+
+    fn convert(x: Self::ScalarField) -> CurveScalar<Self> {
+        CurveScalar(x)
+    }
+
+    fn is_safe_curve() -> bool {
+        // Added additional check to prevent using vulnerabilties in case a discriminant is equal to 0.
+        (Self::A.cube().double().double() + Self::B.square().triple().triple().triple())
+            .is_nonzero()
+    }
+}
+
+/// A point on a short Weierstrass curve, represented in affine coordinates.
+#[derive(Copy, Clone, Debug, Deserialize, Serialize)]
+pub struct AffinePoint<C: Curve> {
+    pub x: C::BaseField,
+    pub y: C::BaseField,
+    pub zero: bool,
+}
+
+impl<C: Curve> AffinePoint<C> {
+    pub const ZERO: Self = Self {
+        x: C::BaseField::ZERO,
+        y: C::BaseField::ZERO,
+        zero: true,
+    };
+
+    pub fn nonzero(x: C::BaseField, y: C::BaseField) -> Self {
+        let point = Self { x, y, zero: false };
+        debug_assert!(point.is_valid());
+        point
+    }
+
+    pub fn is_valid(&self) -> bool {
+        let Self { x, y, zero } = *self;
+        zero || y.square() == x.cube() + C::A * x + C::B
+    }
+
+    pub fn to_projective(&self) -> ProjectivePoint<C> {
+        let Self { x, y, zero } = *self;
+        let z = if zero {
+            C::BaseField::ZERO
+        } else {
+            C::BaseField::ONE
+        };
+
+        ProjectivePoint { x, y, z }
+    }
+
+    pub fn batch_to_projective(affine_points: &[Self]) -> Vec<ProjectivePoint<C>> {
+        affine_points.iter().map(Self::to_projective).collect()
+    }
+
+    #[must_use]
+    pub fn double(&self) -> Self {
+        let AffinePoint { x: x1, y: y1, zero } = *self;
+
+        if zero {
+            return AffinePoint::ZERO;
+        }
+
+        let double_y = y1.double();
+        let inv_double_y = double_y.inverse(); // (2y)^(-1)
+        let triple_xx = x1.square().triple(); // 3x^2
+        let lambda = (triple_xx + C::A) * inv_double_y;
+        let x3 = lambda.square() - self.x.double();
+        let y3 = lambda * (x1 - x3) - y1;
+
+        Self {
+            x: x3,
+            y: y3,
+            zero: false,
+        }
+    }
+}
+
+impl<C: Curve> PartialEq for AffinePoint<C> {
+    fn eq(&self, other: &Self) -> bool {
+        let AffinePoint {
+            x: x1,
+            y: y1,
+            zero: zero1,
+        } = *self;
+        let AffinePoint {
+            x: x2,
+            y: y2,
+            zero: zero2,
+        } = *other;
+        if zero1 || zero2 {
+            return zero1 == zero2;
+        }
+        x1 == x2 && y1 == y2
+    }
+}
+
+impl<C: Curve> Eq for AffinePoint<C> {}
+
+impl<C: Curve> Hash for AffinePoint<C> {
+    fn hash<H: Hasher>(&self, state: &mut H) {
+        if self.zero {
+            self.zero.hash(state);
+        } else {
+            self.x.hash(state);
+            self.y.hash(state);
+        }
+    }
+}
+
+/// A point on a short Weierstrass curve, represented in projective coordinates.
+#[derive(Copy, Clone, Debug)]
+pub struct ProjectivePoint<C: Curve> {
+    pub x: C::BaseField,
+    pub y: C::BaseField,
+    pub z: C::BaseField,
+}
+
+impl<C: Curve> ProjectivePoint<C> {
+    pub const ZERO: Self = Self {
+        x: C::BaseField::ZERO,
+        y: C::BaseField::ONE,
+        z: C::BaseField::ZERO,
+    };
+
+    pub fn nonzero(x: C::BaseField, y: C::BaseField, z: C::BaseField) -> Self {
+        let point = Self { x, y, z };
+        debug_assert!(point.is_valid());
+        point
+    }
+
+    pub fn is_valid(&self) -> bool {
+        let Self { x, y, z } = *self;
+        z.is_zero() || y.square() * z == x.cube() + C::A * x * z.square() + C::B * z.cube()
+    }
+
+    pub fn to_affine(&self) -> AffinePoint<C> {
+        let Self { x, y, z } = *self;
+        if z == C::BaseField::ZERO {
+            AffinePoint::ZERO
+        } else {
+            let z_inv = z.inverse();
+            AffinePoint::nonzero(x * z_inv, y * z_inv)
+        }
+    }
+
+    pub fn batch_to_affine(proj_points: &[Self]) -> Vec<AffinePoint<C>> {
+        let n = proj_points.len();
+        let zs: Vec<C::BaseField> = proj_points.iter().map(|pp| pp.z).collect();
+        let z_invs = C::BaseField::batch_multiplicative_inverse(&zs);
+
+        let mut result = Vec::with_capacity(n);
+        for i in 0..n {
+            let Self { x, y, z } = proj_points[i];
+            result.push(if z == C::BaseField::ZERO {
+                AffinePoint::ZERO
+            } else {
+                let z_inv = z_invs[i];
+                AffinePoint::nonzero(x * z_inv, y * z_inv)
+            });
+        }
+        result
+    }
+
+    // From https://www.hyperelliptic.org/EFD/g1p/data/shortw/projective/doubling/dbl-2007-bl
+    #[must_use]
+    pub fn double(&self) -> Self {
+        let Self { x, y, z } = *self;
+        if z == C::BaseField::ZERO {
+            return ProjectivePoint::ZERO;
+        }
+
+        let xx = x.square();
+        let zz = z.square();
+        let mut w = xx.triple();
+        if C::A.is_nonzero() {
+            w += C::A * zz;
+        }
+        let s = y.double() * z;
+        let r = y * s;
+        let rr = r.square();
+        let b = (x + r).square() - (xx + rr);
+        let h = w.square() - b.double();
+        let x3 = h * s;
+        let y3 = w * (b - h) - rr.double();
+        let z3 = s.cube();
+        Self {
+            x: x3,
+            y: y3,
+            z: z3,
+        }
+    }
+
+    pub fn add_slices(a: &[Self], b: &[Self]) -> Vec<Self> {
+        assert_eq!(a.len(), b.len());
+        a.iter()
+            .zip(b.iter())
+            .map(|(&a_i, &b_i)| a_i + b_i)
+            .collect()
+    }
+
+    #[must_use]
+    pub fn neg(&self) -> Self {
+        Self {
+            x: self.x,
+            y: -self.y,
+            z: self.z,
+        }
+    }
+}
+
+impl<C: Curve> PartialEq for ProjectivePoint<C> {
+    fn eq(&self, other: &Self) -> bool {
+        let ProjectivePoint {
+            x: x1,
+            y: y1,
+            z: z1,
+        } = *self;
+        let ProjectivePoint {
+            x: x2,
+            y: y2,
+            z: z2,
+        } = *other;
+        if z1 == C::BaseField::ZERO || z2 == C::BaseField::ZERO {
+            return z1 == z2;
+        }
+
+        // We want to compare (x1/z1, y1/z1) == (x2/z2, y2/z2).
+        // But to avoid field division, it is better to compare (x1*z2, y1*z2) == (x2*z1, y2*z1).
+        x1 * z2 == x2 * z1 && y1 * z2 == y2 * z1
+    }
+}
+
+impl<C: Curve> Eq for ProjectivePoint<C> {}
+
+impl<C: Curve> Neg for AffinePoint<C> {
+    type Output = AffinePoint<C>;
+
+    fn neg(self) -> Self::Output {
+        let AffinePoint { x, y, zero } = self;
+        AffinePoint { x, y: -y, zero }
+    }
+}
+
+impl<C: Curve> Neg for ProjectivePoint<C> {
+    type Output = ProjectivePoint<C>;
+
+    fn neg(self) -> Self::Output {
+        let ProjectivePoint { x, y, z } = self;
+        ProjectivePoint { x, y: -y, z }
+    }
+}
+
+pub fn base_to_scalar<C: Curve>(x: C::BaseField) -> C::ScalarField {
+    C::ScalarField::from_noncanonical_biguint(x.to_canonical_biguint())
+}
+
+pub fn scalar_to_base<C: Curve>(x: C::ScalarField) -> C::BaseField {
+    C::BaseField::from_noncanonical_biguint(x.to_canonical_biguint())
+}

src/curve/ecdsa.rs

@@ -0,0 +1,84 @@
+use plonky2::field::types::{Field, Sample};
+use serde::{Deserialize, Serialize};
+
+use crate::curve::curve_msm::msm_parallel;
+use crate::curve::curve_types::{base_to_scalar, AffinePoint, Curve, CurveScalar};
+
+#[derive(Copy, Clone, Debug, Deserialize, Eq, Hash, PartialEq, Serialize)]
+pub struct ECDSASignature<C: Curve> {
+    pub r: C::ScalarField,
+    pub s: C::ScalarField,
+}
+
+#[derive(Copy, Clone, Debug, Deserialize, Eq, Hash, PartialEq, Serialize)]
+pub struct ECDSASecretKey<C: Curve>(pub C::ScalarField);
+
+impl<C: Curve> ECDSASecretKey<C> {
+    pub fn to_public(&self) -> ECDSAPublicKey<C> {
+        ECDSAPublicKey((CurveScalar(self.0) * C::GENERATOR_PROJECTIVE).to_affine())
+    }
+}
+
+#[derive(Copy, Clone, Debug, Deserialize, Eq, Hash, PartialEq, Serialize)]
+pub struct ECDSAPublicKey<C: Curve>(pub AffinePoint<C>);
+
+pub fn sign_message<C: Curve>(msg: C::ScalarField, sk: ECDSASecretKey<C>) -> ECDSASignature<C> {
+    let (k, rr) = {
+        let mut k = C::ScalarField::rand();
+        let mut rr = (CurveScalar(k) * C::GENERATOR_PROJECTIVE).to_affine();
+        while rr.x == C::BaseField::ZERO {
+            k = C::ScalarField::rand();
+            rr = (CurveScalar(k) * C::GENERATOR_PROJECTIVE).to_affine();
+        }
+        (k, rr)
+    };
+    let r = base_to_scalar::<C>(rr.x);
+
+    let s = k.inverse() * (msg + r * sk.0);
+
+    ECDSASignature { r, s }
+}
+
+pub fn verify_message<C: Curve>(
+    msg: C::ScalarField,
+    sig: ECDSASignature<C>,
+    pk: ECDSAPublicKey<C>,
+) -> bool {
+    let ECDSASignature { r, s } = sig;
+
+    assert!(pk.0.is_valid());
+
+    let c = s.inverse();
+    let u1 = msg * c;
+    let u2 = r * c;
+
+    let g = C::GENERATOR_PROJECTIVE;
+    let w = 5; // Experimentally fastest
+    let point_proj = msm_parallel(&[u1, u2], &[g, pk.0.to_projective()], w);
+    let point = point_proj.to_affine();
+
+    let x = base_to_scalar::<C>(point.x);
+    r == x
+}
+
+#[cfg(test)]
+mod tests {
+    use plonky2::field::secp256k1_scalar::Secp256K1Scalar;
+    use plonky2::field::types::Sample;
+
+    use crate::curve::ecdsa::{sign_message, verify_message, ECDSASecretKey};
+    use crate::curve::secp256k1::Secp256K1;
+
+    #[test]
+    fn test_ecdsa_native() {
+        type C = Secp256K1;
+
+        let msg = Secp256K1Scalar::rand();
+        let sk = ECDSASecretKey::<C>(Secp256K1Scalar::rand());
+        let pk = sk.to_public();
+
+        let sig = sign_message(msg, sk);
+        let result = verify_message(msg, sig, pk);
+        assert!(result);
+    }
+}

src/curve/glv.rs

@@ -0,0 +1,140 @@
+use num::rational::Ratio;
+use num::BigUint;
+use plonky2::field::secp256k1_base::Secp256K1Base;
+use plonky2::field::secp256k1_scalar::Secp256K1Scalar;
+use plonky2::field::types::{Field, PrimeField};
+
+use crate::curve::curve_msm::msm_parallel;
+use crate::curve::curve_types::{AffinePoint, ProjectivePoint};
+use crate::curve::secp256k1::Secp256K1;
+
+pub const GLV_BETA: Secp256K1Base = Secp256K1Base([
+    13923278643952681454,
+    11308619431505398165,
+    7954561588662645993,
+    8856726876819556112,
+]);
+
+pub const GLV_S: Secp256K1Scalar = Secp256K1Scalar([
+    16069571880186789234,
+    1310022930574435960,
+    11900229862571533402,
+    6008836872998760672,
+]);
+
+const A1: Secp256K1Scalar = Secp256K1Scalar([16747920425669159701, 3496713202691238861, 0, 0]);
+
+const MINUS_B1: Secp256K1Scalar =
+    Secp256K1Scalar([8022177200260244675, 16448129721693014056, 0, 0]);
+
+const A2: Secp256K1Scalar = Secp256K1Scalar([6323353552219852760, 1498098850674701302, 1, 0]);
+
+const B2: Secp256K1Scalar = Secp256K1Scalar([16747920425669159701, 3496713202691238861, 0, 0]);
+
+/// Algorithm 15.41 in Handbook of Elliptic and Hyperelliptic Curve Cryptography.
+/// Decompose a scalar `k` into two small scalars `k1, k2` with `|k1|, |k2| < √p` that satisfy
+/// `k1 + s * k2 = k`.
+/// Returns `(|k1|, |k2|, k1 < 0, k2 < 0)`.
+pub fn decompose_secp256k1_scalar(
+    k: Secp256K1Scalar,
+) -> (Secp256K1Scalar, Secp256K1Scalar, bool, bool) {
+    let p = Secp256K1Scalar::order();
+    let c1_biguint = Ratio::new(
+        B2.to_canonical_biguint() * k.to_canonical_biguint(),
+        p.clone(),
+    )
+    .round()
+    .to_integer();
+    let c1 = Secp256K1Scalar::from_noncanonical_biguint(c1_biguint);
+    let c2_biguint = Ratio::new(
+        MINUS_B1.to_canonical_biguint() * k.to_canonical_biguint(),
+        p.clone(),
+    )
+    .round()
+    .to_integer();
+    let c2 = Secp256K1Scalar::from_noncanonical_biguint(c2_biguint);
+
+    let k1_raw = k - c1 * A1 - c2 * A2;
+    let k2_raw = c1 * MINUS_B1 - c2 * B2;
+    debug_assert!(k1_raw + GLV_S * k2_raw == k);
+
+    let two = BigUint::from_slice(&[2]);
+    let k1_neg = k1_raw.to_canonical_biguint() > p.clone() / two.clone();
+    let k1 = if k1_neg {
+        Secp256K1Scalar::from_noncanonical_biguint(p.clone() - k1_raw.to_canonical_biguint())
+    } else {
+        k1_raw
+    };
+    let k2_neg = k2_raw.to_canonical_biguint() > p.clone() / two;
+    let k2 = if k2_neg {
+        Secp256K1Scalar::from_noncanonical_biguint(p - k2_raw.to_canonical_biguint())
+    } else {
+        k2_raw
+    };
+
+    (k1, k2, k1_neg, k2_neg)
+}
+
+/// See Section 15.2.1 in Handbook of Elliptic and Hyperelliptic Curve Cryptography.
+/// GLV scalar multiplication `k * P = k1 * P + k2 * psi(P)`, where `k = k1 + s * k2` is the
+/// decomposition computed in `decompose_secp256k1_scalar(k)` and `psi` is the Secp256k1
+/// endomorphism `psi: (x, y) |-> (beta * x, y)` equivalent to scalar multiplication by `s`.
+pub fn glv_mul(p: ProjectivePoint<Secp256K1>, k: Secp256K1Scalar) -> ProjectivePoint<Secp256K1> {
+    let (k1, k2, k1_neg, k2_neg) = decompose_secp256k1_scalar(k);
+
+    let p_affine = p.to_affine();
+    let sp = AffinePoint::<Secp256K1> {
+        x: p_affine.x * GLV_BETA,
+        y: p_affine.y,
+        zero: p_affine.zero,
+    };
+
+    let first = if k1_neg { p.neg() } else { p };
+    let second = if k2_neg {
+        sp.to_projective().neg()
+    } else {
+        sp.to_projective()
+    };
+
+    msm_parallel(&[k1, k2], &[first, second], 5)
+}
+
+#[cfg(test)]
+mod tests {
+    use anyhow::Result;
+    use plonky2::field::secp256k1_scalar::Secp256K1Scalar;
+    use plonky2::field::types::{Field, Sample};
+
+    use crate::curve::curve_types::{Curve, CurveScalar};
+    use crate::curve::glv::{decompose_secp256k1_scalar, glv_mul, GLV_S};
+    use crate::curve::secp256k1::Secp256K1;
+
+    #[test]
+    fn test_glv_decompose() -> Result<()> {
+        let k = Secp256K1Scalar::rand();
+        let (k1, k2, k1_neg, k2_neg) = decompose_secp256k1_scalar(k);
+        let one = Secp256K1Scalar::ONE;
+        let m1 = if k1_neg { -one } else { one };
+        let m2 = if k2_neg { -one } else { one };
+
+        assert!(k1 * m1 + GLV_S * k2 * m2 == k);
+
+        Ok(())
+    }
+
+    #[test]
+    fn test_glv_mul() -> Result<()> {
+        for _ in 0..20 {
+            let k = Secp256K1Scalar::rand();
+
+            let p = CurveScalar(Secp256K1Scalar::rand()) * Secp256K1::GENERATOR_PROJECTIVE;
+
+            let kp = CurveScalar(k) * p;
+            let glv = glv_mul(p, k);
+
+            assert!(kp == glv);
+        }
+
+        Ok(())
+    }
+}

src/curve/mod.rs

@@ -0,0 +1,8 @@
+pub mod curve_adds;
+pub mod curve_msm;
+pub mod curve_multiplication;
+pub mod curve_summation;
+pub mod curve_types;
+pub mod ecdsa;
+pub mod glv;
+pub mod secp256k1;

src/curve/secp256k1.rs

@@ -0,0 +1,100 @@
+use plonky2::field::secp256k1_base::Secp256K1Base;
+use plonky2::field::secp256k1_scalar::Secp256K1Scalar;
+use plonky2::field::types::Field;
+use serde::{Deserialize, Serialize};
+
+use crate::curve::curve_types::{AffinePoint, Curve};
+
+#[derive(Debug, Copy, Clone, Deserialize, Eq, Hash, PartialEq, Serialize)]
+pub struct Secp256K1;
+
+impl Curve for Secp256K1 {
+    type BaseField = Secp256K1Base;
+    type ScalarField = Secp256K1Scalar;
+
+    const A: Secp256K1Base = Secp256K1Base::ZERO;
+    const B: Secp256K1Base = Secp256K1Base([7, 0, 0, 0]);
+    const GENERATOR_AFFINE: AffinePoint<Self> = AffinePoint {
+        x: SECP256K1_GENERATOR_X,
+        y: SECP256K1_GENERATOR_Y,
+        zero: false,
+    };
+}
+
+// 55066263022277343669578718895168534326250603453777594175500187360389116729240
+const SECP256K1_GENERATOR_X: Secp256K1Base = Secp256K1Base([
+    0x59F2815B16F81798,
+    0x029BFCDB2DCE28D9,
+    0x55A06295CE870B07,
+    0x79BE667EF9DCBBAC,
+]);
+
+/// 32670510020758816978083085130507043184471273380659243275938904335757337482424
+const SECP256K1_GENERATOR_Y: Secp256K1Base = Secp256K1Base([
+    0x9C47D08FFB10D4B8,
+    0xFD17B448A6855419,
+    0x5DA4FBFC0E1108A8,
+    0x483ADA7726A3C465,
+]);
+
+#[cfg(test)]
+mod tests {
+    use num::BigUint;
+    use plonky2::field::secp256k1_scalar::Secp256K1Scalar;
+    use plonky2::field::types::{Field, PrimeField};
+
+    use crate::curve::curve_types::{AffinePoint, Curve, ProjectivePoint};
+    use crate::curve::secp256k1::Secp256K1;
+
+    #[test]
+    fn test_generator() {
+        let g = Secp256K1::GENERATOR_AFFINE;
+        assert!(g.is_valid());
+
+        let neg_g = AffinePoint::<Secp256K1> {
+            x: g.x,
+            y: -g.y,
+            zero: g.zero,
+        };
+        assert!(neg_g.is_valid());
+    }
+
+    #[test]
+    fn test_naive_multiplication() {
+        let g = Secp256K1::GENERATOR_PROJECTIVE;
+        let ten = Secp256K1Scalar::from_canonical_u64(10);
+        let product = mul_naive(ten, g);
+        let sum = g + g + g + g + g + g + g + g + g + g;
+        assert_eq!(product, sum);
+    }
+
+    #[test]
+    fn test_g1_multiplication() {
+        let lhs = Secp256K1Scalar::from_noncanonical_biguint(BigUint::from_slice(&[
+            1111, 2222, 3333, 4444, 5555, 6666, 7777, 8888,
+        ]));
+        assert_eq!(
+            Secp256K1::convert(lhs) * Secp256K1::GENERATOR_PROJECTIVE,
+            mul_naive(lhs, Secp256K1::GENERATOR_PROJECTIVE)
+        );
+    }
+
+    /// A simple, somewhat inefficient implementation of multiplication which is used as a reference
+    /// for correctness.
+    fn mul_naive(
+        lhs: Secp256K1Scalar,
+        rhs: ProjectivePoint<Secp256K1>,
+    ) -> ProjectivePoint<Secp256K1> {
+        let mut g = rhs;
+        let mut sum = ProjectivePoint::ZERO;
+        for limb in lhs.to_canonical_biguint().to_u64_digits().iter() {
+            for j in 0..64 {
+                if (limb >> j & 1u64) != 0u64 {
+                    sum = sum + g;
+                }
+                g = g.double();
+            }
+        }
+        sum
+    }
+}