Experimenting with elliptic curve implementation in rust

Recently, I have been following the awesome tutorial describing the internal details of zkSNARK system PLONK called PLONK by hand by Joshua Fitzgerald (Metastate team). The rust implementation that closely follows this is given by Adrià Massanet in PLONK by fingers.

For brevity, for rest of this document, PLONK by hand and PLONK by fingers will be abbreviated as PBH and PBF respectively. While PBF is a bigger implementation of the zkSNARK protocol described in PBH, what we concentrate on in this document is how in rust implementation of elliptic curves go.

The purpose of this document remains me describing as a non-cryptography expert the interpretations I draw from what I see. Things might be hand-wavy and could be wrong, but the point remains for this to be a good entry point for others who like to venture into similar territories, which may be slightly wrong at first. May this be a quick disclaimer for readers :P .

What are elliptic curves and how do they feature here

Elliptic curves are an important cryptographic primitive in modern times. What makes them valuable in where they are used are generally two properties (one optional, but used in PLONK).

1. Discrete Log Problem: Points on elliptic curves form a group
   $G$. Further, this group has property of being finite and cyclic. What that entails is that there exists some element
   $g\in G$ (the generator) that all other elements in
   $G$ can be constructed by raising
   $g$ to some power
   $x$ as
   $g^x$. However once you have
   $g$ and
   $g^x$, there exists no efficient "classical computer" algorithm that could deduce
   $x$.
2. EC Pairing(Optional): There exists some bilinear map (a function that takes in two elements from different groups (for EC)
   $G_1$ and
   $G_2$ and presents an element of a third group
   $G_T$)
   $e:G_1\times G_2\to G_T$ such that following property holds:
   $$e(aP,bQ)=e(P,Q{)}^{ab}=e(aP,Q{)}^b=e(P,aQ{)}^b=e(P,abQ)=e(abP,Q)$$ where
   $P$ and
   $Q$ are elements in
   $G_1$ and
   $G_2$ respectively, while
   $a$ and
   $b$ are integers. Why this is useful can be seen later.

See more properties of bilinear maps in appendix here.

Onto the first trait: FiniteField

For coding up the elliptic curve, we need to describe a field. This should be a trait since we are going to define functions that a concrete implementation should provide.

Setup a new project using cargo init --lib pbf, make lib.rs under src/, add

pub mod elliptic_curve

Create elliptic_curve.rs under src/ and add:

pub trait FiniteField {
    /// Order or cardinality of the field
    fn order() -> u64;

    /// Gives the zero element and one element
    /// the multiplicative and additive identities
    fn one() -> Self;
    fn zero() -> Self;

    /// Check whether element is identity
    /// element
    fn is_one(&self) -> bool;
    fn is_zero(&self) -> bool;
}

Running cargo build works for now. Good start! Let's start adding some real function that we will use in real life w.r.t. Finite Fields.

Adding

pub trait FiniteField {
    ...
    /// Calculates the inverse of a non-zeo
    /// element. Returns `None` if is_zero
    fn inverse(&self) -> Option<Self>;
}

makes the compiler unhappy and it returns:

error[E0277]: the size for values of type `Self` cannot be known at compilation time
   --> src/elliptic_curve.rs:23:26
    |
23  |     fn inverse(&self) -> Option<Self>;
    |                          ^^^^^^^^^^^^ doesn't have a size known at compile-time
    |
note: required by a bound in `Option`

What it means is to be able to wrap the type T in Option<T>, the amount of memory required by T should be known at compile time. If we think closely, that is true for any finite field we may implement because finite fields are operated on modulo some large prime

$p$ which is constant for the field, making the memory footprint of any field element known at compile time (as at maximum it may take as much memory as required to represent

$p-1$).

So, we add:

pub trait FiniteField:
    Sized
{
    ...
}

telling that whatever implements FiniteField has to have the Sized trait. This solves current compilation issues. We add additional operations which we will use later such as mathematical operation like negation, addition, multiplication. We also add operations for raising element to some power, finding inverse etc.

use std::{
    fmt::{Debug, Display},
    ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub, SubAssign},
};

pub trait FiniteField:
    Sized               // Size known at compile time
    + From<u64>         // Can create an element from u64
    + Debug + Display   // For debugging and printing
    + Copy              // Easy to copy during assignment
    // Mathemetical Operations ----
    + Neg<Output=Self>
    + Add<Self, Output=Self> + AddAssign<Self> // For + and +=
    + Sub<Self, Output=Self> + SubAssign<Self> // For - and -=
    + Mul<Self, Output=Self> + MulAssign<Self> // For * and *=
    + Div<Self, Output=Option<Self>>
    // Mathemetical Operations ----
{
    /// Order or cardinality of the field
    fn order() -> u64;

    /// Gives the zero element and one element
    /// the multiplicative and additive identities
    fn one() -> Self;
    fn zero() -> Self;

    /// Check whether element is identity
    /// element
    fn is_one(&self) -> bool;
    fn is_zero(&self) -> bool;

    /// Calculates the inverse of a non-zeo
    /// element. Returns `None` if is_zero
    fn inverse(&self) -> Option<Self>;

    /// Raise element to some power
    fn pow(&self, exp: u64) -> Self;

    /// Export to u64
    fn as_u64(&self) -> u64;
}

When defining for traits such as std::ops::Add, the binary operation will have a right hand side that you like to add with. Such type is the generic parameter that needs to be given (defaults to Self). While generally the other element that needs to be added is of same type, in certain instances such an operation may not make sense. For example, only Duration should be added to current time clock; a clock may not be added to a clock.

During instances where such weird behavior is not seen, Add<Self, Output=Self> can just be written as Add. The Output is an associated type here. What it means is that only one implementation of trait Add<rhstype> is expected for some rhstype and no two traits Add<T, Output=K>, Add<T, Output=P> are allowed at the same time.

The trait From<u64> is interesting here, it allows creation of the object from some u64 type, the function as_u64 is the complementary exporting function here.

A note on blanket implementation: If you hover over From<u64>, you will see the following:

Used to do value-to-value conversions while consuming the input value. It is the reciprocal of Into.

One should always prefer implementing From over Into because implementing From automatically provides one with an implementation of Into thanks to the blanket implementation in the standard library.
...

Another trait exists pub trait Into<T>: Sized in rust core library. With From implemented for a type, the following call is satisfied

let a: FiniteField = FiniteField::from(41u64);

With Into<T> implemented, the following call is satisfied:

let a: FiniteField = 41u64.into()

However, the default implementation of Into<T> is provided by rust if From<U> is defined. So while implementing From<u64> makes both the above possible, with just Into<T>, the former may not be possible.

The two curves

$G_1$ and

$G_2$ and pairing function

$e$

We create traits for points on both the curves

$G_1$ and

$G_2$. Since most operation is going to happen in

$G_1$ and we use

$G_2$ just to allow for pairing, we define the following traits:

pub trait G1Point:
    Copy                // Representation small enough for efficient copy
    + Debug + Display   // For debugging and printing
    + PartialEq         // Allow for equality testing but may not be transitively equal
    + Hash              // A hashable type     
    + Ord               // Total order 
    // Mathemetical Operations ----
    + Neg<Output=Self>
    + Sub<Output=Self>
    + Add<Output=Self>
    + Mul<Self::SubField, Output = Self>
    // Mathemetical Operations ----
{
    type Field: FiniteField;
    type SubField: FiniteField;

    /// Constructor
    fn new(x: Self::Field, y: Self::Field) -> Self;

    /// Provide generator element
    fn generator() -> Self;

    /// Give the identity element
    fn identity() -> Self;
    /// Check wheter current point is an identity point
    fn is_identity(&self) -> bool;

    /// Provide the x and y coordinates of the point
    fn x(&self) -> &Self::Field;
    fn y(&self) -> &Self::Field;
}

Two associated types exist for G1Point namely Field and SubField. (TODO: Describe why).

Similarly, we describe:

pub trait G2Point:
    Copy                // Representation small enough for efficient copy
    + Debug + Display   // For debugging and printing
    + PartialEq         // Allow for equality testing but may not be transitively equal
    // Mathemetical Operations ----
    + Neg<Output=Self>
    + Sub<Output=Self>
    + Add<Output=Self>
    + Mul<Self::SubField, Output = Self>
    // Mathemetical Operations ----
{
    type Field: FiniteField;
    type SubField: FiniteField;

    /// Constructor
    fn new(x: Self::Field, y: Self::Field) -> Self;

    /// Provide generator element
    fn generator() -> Self;

    /// Provide the x and y coordinates of the point
    fn x(&self) -> &Self::Field;
    fn y(&self) -> &Self::Field;

    /// Get the embedding degree of the curve
    fn embedding_degree() -> u64;
}

For pairing, we define the target groupt

$G_T$ and a pairing trait as follows:

pub trait GTPoint: Display + Copy + PartialEq + Mul {
    type S;
    fn pow(&self, n: Self::S) -> Self;
}

pub trait Pairing {
    type G1: G1Point;
    type G2: G2Point;
    type GT: GTPoint;

    fn pairing(p: Self::G1, q: Self::G2) -> Self::GT;
}

Implementing a concrete field: the u64 field

Much of the implementation for u64 field is simple to understand. We will build our actual curves

$G_1$,

$G_2$ and

$G_T$ on top of this concrete implementation later. Create u64_field.rs under src/

use std::{
    fmt::{self, Display},
    ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub, SubAssign},
};

use crate::{elliptic_curve::FiniteField, extended_euclid};

#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, PartialOrd, Ord)]
pub struct U64Field<const P: u64>(u64);

impl<const P: u64> U64Field<P> {
    fn new(n: u64) -> Self {
        Self(n % P)
    }
}

impl<const P: u64> FiniteField for U64Field<P> {
    fn order() -> u64 {
        P
    }

    fn one() -> Self {
        Self::from(1u64)
    }
    fn zero() -> Self {
        Self::from(0u64)
    }

    fn is_one(&self) -> bool {
        self.0 == 1
    }
    fn is_zero(&self) -> bool {
        self.0 == 0
    }

    fn inverse(&self) -> Option<Self> {
        let (gcd, x, _) = extended_euclid::eea_recursive(self.0 as i64, P as i64);
        if gcd == 1 {
            Some(Self::new(x as u64))
        } else {
            None
        }
    }

    fn pow(&self, mut exp: u64) -> Self {
        let mut result = Self::one();
        let mut base = *self;

        while exp > 0 {
            if exp % 2 == 1 {
                result *= base;
            }
            exp >>= 1;
            base *= base;
        }

        result
    }

    fn as_u64(&self) -> u64 {
        self.0
    }
}

impl<const P: u64> From<u64> for U64Field<P> {
    fn from(n: u64) -> Self {
        Self::new(n)
    }
}

impl<const P: u64> Display for U64Field<P> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "({})", self.0)
    }
}

impl<const P: u64> Add for U64Field<P> {
    type Output = Self;
    fn add(self, rhs: Self) -> Self::Output {
        Self::Output::new(self.0 + rhs.0)
    }
}

impl<const P: u64> AddAssign for U64Field<P> {
    fn add_assign(&mut self, rhs: Self) {
        self.0 = self.0 + rhs.0
    }
}

impl<const P: u64> Neg for U64Field<P> {
    type Output = Self;
    fn neg(self) -> Self::Output {
        Self::Output::new(P - self.0)
    }
}

impl<const P: u64> Sub for U64Field<P> {
    type Output = Self;
    fn sub(self, rhs: Self) -> Self::Output {
        self + -rhs
    }
}

impl<const P: u64> SubAssign for U64Field<P> {
    fn sub_assign(&mut self, rhs: Self) {
        *self += -rhs;
    }
}

impl<const P: u64> Mul for U64Field<P> {
    type Output = Self;
    fn mul(self, rhs: Self) -> Self::Output {
        Self::new(self.0 * rhs.0)
    }
}

impl<const P: u64> MulAssign for U64Field<P> {
    fn mul_assign(&mut self, rhs: Self) {
        self.0 = (self.0 * rhs.0) % P
    }
}

impl<const P: u64> Div for U64Field<P> {
    type Output = Option<Self>;
    fn div(self, rhs: Self) -> Self::Output {
        rhs.inverse().map(|e| self * e)
    }
}

Finding the inverse: Extended euclidean algorithm

What does it mean for an element to have an inverse in our u64 field? We deem this property to hold true for any non-zero element

$a$ in the field:

$$a\times a^{-1}=a^{-1}\times a=1$$ where

$1$ is the multiplicative identity. Of course, all of this is done mod

$p$.

To find the inverse of some element

$a$, we resort to first find coefficients that satisfy Bézout's identity for which we need to deploy extended euclidean algorithm for element

$a$ and field prime

$p$.

Bézout's identity: Let

$a$ and

$b$ be integers with greatest common divisor

$d$. Then there exist integers

$x$ and

$y$ such that

$ax+by=d$. Moreover, the integers of the form

$az+bt$ are exactly the multiples of

$d$.

This excellent video describes how extended euclidean algorithm can be used to get these coefficients

$x$ and

$y$.

Once we have the coefficients

$x$ and

$y$, we see using following (assuming

$p$ to be prime and hence

$a$ and

$p$ are coprimes):

$$\begin{gathered} xa+yp=1 \\ (xa+yp)(\mathrm{mod}p)\equiv 1(\mathrm{mod}p) \\ xa(\mathrm{mod}p)\equiv 1(\mathrm{mod}p) \\ x=a^{-1} \end{gathered}$$ Hence we see

$x$ is inverse of

$a$ and we don't really need

$y$ for calculating the inverse.

Let's begin by implementing this algorithm.

Base case: Let's say we are asked to find

$x$ and

$y$ for

$(a,0)$ where

$d=a$, the following is seen:

$a\times 1+b\times 0=a$. Any GCD calculation always ends up at this base case for e.g.:

$gcd(22,6)=gcd(6,4)=gcd(4,2)=gcd(2,0)$. At the end, 2 is the gcd,

$2\times 1+0\times 0=2$

Recursive case: In case

$b\ne 0$, we know if we decompose

$(a,b)$ as

$a=bq+r$,

$gcd(a,b)=gcd(b,r)=g$. Let's say there exists some "extended gcd" algorithm that replies a 3-tuple containing

$(g,x,y)$ for some given

$(a,b)$;

$$\begin{gathered} g=xb+yr \\ g=xb+y(a-bq) \\ g=xb+ya-ybq \\ g=ya+(x-yq)b \end{gathered}$$

Hence, the algorithm can be written as (in extended_euclid.rs under src/):

pub fn eea_recursive(a: i64, b: i64) -> (i64, i64, i64) {
    if b == 0 {
        return (a, 1, 0)
    }
    let (gcd, x, y) = eea_recursive(b, a % b);
    (gcd, y, x-(a/b)*y)
}

Further optimizations can be done by making it iterative and not recursive available here in PBF.

Implementing the

$G_1$ curve

We implement the first curve in PBH via what we will now call as F101G1PointAffine representations. As the name suggests, we will use u64 field over prime 101 as coordinate domain, and present it in "affine" form i.e. have an explicit point at inifinity.

Create f101_g1.rs under src/ with following:

use std::{
    fmt::Display,
    ops::{Add, Neg, Mul, Sub}, error::Error,
};
use thiserror::Error;

use crate::{
    elliptic_curve::{FiniteField, G1Point},
    u64_field::U64Field,
};

// Type alias
type F101 = U64Field<101>;

/// A point in the $y^2+x^3+3$ curve, on the $\mathbb{F}_{101}$ field.
/// The generator $g=(1,2)$ generates a subgroup of order 17: $17g=g$
#[derive(Debug, Copy, Clone, PartialEq, Eq, Hash, PartialOrd, Ord)]
pub struct F101G1PointAffine {
    pub x: F101,
    pub y: F101,
    pub is_infinity: bool,
}

#[derive(Error, Debug)]
pub enum G1Error {
    #[error("point not on curve")]
    NotOnCurve
}

impl G1Point for F101G1PointAffine {
    type Field = F101;
    type SubField = F101;

    fn new(x: Self::Field, y: Self::Field) -> Result<Self, Box<dyn Error>> {
        let candidate_point = Self {
            x,
            y,
            is_infinity: false,
        };
        if !candidate_point.in_curve() {
            return Err(Box::new(G1Error::NotOnCurve))
        }
        Ok(candidate_point)
    }

    /// Checks if the point is on the curve. i.e. the equation holds
    /// $y^2 = x^3 + 3$. Requires pow() to provide fresh elements when
    /// called.
    fn in_curve(&self) -> bool {
        self.y.pow(2) == self.x.pow(3) + F101::new(3u64)
    }

    /// Affine representation has point on infinity as identity element
    fn is_identity(&self) -> bool {
        self.is_infinity
    }

    /// Returns the point at infinity
    fn identity() -> Self {
        Self {
            x: Self::Field::zero(),
            y: Self::Field::zero(),
            is_infinity: true,
        }
    }

    /// Described generator in PLONK by hand
    fn generator() -> Self {
        Self {
            x: F101::new(1u64),
            y: F101::new(2u64),
            is_infinity: false,
        }
    }

    /// Returns the size of the subgroup generated by generator $g=(1,2)$
    fn generator_subgroup_size() -> Self::Field {
        F101::new(17u64)
    }

    fn x(&self) -> &Self::Field {
        &self.x
    }
    fn y(&self) -> &Self::Field {
        &self.y
    }
}

impl Display for F101G1PointAffine {
    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
        if self.is_infinity {
            write!(f, "infinite")
        } else {
            write!(f, "({},{})", self.x, self.y)
        }
    }
}

impl Neg for F101G1PointAffine {
    type Output = F101G1PointAffine;

    /// On an elliptic curve, the point's negative is
    /// mirror point along x-axis. Since elliptic
    /// curves are symmetric along x-axis, (x, -y) is
    /// negative point for all non infinity points (x, y)
    fn neg(self) -> Self::Output {
        if self.is_infinity {
            self
        } else {
            F101G1PointAffine::new(self.x, -self.y).unwrap()
        }
    }
}

impl Add for F101G1PointAffine {
    type Output = F101G1PointAffine;
    fn add(self, rhs: F101G1PointAffine) -> Self {
        // Anything added to infinity is itself
        if self.is_infinity {
            return rhs;
        } else if rhs.is_infinity {
            return self;
        }

        // Anything added to its negative is
        // additive identity element
        if self == -rhs {
            return F101G1PointAffine::identity()
        }

        // Actual point addition, refer:
        // https://en.wikipedia.org/wiki/Elliptic_curve_point_multiplication#Point_addition
        // lambda is never used with its init value (0) ever
        let mut lambda = F101::new(0);

        if self != rhs {
            // Normal case
            // Unwrap safety: rhs.x and self.x can not be same
            // since then rhs.y == self.y also, leading to limiting case below
            // or rhs.y == -self.y leading to case handled above already
            lambda = ((rhs.y - self.y) / (rhs.x - self.x)).unwrap()
        } else {
            // Limiting case
            lambda = (F101::new(3) * self.x.pow(2) / (F101::new(2) * self.y)).unwrap()
        }

        let x = lambda.pow(2) - self.x - rhs.x;
        let y = lambda * (self.x - x) - self.y;

        Self {
            x,
            y,
            is_infinity: false,
        }
    }
}

impl Sub for F101G1PointAffine {
    type Output = F101G1PointAffine;
    fn sub(self, rhs: Self) -> Self::Output {
        self + rhs.neg()
    }
}

// Scalar multiplication
impl Mul<F101> for F101G1PointAffine {
    type Output = F101G1PointAffine;
    fn mul(self, rhs: F101) -> Self::Output {
        let mut scalar = rhs.as_u64();

        if scalar == 0 || self.is_identity() {
            return F101G1PointAffine::identity();
        }

        let mut result = F101G1PointAffine::identity();
        let mut base = self;

        while scalar > 0 {
            if scalar % 2 == 1 {
                result = result + base
            }
            scalar >>= 1;
            base = base + base;
        }

        result
    }
}

Appendix

Properties of a bilinear maps (ChatGPT)

Bilinear maps possess several important properties. Let's consider a bilinear map

$B$:

$V$ ×

$W$ →

$X$, where V, W, and X are vector spaces over the same field. The properties of a bilinear map are as follows:

Linearity in the first argument: For any fixed w in W, the function v ↦ B(v, w) is a linear transformation from V to X. This means that it satisfies the following properties:

B(v1 + v2, w) = B(v1, w) + B(v2, w) (additivity) B(cv, w) = cB(v, w) (homogeneity), where v, v1, v2 are vectors in V, c is a scalar, and w is a fixed vector in W. Linearity in the second argument: For any fixed v in V, the function w ↦ B(v, w) is a linear transformation from W to X. This implies:

B(v, w1 + w2) = B(v, w1) + B(v, w2) (additivity) B(v, cw) = cB(v, w) (homogeneity), where w, w1, w2 are vectors in W, c is a scalar, and v is a fixed vector in V. Compatibility with scalar multiplication: The bilinear map satisfies the property:

B(cv, w) = B(v, cw) = cB(v, w), where v is a vector in V, w is a vector in W, and c is a scalar. Distributive property: Bilinear maps distribute over vector addition in both arguments:

B(v1 + v2, w1 + w2) = B(v1, w1) + B(v1, w2) + B(v2, w1) + B(v2, w2), where v1, v2 are vectors in V, w1, w2 are vectors in W. These properties ensure that the bilinear map behaves linearly with respect to each argument and respects the operations of vector addition and scalar multiplication.

Additionally, there are other properties specific to certain types of bilinear maps, such as symmetric bilinear maps (satisfying B(v, w) = B(w, v)) or alternating bilinear maps (satisfying B(v, v) = 0 for all v). These properties depend on the specific context and requirements of the application.

Bilinear maps are not necessarily one-to-one (injective) or onto (surjective). The properties of injectivity and surjectivity depend on the specific bilinear map and the vector spaces involved.

Injectivity: A bilinear map B: V × W → X is injective if distinct pairs of vectors in V and W always map to distinct elements in X. In other words, if B(v1, w1) = B(v2, w2), then it implies that (v1, w1) = (v2, w2). However, in general, bilinear maps can have non-trivial kernels, meaning that different pairs of vectors may map to the same element in X.

Surjectivity: A bilinear map B: V × W → X is surjective if every element in the target space X can be obtained as the result of applying the map to some pair of vectors from V and W. Surjectivity is also dependent on the specific bilinear map and the vector spaces involved. Not all bilinear maps are surjective onto their entire target space.

It's worth noting that when bilinear maps are used in the context of pairing-based cryptography, the bilinear maps involved in cryptographic pairings are carefully constructed to possess desired properties and security characteristics. These constructions often involve selecting appropriate elliptic curves and groups to ensure the desired properties of bilinearity, non-degeneracy, and computability, while also taking into account security considerations.