# EPF Update 10&11 - In Math we Trust

After understanding the entire flow of the Verkle Cryptography API, and how it's supposed to work in Nim, for the Nimbus client.

I raised an [issue](https://github.com/mratsim/constantine/issues/275) soon after in [Constantine](https://github.com/mratsim/constantine/), and broke down my work into 15 objective tasks.

They were soon after listed in my [PR](https://github.com/mratsim/constantine/pull/276). The tasks were as follows:

Fixes #[275](https://github.com/mratsim/constantine/issues/275)

- [X] Transcript generation using Fiat Shamir style
- [X] Barycentric form using precompute optimisation over domain
- [X] Common Reference String generation
- [X] Common util funs for ipa
- [X] IPA functions for prover
- [X] IPA functions for verifier
- [X] Functions for creating multiproofs
- [X] Functions for checking multiproofs
- [X] Random element generator for Pedersen
- [X] Assert for no duplicate elements for random elem gen
- [ ] Individual tests
- [ ] Integrated tests
- [X] Functions for Pedersen Commitments
- [X] Adding more comments for explanation/ more resource links
- [X] Refactoring to generics

This to-do is updated to 7/10/2023.

## Barycentric form

For the past few weeks I have been racking my head on how to implement the precompute approach for Lagrangian interpolation of the polynomials. However, soon after I figured out, thanks to this article.

The followting is an excerpt from Kev's article.

### The Lagrange Polynomial

We briefly restate the formula for a lagrange polynomial:

$$
\mathcal{L_i}(X) = \prod_{j \neq i, j = 0}\frac{X -x_j}{x_i - x_j}
$$

> The i'th lagrange polynomial evaluated at $x_i$ is 1 and 0 everywhere else **on the domain**

## First form of the barycentric interpolation formula

We introduce the polynomial $A(X) = (X - x_0)(X - x_1)...(X-x_n)$.

We also introduce the derivative of $A'(X) = \sum_{j=0}^{d-1}\prod_{i \neq j}(X - x_i)$ .

> You can derive this yourself by generalising the product rule: <https://en.wikipedia.org/wiki/Product_rule#Product_of_more_than_two_factors>

In general this derivative does not have a succint/sparse form. We do have a succint form if the domain is the roots of unity!

Now note that $A'(x_j) = \prod_{i=0,i \neq j}(x_j - x_i)$

> If we plug in $x_k$ into $A'(X)$ all the terms with $X - x_k$ will vanish, this is why the sum disappears into a single product.

We can use $A$ and $A'$ to re-define our lagrange polynomial as :

$$
\mathcal{L_i}(X) = \frac{A(X)}{A'(x_i) (X - x_i)}
$$

>Looking at the original lagrange formula, $A'(x_i)$ is the denominator and $\frac{A(X)}{X - x_i}$ is the numerator.

The first barycentric form for a polynomial $f(X)$ can now be defined as :

$$
f(X) = \sum_{i=0}^{d-1}{\frac{A(X)}{A'(x_i) (X - x_i)} f_i}
$$

#### Remarks

- $A(X)$ is not dependent on the values of $f_i$ and so can be brought out of the summation.
- $A'(X)$ is only dependent on the domain, so it can be precomputed, along with $A(X)$

The Nim equivalent is stated [here](https://github.com/mratsim/constantine/blob/3292cd9a6c6b3e47da83d392fb29dc5c74e439ac/constantine/commitments/ipa/barycentric_form.nim).

One of the main portions of the code is as follows:

func newPrecomputedWeights* [PrecomputedWeightsObj] (res: var PrecomputedWeights)=
 var midpoint: uint64 = DOMAIN
 var barycentricWeightsInst {.noInit.} : array[(midpoint*2), EC_P_Fr]
 
 for i in uint64(0)..midpoint:
  var weights {.noInit.}: EC_P_Fr
  weights.barycentricWeights(i)

  var inverseWeights {.noInit.}: FF
  inverseWeights.inv(weights)

  barycentricWeightsInst[i] = weights
  barycentricWeightsInst[i+midpoint] = inverseWeights

  midpoint = DOMAIN - 1
  var invertedDomain: array[(midpoint*2), EC_P_Fr] 

  for i in uint64(0)..DOMAIN:
   var k: EC_P_Fr = cast[EC_P_Fr](i)

   k.inv(k)

   var neg_k : EC_P_Fr

   var zero : FF

   zero.setZero()

   neg_k.diff(zero, k)

   invertedDomain[i-1] = k
   invertedDomain[(i-1) + midpoint] = neg_k
   
   res.barycentricWeights = barycentricWeightsInst
   res.invertedDomain = invertedDomain

Here the precomputed weights are $A(X)$ and $A'(X)$ as follows. The type defined to store them is as follows :

type 
 PrecomputedWeights* = object
  #This stores A'(x_i) and 1/A'(x_i)
  barycentricWeights: openArray[EC_P_Fr]
  #This stores 1/k and -1/k for k in [0, 255]
  invertedDomain: openArray[EC_P_Fr]

The theoretical equivalent of this is here-
For our case where the domain is the discrete interval $[0, 255]$

We only need to precompute $\frac{1}{x_i}$ for $x_i \in [-255, 255]$. This is 510 values, so we would store 510 * 32 = 16Kb.

**How would we lookup these values?**

First we imagine that we have stored the values in an array:

$[\frac{1}{1}, \frac{1}{2}, \frac{1}{3}, \frac{1}{4}... \frac{1}{255},\frac{1}{-1},\frac{1}{-2},...\frac{1}{-255}]$

We first note that we can easily get from $\frac{1}{k}$ to $\frac{1}{-k}$ in the array by jumping forward 255 indices. Our strategy will be to find $\frac{1}{k}$ then jump to $\frac{1}{k}$ if we need to.

**Example**

We want to compute $\frac{1}{0 - 255}$.

- Compute the $abs(0-255) = 255 = i$

In practice, we can use an if statement to check whether 255 or 0 is larger, and subtract accordingly.

- Note that $\frac{1}{i}$ is at index $i-1$
- Since our original computation was $0 - 255$ which is negative, we need to get the element at index $i - 1 + 255$

### Precomputing $\frac{A'(x_m)}{A'(x_i)}$

With the roots of unity, we did not need this optimisation as $\frac{A'(x_m)}{A'(x_i)}$ equaled $\frac{\omega^i}{\omega^m}$ which was trivial to fetch from the domain.

For this, we will need to store precomputed values, if we want to efficiently compute $q_m$ in $O(d)$ steps, and to also avoid inversions.

(Dankrad): We precompute $A'(x_i)$ and $\frac{1}{A'(x_i)}$. Given that we have 256 points in the domain. This will cost us 256 *2* 32 bytes = 16kB.

**How would I lookup these values?**

Similar to the previous optimisation, we store $A'(x_i)$ in an array.

$[A'(0), A'(1), A'(2), A'(3)... A'(255),\frac{1}{A'(0)},\frac{1}{A'(1)},\frac{1}{A'(2)},...\frac{1}{A'(255)}]$

**Example**

We want to compute $\frac{A'(0)}{A'(5)}$

- We can fetch $A'(0)$ by looking up the element at index $0$ in the array.
- We can fetch $\frac{1}{A'(5)}$ by looking up the element at index 5, then jumping forward 256 positions.

This is an optimisation that is being followed only in **Go** implementation apart from the Nim (Constantine) implementation, which is newly being developed.

## Other important Math functions for Transcript generation and Inner Product Arguments.

Most of these functions exist in the `common_utils.nim` file in my Constantine PR.

One of the very important tasks for me in here was writing a function to **Derive the generator element for Pedersen Commitment**.

These are the set of dependencies

## Seed

The current seed being used is: `eth_verkle_oct_2021`

## Rust implementation

 fn generate_random_elements(num_required_points: usize) -> Vec<EdwardsAffine> {
    use bandersnatch::Fq;
    use sha2::{Digest, Sha256};

    let mut hasher = Sha256::new();
    hasher.update(b"eth_verkle_oct_2021");
    let bytes = hasher.finalize().to_vec();

    let u = bandersnatch::Fq::from_be_bytes_mod_order(&bytes);
    let choose_largest = false;

    (0..)
        .into_iter()
        .map(|i| Fq::from(i as u128) + u)
        .map(|x| EdwardsAffine::get_point_from_x(x, choose_largest))
        .filter_map(|point| point)
        .filter(|point| point.is_in_correct_subgroup_assuming_on_curve())
        .take(num_required_points)
        .collect()
}

For Nim we use the `Projective` coordinates, as that's what we've been using in order to make the library comparitively as fast as Go. Here is how the code looks like -

func generate_random_elements* [FF]
    (points: var  openArray[FF], 
     num_points: uint64)  =

    var incrementer: uint64 = 0

    while uint64(len(points)) !=  num_points:

        var digest : sha256
        digest.init()
        digest.update(seed)

        var b {.noInit.} : array[8, byte]
        digest.update(b)

        var hash {.noInit.} : array[32, byte]

        digest.finish(hash)

        var x {.noInit.}:  FF

        x.deserialize(hash)
        doAssert(cttCodecEcc_Success)
        ## Deserialisation success !
        incrementer=incrementer+1

        var x_as_Bytes {.noInit.} : array[32, byte]
        x_as_Bytes.serialize(x)
        doAssert(cttCodecEcc_Success)
        ## Serialisation Success !

        var point_found {.noInit.} : EC_P
        point_found.deserialize(x_as_Bytes)

        doAssert(cttCodecEcc_Success)
        points[incrementer] = point_found

## For the next article

Sorry for writing articles biweekly, I just get lost into implementation sometimes, next update will be about how I'll be `implementing splitting and folding of scalars` while performing Pedersen Commitments inside the Inner Product Argument setting.

I'll be also updating about how I'll be creating IPA proofs, creating a Multiproof using those IPA proofs and finally calling those proofs in a Multiproof Verifying Mechanism.

Till then, best wishes :+1: