# How to do Circle STARK math in Bitcoin?

In this document, we explain how the field operations related to Circle STARK are done in the [Bitcoin Wildlife Sanctuary's implementation](https://github.com/Bitcoin-Wildlife-Sanctuary/rust-bitcoin-m31) in the Bitcoin script. This implementation was a fork from the [BitVM implementation](https://github.com/BitVM/rust-bitcoin-m31-or-babybear) of M31 and BabyBear that focuses more on M31-specific optimization (such as two-layer Karatsuba for QM31).

There are three fields being implemented:
- M31: The prime field under the prime modulus $2^{31} - 1$. This is where the FFT in Circle STARK is defined over. The trace in STARK is also encoded as M31 elements.
- CM31: A degree-2 extension field of M31, under the irreducible polynomial $x^2 + 1$. This is in essence the complex field $(a + ib)$ where $i$ is the imaginary unit.
- QM31: A degree-2 extension field of CM31, under another irreducible polynomial $y^2 - 2 - i$. One can write an QM31 element as $(a + ib) + j (c + id)$, where $a, b, c, d$ are M31 elements. Probabilistic tests in Circle STARK are done in this field because it is big enough to offer a sufficient security level for these tests.

As follows, we explain the algorithm that we use to perform field operations (focusing on addition and multiplication). For readibility, we would not present the Bitcoin script. Readers can find the corresponding scripts in the [GitHub repository](https://github.com/Bitcoin-Wildlife-Sanctuary/rust-bitcoin-m31). Instead, we present pseudocode. All the pseudocode is running over the Bitcoin integers.

## M31

In the following, we use $q = 2^{31} - 1$ to represent the prime.

To add two numbers $a$ and $b$, the script works as follows:
- Compute $b' \leftarrow b - q$. This makes $b'$ a negative number.
- Compute $c' \leftarrow a + b'$.
- If $c' < 0$, $c \leftarrow c' + q$. Otherwise, $c \leftarrow c'$.
- Return $c$ as the sum.

The reason for subtracting $q$ from $b$ before adding it with $a$ is to make sure $c'$ stays within the Bitcoin integer limit.

To multiply two numbers $a$ and $b$, we use an algorithm that incorporates a small lookup table. Here are the concrete algorithm:
- Bit-decompose $b$ into bits. Let $d_1, d_2, ..., d_{31}$ represent the 31 bits of $b$, with $d_1$ be the most significant bit.
- Compute a three-element lookup table:
* $(a \mod{q}) - q$
* $(2a \mod{q}) - q$
* $(3a \mod{q}) - q$
- Initialize an aggregator variable $s \leftarrow 0$.
- If $d_{1} = 1$:
* $s' \leftarrow s + ((a \mod{q}) - q)$.
* If $s' < 0$, $s \leftarrow s' + q$. Otherwise, $s \leftarrow s'$.
- For each pair $(d_2, d_3), (d_4, d_5), ..., (d_{30}, d_{31})$:
* Double $s$ twice using the addition algorithm.
* Let the first bit in a pair be $d_l$ and the second bit be $d_r$.
* If $d_l = 1$ and $d_r = 1$:
- $s' \leftarrow s + ((3a \mod{q}) - q)$.
- If $s' < 0$, $s \leftarrow s' + q$. Otherwise, $s \leftarrow s'$.
* Else, if $d_l = 1$ and $d_r = 0$:
- $s' \leftarrow s + ((2a \mod{q}) - q)$.
- If $s' < 0$, $s \leftarrow s' + q$. Otherwise, $s \leftarrow s'$.
* Else, if $d_l = 0$ and $d_r = 1$:
- $s' \leftarrow s + ((a \mod{q}) - q)$.
- If $s' < 0$, $s \leftarrow s' + q$. Otherwise, $s \leftarrow s'$.
- Return $s$ as the product.

The lookup table enables the algorithm to perform fewer than 31 additions. Note that the lookup table includes elements that have already been subtracted by $q$. This is to save the step $b' \leftarrow b - q$ in the addition algorithm.

## CM31

We now move our attention into CM31. Addition is trivial as it is just adding the real part and the imaginary part separately. So, we will now focus on multiplication.

We can represent the two CM31 elements as $a_1 + ib_1$ and $a_2 + ib_2$. The product is:

$$(a_1 a_2 - b_1 b_2) + i (a_1 b_2 + a_2 b_1)$$

The naive way is to compute the four products $a_1 a_2$, $b_1 b_2$, $a_1 b_2$, $a_2 b_1$ individually, but one can do so more efficiently through Karatsuba, the algorithm of which is as follows:
- Compute $p_1 = a_1 a_2$ and $p_2 = b_1 b_2$.
- Compute $t_1 = a_1 + b_1$ and $t_2 = a_2 + b_2$.
- Compute $p_3 = t_1 t_2 = a_1 a_2 + a_1 b_2 + a_2 b_1 + b_1 b_2$.
- Compute $a_3 = p_1 - p_2$.
- Compute $b_3 = p_3 - p_1 - p_2$.

One can see that this algorithm only uses three M31 multiplications.

## QM31

Similarly, addition in QM31 is trivial as it is adding the first element and the second element, both of which are CM31, separately. Here, we describe the multiplication algorithm, which is slightly more complicated as we need to multiply the square of the nonresidue.

Let us represent a QM31 element as $c_1 + jd_1$ and $c_2 + jd_2$ where $c_1$, $d_1$, $c_2$, $d_2$ are all CM31 elements and $j$ is the unit for the second element. Note that the product would look as follows:

$$(c_1 c_2 + j^2 d_1 d_2) + j(c_1 d_2 + c_2 d_1) $$

Here $j^2 = 2 + i$.

The algorithm is as follows.
- Compute $p_1 = c_1 c_2$ and $p_2 = d_1 d_2$.
- Compute $t_1 = c_1 + d_1$ and $t_2 = c_2 + d_2$.
- Compute $p_3 = t_1 t_2 = c_1 c_2 + c_1 d_2 + c_2 d_1 + d_1 d_2$.
- Compute $c_3 = c_1 c_2 + (2 + i) d_1 d_2$.
- Compute $d_3 = p_3 - p_1 - p_2$.

One can see that it uses three CM31 multiplications, which also means that QM31 multiplication in total takes 9 M31 multiplications.

It is useful to note that since $c_1 c_2$ and $d_1 d_2$ are in CM31, so they already have an imaginary part. Let $c_1 c_2 = \mathbf{Re}(c_1 c_2) + i\cdot \mathbf{Im}(c_1 c_2)$ and $d_1 d_2 = \mathbf{Re}(d_1 d_2) + i\cdot \mathbf{Im}(d_1 d_2)$.