# Pippenger's Algorithm

> AKA "Bucket" method

## Motivation

Improve upon [Double-and-add](https://fractalyze.gitbook.io/intro/~/revisions/GBuVDfZJZBNDtoleSgnx/primitives/abstract-algebra/elliptic-curve/scalar-multiplication/double-and-add) by reducing the total number of double and add operations needed for [MSM](https://fractalyze.gitbook.io/intro/~/revisions/GBuVDfZJZBNDtoleSgnx/primitives/abstract-algebra/elliptic-curve/msm). 

## Algorithm Process

### Introduction

**Multi-Scalar Multiplication (MSM)** refers to the efficient computation of expressions in the following form:

$$
S = \sum\_{i=1}^{n} k\_i P\_i
$$

* $$k\_i$$: the $$i$$-th $$\lambda$$-bit scalar
* $$P\_i$$: the $$i$$-th point on an elliptic curve

**Pippenger's algorithm** speeds up MSM by dividing each scalar $$k\_i$$ into multiple $$s$$**-bit chunks**, and then grouping and processing scalars by their corresponding bit positions (windows).

### 1. Scalar Decomposition

First, we decompose the scalars $$k\_i$$:

$$
\begin{align\*}
k\_i &= k\_{i,1} + 2^{s}k\_{i,2} + ... + 2^{(\lceil\frac{\lambda}{s}\rceil -1)s}k\_{i,\lceil\frac{\lambda}{s}\rceil } \\
k\_i &= \sum\_{j=1}^{\lceil\frac{\lambda}{s}\rceil }2^{(j - 1)s} k\_{i,j}
\end{align\*}
$$

* $$s$$ : the number of bits per window
* $$\lceil\frac{\lambda}{s}\rceil$$ : the total number of windows

Through this definition, we can rewrite the original MSM summation as follows:

$$
\begin{align\*}
S &= \sum\_{i=1}^{n} k\_i P\_i \\
&= \sum\_{i=1}^{n} \sum\_{j=1}^{\lceil\frac{\lambda}{s}\rceil } 2^{(j-1)s}k\_{i,j} P\_i \\
&= \sum\_{j=1}^{\lceil\frac{\lambda}{s}\rceil } 2^{(j-1)s} \left( \sum\_{i=1}^{n} k\_{i,j} P\_i \right) \\
&= \sum\_{j=1}^{\lceil\frac{\lambda}{s}\rceil } 2^{(j-1)s} W\_j
\end{align\*}
$$

Note that $$W\_j$$ is the sum for the $$j$$-th window.

### 2. Bucket Accumulation 

Given our $$W\_j$$, we denote buckets $$B$$ as the following:

$$
B\_{\ell, j} = \sum\_{i: k\_{ij} = \ell} P\_i
$$

Or in other words, $$B\_{\ell, j}$$ refers to the $$\ell$$ bucket in window $$W\_j$$ consisting of the summation of points $$P\_i$$ with the same $$k\_{i,j}$$.

Note that $$k\_{i,j}$$ must be within the range $$\[0, 2^w - 1]$$; however, if $$k\_{i,j}$$ is 0, as $$k\_{i,j}$$ is the scalar multiplier, the result of the bucket will be 0. Therefore, the practical range of $$k\_{i,j}$$ is $$\[1, 2^w - 1]$$.

### 3. Bucket Reduction

Thus, in total, a window $$W\_j$$ will be re-structured as follows:

$$
W\_j = \sum\_{i=1}^{n} k\_{ij}P\_i = \sum\_{\ell=1}^{2^s - 1} \ell B\_{\ell, j}
$$

Note that the computation of a singular $$\ell B\_{\ell,j}$$ is computed through a running sum of recursive summations and not through further scalar multiplications.

#### Complexity Cost:

* Original: $$(2^s - 1) \times \mathsf{ScalarMul} + (2^s - 2) \times \mathsf{PointAdd}$$
* Bucketed: $$(2^{s+1} - 2) \times \mathsf{PointAdd}$$

#### C++ Implementation

```cpp
Bucket ReduceBucketPoints(const std::vector<Bucket>& buckets) {
  Bucket res = Bucket::Zero();
  Bucket sum = Bucket::Zero();
  for (size_t i = buckets.size() - 1; i != SIZE_MAX; --i) {
    sum += buckets[i];
    res += sum;
  }
  return res;
}
```

#### Example

$$
W\_1 = B\_{1,1} + 2B\_{2,1} + 3B\_{3,1} + 4B\_{4,1}
$$

```
// Round 1
sum = B_(4,1);
res = B_(4,1);
// Round 2
sum = B_(3,1) + B_(4,1);
res = B_(4,1) + [B_(3,1) + B_(4,1)];
// Round 3
sum = B_(2,1) + B_(3,1) + B_(4,1);
res = B_(4,1) + [B_(3,1) + B_(4,1)] + [B_(2,1) + B_(3,1) + B_(4,1)];
// Round 3
sum = B_(1,1) + B_(2,1) + B_(3,1) + B_(4,1);
res = B_(4,1) + [B_(3,1) + B_(4,1)] + [B_(2,1) + B_(3,1) + B_(4,1)] + [B_(1,1) + B_(2,1) + B_(3,1) + B_(4,1)];

// Total:
res = B_(1,1) + 2*B_(2,1) + 3*B_(3,1) + 4*B_(4,1);
```

### 4.  Window Reduction - Final MSM Result

With the finished individual window computations, we return to the total MSM summation and add all windows together:

$$
S= \sum\_{j=1}^{\lceil\frac{\lambda}{s}\rceil } 2^{(j-1)s} W\_j
$$

In practice, the computation is calculated recursively from $$j = \left\lceil \frac{\lambda}{s} \right\rceil$$ to $$j = 1$$ like so:

$$
T\_j = \begin{cases}
0 & j = \left\lceil \frac{\lambda}{s} \right\rceil \\
W\_j + 2 ^sT\_{j+1} & j \ge 1 \\
\end{cases}
$$

#### Complexity Cost

$$
\left\lceil \frac{\lambda}{s} \right\rceil(s \times \mathsf{PointDouble} + \mathsf{PointAdd})
$$

#### C++ Implementation

```cpp
Bucket ReduceWindows(uint32_t lambda, uint32_t s, const std::vector<Bucket>& buckets) {
  Bucket res = Bucket::Zero();
  uint32_t window_num = static_cast<uint32_t>(std::ceil(lambda / s));
  for (uint32_t j = window_num - 1; j != UINT32_MAX; --i) {
    for (uint32_t i = 0; i < s; ++i) {
      res = res * 2;
    }
    res += buckets[j];
  }
  return res;
}
```

## Complexity Analysis

* To compute $$W\_j$$, $$n \times \mathsf{PointAdd}$$ is required: $$\left\lceil \frac{\lambda}{s} \right\rceil(n \times \mathsf{PointAdd})$$
* Each $$W\_j$$ needs $$(2^{s+1} - 2) \times \mathsf{PointAdd}$$: $$\left\lceil \frac{\lambda}{s} \right\rceil((2^{s+1} - 2) \times \mathsf{PointAdd})$$
* Final $$Q$$ needs: $$\left\lceil \frac{\lambda}{s} \right\rceil(s \times \mathsf{PointDouble} + \mathsf{PointAdd})$$

### Total complexity

$$
\left\lceil \frac{\lambda}{s} \right\rceil(s \times \mathsf{PointDouble} + (n + 2^{s+1} - 1) \times \mathsf{PointAdd}) \approx \\\lambda \times \mathsf{PointDouble} + \left\lceil \frac{\lambda}{s} \right\rceil(n + 2^{s+1} - 1)\times \mathsf{PointAdd}
$$

### Optimization Table (When , $$n = 2^{20}$$)

<table><thead><tr><th width="298.8214111328125">s</th><th>PointAdd Count</th></tr></thead><tbody><tr><td>2</td><td>134,218,624</td></tr><tr><td>4</td><td>67,110,848</td></tr><tr><td>8</td><td>33,570,784</td></tr><tr><td><mark style="color:green;">16</mark></td><td><mark style="color:green;">18,874,352</mark></td></tr><tr><td>32</td><td>68,727,865,336</td></tr><tr><td>64</td><td>147,573,952,589,680,607,228</td></tr><tr><td>128</td><td>1,361,129,467,683,753,853,853,498,429,727,074,942,974</td></tr></tbody></table>

> 🔍 This clearly shows that selecting an optimal $$s$$ value is critical for performance.

## Example

Let's try to run Pippenger's on the following simple example:

$$
\begin{align\*}
S &= k\_1P\_1 + k\_2P\_2 + k\_3P\_3\\
&=12P\_1 + 9P\_2 + 13P\_3
\end{align\*}
$$

### 1. Scalar Decomposition

Let's define our bits per window $$s$$ to be 2. In this case, since the maximum amount of bits needed to constrain the values 12, 9, and 13 is 4, we will have $$\lceil\frac{\lambda}{s}\rceil  =\mathsf{max\_bits} / s = 4/2 = 2$$ windows.

Note that in binary we have:

$$
k\_1 =12 = \underbrace{0}*{2^0}\underbrace{0}*{2^1}\underbrace{1}*{2^2}\underbrace{1}*{2^3}\\
k\_2 = 9 = \underbrace{1}*{2^0}\underbrace{0}*{2^1}\underbrace{0}*{2^2}\underbrace{1}*{2^3} \\
k\_3 =13 = \underbrace{1}*{2^0}\underbrace{0}*{2^1}\underbrace{1}*{2^2}\underbrace{1}*{2^3}
$$

Decomposing our $$k\_0$$, $$k\_1$$, and $$k\_2$$ results in the following:

$$
\begin{align\*}
k\_i &= k\_{i,1} + 2^{s}k\_{i,2} + ... + 2^{(\lceil\frac{\lambda}{s}\rceil -1)s}k\_{i,\lceil\frac{\lambda}{s}\rceil } \\

\end{align\*}
$$

$$
\begin{align\*}
\begin{array}{rlrlrl}
k\_1 = 12 &= \overbrace{\underbrace{0}*{2^0}\underbrace{0}*{2^1}}^{\text{window 1}}\overbrace{\underbrace{1}*{2^2}\underbrace{1}*{2^3}}^{\text{window 2}}
& \quad
k\_2 = 9 &= \overbrace{\underbrace{1}*{2^0}\underbrace{0}*{2^1}}^{\text{window 1}}\overbrace{\underbrace{0}*{2^2}\underbrace{1}*{2^3}}^{\text{window 2}}
& \quad
k\_3 = 13 &= \overbrace{\underbrace{1}*{2^0}\underbrace{0}*{2^1}}^{\text{window 1}}\overbrace{\underbrace{1}*{2^2}\underbrace{1}*{2^3}}^{\text{window 2}} \\\[1ex]

&= k\_{1,1} + 2^sk\_{1,2}
& \quad &= k\_{2,1} + 2^sk\_{2,2}
& \quad &= k\_{3,1} + 2^sk\_{3,2} \\\[1ex]

&= 0 + 3(2^2)
& \quad &= 1 + 2(2^2)
& \quad &= 1 + 3(2^2) \\\[1ex]
\end{array}
\end{align\*}
$$

### 2+3. Bucket Accumulation + Reduction

$$
W\_j = \sum\_{i=1}^{n} k\_{ij}P\_i = \sum\_{\ell=1}^{2^s - 1} \ell B\_{\ell, j}
$$

Given our $$s$$ of 2, we now have buckets of $$\[B\_{1,j}, B\_{2,j}, B\_{3,j}]$$ for a given window $$W\_j$$. Our windows are now re-structured to buckets like so:

$$
\begin{align\*}
W\_1 &= \overbrace{1}^{k\_{2,1}}P\_2 + \overbrace{1}^{k\_{3,1}}P\_3 = 1(P\_2+P\_3) = B\_{1,1}\\
W\_2 &= \overbrace{3}^{k\_{1,2}}P\_1 + \overbrace{2}^{k\_{2,2}}P\_2 + \overbrace{3}^{k\_{3,2}}P\_3 = 2(P\_2) + 3(P\_1 + P\_3) = 2B\_{2,2} + 3B\_{3,2}
\end{align\*}
$$

From bottom up, in our **bucket accumulation** stage we compute the buckets:

| Bucket       | Computation     |
| ------------ | --------------- |
| $$B\_{1,1}$$ | $$P\_2 + P\_3$$ |
| $$B\_{2,2}$$ | $$P\_2$$        |
| $$B\_{3,2}$$ | $$P\_1+P\_3$$   |

...and in our **bucket reduction** stage we compute the windows:

$$
\begin{align\*}
W\_1 &= B\_{1,1}\\
W\_2 &= (B\_{2,2} + B\_{2,2}) + (B\_{3,2} + B\_{3,2} + B\_{3,2})
\end{align\*}
$$

### 4.  Window Reduction - Final MSM Result

Putting all the windows together for the MSM sum, we get the following:

$$
\begin{align\*}
S &=\sum\_{j=1}^{\lceil\frac{\lambda}{s}\rceil } 2^{(j-1)s} W\_j \\
&= W\_1 + 2^2W\_2 \\
\end{align\*}
$$

## References

* <https://hackmd.io/lNOsNGikQgO0hLjYvH6HJA>
* <https://eprint.iacr.org/2022/1321.pdf>

> Written by [Ryan Kim](https://app.gitbook.com/u/cPk8gft4tSd0Obi6ARBfoQ16SqG2 "mention") & [Ashley Jeong](https://app.gitbook.com/u/PNJ4Qqiz7kSxSs58UIaauk95AO43 "mention") of Fractalyze