If you are a **Snarker** you must have already heard of **KZG** commitment, if you heard of **KZG** you must know **Pairings**. This is what we will talk about.

<br />

# What we do not have

- group theory, field theory and homomorphism

basic concepts will not be covered here, you can refer any books on **abstract algebraic**

- divisors

basic concepts will not be covered here, [Pairing for Beginners](https://static1.squarespace.com/static/5fdbb09f31d71c1227082339/t/5ff394720493bd28278889c6/1609798774687/PairingsForBeginners.pdf) is enough for pairings, if you want to know more about that you might need to dive into any books on **algebraic geometry**

- structure of elliptic curve over finite field and its arithmetics (scalar multiplication)

theories and algorithms will not be covered here, you can refer [Guide to Elliptic Curve Cryptography](http://tomlr.free.fr/Math%E9matiques/Math%20Complete/Cryptography/Guide%20to%20Elliptic%20Curve%20Cryptography%20-%20D.%20Hankerson,%20A.%20Menezes,%20S.%20Vanstone.pdf)

- hash to curve

will be filled in later on

- non-affine coordinate

will be filled in later on

- advanced scalar multiplication algorithms GLV/GLS

will be filled in later on

<br />

---

# What we have

This artical mainly focus on variant pairings:
- [Weil Pairing](#Weil-Pairing)
- [Tate Pairing](#Tate-Pairing)
- [Ate Pairing](#Ate-Pairing)
- [Optimal Pairing](#Optimal-Pairing)

and their practical implementations. Besides, we include some neccesary tricks in implementation, especially:
- [Tower Fields](#Tower-Fields)
- [Frobenius Map](#Magic-Power-of-Frobenius-Map)
- [Twist](#Magic-Power-of-Twist)
- [Cyclotomic Group](#Arithmetics-over-Cyclotomic-Group)

---

# About Code

- [python implementation](https://github.com/PayneJoe/crypto_research/tree/main/ecc/src/pairings)

Mainly focus on pairing computations, including **Miller Loop** and **Final Exponentiation**. It has already been done.

Arithmetics of finite field and elliptic curve not been cracked yet, just use **Galois Field** and **Elliptic Curve** of sagemath instead.

- [rust implementation](https://github.com/PayneJoe/crypto_research/tree/main/ecc/src/pairings)

All written with dirty hands, from **Bigint** Arithmetics to **Finite Field** arithmetics to **Elliptic Curve** arithmetics, to **Pairings**. Much testation work to be done ...

<br />

---

# Public Information

<br />

- Modulus of base prime field (characteristic) $F_p$ with 381-bits:
$$
p = 4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129015664037894272559787
$$

<br />

- Embedding degree, or the degree of full extension field $F_{p^k}$:
$$
k = 12
$$

<br />

- Elliptic Curve (additive group) over **base prime field** $F_p$:
$$
\mathbb{G}_1/E(F_p): y^2 = x^3 + 4
$$

<br />

- Elliptic Curve (additive group) over extension field $F_{p^k}$:
$$
\mathbb{G}_2/E(F_{p^k}): y^2 = x^3 + 4
$$

<br />

- Largest prime factor of $|E(F_p)|$ with 255-bits:
$$
r = 52435875175126190479447740508185965837690552500527637822603658699938581184513
$$

<br />

- Trace of Frobenius:
$$
t = p + 1 - |E(F_p)| = -15132376222941642751
$$

<br />

- Parameter for BLS12 Pairing-family:
$$
x = -15132376222941642752
$$
while:
$$
\begin{aligned}
r(x) &= x^4 - x^2 + 1 \\
p(x) &= (x - 1)^2 \cdot r(x) \cdot \frac{1}{3} + x\\
t(x) &= x + 1 \\
\end{aligned}
$$

<br />

- Target (multiplicative) group with order $r$ defined over $F_{p^k}$:
$$
\mathbb{G}_T: F_{p^k}^{\times}[r]
$$

<br />

----

# Field Arithmetics Underneath

<br />

BLS12-381 curve is defined as below:
$$
E(F_{p^{12}}): y^2 = x^3 + 4
$$
where $F_{p^{12}} = F_p[v] / X^{12} - 2 X^6 + 2$. Where does the field come from?

<br />

## Tower Fields

### Definition

<br />

As we know $G_T$ within BLS12-381 pairing is a $r$-torsion subgroup of multiplicative group defined over $F_{p^{12}}$, usually denoted as $F_{p^{12}}^{\times}[r]$.

<br />

How is $F_{p^{12}}$ constructed? **tower fields** comes:

$$
F_p[u] \xrightarrow{X^2 - \alpha} F_{p^2}[v] \xrightarrow{X^3 - \beta} F_{p^6}[w] \xrightarrow{X^2 - \gamma} F_{p^{12}}
$$

<br />

where these constant parameters of extension field modulus are (specially for BLS12-381):
$$
\begin{aligned}
\alpha &= -1 \in F_{p} \\
\beta &= u + 1 = \sqrt{\alpha} + 1 = \sqrt{-1} + 1 \in F_{p^2} \\
\gamma &= v = \sqrt[3]{\beta} = \sqrt[3]{u + 1} = \sqrt[3]{\sqrt{\alpha} + 1} = \sqrt[3]{\sqrt{-1} + 1} \in F_{p^6}
\end{aligned}
$$

<br />

:::info
$\bigstar \bigstar \bigstar$

- $X^2 - \gamma = X^2 - v$ is irreducible in $F_{p^6}$

- since $v$ is one root of $X^3 - \beta$, then we have $X^6 - \beta$ is irreducible in $F_{p^2}$

- since $\beta - 1$ is one root of $X^2 - \alpha$, then we have $(X^6 - 1)^2 - \alpha = X^{12} - 2 X^6 + 2$ is irreducible in $F_p$

- therefore we have $F_{p^{12}} = F_p[v] / X^{12} - 2 X^6 + 2$
:::

<br />

That's to say, arithmetics over $F_{p^{12}}$ can be accomplished through $F_{p^6}$, while arithmetics over $F_{p^6}$ can be accomplished through $F_{p^2}$, and $F_{p^2}$ can be accomplished through $F_p$. This what we called **tower fields**.

<br />

You may have already noticed that $F_{p} \longrightarrow F_{p^2}$ and $F_{p^6} \longrightarrow F_{p^{12}}$ are both **quadratic extension**, $F_{p^2} \longrightarrow F_{p^6}$ is a **cubic extension**. So **quadratic extension** and **cubic extension** plays an important role within arithmetics over extension fields.

<br />

### Arithmetics for Quadratic Extension

This part you can directly refer 5.2.1 of [Guide to Pairing-based Cryptography](https://www.math.u-bordeaux.fr/~damienrobert/csi/book/pairings.pdf).

<br />

### Arithmetics for Cubic Extension

This part you can directly refer 5.2.2 of [Guide to Pairing-based Cryptography](https://www.math.u-bordeaux.fr/~damienrobert/csi/book/pairings.pdf).

<br />

### Arithmetics over Cyclotomic Group

Cyclotomic group plays an important role in **Final Expoentiation** especially for **Tate/Ate Pairings**. It has special arithmetics, especially **squaring** and **exponentiation**.

<br />

Assume $\alpha \in F_{p^3}, p = q^2$, $F_q$ denotes **base prime field**, how to compute square of $\alpha$?

<br />

First of all, we need a way to represent $\alpha$, for example:
$$
\begin{aligned}
F_{p^3} &= F_{(q^2)^3} = F_{q^2}[v] / v^3 - u \\
F_{q^2} &= F_q[u] / u^2 - \xi \\
\end{aligned}
$$

<br />

Assume :
$$
\alpha = a + b \cdot v + c \cdot v^2
$$

how about $\alpha^2$ ?:
$$
\begin{aligned}
\alpha^2 &= (a + bv + cv^2)^2 \\
&= a^2 + b^2 v^2 + c^2 v^4 + 2(a \cdot bv + a \cdot cv^2 + bv \cdot cv^2) \\
\end{aligned}
$$

since $v^3 - u = 0$, we can proceed forward:
$$
\begin{aligned}
\alpha^2 &= (a^2 + 2 bc \cdot u) + (c^2 \cdot u + 2 ab) \cdot v + (b^2 + 2 ac) \cdot v^2 \\
&= A + B \cdot v + C \cdot v^2
\end{aligned}
$$

<br />

So there are at least 3 squares in $F_{q^2}$ (they are $a^2, b^2, c^2$ respectively), and 5 multiplications in $F_{q^2}$ (they are $a \cdot b, b \cdot c, a \cdot c, bc \cdot u, c^2 \cdot u$ respectively).

<br />

Multiplications over $F_{q^2}$ might be a little more expensive. Is there any improvements on this? YES

<br />

#### Squaring Friendly Field

<br />

If $6 | n$, and $p$ is a large prime characteristic, the order of multiplicative group $F_{p^n}^{\times}$ can be represented multiplication of multiple cyclotomic polymomials:
$$
p^n - 1 = \prod_{i = 1}^{i | n} \varPhi_i(p)
$$

then we call $F_{p^n}$ is a **Squaring Friendly Field**.

<br />

For example:
$$
p^{12} - 1 = \varPhi_1(p) \cdot \varPhi_2(p) \cdot \varPhi_3(p) \cdot \varPhi_4(p) \cdot \varPhi_6(p) \cdot \varPhi_{12}(p)
$$

where:
$$
\begin{aligned}
\varPhi_1(p) &= p - 1 \\
\varPhi_2(p) &= p + 1 \\
\varPhi_3(p) &= (p - \zeta_3^1) \cdot (p - \zeta_3^2) = p^2 + p + 1, \zeta_3^3 = 1 \\
\varPhi_4(p) &= (p - \zeta_4^1) \cdot (p - \zeta_4^3) = p^2 + 1, \zeta_4^4 = 1 \\
\varPhi_6(p) &= (p - \zeta_6^1) \cdot (p - \zeta_6^5) = p^2 - p + 1, \zeta_6^6 = 1 \\
\varPhi_{12}(p) &= (p - \zeta_{12}^1) \cdot (p - \zeta_{12}^5) \cdot (p - \zeta_{12}^7) \cdot (p - \zeta_{12}^{11}) = p^4 - p^2 + 1, \zeta_{12}^{12} = 1 \\
\end{aligned}
$$

<br />

another words, the order of group $F_{p^{12}}^{\times}$ can be factored into:
$$
|F_{p^{12}}^{\times}| = p^{12} - 1 = (p - 1) \cdot (p + 1) \cdot (p^2 + p + 1) \cdot (p^2 + 1) \cdot (p^2 - p + 1) \cdot (p^4 - p^2 + 1)
$$

<br />

So there must be a (r-torsion) subgroup $\mathbb{G}_{\varPhi_{12}(p)}$ of order $r = p^4 - p^2 + 1$ within multiplicative group $F_{p^{12}}^{\times}$.

<br />

In pairing cryptography (if embedding degree is 12), computations are focus on arithmetics over $\mathbb{G}_{\varPhi_{12}(p)}$. There is a great property upon this:
$$
\begin{aligned}
&\alpha \in \mathbb{G}_{\varPhi_{12}(p)} \subset F_{p^{12}}^{\times} \\
&\Downarrow \\
&\alpha^{\varPhi_{12}(p)} = \alpha^{p^4 - p^2 + 1} = 1
\end{aligned}
$$

<br />

#### Fast Squaring

<br />

Back to squaring arithmetic of $\alpha \in F_{(q^2)^3} = F_{q^2}[v] / v^3 - u$:
$$
\begin{aligned}
\alpha = a + b \cdot v + c \cdot v^2, \{a, b, c\} \in F_{q^2} \\
\end{aligned}
$$

according to above formula, the squaring is:
$$
\begin{aligned}
\alpha^2 &= (a^2 + 2 bc \cdot u) + (c^2 \cdot u + 2 ab) \cdot v + (b^2 + 2 ac) \cdot v^2 \\
&= A + B \cdot v + C \cdot v^2
\end{aligned}
$$

The problem is how to compute $\{a \cdot b, b \cdot c, a \cdot c\}$ efficiently?

<br />

Since $\alpha \in \mathbb{G}_{\varPhi_6(q)}$, then we have:
$$
\begin{aligned}
\alpha^{\varPhi_6(q)} = \alpha^{q^2 - q + 1} = 1 \\
\Longrightarrow \alpha^{q^2} \cdot \alpha = \alpha^q \\
\end{aligned}
$$

where:
$$
\begin{aligned}
\alpha^q &= (a + b \cdot v + c \cdot v^2)^q \\
&= a^q + b^q \cdot v^q + c^q \cdot (v^q)^2 \\
\end{aligned}
$$

How to compute individual items such as $a^q, b^q, c^q, v^q, v^{2q}$ effectively? Let's get it done piece by piece:

- <span style='color:red'>$a^q, b^q, c^q = ?$</span>

Put the conclusion first: $a^q = \bar{a}, b^q = \bar{b}, c^q = \bar{c}$. The above formula can be simplified into:
$$
\begin{aligned}
\alpha^q = \bar{a} + \bar{b} \cdot v^q + \bar{c} \cdot (v^q)^2
\end{aligned}
$$

<br />

<br />

- <span style="color:red">$v^q = ?$</span>

since $v^3 = u, u^2 = \xi, u \in F_{q^2}, \xi \in F_q$, then we have
$$
v^q = v^{3 \cdot (q - 1)/3} \cdot v= u^{\frac{q - 1}{3}} \cdot v = \xi^{\frac{q - 1}{6}} \cdot v= \zeta_6 \cdot v
$$
where $\zeta_6$ is a **primitive 6-th root of unity** in $F_q$.

<br />

*More properties of primitive 6-th root of unity in $F_q$ can see this* :-1:
:::info
- $1 + \zeta_6^2 = \zeta_6$
- $\zeta_6^2 + \zeta_6^4 = -1$
- $\zeta_6^4 + \zeta_6 = 0$
:::

<br />

All together we have:
$$
\alpha^q = \bar{a} + (\bar{b} \cdot \zeta_6) \cdot v + (\bar{c} \cdot \zeta_6^2) \cdot v^2
$$

<br />

Appling Frobenius map we have:
$$
\begin{aligned}
\alpha^{q^2} &= \bar{a}^q + (\bar{b}^q \cdot \zeta_6^q) \cdot v^q + (\bar{c}^q \cdot \zeta_6^{2q}) \cdot v^{2q} \\
&= a + (b \cdot \zeta_6^{q + 1}) \cdot v + (c \cdot \zeta_6^{2q + 2}) \cdot v^2 \\
&= a + (b \cdot \zeta_6^2) \cdot v + (c \cdot \zeta_6^4) \cdot v^2
\end{aligned}
$$

<br />

Apply the property of Squaring Friendly Field we have:
$$
\begin{aligned}
&\alpha^{q^2} \cdot \alpha = \alpha^q \\
&\Updownarrow \\
&(a + (b \cdot \zeta_6^2) \cdot v + (c \cdot \zeta_6^4) \cdot v^2) \cdot (a + b \cdot v + c \cdot v^2) = \bar{a} + (\bar{b} \cdot \zeta_6) \cdot v + (\bar{c} \cdot \zeta_6^2) \cdot v^2 \\
\end{aligned}
$$

expand the equation:
$$
\begin{aligned}
\bar{a} &= (a^2 + bc \zeta_6^2 \cdot u + bc \zeta_6^4 \cdot u) \\
\bar{b} \cdot \zeta_6 &= (ab + ab \zeta_6^2 + c^2 \zeta_6^4 \cdot u) \\
\bar{c} \cdot \zeta_6^2 &= (ac + ac \zeta_6^4 + b^2 \zeta_6^2) \\
\end{aligned}
$$

So the THREE muliplications will be transformed into:
$$
\begin{aligned}
b \cdot c &= \frac{\bar{a} - a^2}{(\zeta_6^2 + \zeta_6^4) \cdot u} = \frac{\bar{a} - a^2}{-u} \\
a \cdot b &= \frac{(\bar{b} + c^2 \cdot u) \cdot \zeta_6}{1 + \zeta_6^2} = \frac{(\bar{b} + c^2 \cdot u) \cdot \zeta_6}{\zeta_6} = c^2 \cdot u + \bar{b} \\
a \cdot c &= \frac{(\bar{c} - b^2) \cdot \zeta_6^2}{1 + \zeta_6^4} = \frac{(\bar{c} - b^2) \cdot \zeta_6^2}{-\zeta_6^2} = b^2 - \bar{c} \\
\end{aligned}
$$

Finally:
$$
\begin{aligned}
A &= a^2 + 2 bc \cdot u = 3a^2 - 2 \bar{a} \\
B &= c^2 \cdot u + 2ab = 3c^2 \cdot u + 2 \bar{b} \\
C &= b^2 + 2 ac = 3 b^2 - 2 \bar{c} \\
\\
\Longrightarrow \alpha^2 &= (3a^2 - \bar{a}) + (3c^2 \cdot u - 2 \bar{b}) \cdot v + (3 b^2 - 2 \bar{c}) \cdot v^2
\end{aligned}
$$

<br />

### Intuition of Field Transition

In order to have intuitive sense, let's briefly introduce the <span style="color:red">field transitions</span> through **Tate/Ate pairings**.

Assume that $P \in \mathbb{G}_1$ defined over $F_p$, while $Q \in \mathbb{G}_2$ defined over $F_{p^{12}}$, actually coordinates of $Q$ must be defined over some subfield of $F_{p^{12}}$ (we will prove for this in later following sections), such as we say $x_Q \in F_{p^6}, x_P \in F_p$:

- Miller Loop

- Double-Add

the line function $f_{r, P}$ still lies where they originally stayed, such as $F_p$:


- Evaluation Line Function

for instance, the single line function:
$$
l_{r, P}(Q) = y_Q - y_T - \alpha \cdot (x_Q - x_T)
$$

the final evaluation result $f_{r, P}(Q)$ will not in $[r]F_{p^{12}}^{\times}$ anymore.


- Final Exponentiation

- Easy Part

push the **Mill Loop** result $f_{r, P}(Q)$ into a special field, what we called **Cyclotomic Group**, $F_{\varPhi_{12}}$:


- Hard Part

push into target field $F_{p^{12}}^{\times}[r]$:


<br />

## Magic Power of Twist

### Why Twist

<br />

Though we have tower fields to represent $F_{p^{12}}$, unfortunately computation over $F_{p^{12}}$ is still expensive. Twist is an advanced trick to solve this problem.

<br />

We can map elements of $F_{p^{12}}$ into a much lower extension degree $F_{p^2}$ with a **sextic-twist**(twist degree is 6), then operations over $F_{p^{12}}$ will be replaced to operations over $F_{p^2}$:
$$
\varphi: F_{p^{12}} \longmapsto F_{p^2}
$$
But how to achieve this?

<br />

### Sextic Twist

<br />

One elliptic curve, denoted as $E_{12}$, defined over high extension degree field $F_{p^{12}}$:
$$
E_{12}: y^2 = x^3 + b
$$

<br />

Another elliptic curve, denoted as $E_2'$, defined over low extension degree field $F_{p^2}$, and has a twist or **isomorphism** relationship with $E_{12}$:
$$
E_2': y'^2 = x'^3 + b \cdot \xi
$$
where $\xi \in F_{p^2}$, $\xi$ is both **non-quadratic** and **non-cubic** residual, namely to say:
$$
\sqrt{\xi} \in F_{p^4}, \sqrt[3]{\xi} \in F_{p^6}
$$

therefore we must have:
$$
\sqrt[6]{\xi} = w \in F_{p^{12}}
$$

<br />

But how to choose a proper $\xi$? It seems like $\xi$ just extends from $F_{p^2}$ to $F_{p^{12}}$. Fortunately, $\beta$ is the very just one!

<br />

:::info
$\bigstar \bigstar \bigstar \bigstar$

According to above tower fields, we can easily have:
$$
F_{p^{2}}[x] \xrightarrow{X^6 - \beta} F_{p^{12}}
$$
:::

<br />

So the twisted curve will be :-1: :
$$
E_2': y'^2 = x'^3 + b \cdot \beta
$$

Is this what we want? After doing some testation you may find that the order of twisted curve maybe is not what we want, it's possible that:
$$
|E_2'| \mod r \ne 0
$$
$r$ is not a prime factor of $|E_2'|$. But what we want is map elements (points) of $E_{12} = E(F_{p^{12}})$ into $E_2' = E'(F_{p^2})$, and this is a one-to-one isomorphism:
$$
\varphi: E(F_{p^{12}}) \longmapsto E(F_{p^2})'
$$
so if $r \nmid |E_2'|$, there's not a cooresponding subgroup with order $r$ in $E_2'$ as $E_{12}$. So must be careful when you are finding a twisted curve (group)!

<br />

*If $\beta$ is not proper, the $\beta^5$ must be the very proper ONE*([other papers](http://indigo.ie/%7Emscott/twists.pdf) also mentioned this), you can try it! So the twisted curve is:
$$
\begin{aligned}
&E_2': y'^2 = x'^3 + b \cdot \xi \\
&\Updownarrow \\
&E_{12}: (\frac{y'}{\sqrt[2]{\xi}})^2 = (\frac{x'}{\sqrt[3]{\xi}})^3 + b
\end{aligned}
$$
where <span style="color:red">$\xi = \beta^5 \in F_{p^2}$ or $\xi = \beta \in F_{p^2}$</span> and $\sqrt[2]{\xi} \in F_{p^4}, \sqrt[3]{\xi} \in F_{p^6}$.

<br />

### Sextic Twist Map

<br />

Last thing on this section is element mapping between these two groups $E_{12}$ and $E_2'$:

- **Twist Operation**

map element of $E_{12}$ to $E_2'$:
$$
\varphi: (x, y) \longmapsto (x \cdot \sqrt[3]{\xi}, y \cdot \sqrt[2]{\xi})
$$

- **Untwist Operation**

map element of $E_2'$ to $E_{12}$:
$$
\varphi^{-1}: (x', y') \longmapsto (\frac{x'}{\sqrt[3]{\xi}}, \frac{y'}{\sqrt[2]{\xi}})
$$

<br />

The cost of conversion is quite cheap, especially after we have the constant ones $\sqrt[2]{\xi}, \sqrt[3]{\xi}$ **precomputed in advance**. Last but not least , we just need to do a few times of untwist operation when neccessary.

<br />

## Magic Power of Frobenius Map

Frobenius Map plays an essential role in arithemtics of high extension fields, especially during **Final Expoentiation** of **Tate/Ate Pairings**.

<br />

### Frobenius Map over $F_{p^2}$

Assume:
$$
\begin{aligned}
F_{p^2} &= F_p[u] / X^2 - \alpha \\
a &= (a_0 + a_1 u) \in F_{p^2} \\
\end{aligned}
$$

where $u^2 = \alpha = -1$ and $a_0, a_1 \in F_p$, we have $a_0^p = a_0, a_1^p = a_1$.

<br />

Then:
$$
\begin{aligned}
a^p &= a_0^p + a_1^p \cdot u^p \\
&= a_0 + a_1 \cdot \alpha^{(p - 1)/2} \cdot u
\end{aligned}
$$

since $(p - 1) / 2$ must be an odd number, so we have:
$$
\begin{aligned}
a^p &= a_0 - a_1 \cdot u \\
a^{p^2} &= a_0 + a_1 \cdot u \\
... \\
a^{p^d} &= a_0 + (-1)^d a_1 \cdot u
\end{aligned}
$$

<br />

**For conclusion, Frobenius Map $\Phi_d(a) = a^{p^d} = a \in F_{p^2}$ if and only if $2 | d$.**

<br />

### Frobenius Map over $F_{p^6}$

Assume:
$$
\begin{aligned}
F_{p^6} &= F_{p^2}[v] / X^3 - \beta \\
a &= (a_0 + a_1 v + a_2 v^2) \in F_{p^6} \\
\end{aligned}
$$

where $v^3 = \beta = u + 1$ and $\beta, a_0, a_1, a_2 \in F_{p^2}$, we have $a_0^p = \bar{a}_0, a_1^p = \bar{a}_1, a_2^p = \bar{a}_2$, and
$$
\begin{aligned}
v^p &= \beta^{(p - 1)/3} \cdot v \\
(v^p)^2 &= \beta^{2 \cdot \frac{p - 1}{3}} \cdot v^2 \\
\beta \cdot \bar{\beta} &= 1 - u^2 = N_2(\beta) \in F_p \\
\\
\end{aligned}
$$

<br />

Then:
$$
\begin{aligned}
a^p &= a_0^p + a_1^p \cdot v^p + a_2^p \cdot (v^p)^2 \\
&= \bar{a}_0 + \bar{a}_1 \cdot \beta^{(p - 1)/3} \cdot v + \bar{a}_2 \cdot \beta^{2 \cdot \frac{p - 1}{3}} \cdot v^2
\\
\\
a^{p^2} &= a_0 + a_1 \cdot \bar{\beta}^{(p - 1)/3} \cdot \beta^{(p - 1)/3} \cdot v + a_2 \cdot \bar{\beta}^{2 \cdot \frac{p - 1}{3}} \cdot \beta^{2 \cdot \frac{p - 1}{3}} \cdot v^2 \\
\\
...\\
\\
a^{p^d} &= C_d(a_0) + C_d(a_1) \cdot N_d(\beta)^{\frac{p - 1}{3}} \cdot v + C_d(a_2) \cdot N_d(\beta)^{2 \cdot \frac{p - 1}{3}} \cdot v^2 \\
\end{aligned}
$$
where $C_d(a_i)$ denotes conjugate $d$ times on $a_i$, and $N_d(\beta)$ denotes norm $d$ times on $\beta$.

<br />

Two aspects need to be considered:

- for $C_d(a_i)$

We can easily have $C_2(a_i) = a_i, C_4(a_i) = a_i, C_6(a_i) = a_i$.

- for $N_d(\beta)$

Since $N_2(\beta) = (1 - u^2) \in F_p$, then $N_4(\beta) = (1 - u^2)^2, N_6(\beta) = (1 - u^2)^3 \in F_p$, so we have:
$$
\begin{aligned}
a^{p^6} &= C_d(a_0) + C_d(a_1) \cdot (1 - u^2)^{p - 1} \cdot v + C_d(a_2) \cdot ((1 - u^2)^{p - 1})^2 \cdot v^2 \\
&= a_0 + a_1 \cdot v + a_2 \cdot v^2 \\
&= a
\end{aligned}
$$

<br />

**For conclusion, $\Phi_d(a) = a^{p^d} = a \in F_{p^6}$ if and only if $6 | d$.**

<br />

### Frobenius Map over $F_{p^{12}}$

Assume:
$$
\begin{aligned}
F_{p^{12}} &= F_{p^6}[w] / X^2 - v \\
a &= (a_0 + a_1 w) \in F_{p^{12}} \\
\end{aligned}
$$

where $w^2 = v, v^3 = u + 1, u^2 = -1$ and $w, a_0, a_1 \in F_{p^6}$.

<br />

Similarily,
$$
\begin{aligned}
a^p &= a_0^p + a_1^p \cdot w^p \\
&= \Phi_1(a_0) + \Phi_1(a_1) \cdot v^{(p - 1)/2} \cdot w \\
\\
a^{p^2} &= \Phi_2(a_0) + \Phi_2(a_1) \cdot (v^{p + 1})^{(p - 1)/2} \cdot w \\
\\
...\\
a^{p^d} &= \Phi_d(a_0) + \Phi_d(a_1) \cdot (v^{p^d + p^{d - 1} + ... + 1})^{(p - 1)/2} \cdot w \\
&= \Phi_d(a_0) + \Phi_d(a_1) \cdot (v^{\frac{p^{d} - 1}{p - 1}})^{(p - 1)/2} \cdot w \\
&= \Phi_d(a_0) + \Phi_d(a_1) \cdot v^{(p^d - 1)/2} \cdot w \\
\end{aligned}
$$

when $d = 12$, since $6 | d$, then we have $\Phi_d(a_i) = a_i$. since $\frac{p^{12} - 1}{2} = (p^6 - 1) \cdot \frac{p^6 + 1}{2}$, and $v \in F_{p^6}$, then we have:
$$
(v^{p^6 - 1})^{\frac{p^6 + 1}{2}} = 1
$$

therefore:
$$
\begin{aligned}
a^{p^{12}}&= a_0 + a_1 \cdot w \\
&= a \\
\end{aligned}
$$

**For conclusion, $\Phi_d(a) = a^{p^d} = a \in F_{p^{12}}$ if and only if $12 | d$.**

<br />

### Frobenius Map and Conjunction

There is a quadratic extension:
$$
F_{q^2} = F_{q}[u] / X^2 - \alpha
$$
where $q = p^m$, assume $a = a_0 + a_1 \cdot u \in F_{q^2} = F_{p^{2m}}$, where $a_0, a_1 \in F_{q} = F_{p^m}$. If we want to do m-times **Frobenius Map** on $a$, then:
$$
\begin{aligned}
a^{p^m} &= a_0^{p^m} + a_1^{p^m} \cdot u^{p^m} \\
&= a_0 + a_1 \cdot \alpha^{(p^m - 1)/2} \cdot u\\
\end{aligned}
$$
since $\alpha$ is non-quadratic residual, namely $\alpha^{(p^m - 1)/2} = -1$, therefore we have:
$$
a^q = a^{p^m} = \bar{a}
$$
It's free cost!

<br />

For instance $(F_{p^{12}})^{p^6} = \bar{F}_{p^{12}}, (F_{p^{6}})^{p^3} = \bar{F}_{p^{6}}, (F_{p^{2}})^{p} = \bar{F}_{p^2}, ...$.

----

# Curve Arithmetics Underneath

This is where **Scalar Multiplication** mainly works. There are two curves (groups) we need to instantiate, $\mathbb{G}_1$ and $\mathbb{G}_2$, in curve BLS12-381:
$$
\begin{aligned}
\mathbb{G}_1 &= E(F_p)[r] \\
\mathbb{G}_2 &= E'(F_{p^2})[r]\\
\end{aligned}
$$
where $\mathbb{G}_1$ is the $r$-torsion curve (subgroup) defined over **Base Prime Field** $F_p$, $\mathbb{G}_2$ is the $r$-torsion **twisted curve** (subgroup) defined over **sextic-twisted field** (relative to $F_{p^{12}}$) $F_{p^{12 / 6}} = F_{p^2}$.

<br />

## Arithmetics of Curve $\mathbb{G}_1$

$$
E(F_p): y^2 = x^3 + 4
$$
where $x, y \in F_p$.

<br />

While curve $\mathbb{G}_1$ also defined over prime field, it's just a $r$-torsion subgroup of $E(F_p)$. Arithmetics (Scalar Multiplication) is the same, nothing special.

<br />

## Arithmetics of Curve $\mathbb{G}_2$

$$
E'(F_{p^2}): y^2 = x^3 + 4 \cdot \beta
$$
where $x, y, \beta \in F_{p^2}, w = \sqrt[6]{\beta} \in F_{p^{12}}$.

<br />

Similarly $\mathbb{G}_2$ is just $r$-torsion subgroup of $E'(F_{p^2})$, Arithmetics (Scalar Multiplication) is over $F_{p^2}$, not $F_p$ anymore. Arithmetics over $F_{p^2}$ has already illustrated above.

<br />

## Arithmetics of Target Group $\mathbb{G}_T$

The target group $\mathbb{G}_T$ is actually not a curve (additional group), but a $r$-torsion subgroup of multiplicative group, represented as $F_{p^k}^{\times}[r]$. So its arithmetics is same as $F_{p^k}^{\times}$.

<br />

----

# Evolution of Pairings

## Weil Reciprocity

$g$ and $f$ are two zero (divisor) functions defined over elliptic curve over finite field, $f, g \in K(E)$, whose divisors are **disjoints each other**, $supp((f)) \land supp((g)) = \emptyset$. Then we have:

$$
g((f)) \equiv f((g))
$$

where $(f)$ denotes divisors of functions $f$, $g((f))$ denotes evaluation of function $g$ on divisor $(f)$, similarly as $f((g))$.

<br />

If we relax the constraints, if $supp((f)) \land supp((g)) \ne \emptyset$, then a more general **Weil Reciprocity** comes:
$$
g((f)) \equiv \epsilon((f), (g)) \cdot f((g))
$$
where $\epsilon((f), (g)) = 1$ if $(f)$ and $(g)$ are disjoint, otherwise $\epsilon((f), (g)) = -1$.

<br />

:::success
Details of general definition of Weil Reciprocity, you can refer THEOREM 3.9 of [Guide to Pairing-based Cryptography](https://www.math.u-bordeaux.fr/~damienrobert/csi/book/book.pdf).
:::

<br />

What's the significance of Weil Reciprocity? It directly gives birth to **Weil Pairing**.

<br />

## Weil Pairing

Assume there are two **linearly disjoint points** of $r$-torsion subgroup defined on curve, $P, Q \in E[r], P \ne \lambda Q$. Upon that, assume $(f) = r \cdot D_P$, and $D_P \equiv (P) - (\mathcal{O})$, similarly $(g) = r \cdot D_Q, D_Q \equiv (Q) - (\mathcal{O})$. Satisfing $Supp(D_P) \land Supp(D_Q) = \emptyset$.

<br />

Then we have:
$$
\begin{aligned}
g_{rD_Q}(r \cdot D_P) &\equiv f_{rD_P}(r \cdot D_Q) \\
&\Downarrow \\
g_{rD_Q}(D_P)^r &\equiv f_{rD_P}(D_Q)^r \\
&\Downarrow \\
(\frac{g_{rD_Q}(D_P)}{f_{rD_P}(D_Q)})^r &\equiv 1 \\
\end{aligned}
$$

<br />

So does **Weil Pairing** come:
$$
\frac{f_{rD_Q}(D_P)}{f_{rD_P}(D_Q)} = \mu_r \in F_{p^k}^{\times}[r]
$$
where $(f_{rD_Q}) = rD_Q, (f_{rD_P}) = rD_P$, and $\mu_r$ is $r$th unity root of multiplicative group $F_{p^k}^{\times}$, namely $\mu_r^r \equiv 1 \mod p^k - 1$.

<br />

#### How to choose proper $D_P$ and $D_Q$

Theoretically we have to choose proper divisors $D_P$ and $D_Q$, you may wonder whether the result $\mu_r$ would be different or not if we choose different $D_P$ and $D_Q$?

Actually the result value of **Weil Pairing** $\mu_r$ is independent of $D_P$ and $D_Q$. Let's prove it!

<br />

Assume $D_{P1}$ and $D_{P2}$ are both equivalent with $(P) - (\mathcal{O})$, then there must exists a divisor $(t)$ satisfing $D_{P1} = D_{P2} + (t)$, then:
$$
\frac{f_{rD_Q}(D_{P1})}{f_{rD_{P1}}(D_Q)} = \frac{f_{rD_Q}(D_{P2}) \cdot f_{rD_Q}((t))}{f_{rD_{P2}}(D_Q) \cdot f_{(t)}(D_Q)^r}
$$

According to **Weil Reciprocity**, since $supp((t)) \land supp(rD_Q) = \emptyset$, we have $f_{rD_Q}((t)) \equiv f_{(t)}(rD_Q) = f_{(t)}(D_Q)^r$. Therefore:
$$
\frac{f_{rD_Q}(D_{P1})}{f_{rD_{P1}}(D_Q)} = \frac{f_{rD_Q}(D_{P2})}{f_{rD_{P2}}(D_Q)}
$$

<br />

Since then we choose the simplest divisors, $D_P = (P) - (\mathcal{O}), D_Q = (Q) - (\mathcal{O})$. According to the general **Weil Reciprocity**, we have:
$$
(-1)^r \cdot \frac{f_{r((Q) - (\mathcal{O}))}((P) - (\mathcal{O}))}{f_{r((P) - (\mathcal{O}))}((Q) - (\mathcal{O}))} = \mu_r \in F_{p^k}^{\times}[r]
$$

<br />

#### How to evaluate divisors $D_P$ and $D_Q$

Evaluation on $D_P = (P) - (\mathcal{O})$ can be simplified further:

$$
f_{rD_Q}(D_P) \equiv f_{rD_Q}(P)
$$

if and only if $P$ and $Q$ are **linearly disjoint** (independent) with each other, $P \ne \lambda Q, \lambda \le r$.

Notice that they are equivalent $\equiv$ not $=$, another words the value maybe different, but it does not have effects on the result of **Weil Pairing** $\mu_r$.

<br />

Therefore **Weil Pairing** becomes:
$$
(-1)^r \cdot \frac{f_{r((Q) - (\mathcal{O}))}(P)}{f_{r((P) - (\mathcal{O}))}(Q)} = \mu_r \in F_{p^k}^{\times}[r]
$$

<br />

**Miller Loop** makes evaluation of function $f_{r((Q) - (\mathcal{O}))}(P)$ practical, obviousely **Weil Pairing** is symmetric, it actually runs two times of **Miller Loop**. Seems like not economic? Actually single one is enough, this is what **Tate Pairing** does.

<br />

#### Algorithm

Directly come from Algorithm 3.2 of [Guide to Pairing-based Cryptography](https://www.math.u-bordeaux.fr/~damienrobert/csi/book/book.pdf):


<br />

## Tate Pairing

You may wonder what does $f_{r((P) - (\mathcal{O}))}(Q)$ look like? Since $P, Q \in E_{F_{p^k}}[r]$, then $P_x, P_y, Q_x, Q_y \in F_{p^k}$, so $f_{r((P) - (\mathcal{O}))}(Q) \in F_{p^k}$. Appling property of coset, **Tate Pairing** comes:
$$
f_{r((P) - (\mathcal{O}))}(Q)^{\frac{p^k - 1}{r}} \equiv \mu_r \in F_{p^k}^{\times}[r]
$$

It has two steps, **Miller Loop** and **Final Exponentiation**. This is also where **Final Exponentiation** comes from.

<br />

#### Formal Definition

Actually **Tate Pairing** has a more formal definition:
$$
e_{T, r}(P, Q): E[r](F_{p^k}) \times E(F_{p^k}) / rE(F_{p^k}) \rightarrow \mu_r
$$
specially $P \in E[r](F_{p^k}), Q \in E(F_{p^k}) / rE(F_{p^k})$, **$Q$ is not element of $r$-torsion subgroup**, it does not lie on $E[r](F_{p^k})$ as $P$ anymore.

It is actually an element of quotient group, any element in $E(F_{p^k})$ which is **linearly disjoint** with $P$.

Since then, how about the evaluation of function $f_{r, P}(Q)$ (result of **Miller Loop**)? It must be also an element of quotient group:

$$
f_{r, P}(Q) \in F_{p^k}^{\times} / r F_{p^k}^{\times}
$$

<br />

Seems like **Tate Pairing** is more general (more relaxed constraints) than **Weil Pairing**, right?

<br />

:::success
Since $P \ne \lambda Q, \lambda \le r$, usually for the convenience of computation(utilization of Frobenius Automorphism) we let $P \in \mathbb{G}_1 = \pi[1]$, namely $\pi(P) = 1 \cdot P$, $\mathbb{G}_1$ is so-called **Base Group**. While $Q \in \mathbb{G}_2 = \pi[p]$, namely $\pi(Q) = [p] Q$, $\mathbb{G}_2$ is so-called **Trace-zero Group**.
:::

<br />

#### Algorithm

Directly come from Algorithm 3.3 of [Guide to Pairing-based Cryptography](https://www.math.u-bordeaux.fr/~damienrobert/csi/book/book.pdf):


<br />

#### Miller Loop

You may notice that, lengths of **Miller Loop** of **Weil Pairing** and **Tate Pairing** are both $\log(r)$ (bit length of $r$).

*Theoretically **Tate Pairing** is practically good enough, more research are focus on practical implementation of pairings later on.* Specially on shorten **Miller Loop** and more efficient **Final Expoentiation**.

Is there any shorter **Miller Loop** algorithms? Yes, but we have to dive into multiplicative group $F_{p^k}^{\times}$ and explore its structure a bit more.

<br />

## Ate Pairing

In **Ate Pairing**, point $Q$ is limited within $\mathbb{G}_2 = \pi[p]$, and point $P$ is limited within $\mathbb{G}_1 = \pi[1]$. This constraints can greatly accelerate pairing computation.

#### Great Property of Miller Algorithm

$$
\begin{aligned}
f_{a + b}, Q(P) &= f_{a, Q}(P) \cdot f_{b, Q}(P) \cdot \frac{l_{[a]Q,[b]Q}}{v_{[a + b]Q}} \\
f_{a \cdot b, Q}(P) &= f_{b, Q}^{a}(P) \cdot f_{a, [b]Q}(P) \\
\end{aligned}
$$

<br />

#### Shorter Miller Loop

$r | p^k - 1$, assume $u = (p^k - 1) // r$, then we have:
$$
f_{p^k, Q}(P) = f_{p^k - 1, Q}(P) \cdot f_{1, Q}(P) \cdot \frac{l_{[p^k - 1]Q,[1]Q}(P)}{v_{[p^k - 1 + 1]Q}(P)}
$$

<br />

Since $Q \in \mathbb{G}_2$, then $[r] Q = \mathcal{O}, [p^k - 1]Q = \mathcal{O}$, we have $l_{[p^k - 1]Q,[1]Q} = v_{[p^k - 1 + 1]Q} = (Q) + (-Q) - 2(\mathcal{O})$, therefore:
$$
f_{p^k, Q}(P) = f_{p^k - 1, Q}(P) = f_{r, Q}(P)^u
$$

It seems like we can replace $r$ into $p$, so it is uncessary since $p > r$.

<br />

If $\lambda \equiv p \mod r$ and $r | p^k - 1$, then $r | \lambda^k - 1$, assume $m = (\lambda^k - 1) // r$, similarly we have:
$$
f_{\lambda^k, Q}(P) = f_{\lambda^k - 1, Q}(P) = f_{r, Q}(P)^m
$$

<br />

According to above great property of **Miller Loop** we have:
$$
\begin{aligned}
f_{\lambda^k, Q} &= f_{\lambda, Q}^{\lambda^{k - 1}} \cdot f_{\lambda^{k - 1}, [\lambda]Q} \\
&= f_{\lambda, Q}^{\lambda^{k - 1}} \cdot f_{\lambda, [\lambda]Q}^{\lambda^{k - 2}} \cdot f_{\lambda^{k - 2}, [\lambda^2]Q} \\
&= f_{\lambda, Q}^{\lambda^{k - 1}} \cdot f_{\lambda, [\lambda]Q}^{\lambda^{k - 2}} \cdot f_{\lambda, [\lambda^2]Q}^{\lambda^{k - 3}} \cdot f_{\lambda^{k - 3}, [\lambda^3]Q} \\
&= ... \\
&= f_{\lambda, Q}^{\lambda^{k - 1}} \cdot f_{\lambda, Q}^{\lambda^{k - 2} \cdot p} \cdot f_{\lambda, Q}^{\lambda^{k - 3} \cdot p^2} \cdot ... \cdot f_{\lambda, Q}^{1 \cdot p^{k - 1}} \\
&= f_{\lambda, Q}^{\sum_{i = 0}^{k - 1} (k - 1 - i) \cdot p^i}
\end{aligned}
$$
let $c = \sum_{i = 0}^{k - 1} (k - 1 - i) \cdot p^i$, we have:
$$
f_{\lambda^k, Q}(P) = f_{\lambda, Q}(P)^c = f_{r, Q}(P)^m
$$

since $\lambda < r$, we can absolutely replace $r$ into $\lambda$, but how to get $\lambda$?

<br />

#### Trace of Frobenius

According to **Hessen Bound**, we have:
$$
|p + 1 - t| = \#E(F_p)
$$
let $T = t - 1$ where $t$ is the Trace of Frobenius, since $r | \#E(F_p)$, then $\lambda = T \equiv p \mod r$.

Finally we have:
$$
\begin{aligned}
a_{T}(P, Q) &= f_{T, Q}(P)^{\frac{p^k - 1}{r}} \\
(f_{T, Q}(P)^c)^{\frac{p^k - 1}{r}} &= (f_{r, Q}(P)^m)^{\frac{p^k - 1}{r}} \\
\Downarrow \\
a_{T}(P, Q)^c &= e_{T, r}(P, Q)^m \\
\end{aligned}
$$

Obviousely we can see the strong relationship between **Tate Pairing** $e_{T, r}(P, Q)$ and **Ate Pairing** $a_{T}(P, Q)$. You may notice that the values of pairings are different, but uniqueness remained.

<br />

What we actually do in **Ate Pairing** is finding a multiple of $r$, $T$ is what we just found, $r | T^k - 1$. But, is $\log{T}$ the shortest length of **Miller Loop**? Maybe (or not)

<br />

## Optimal Ate Pairing


<br />

Appling **Ate Pairing** above, trivally replacing $T$ into $p$, we have:
$$
f_{p, Q}^{x2} = f_{r, Q}^u
$$
where $u = (p^k - 1) // r$, and $x2 = k \cdot p^{k - 1}$.

<br />

According to **Ate Pairing** above, we have:
$$
f_{\lambda, Q}^{x1} = f_{r, Q}^m
$$
where $m = (\lambda^k - 1) // r$, and $x1 = \sum_{i = 0}^{k - 1} (k - 1 - i) \cdot p^i$.

<br />

Upon this, we have:
$$
\begin{aligned}
f_{p, Q}^{x2 \cdot m} &= f_{T, Q}^{x1 \cdot u} \\
\Longrightarrow f_{\lambda, Q}^{x1} &= f_{p, Q}^{\frac{x2 \cdot m}{u}} \\
&= f_{p, Q}^{\frac{k \cdot p^{k - 1} \cdot m}{(p^k - 1) // r}}
\end{aligned}
$$

<br />

Therefore, $x1$ powers of **Ate Pairing**:
$$
a_{\lambda}(P, Q)^{x1} = (f_{\lambda, Q}(P)^{x1})^{\frac{p^k - 1}{r}}
= (f_{p, Q}(P)^{\frac{k \cdot p^{k - 1} \cdot m}{(p^k - 1) // r}})^{\frac{p^k - 1}{r}}$$

Let's leave it alone for a while.

<br />

We generalize the $\lambda = \sum_{i = 0}^l c_i \cdot p^i, l < k$, and $\lambda = m \cdot r$. Same as **Ate Pairing**, finding a multiple of $r$. According to properties of **Miller Loop**, we have:

$$
\begin{aligned}
f_{\lambda, Q} &= f_{\sum_{i = 0}^l c_i \cdot p^i, Q} \\
&= \prod_{i = 0}^l f_{c_i \cdot p^i, Q} \cdot \prod_{j = 0}^{l - 1} \frac{l_{[s_{j + 1}]Q,[c_j \cdot q^j]Q}}{v_{[s_j]Q}} \\
&= \prod_{i = 0}^l f_{p^i, Q}^{c_i} \cdot (\prod_{i = 0}^l f_{c_i, Q}^{p^i} \cdot \prod_{j = 0}^{l - 1} \frac{l_{[s_{j + 1}]Q,[c_j \cdot q^j]Q}}{v_{[s_j]Q}}) \\
\end{aligned}
$$
where applying properties of **Miller Loop**, we have:
$$
\prod_{i = 0}^l f_{p^i, Q}^{c_i} = \prod_{i = 0}^l f_{p, Q}^{i \cdot c_i \cdot p^{i - 1}} = f_{p, Q}^{\sum_{i = 0}^l i \cdot c_i \cdot p^{i - 1}}
$$

<br />

Then **Ate Pairing** is transformed into:
$$
\begin{aligned}
a_{\lambda}(P, Q) &= f_{\lambda, Q}(P)^{\frac{p^k - 1}{r}} \\
&= (f_{p, Q}(P)^{\sum_{i = 0}^l i \cdot c_i \cdot p^{i - 1}})^{\frac{p^k - 1}{r}} \cdot (\prod_{i = 0}^l f_{c_i, Q}(P)^{p^i} \cdot \prod_{j = 0}^{l - 1} \frac{l_{[s_{j + 1}]Q,[c_j \cdot q^j]Q}(P)}{v_{[s_j]Q}(P)})^{\frac{p^k - 1}{r}} \\
&= (f_{p, Q}(P)^{\frac{k \cdot p^{k - 1} \cdot m}{(p^k - 1) // r}})^{\frac{p^k - 1}{r}} \\
\end{aligned}
$$

Obviousely **Ate Pairing** $a_{\lambda}(P, Q)$ is divided into two parts, since the left part $(f_{p, Q}^{\sum_{i = 0}^l i \cdot c_i \cdot p^{i - 1}})^{\frac{p^k - 1}{r}}$ is a trival **Ate Pairing** (length of **Miller Loop** is $\log{p}$), so **the right part $(\prod_{i = 0}^l f_{c_i, Q}^{p^i} \cdot \prod_{j = 0}^{l - 1} \frac{l_{[s_{j + 1}]Q,[c_j \cdot q^j]Q}}{v_{[s_j]Q}})^{\frac{p^k - 1}{r}}$ must be an Ate Pairing too if and only if $\frac{k \cdot p^{k - 1} \cdot m}{(p^k - 1) // r} \ne \sum_{i = 0}^l i \cdot c_i \cdot p^{i - 1}$**.

<br />

So the **Optimal Ate Pairing** comes:
$$
a_{[c_0, c_1,...,c_l]}(P, Q) = (\prod_{i = 0}^l f_{c_i, Q}(P)^{p^i} \cdot \prod_{j = 0}^{l - 1} \frac{l_{[s_{j + 1}]Q,[c_j \cdot q^j]Q}(P)}{v_{[s_j]Q}(P)})^{\frac{p^k - 1}{r}}
$$
where the coefficients $c_i$ are as small as possible.

<br />

No need to worry about powers $p^i$, it's almost free cost after appling few times **Frobenius Map**. What **Optimal Ate Pairing** need to do is concurrently calculuating $f_{c_i, Q}(P)$ and $l_{[s_{j + 1}]Q,[c_j \cdot q^j]Q}(P)$, the length of **Miller Loop** is $max(\log{c_i})$.

<br />

The only problem is how to find these coefficients $c_i$? It is actually a **Lattice Problem**, more details you can refer [Optimal Pairings](https://eprint.iacr.org/2008/096.pdf).

<br />

---

# Python Implementation

## Instantiation of Curve BLS12-381

### Trace/Unti-trace Map

Trace Map:
$$
Tr(Q) = \sum_{i = 0}^k (Q_x^{p^i}, Q_y^{p^i})
$$
where $k$ is the full extension degree. It maps any $r$-torsion points of $E(F_{p^k})[r]$ into $\mathbb{G}_1 = E(F_p)[r]$.

<br />

Untri-trace Map:
$$
UnTr(Q) = [k] Q - Tr(Q)
$$
It maps any $r$-torsion points of $E(F_{p^k})[r]$ into $\mathbb{G}_2$, whose **Trace Map** result is $\mathcal{O}$.

<br />

Since $\mathbb{G}_1$ is defined over $F_p$, so $Q \in \mathbb{G}_1, Tr(Q) = [d] Q$. While $Q \in \mathbb{G}_2, Tr(Q) = \mathcal{O}$, this is where **Trace-zero subgroup** come from.

<br />

def anti_trace_map(point, d, p, E):
    return d * point - trace_map(point, d, p, E)

def trace_map(point, d, p, E):
    result = point
    point_t = point
    for i in range(1, d):
        point_x, point_y = list(point_t)[0], list(point_t)[1]
        point_t = E(point_x ** p, point_y ** p)
        result = result + point_t
    return result

<br />

### Finite Field Conversion

## map element of Fp2 into Fp12
def into_Fp12(e_fp2, beta, F, gen):
    a = beta.polynomial().list()
    if len(a) == 1 :
        a = a + [0]
    e = e_fp2.polynomial().list()
    if len(e) == 1:
        e = e + [0]
    return F(e[0]) + F(e[1]) * (gen ** 6 - F(a[0])) / F(a[1])

## map elements of Fp12 into Fp2 with critical conditions
def into_Fp2(e_fp12, F, gen):
    coef = e_fp12.polynomial().list()
    zero_coeff = [1 for i in range(12) if ((len(coef) > i) and (i != 0) and (i != 6) and (F(coef[i]) == F(0)))]
    assert(reduce(mul, zero_coeff) == 1)
    
    return (F(coef[0]) + F(coef[6])) + gen * F(coef[6])

## map elements of Fp12_t into Fp12
def Fp12_t_into_Fp12(e_fp12_t, F, gen):
    coef = list(e_fp12_t)
    result = []
    for i in range(len(coef)):
        result.append([(F(c) * (((gen ** 6) - F(1)) ** j) * (gen ** i)) for j, c in enumerate(coef[i].polynomial().list())])
    
    return reduce(add, sum(result, [])) 

<br />

### Twist and Untwist

def untwist(x, y, t_x, t_y):
    return x / t_x, y / t_y

def twist(x, y, t_x, t_y):
    return x * t_x, y * t_y

<br />

### Definition of $\mathbb{G}_1$

$\mathbb{G}_1$ denotes curve defined **base prime field**, namely $E(F_p)[r]$

p = 4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129015664037894272559787
q = 52435875175126190479447740508185965837690552500527637822603658699938581184513

A = 0
B = 4

## base prime field
Fp = GF(p)

## E1 over base prime field, map any point on Efp into the q-torsion subgroup
Efp = EllipticCurve(Fp, [A, B])
r_E = Efp.order()
cofactor_E1 = r_E // q
# g_E1 = Efp(0)
# while g_E1 == Efp(0):
#     a = Efp.random_element()
#     g_E1 = cofactor * a
g_E1 = Efp(
    2262513090815062280530798313005799329941626325687549893214867945091568948276660786250917700289878433394123885724147,
    3165530325623507257754644679249908411459467330345960501615736676710739703656949057125324800107717061311272030899084
)
assert(q * g_E1 == Efp(0))
## trace map on E1 is trival, stays on E1
assert(trace_map(g_E1, 12, p, Efp) == 12 * g_E1)

print('\n ##################################### Curve G1: \n cofactor = {}, \n generator = {}, \n order = {} \n'.format(cofactor_E1, g_E1, r_E))

<br />

### Definition of $\mathbb{G}_2'$

$\mathbb{G}_2'$ denotes curve defined over field $F_{p^{k/d}}$, namely $E'(F_{p^{k/d}})[r]$, who is $d$-twisted with $E(F_{p^k})[r]$. In BLS12-381 (sextic-twist), $\mathbb{G}_2' = E'(F_{p^{2}})[r]$.

########## Fp2 = Fp[X] / X^2 - alpha
## alpha = -1
d = 2
alpha = Fp(-1)
X = Fp['X'].gen()
pol2 = X ** d - alpha
assert(pol2.is_irreducible() == True)
Fp2 = GF(p ** d, 'u', modulus = pol2)
u = Fp2.gen()

## Fp12 = Fp2[X] / X^6 - beta
d = 6
beta = u + 1
XX = Fp2['XX'].gen()
pol12 = XX ** d - beta
assert(pol2.is_irreducible() == True)
beta_t = beta 
Efp2_t = EllipticCurve(Fp2, [A, B * beta_t])
## find the proper twisted curve, who has a q-torsion subgroup which is isomorphism with Efpk's one
if Efp2_t.order() % q != 0:
    beta_t = beta ** 5
    Efp2_t = EllipticCurve(Fp2, [A, B * beta_t])

<br />

### Definition of $\mathbb{G}_{12}'$

$\mathbb{G}_{12}'$ denotes twisted curve defined over $F_{p^k}$, namely $E'(F_{p^k})[r]$.

## twist curve E' over Fp12
Fp12_t = Fp2.extension(pol12, 'x')
Efp12_t = Efp2_t.change_ring(Fp12_t)
print('\n Twist curve E defined over Fp12: {}\n'.format(Efp12_t))

<br />

### Definition of $\mathbb{G}_{12}$

$\mathbb{G}_{12}$ denotes curve defined over $F_{p^k}$, namely $E(F_{p^k})[r]$.

## Fp12 = Fp[X] / X^12 - 2X^6 + 2
Fp12 = GF(p ** 12, 'w', modulus = X ** 12 - 2 * (X ** 6) + 2)
w = Fp12.gen()

## constant parameters of twist/untwist 
beta_t_x = w ** 2
beta_t_y = w ** 3

## make sure g_E2 is in the q-torsion subgroup on Efp2_t
r_E2_t = Efp2_t.order()
cofactor_E2_t = r_E2_t // q
# g_E2 = Efp2_t(0)
# while g_E2 == Efp2_t(0):
#     b = Efp2_t.random_element()
#     g_E2 = cofactor_E2_t * b
g_E2 = Efp2_t([
    [
        1265792444950586559339325656560420460408530841056393412024045461464508512562612331578200132635472221512040207420018,
        12405554917932443612178266677500354121343140278261928092817953758979290953103361135966895680930226449483176258412
    ],
    [
        3186142311182140170664472972219788815967440631281796388401764195993124196896119214281909067240924132200570679195848,
        1062539859838502367600126754068373748370820338894390252225574631210227991825937548921368149527995055326277175720251
    ],
])
assert(q * g_E2 == Efp2_t(0))
print('\n #################################### Curve G2: \n cofactor = {}, \n generator = {}, \n order = {} \n'.format(cofactor_E2_t, g_E2, r_E2_t))

## make sure g_E2 is in Fp12 first, uniform the field before untwist
Efp12 = Efp.change_ring(Fp12)
g_E12 = into_E12(g_E2, beta, Fp, w, beta_t_x, beta_t_y, Efp12)

## For the convenience of do Frobenius Map within Fp2, namely (x^p, y^p)
## traditionaly need 3 steps:
## 1. untwist (x, y) to (x', y'), (x', y') = (x / beta_t_x, y / beta_t_y)
## 2. do Frobenius Map within Fp12, (x'^p, y'^p) = (x^p / beta_t_x^p, y^p / beta_t_y^p)
## 3. twist back to (x, y), (x, y) = (x'^p * beta_t_x, y'^p * beta_t_y) = (x^p / beta_t_x^{p - 1}, y^p / beta_t_y^{p - 1})
## 
## Someone may wonder why wouldn't we do Frobenius Map within Fp2 directly? 
## Since one time of Frobenius Map within Fp2, phi(P), may skip out of G2, though P belongs to G2, 
## so we must do it within the FULL EXTENSION Fp12.
##
## Caching beta_t_x^{-(p - 1)} or beta_t_y^{-(p - 1)} would be much preferable
## 
twist_frob_x = into_Fp2(1 / (beta_t_x ** (p - 1)), Fp, u)
twist_frob_y = into_Fp2(1 / (beta_t_y ** (p - 1)), Fp, u)
print('\n Twist parameters: cubic_root(beta_t)^-1 = {}, sqrt(beta_t)^-1 = {} \n'.format(beta_t_x, beta_t_y))
print('\n Twist parameters for Frobenius Map within Fp2: \n cubic_root(beta_t)^-(p - 1) = {}, \n sqrt(beta_t)^-(p - 1) = {} \n'.format(
    twist_frob_x, twist_frob_y
))

print('\n ==================================== DEBUG ====================================\n ')
## make sure g_E12 is in the zero-trace subgroup of q-torsion
assert(q * g_E12 == Efp12(0))
assert(trace_map(g_E12, 12, p, Efp12) == Efp12(0))
print('\n #### UNTWIST: Point of E2 \n {} \n is mapped into E12 \n {} \n successfully! \n'.format(g_E2, g_E12))

## make sure it can be twisted back
x, y = twist(list(g_E12)[0], list(g_E12)[1], beta_t_x, beta_t_y)
x, y = (into_Fp2(x, Fp, u), into_Fp2(y, Fp, u))
assert(Efp2_t(x, y) == g_E2)
print('\n #### TWIST: Point of E12 \n {} \n is mapped into E2 \n {} \n successfully! \n'.format(g_E12, Efp2_t(x, y)))

<br />

## Weil Pairing

### Evaluation of Double-line Function

## evaluation of double line divisor function
## arithmetics on fields, not on multiplicative group
def double_line(line_point, eval_point, E, phi, reverse = False):
    ######################## arithemtic on finite field of line_point
    ## lambda = 3x^2 / 2y
    (x_L, y_L) = (list(line_point)[0], list(line_point)[1])
    (x_E, y_E) = (list(eval_point)[0], list(eval_point)[1])
    alpha = (3 * x_L^2) / (2 * y_L)
    x_2L = alpha * alpha - 2 * x_L
    y_2L = -y_L - alpha * (x_2L - x_L)
    
    ######################## arithmetic on mixed finite field
    ## x_E, y_E \in F2
    ## y_L, x_L, alpha, x_2L \in F1
    if reverse:
        ## evaluation of slop line l_2T
        e_1 = phi(y_E) - y_L - alpha * (phi(x_E) - x_L)
        ## evaluation of vertical line v_2T
        e_2 = phi(x_E) - x_2L
    else:
        ## evaluation of slop line l_2T
        e_1 = y_E - phi(y_L) - phi(alpha) * (x_E - phi(x_L))
        ## evaluation of vertical line v_2T
        e_2 = x_E - phi(x_2L)

    return E(x_2L, y_2L), e_1, e_2

<br />

### Evaluation of Add-line Function

## evaluation of add line divisor function
## arithmetics on fields, not on multiplicative group
def add_line(line_left_point, line_right_point, eval_point, E, phi, reverse = False):
    ######################## arithemtic on finite field of line_point
    ## lambda = (y2 - y1) / (x2 - x1)
    (x_L, y_L) = (list(line_left_point)[0], list(line_left_point)[1])
    (x_R, y_R) = (list(line_right_point)[0], list(line_right_point)[1])
    (x_E, y_E) = (list(eval_point)[0], list(eval_point)[1])
    alpha = (y_L - y_R) / (x_L - x_R)
    x_LR = alpha * alpha - x_L - x_R
    y_LR = -y_L - alpha * (x_LR - x_L)
    
    ######################## arithmetic on mixed finite field
    ## x_E, y_E \in F2
    ## y_L, x_L, alpha, x_LR \in F1
    if reverse:
        ## evaluation of slop line l_{T + P}
        e_1 = phi(y_E) - y_L - alpha * (phi(x_E) - x_L)
        ## evaluation of vertical line v_{T + P}
        e_2 = phi(x_E) - x_LR
    else:
        ## evaluation of slop line l_{T + P}
        e_1 = y_E - phi(y_L) - phi(alpha) * (x_E - phi(x_L))
        ## evaluation of vertical line v_{T + P}
        e_2 = x_E - phi(x_LR)
    
    return E(x_LR, y_LR), e_1, e_2

<br />

### Miller Loop

## Miller Loop of Weil Pairing
def MillerLoop(P, Q, G, q, phi, reverse = False):
    T = P
    f1 = 1
    f2 = 1
    e_bits = [int(i) for i in bin(q)[2:]]
    ## last bit cannot be evaluated, since the slope would be a vertical line
    for i in range(1, len(e_bits)):
        if (i == len(e_bits) - 1) and (e_bits[i] == 0):
            f1 = f1 * (list(Q)[0] - list(T)[0])
            T = 2 * T
            break
        T, e_1, e_2 = double_line(T, Q, G, phi, reverse)
        f1, f2 = (f1 * f1 * e_1, f2 * f2 * e_2)
        if (i == len(e_bits) - 1) and (e_bits[i] == 1):
            f1 = f1 * (list(Q)[0] - list(T)[0])
            T = T + P
            break
        if e_bits[i] == 1:
            T, e_1, e_2 = add_line(T, P, Q, G, phi, reverse)
            f1, f2 = (f1 * e_1, f2 * e_2)
    assert(T == G(0))
    
    return f1 / f2

<br />

### Testation of Weil Pairing

## Weil Pairing Entry
def WeilPairing(P, Qx, G1, G12, q, phi):
    t0 = time.perf_counter()
    f_rP_Q = MillerLoop(P, Qx, G1, q, phi, False)
    t1 = time.perf_counter()
    f_rQ_P = MillerLoop(Qx, P, G12, q, phi, True)
    t2 = time.perf_counter()
    mu_r = ((-1) ** q) * (f_rP_Q / f_rQ_P)
    print('\n ##[Weil Pairing] Time consuming: t[f(P, Qx)] = {:.3f},  t[f(Qx, P)] = {:.3f}'.format(t1 - t0, t2 - t1))
    
    return mu_r


G1, G2_t, G12, G12_t = (Efp, Efp2_t, Efp12, Efp12_t)
C1, C2 = (cofactor_E1, cofactor_E2_t)

## make sure they are in G1 and G2_t repectively
P, Q = (C1 * G1.random_element(), C2 * G2_t.random_element())
assert(q * P == G1(0))
assert(q * Q == G2_t(0))

## untwist from E2_t to E12: Q -> Qx
Qx = into_E12(Q, beta, Fp, w, beta_t_x, beta_t_y, G12)
assert(q * Qx == G12(0))
assert(trace_map(Qx, 12, p, G12) == G12(0))

####################################### Weil Pairing Testation 
## P is defined over E(Fp), Qx is defined over E(Fpk)
## phi maps Fp to Fp12
phi = Hom(Fp, Fp12)(Fp.gen().minpoly().roots(Fp12)[0][0])
assert(P.curve() is not Qx.curve())
mu_r_weil = WeilPairing(P, Qx, G1, G12, q, phi)
## make sure pairing result is in q-torsion subgroup
assert(mu_r_weil ** q == Fp12(1))
#######################################

<br />

Output:

## Time consuming: t[f(P, Qx)] = 0.060,  t[f(Qx, P)] = 0.095

<br />

Obviousely time cost of $f_{r, Q}(P)$ is much more than that of $f_{r, P}(Q)$, since $P$ is defined over **Base Prime Field**, $P \in E(F_p)$, while $Q$ is defined over **Full Extension Field**, $Q \in E(F_{p^k})$.

- Double-add on $F_{p^k}$ is more expensive than on $F_p$
- Function evaluation is absolutely defined over $F_{p^k}$, so this part would be almost equal

<br />

## Tate Pairing

Actually in **Tate Pairing** the vertical line evaluation can be ommited due to the **Final Exponentiation**. Let's prove that!

<br />

Recall twist/ untwist operation:
$$
\begin{aligned}
\varphi: (x', y') \mapsto (x, y) \\
\\
\Longrightarrow \begin{cases}
x = x' \cdot w^2 \\
y = y' \cdot w^3 \\
\end{cases}
\end{aligned}
$$
where $x', y' \in F_{p^2}$, $x, y \in F_{p^{12}}, w \in F_{p^{12}}, w^2 \in F_{p^6}, w^3 \in F_{p^4}$.

<br />

According to definition of embedding degree, $k = 12$ is the **minimal value** satisfying $r | p^k - 1$, namely $q \nmid p^2 - 1, q \nmid p^4 - 1, q \nmid p^6 - 1$, so we must have $p^2 - 1 | (p^{12} - 1) / q, p^4 - 1 | (p^{12} - 1) / q, p^6 - 1 | (p^{12} - 1) / q$.

<br />

Also since $x'^{p^2 - 1} = 1, (w^2)^{p^6 - 1} = 1, (w^3)^{p^4 - 1} = 1$, assuming $(p^{12} - 1) / q = c_1 \cdot (p^2 - 1) = c_2 \cdot (p^6 - 1)$, then we must have $x^{(p^{12} - 1) / q} = (x'^{p^2 - 1})^{c_1} \cdot ((w^2)^{p^6 - 1})^{c_2} = 1$.

<br />

Before untwisting $Q \in \mathbb{G}_2' = E'(F_{p^2})[r]$, after untwisting $Q \in \mathbb{G}_{12} = E(F_{p^{12}})[r]$. The vertical line funcion $x - x_T$, the evaluation would be $x_Q - x_T$, where $x_Q$ is untwisted value and $x_T \in F_p$. Finaly we have $(x_Q - x_T)^{(p^{12} - 1) / q} \equiv 1$.

<br />

### Optimized Evaluation of Double-line Function

## evaluation of double line divisor function
## arithmetics on fields, not on multiplicative group
def double_line(line_point, eval_point, E, phi, reverse = False):
    ######################## arithemtic on finite field of line_point
    ## lambda = 3x^2 / 2y
    (x_L, y_L) = (list(line_point)[0], list(line_point)[1])
    (x_E, y_E) = (list(eval_point)[0], list(eval_point)[1])
    alpha = (3 * x_L^2) / (2 * y_L)
    x_2L = alpha * alpha - 2 * x_L
    y_2L = -y_L - alpha * (x_2L - x_L)
    
    ######################## arithmetic on mixed finite field
    ## x_E, y_E \in F2
    ## y_L, x_L, alpha, x_2L \in F1
    if reverse:
        ## evaluation of slop line l_2T
        e_1 = phi(y_E) - y_L - alpha * (phi(x_E) - x_L)
        # ## evaluation of vertical line v_2T
        # e_2 = phi(x_E) - x_2L
    else:
        ## evaluation of slop line l_2T
        e_1 = y_E - phi(y_L) - phi(alpha) * (x_E - phi(x_L))
        # ## evaluation of vertical line v_2T
        # e_2 = x_E - phi(x_2L)

    return E(x_2L, y_2L), e_1

<br />

### Optimized Evaluation of Add-line Function

## evaluation of add line divisor function
## arithmetics on fields, not on multiplicative group
def add_line(line_left_point, line_right_point, eval_point, E, phi, reverse = False):
    ######################## arithemtic on finite field of line_point
    ## lambda = (y2 - y1) / (x2 - x1)
    (x_L, y_L) = (list(line_left_point)[0], list(line_left_point)[1])
    (x_R, y_R) = (list(line_right_point)[0], list(line_right_point)[1])
    (x_E, y_E) = (list(eval_point)[0], list(eval_point)[1])
    alpha = (y_L - y_R) / (x_L - x_R)
    x_LR = alpha * alpha - x_L - x_R
    y_LR = -y_L - alpha * (x_LR - x_L)
    
    ######################## arithmetic on mixed finite field
    ## x_E, y_E \in F2
    ## y_L, x_L, alpha, x_LR \in F1
    if reverse:
        ## evaluation of slop line l_{T + P}
        e_1 = phi(y_E) - y_L - alpha * (phi(x_E) - x_L)
        # ## evaluation of vertical line v_{T + P}
        # e_2 = phi(x_E) - x_LR
    else:
        ## evaluation of slop line l_{T + P}
        e_1 = y_E - phi(y_L) - phi(alpha) * (x_E - phi(x_L))
        # ## evaluation of vertical line v_{T + P}
        # e_2 = x_E - phi(x_LR)
    
    return E(x_LR, y_LR), e_1

### Optimized Miller Loop

## General Miller Loop Entry
def MillerLoop(P, Q, G, q, phi, reverse = False):
    T = P
    f1 = 1
    f2 = 1
    e_bits = [int(i) for i in bin(q)[2:]]
    
    print('Miller Loop Length: {}'.format(len(e_bits)))
    
    ## last bit cannot be evaluated, since the slope would be a vertical line
    for i in range(1, len(e_bits)):
        if (i == len(e_bits) - 1) and (e_bits[i] == 0):
            f1 = f1 * (list(Q)[0] - list(T)[0])
            T = 2 * T
            break
        T, e_1 = double_line(T, Q, G, phi, reverse)
        f1 = f1 * f1 * e_1
        if (i == len(e_bits) - 1) and (e_bits[i] == 1):
            f1 = f1 * (list(Q)[0] - list(T)[0])
            T = T + P
            break
        if e_bits[i] == 1:
            T, e_1 = add_line(T, P, Q, G, phi, reverse)
            f1 = f1 * e_1
    assert(T == G(0))
    
    return f1

<br />

### Easy-part of Final Exponentiation

For illustration convenience, we does not use **Frobenius Map** trick here, just directly use time-consuming trivial power. Actually it's almost free cost after using **Frobenius Map**.

<br />

## trival implementation of easy part, Frobenius not used here actually
## exp = (p^6 - 1) * (p^2 + 1)
## 2 * Frobenius + 2 * Mul + 1 * Inv
def easy_part(f):
    ff = f
    ## 1 * Frobenius
    t0 = f ** (p ** 6)
    ## 1 * Inv
    t1 = 1 / f
    ## 1 * Mul
    f = t0 * t1
    ## 1 * Frobenius
    t0 = f ** (p ** 2)
    ## 1 * Mul
    f = t0 * f
    
    actual = ff ** (((p ** 6) - 1) * ((p ** 2) + 1))
    assert(actual == f)
    
    return f

<br />

### Hard-part of Final Exponentiation

Same as above, we does not use **Frobenius Map** here.

<br />

As we know, the hard part is arithmetics on **Cyclotomic Subgroup**, namely $F_{\varPhi_{12}}^{\times}$. According to [On the Computation of the Optimal Ate Pairing at the 192-bit Security Level](https://eprint.iacr.org/2016/130.pdf), the power of hard part is not $\frac{p^4 - p^2 + 1}{r}$, but **three times of that**:
$$
f^{3 \cdot \frac{p^4 - p^2 + 1}{r}} = f^{\lambda_0 + \lambda_1 \cdot p + \lambda_2 \cdot p^2 + \lambda_3 \cdot p^3}
$$
where:
$$
\begin{aligned}
\lambda_3 &= x^2 - 2 x + 1 \\
\lambda_2 &= \lambda_3 \cdot x \\
\lambda_1 &= \lambda_2 \cdot x - \lambda_3 \\
\lambda_0 &= \lambda_1 \cdot x + 3 \\
\end{aligned}
$$
In conclusion :
$$
e_{T, r}(P, Q) = (f_{r, P}(Q)^{\frac{p^k - 1}{r}})^3 \ne f_{r, P}(Q)^{\frac{p^k - 1}{r}}
$$

<br />

## reference from Algorithm 1 of "On the Computation of the Optimal Ate Pairing at the 192-bit Security Level"
## trival implementation of hard part, Frobenius not used here actually
## exp = (p^4 - p^2 + 1) / r
def hard_part(f, u, p, q):
    ## 1 * Sqr + 1 * Inv
    t0 = 1 / (f * f)
    ## 1 * Pow
    t5 = f ** u
    ## 1 * Sqr
    t1 = t5 * t5
    ## 1 * Mul
    t3 = t0 * t5
    
    ## 1 * Pow
    t0 = t3 ** u
    
    ## 1 * Pow
    t2 = t0 ** u
    
    ## 1 * Pow
    t4 = t2 ** u
    
    ## 1 * Mul
    t4 = t1 * t4
    ## 1 * Pow
    t1 = t4 ** u
    ## 1 * Inv
    t3 = 1 / t3
    ## 1 * Mul
    t1 = t3 * t1
    ## 1 * Mul
    t1 = t1 * f # f^\lambda_0
    
    # 1 * Inv
    t3 = 1 / f
    ## 1 * Mul
    t0 = t0 * f
    ## 1 * Frobenius
    t0 = t0 ** (p ** 3) # f^\lambda_3
    
    ## 1 * Mul
    t4 = t3 * t4
    ## 1 * Frobenius
    t4 = t4 ** p # f^\lambda_1
    
    ## 1 * Mul
    t5 = t2 * t5
    ## 1 * Frobenius
    t5 = t5 ** (p ** 2) # f^\lambda_2
    
    ## 3 * Mul
    t5 = t5 * t0
    t5 = t5 * t4
    t5 = t5 * t1
    
    ## third power of actual pairing result
    actual = f ** (((p ** 4) - (p ** 2) + 1) // q)
    assert(t5 == actual ** 3)
    assert(t5 ** q == 1)
    
    return t5

<br />

### Final Exponentiation

## Final Exponentiation Entry
def FinalExponentiation(f, p, k, q, u, trivial = True):
    if trivial:
        mu_r = f ** (((p ** k) - 1) // q)
    else:
        t0 = time.perf_counter()
        f = easy_part(f)
        t1 = time.perf_counter()
        mu_r = hard_part(f, u, p, q)
        t2 = time.perf_counter()
        print('\n     ##[Hard Part of Tate Pairing] Time consuming: t[easy] = {:.3f},  t[hard] = {:.3f}'.format(t1 - t0, t2 - t1))
    return mu_r

<br />

### Testation of Tate Pairing

## Tate Pairing Entry
def TatePairing(P, Qx, G1, q, phi, p, k, u, trivial = True):
    t0 = time.perf_counter()
    f = MillerLoop(P, Qx, G1, q, phi, False)
    t1 = time.perf_counter()
    mu_r = FinalExponentiation(f, p, k, q, u, trivial)
    t2 = time.perf_counter()
    print('\n ##[Tate Pairing] Time consuming: t[f(P, Qx)] = {:.3f},  t[exp] = {:.3f}'.format(t1 - t0, t2 - t1))
    
    return mu_r


G1, G2_t, G12, G12_t = (Efp, Efp2_t, Efp12, Efp12_t)
C1, C2 = (cofactor_E1, cofactor_E2_t)

## make sure they are in G1 and G2_t repectively
P, Q = (C1 * G1.random_element(), C2 * G2_t.random_element())
assert(q * P == G1(0))
assert(q * Q == G2_t(0))

## untwist from E2_t to E12: Q -> Qx
Qx = into_E12(Q, beta, Fp, w, beta_t_x, beta_t_y, G12)
assert(q * Qx == G12(0))
assert(trace_map(Qx, 12, p, G12) == G12(0))

####################################### Trivial Tate Pairing Testation
mu_r_tate_1 = TatePairing(P, Qx, G1, q, phi, p, k, True)
assert(mu_r_tate ** q == Fp12(1))
#######################################

####################################### parameter for p(x), q(x), and t(x)
x = -15132376222941642752
t = x + 1
## p = ((x - 1)^2 * (x^4 - x^2 + 1)) / 3 + x
assert((pow((x - 1), 2) * (pow(x, 4) - pow(x, 2) + 1)) // 3 + x == p)
## q = x^4 - x^2 + 1
assert(pow(x, 4) - pow(x, 2) + 1 == q)
## t = x + 1
assert(abs(p + 1 - t) == Efp.order())

####################################### Nontrivial Tate Pairing Testation
mu_r_tate_2 = TatePairing(P, Qx, G1, q, phi, p, k, x, False)
assert(mu_r_tate ** q == Fp12(1))

## The hard part is 3rd power of pairing
assert(mu_r_tate_1 ** 3 == mu_r_tate_2)

<br />

The running output:

Miller Loop Length: 255

 ##[Tate Pairing] Time consuming: t[f(P, Qx)] = 0.039,  t[exp] = 0.079
Miller Loop Length: 255

     ##[Hard Part of Tate Pairing] Time consuming: t[easy] = 0.114,  t[hard] = 0.082

 ##[Tate Pairing] Time consuming: t[f(P, Qx)] = 0.051,  t[exp] = 0.195

After applying **Frobenius Map**, the time cost of final exponentiation would greately reduced.

<br />

## Ate Pairing

### Miller Loop

In **Ate Pairing**, since $[r] P = \mathcal{O}$, $l_{[r - 1]P, P}$ actually is a vertical line, the last step of **Miller Loop** cannot evaluated directly, so we used a specific manner to deal with it.

But in **Ate Pairing**, $[T] P \ne \mathcal{O}$ which is far away from $\mathcal{O}$, no need to worry $l_{[r - 1]P, P}$, so we will strip that specific manner used in **Tate Pairing**.

<br />

## General Miller Loop Entry
def MillerLoop(P, Qx, Qy, G, q, phi, reverse = False):
    ## if power q is negative or not
    P = P if q > 0 else -P
    q = q if q > 0 else -q
    
    T = P
    f1 = 1
    e_bits = [int(i) for i in bin(q)[2:]]

    print('Miller Loop Length: {}'.format(len(e_bits)))
    
    for i in range(1, len(e_bits)):
        ##### strip this specific manner used in Tate Pairing
        # if (i == len(e_bits) - 1) and (e_bits[i] == 0):
        #     f1 = f1 * (list(Q)[0] - list(T)[0])
        #     T = 2 * T
        #     break
        T, e_1 = double_line(T, Qx, Qy, G, phi, reverse)
        f1 = f1 * f1 * e_1
        ##### strip this specific manner used in Tate Pairing
        # if (i == len(e_bits) - 1) and (e_bits[i] == 1):
        #     f1 = f1 * (list(Q)[0] - list(T)[0])
        #     T = T + P
        #     break
        if e_bits[i] == 1:
            T, e_1 = add_line(T, P, Qx, Qy, G, phi, reverse)
            f1 = f1 * e_1
    
    return f1

<br />

### Testation of Ate Pairing

Notice that in curve BLS12-381, the parameter $x$ for polynomials $p(x), q(x), t(x)$ is a **negative one**:

$$
\begin{aligned}
q(x) &= x^4 - x^2 + 1 \\
p(x) &= (x - 1)^2 \cdot q(x) \cdot \frac{1}{3} + x \\
t(x) &= x + 1
\end{aligned}
$$

where $x = -15132376222941642752$.

<br />

Therefore we must deal with it properly in **Miller Loop** before looping.

<br />

## Ate Pairing Entry
def AtePairing(P, Qx, Qy, G1, q, phi, p, k, u, T):
    t0 = time.perf_counter()
    f = MillerLoop(P, Qx, Qy, G1, T, phi, False)
    t1 = time.perf_counter()
    mu_r = FinalExponentiation(f, p, k, q, u)
    t2 = time.perf_counter()
    print('\n ##[Ate Pairing] Time consuming: t[f(P, Qx)] = {:.3f},  t[exp] = {:.3f}'.format(t1 - t0, t2 - t1))
    
    return mu_r

G1, G2_t, G12, G12_t = (Efp, Efp2_t, Efp12, Efp12_t)
C1, C2 = (cofactor_E1, cofactor_E2_t)

## make sure they are in G1 and G2_t repectively
P, Q = (C1 * G1.random_element(), C2 * G2_t.random_element())
assert(q * P == G1(0))
assert(q * Q == G2_t(0))

## untwist from E2_t to E12: Q -> Qx
# Qx = into_E12(Q, beta, Fp, w, beta_t_x, beta_t_y, G12)
# assert(q * Qx == G12(0))
# assert(trace_map(Qx, 12, p, G12) == G12(0))
# Px = into_E2_t(P, Fp, u, beta_t_x, beta_t_y, G2_t)
Px_t, Py_t = twist(P.xy()[0], P.xy()[1], beta_t_x, beta_t_y)

## parameter for p(x), q(x), and t(x)
x = -15132376222941642752
t = x + 1
## p = ((x - 1)^2 * (x^4 - x^2 + 1)) / 3 + x
assert((pow((x - 1), 2) * (pow(x, 4) - pow(x, 2) + 1)) // 3 + x == p)
## q = x^4 - x^2 + 1
assert(pow(x, 4) - pow(x, 2) + 1 == q)
## t = x + 1
assert(abs(p + 1 - t) == Efp.order())

## p \equiv T \mod q
T = t - 1
####################################### Ate Pairing Testation
phi = Hom(Fp2, Fp12)(Fp2.gen().minpoly().roots(Fp12)[0][0])
mu_r_ate = AtePairing(Q, Px_t, Py_t, G2_t, q, phi, p, k, x, T)
assert(mu_r_ate ** q == Fp12(1))

<br />

The running output:

Miller Loop Length: 64

     ##[Hard Part of Ate Pairing] Time consuming: t[easy] = 0.106,  t[hard] = 0.073

 ##[Ate Pairing] Time consuming: t[f(P, Qx)] = 0.025,  t[exp] = 0.179

Obviousely time cost of **Miller Loop** is greatly reduced, since $\log{T}$ is far more less than $\log{q}$ (64 vs 255).

<br />

---

# Rust Implementation

Much testation work need to be done, [code](https://github.com/PayneJoe/crypto_research/tree/main/ecc/src) to be updated...

<br />

---

# References

[1] [A note on twists for pairing friendly curves](http://indigo.ie/%7Emscott/twists.pdf)

[2] [Pairing-Friendly Elliptic Curves of Prime Order](https://eprint.iacr.org/2005/133.pdf)

[3] [Pairing for Beginners](https://static1.squarespace.com/static/5fdbb09f31d71c1227082339/t/5ff394720493bd28278889c6/1609798774687/PairingsForBeginners.pdf)

[4] [Guide to Pairing-based Cryptography](https://www.math.u-bordeaux.fr/~damienrobert/csi/book/book.pdf)

[5] [Faster pairing computations on curves with high-degree twists](https://www.iacr.org/archive/pkc2010/60560227/60560227.pdf)

[6] [Optimal Pairings](https://eprint.iacr.org/2008/096.pdf)

[7] [On the Computation of the Optimal Ate Pairing at the 192-bit Security Level](https://eprint.iacr.org/2016/130.pdf)

[8] [Efficient Final Exponentiation via Cyclotomic Structure for Pairings over Families of Elliptic
Curves](https://eprint.iacr.org/2020/875.pdf)

[9] [Faster Squaring in the Cyclotomic Subgroup of Sixth Degree Extensions](https://eprint.iacr.org/2009/565.pdf)

[10] [A Guide to Plane Algebraic Curves]()