Computing the Optimal Ate Pairing over the BN254 Curve

The post BN254 For The Rest Of Us constitutes the base point for this article, so I highly recommend reading it to get a solid background of the BN254 curve. This post will reuse some parts with an emphasy on the computational aspect of the optimal Ate pairing over the BN254 curve.

For those with a poor background on Pairing-Based Cryptography (PBC) but with a decent background on Elliptic-Curve Cryptoraphy (ECC) I recommend reading the following resources in order:

1. Pairings For Beginners.
2. Pairing-Friendly Elliptic Curves of Prime Order.
3. High-Speed Software Implementation of the Optimal Ate Pairing over Barreto–Naehrig Curves.

Only after (and not before) you feel sufficiently comfortable with both PBC and ECC come back to this post and start to make sense out of it.

The implementation in zkASM can be found here.

The Curve

The BN254 is the elliptic curve

$E$ of the form:

$$y^2=x^3+3$$
defined over a finite field

${\mathbb{F}}_p$, where the prime

$p$ is the following

$254$-bit number:

p = 21888242871839275222246405745257275088696311157297823662689037894645226208583

The number of points

$\mathrm{\#}E({\mathbb{F}}_p)=p+1-t$ happens to also be a prime

$r$ of

$254$ bits:

r = 21888242871839275222246405745257275088548364400416034343698204186575808495617

where t = 147946756881789318990833708069417712967 is known as the trace of Frobenius.

The point

$(1,2)\in E({\mathbb{F}}_p)$ is a point of order

$r$, i.e. a generator of the group of points of the curve. In fact, all points in

$E({\mathbb{F}}_p)\mathrm{\setminus }\{\mathcal{O}\}$ are generators, as Lagrange's Theorem shows.

The embedding degree

$k$ of this curve is

$12$.

For more information, refer to BN254 For The Rest Of Us or the EIPs 196 and 197.

Tower of Extension Fields

To represent

${\mathbb{F}}_{p^{12}}$ we use the following tower of extension fields:
fields:

$$\begin{array}{rl}{\mathbb{F}}_{p^2}& ={\mathbb{F}}_p[u]/(u^2-\beta ),\\ {\mathbb{F}}_{p^6}& ={\mathbb{F}}_{p^2}[v]/(v^3-\xi ),\\ {\mathbb{F}}_{p^{12}}& ={\mathbb{F}}_{p^6}[w]/(w^2-v),\end{array}$$
where

$\beta$ is not a quadratic residue in

${\mathbb{F}}_p$ and

$\xi$ is neither a quadratic residue or a cubic residue in

${\mathbb{F}}_{p^2}$.

In the particular case of the BN254, we have

$\beta =-1$ and

$\xi =9+u$ and therefore:

$$\begin{array}{rl}{\mathbb{F}}_{p^2}& ={\mathbb{F}}_p[u]/(u^2+1),\\ {\mathbb{F}}_{p^6}& ={\mathbb{F}}_{p^2}[v]/(v^3-(9+u)),\\ {\mathbb{F}}_{p^{12}}& ={\mathbb{F}}_{p^6}[w]/(w^2-v),\end{array}$$

from which we can extract the identities

$u^2=-1$,

$w^2=v$,

$v^3=9+u$ and consequently the idenity

$w^6=9+u$.

Using the latter identity, we can also write:

$${\mathbb{F}}_{p^{12}}={\mathbb{F}}_{p^2}[w]/(w^6-(9+u))$$

These two representations of

${\mathbb{F}}_{p^{12}}$ will allow us to apply some optimizations:

1. Writing any element
   $f\in {\mathbb{F}}_{p^{12}}$ as
   $f=g+hw$, with
   $g,h\in {\mathbb{F}}_{p^6}$ implies ^[1]
   $f^{p^6}=\overline{f}=g-hw$. Thus, computing the
   $p^6$-th power of any element in
   ${\mathbb{F}}_{p^{12}}$ is computationally "free".
2. Now, if we write
   $g=g_0+g_{1}v+g_2{v}^2$ and
   $h=h_0+h_{1}v+h_2{v}^2$ where
   $g_0$,
   $g_1$,
   $g_2$,
   $h_0$,
   $h_1$,
   $h_2\in {\mathbb{F}}_{p^2}$ we have:
   
   $$f=g_0+h_{0}w+g_1{w}^2+h_1{w}^3+g_2{w}^4+h_2{w}^5,$$
   which will allow us to compute
   $f^p,f^{p^2}$ and
   $f^{p^3}$ (needed for the final exponentiation) very efficiently.

Twist

BN curves admits a sextic twist

$E^{\prime }$ defined over

${\mathbb{F}}_{p^2}$ of the form:

$$y^{\prime 2}=x^{\prime 3}+3/(9+u)$$
or more specifically:

(y')² = (x')³ + (19485874751759354771024239261021720505790618469301721065564631296452457478373 + 266929791119991161246907387137283842545076965332900288569378510910307636690·u)

This means that pairing computations can be restricted to points

$P\in E({\mathbb{F}}_p)$ and

$Q\in E^{\prime }({\mathbb{F}}_{p^2})$ and avoid computations over

${\mathbb{F}}_{p^{12}}$ entirely.

The number of points of

$E^{\prime }$ is:

$$\mathrm{\#}{E}^{\prime }({\mathbb{F}}_{p^2})=(p+1-t)(p-1+t)=r\cdot c$$
which concrete to:

#E' = 479095176016622842441988045216678740799252316531100822436447802254070093686356349204969212544220033486413271283566945264650845755880805213916963058350733
c = 21888242871839275222246405745257275088844257914179612981679871602714643921549

The mapping

$$\mathrm{\Psi }(x^{\prime },y^{\prime })=(w^2{x}^{\prime },w^3{y}^{\prime })$$
defines a group homomorphism from

$E^{\prime }({\mathbb{F}}_{p^2})$ to

$E({\mathbb{F}}_{p^{12}})$ and we are going to use it to send points from the twisted curve

$E^{\prime }$ to the curve

$E$ (but not the other way around).

The Pairing

The optimal Ate pairing

$e:{\mathbb{G}}_1\times {\mathbb{G}}_2\to {\mathbb{G}}_T$ over the BN254 curve is defined by:

$$e(P,Q)={(f_{6x+2,Q}(P)\cdot {\ell }_{[6x+2]\mathrm{\Psi }(Q),{\pi }_p(\mathrm{\Psi }(Q))}(P)\cdot {\ell }_{[6x+2]\mathrm{\Psi }(Q)+{\pi }_p(\mathrm{\Psi }(Q)),-{\pi }_p^2(\mathrm{\Psi }(Q))}(P))}^{\frac{p^{12}-1}{r}}$$
where:

• ${\mathbb{G}}_1=E({\mathbb{F}}_p)[r]$ is the only subgroup of the
  $r$-torsion of
  $E$ over
  ${\mathbb{F}}_p$.

• ${\mathbb{G}}_2=E^{\prime }({\mathbb{F}}_{p^2})[r]$ is the only subgroup of the
  $r$-torsion of
  $E^{\prime }$ over
  ${\mathbb{F}}_{p^2}$.

• ${\mathbb{G}}_T={\mu }_r$ is the set of
  $r$-th roots of unity over
  ${\mathbb{F}}_{p^{12}}^{\ast }$.

• x = 4965661367192848881 and therefore 6x+2 = 29793968203157093288, which is a number of

  $65$ bits that can be expressed in the base
  $\{-1,0,1\}$ as follows:
  ​​​​2^3 + 2^5 - 2^7 + 2^10 - 2^11 + 2^14 + 2^17 + 2^18 - 2^20 + 2^23 - 2^25 + 2^30 + 2^31 + 2^32 - 2^35 + 2^38 - 2^44 + 2^47 + 2^48 - 2^51 + 2^55 + 2^56 - 2^58 + 2^61 + 2^63 + 2^64
  ​​​​[0, 0, 0, 1, 0, 1, 0, -1, 0, 0, 1, -1, 0, 0, 1, 0, 0, 1, 1, 0, -1, 0, 0, 1, 0, -1, 0, 0, 0, 0, 1, 1, 1, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 1, 1, 0, 0, -1, 0, 0, 0, 1, 1, 0, -1, 0, 0, 1, 0, 1, 1]
  ​​​​0001010-1001-100100110-10010-1000011100-100100000-1001100-1000110-1001011
  

  which has a Hamming weight of 26.

  We prefer base

  $\{-1,0,1\}$ over base
  $\{0,1\}$ because 6x+2 expressed in the latter is:
  ​​​​2^3 + 2^5 + 2^7 + 2^8 + 2^9 + 2^11 + 2^12 + 2^13 + 2^17 + 2^18 + 2^20 + 2^21 + 2^22 + 2^25 + 2^26 + 2^27 + 2^28 + 2^29 + 2^31 + 2^32 + 2^35 + 2^36 + 2^37 + 2^44 + 2^45 + 2^46 + 2^48 + 2^51 + 2^52 + 2^53 + 2^54 + 2^56 + 2^58 + 2^59 + 2^60 + 2^63 + 2^64 
  ​​​​[0,0,0,1,0,1,0,1,1,1,0,1,1,1,0,0,0,1,1,0,1,1,1,0,0,1,1,1,1,1,0,1,1,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,1,0,0,1,1,1,1,0,1,0,1,1,1,0,0,1,1]
  ​​​​00010101110111000110111001111101100111000000111010011110101110011
  

  which has a Hamming weigth of 37.

  I moved from base

  $\{0,1\}$ to base

  $\{-1,0,1\}$ by using the identity^[2]

  $$2^{n+1}-2^m=2^n+2^{n-1}+...+2^m,$$ that holds for any

  $n,m\in \mathbb{N}$ such that

  $n\ge m$. In words, it is unpacking groups of

  $1$'s of any size to a single

  $1$ and

  $-1$ in the binary decomposition.

• $f_{i,Q}$, for any
  $i\in \mathbb{N}$ and
  $Q\in {\mathbb{G}}_2$, is a
  ${\mathbb{F}}_{p^{12}}$ rational function on
  $E$ that is computed using Miller's "double-and-add"-style algorithm:
  
  $$\begin{array}{r}{f}_{i+j,Q}=f_{i,Q}\cdot f_{j,Q}\cdot \frac{{\ell }_{[i]\mathrm{\Psi }(Q),[j]\mathrm{\Psi }(Q)}}{v_{[i+j]\mathrm{\Psi }(Q)}},\\ f_{i+1,Q}=f_{m,Q}\cdot \frac{{\ell }_{[i]\mathrm{\Psi }(Q),\mathrm{\Psi }(Q)}}{v_{[i+1]\mathrm{\Psi }(Q)}},\end{array}$$
  starting with
  $f_{1,Q}=1$.

• ${\ell }_{P_1,P_2}$ denotes the equation of a line passing through
  $P_1$ and
  $P_2$ while
  $v_P$ denotes the equation of a vertical line passing through
  $P$. However, due to the final exponentiation, the value of
  $e(P,Q)$ remains unchanged if we omit the division by the vertical lines, so we will not care about such computation.

• ${\pi }_p^i(x,y)=(x^{p^i},y^{p^i})$ is known as the Frobenius operator.

Therefore, what we want to compute is the following:

bound = [0, 0, 0, 1, 0, 1, 0, -1, 0, 0, 1, -1, 0, 0, 1, 0, 0, 1, 1, 0, -1,
         0, 0, 1, 0, -1, 0, 0, 0, 0, 1, 1, 1, 0, 0, -1, 0, 0, 1, 0, 0, 0,
         0, 0, -1, 0, 0, 1, 1, 0, 0, -1, 0, 0, 0, 1, 1, 0, -1, 0, 0, 1, 0,
         1, 1] 

1101111010011011000001001001000111011101011100001011001111110111111100000111110100101101000001011101000011111001111110101000110

         # This is 6x+2 in base {-1,0,1}, which has a Hamming weight of 26

def e(P,Q):
    assert(subgroup_check_G1(P))
    assert(subgroup_check_G2(Q))
    if ((P == 0) or (Q == 0)) return 1 # Here, 0 is the identity element of E
    
    R = Q
    f = 1
    for i in range(len(bound)-2,-1,-1): # len(bound)=65
        f = f * f * line(twist(R),twist(R),P)
        R = 2 * R
        if bound[i] == 1:
            f = f * line(twist(R),twist(Q),P)
            R = R + Q
        elif bound[i] == -1:
            f = f * line(twist(R),twist(-Q),P)
            R = R - Q
            
    Qp = FrobEndo(Q); # Q'
    f = f * line(twist(R), twist(Qp), P);
    R = R + Qp;
    Qpp = -FrobEndo(Qp); # Q''
    f = f * line(twist(R), twist(Qpp), P);

    return final_exponentiation(f)

In what follows we explain how to compute each of the components of the previous pseudocode snippet.

Subgroup Checks

To summary BN254 For The Rest Of Us, we have:

• $P\in {\mathbb{G}}_1=E({\mathbb{F}}_p)[r]$ if and only if
  $P\in E({\mathbb{F}}_p)$. That is, if we write
  $P=(x,y)$ then we simply need to ensure that the equation
  $y^2=x^3+3$ holds over
  ${\mathbb{F}}_p$.

  Costs:

  $2s+1m+1a$

• A point

  $Q=(x^{\prime },y^{\prime })\in E^{\prime }({\mathbb{F}}_{p^2})$ belongs to
  ${\mathbb{G}}_2=E^{\prime }({\mathbb{F}}_{p^2})[r]$ if and only if
  $\psi (Q)=[6x^2]Q$, where
  $\psi :E^{\prime }({\mathbb{F}}_{p^2})\to E^{\prime }({\mathbb{F}}_{p^2})$ is defined as:
  
  $$\psi (x^{\prime },y^{\prime })=({\gamma }_{1,2}{\overline{x}}^{\prime },{\gamma }_{1,3}{\overline{y}}^{\prime }),$$
  where
  ${\gamma }_{1,2}=(9+u{)}^{\frac{p-1}{3}},{\gamma }_{1,3}=(9+u{)}^{\frac{p-1}{2}}$.

  This was proven in Proposition 3 of 2022/352 and it is clearly much faster than the naive check

  $[r]Q=\mathcal{O}$, since the constants
  ${\gamma }_{1,2},{\gamma }_{1,3}$ can be precomputed.

  TODO:
  A more efficient approach appears in 2022/348, which ensures that

  $Q\in E^{\prime }({\mathbb{F}}_{p^2})$ belongs to

  ${\mathbb{G}}_2$ if and only if:
  

  $$[x+1]Q+\psi ([x]Q)+{\psi }^2([x]Q)={\psi }^3([2x]Q)$$

Costs:

$2\tilde{c}+2\tilde{m}+\text{cost of }[6x^2]P$

Line Equations

Different Points

Given two points

$Q_1=(x_1^{\prime },y_1^{\prime })$ and

$Q_2=(x_2^{\prime },y_2^{\prime })$ in the twisted curve

$E^{\prime }({\mathbb{F}}_{p^2})$ such that

$x_1^{\prime }\ne y_1^{\prime }$ and

$x_2^{\prime }\ne y_2^{\prime }$ and a point

$P=(x,y)$ in the curve

$E({\mathbb{F}}_p)$, the evaluation at

$P$ of the line passing through

$\mathrm{\Psi }(Q_1)$ and

$\mathrm{\Psi }(Q_2)$ is given by:

$${\ell }_{\mathrm{\Psi }(Q_1),\mathrm{\Psi }(Q_2)}(P)=w^2(x_2^{\prime }-x_1^{\prime })y+w^3(y_1^{\prime }-y_2^{\prime })x+w^5(x_1^{\prime }{y}_2^{\prime }-x_2^{\prime }{y}_1^{\prime })$$

Costs:

$3\tilde{a}+2\tilde{m}+4m$

The Same Points

Given a point

$Q=(x^{\prime },y^{\prime })$ in the twisted curve

$E^{\prime }({\mathbb{F}}_{p^2})$ and a point

$P=(x,y)$ in the curve

$E({\mathbb{F}}_p)$, the evaluation at

$P$ of the line passing through

$\mathrm{\Psi }(Q)$ is given by:

$${\ell }_{\mathrm{\Psi }(Q),\mathrm{\Psi }(Q)}(P)=(3x^{\prime 3}-2y^{\prime 2})(9+u)+w^3(2y{y}^{\prime })+w^4(-3x{x}^{\prime 2})$$

Costs:

$1\tilde{a}+2\tilde{s}+2\tilde{m}+10m$

The Last two Lines

Now we explain some tricks for the computation of the lines:

$${\ell }_{[6x+2]\mathrm{\Psi }(Q),{\pi }_p(\mathrm{\Psi }(Q))}(P),{\ell }_{[6x+2]\mathrm{\Psi }(Q)+{\pi }_p(\mathrm{\Psi }(Q)),-{\pi }_p^2(\mathrm{\Psi }(Q))}(P).$$

• For the first line notice that if

  $Q=(x^{\prime },y^{\prime })\in E^{\prime }({\mathbb{F}}_{p^2})$, then:
  
  $${\pi }_p(\mathrm{\Psi }(Q))=((w^2{x}^{\prime }{)}^p,(w^3{y}^{\prime }{)}^p)=(w^2({\gamma }_{1,2}{x}^{\prime p}),w^3({\gamma }_{1,3}{y}^{\prime p}))=\mathrm{\Psi }({\gamma }_{1,2}{x}^{\prime p},{\gamma }_{1,3}{y}^{\prime p})=\mathrm{\Psi }({\gamma }_{1,2}{\overline{x}}^{\prime },{\gamma }_{1,3}{\overline{y}}^{\prime }),$$
  where
  ${\gamma }_{1,2}=(9+u{)}^{\frac{p-1}{3}},{\gamma }_{1,3}=(9+u{)}^{\frac{p-1}{2}}$ are precomputed.

  The conjugate of an element in

  $f=a+bu$ in

  ${\mathbb{F}}_{p^2}$ is

  $\overline{f}=a-bu$.

  Moreover, due to

  $\mathrm{\Psi }$ being an homomorphism we have that
  $[n]\mathrm{\Psi }(Q)=\mathrm{\Psi }([n]Q)$ for all
  $n\in {\mathbb{F}}_r$. Consequently, we compute
  ${\ell }_{\mathrm{\Psi }([6x+2]Q),\mathrm{\Psi }(Q^{\prime })}(P)$ as in the previous section, with
  $Q^{\prime }=({\gamma }_{1,2}{\overline{x}}^{\prime },{\gamma }_{1,3}{\overline{y}}^{\prime })=(x_1,y_1)$ (it's a good exercise to check that
  $Q^{\prime }$ belongs to
  $E$).

  Costs:

  $\underset{Q^{\prime }}{\underbrace{2\tilde{c}+2\tilde{m}}}+{\ell }_{\text{diff}}$

  I don't include the cost of computing

  $[6x+2]Q$, since it is obtained in the Miller loop.

• For the second line we similarly have:
  

  $$-{\pi }_p^2(\mathrm{\Psi }(Q))=-{\pi }_p({\pi }_p(\mathrm{\Psi }(Q)))=-{\pi }_p(\mathrm{\Psi }(x_1,y_1))=\mathrm{\Psi }({\gamma }_{1,2}\overline{x_1},-{\gamma }_{1,3}\overline{y_1})$$
  Again, due to the homomorphic property of
  $\mathrm{\Psi }$, we compute
  ${\ell }_{\mathrm{\Psi }([6x+2]Q+Q^{\prime }),\mathrm{\Psi }(Q^{″})}(P)$ as in the previous section, with
  $Q^{″}=({\gamma }_{1,2}\overline{x_1},-{\gamma }_{1,3}\overline{y_1})$ (again, check that
  $Q^{″}$ belongs to
  $E$).

  Costs:

  $1E_{\text{add}}^{\prime }+\underset{Q^{″}}{\underbrace{2\tilde{c}+2\tilde{m}+2m}}+{\ell }_{\text{diff}}$

Hence, in both cases, the Frobenius operator comes for free.

Final Exponentiation

In the final exponentiation, we raise an element

$f\in {\mathbb{F}}_{p^{12}}$ to the power

$(p^{12}-1)/r$.

First, divide the exponent as follows:

$$\frac{p^{12}-1}{r}=\underset{\text{easy part}}{\underbrace{(p^6-1)(p^2+1)}}\underset{\text{hard part}}{\underbrace{\frac{p^4-p^2+1}{r}}}$$

The Easy Part

Let's start with the easy part of the easy part, when we first attempt to compute

$f^{p^6-1}$. But since,

$f^{p^6}=\overline{f}$, we have

$f^{p^6-1}=\overline{f}\cdot f^{-1}$, just one conjugate, one inversion and one multiplication in

${\mathbb{F}}_{p^{12}}$.

The conjugate of an element in

$f=a+bw$ in

${\mathbb{F}}_{p^{12}}$ is

$\overline{f}=a-bw$. Hence, if we write

$f=g_0+h_{0}w+g_1{w}^2+h_1{w}^3+g_2{w}^4+h_2{w}^5$ we have:

$$\overline{f}=g_0-h_{0}w+g_1{w}^2-h_1{w}^3+g_2{w}^4-h_2{w}^5.$$

To finish the easy part, we first compute

$(f^{(p^6-1)}{)}^{p^2}$ which is the Frobenius operator applied to

$f^{p^6-1}$ and then multiply the result by

$f^{p^6-1}$ to finally obtain

$m:=f^{(p^6-1)(p^2+1)}\in {\mathbb{F}}_{p^{12}}$.

Importantly, after the computation of the easy part, the resulting field element becomes a member of the cyclotomic subgroup

${\mathbb{G}}_{{\varphi }_6(p^2)}={\mathbb{G}}_{{\varphi }_{12}(p)}=\{a\in {\mathbb{F}}_{p^{12}}\mid a^{{\varphi }_{12}(p)}=a^{p^4-p^2+1}=a^{p^6+1}=1\}$, which means that inversion from now on can be computed by simply conjugation. Moreover, when an element

$a$ belongs to

${\mathbb{G}}_{{\varphi }_6(p^2)}$, one can compute

$a^2$ using fewer multiplications and therefore the overall costs of performing exponentiations (see this section).

Costs:

$\underset{f^{p^6-1}}{\underbrace{1C+1I+1M}}+5\tilde{m}+1M$

The Hard Part

Following Section 5 of this 2008/490, we can write:

$$\frac{p^4-p^2+1}{r}={\lambda }_0+{\lambda }_{1}p+{\lambda }_2{p}^2+{\lambda }_3{p}^3,$$
where:

• ${\lambda }_0=-2-18x-30x^2-36x^3$.
• ${\lambda }_1=1-12x-18x^2-36x^3$.
• ${\lambda }_2=1+6x^2$.
• ${\lambda }_3=1$.

So we just need to compute:

$$m^{{\lambda }_0},m^{{\lambda }_{1}p},m^{{\lambda }_2{p}^2},m^{{\lambda }_3{p}^3},$$
and then multiply them all.

In what follows, remember that x = 4965661367192848881, which is a number of

$63$ bits that can be expressed in the base

$\{-1,0,1\}$ as follows:

1 - 2^4 + 2^9 + 2^11 + 2^16 + 2^19 + 2^21 + 2^22 + 2^25 + 2^27 + 2^30 + 2^34 + 2^36 + 2^37 + 2^39 + 2^41 + 2^44 + 2^47 + 2^48 + 2^51 - 2^53 + 2^56 + 2^58 + 2^62
[1,0,0,0,-1,0,0,0,0,1,0,1,0,0,0,0,1,0,0,1,0,1,1,0,0,1,0,1,0,0,1,0,0,0,1,0,1,1,0,1,0,1,0,0,1,0,0,1,1,0,0,1,0,-1,0,0,1,0,1,0,0,0,1]
1000-1000010100001001011001010010001011010100100110010-1001010001

which has a Hamming weight of 24.

We prefer base

$\{-1,0,1\}$ over base

$\{0,1\}$ because x expressed in the latter is:

1 + 2^4 + 2^5 + 2^6 + 2^7 + 2^8 + 2^11 + 2^16 + 2^19 + 2^21 + 2^22 + 2^25 + 2^27 + 2^30 + 2^34 + 2^36 + 2^37 + 2^39 + 2^41 + 2^44 + 2^47 + 2^48 + 2^51 + 2^53 + 2^54 + 2^55 + 2^58 + 2^62
[1,0,0,0,1,1,1,1,1,0,0,1,0,0,0,0,1,0,0,1,0,1,1,0,0,1,0,1,0,0,1,0,0,0,1,0,1,1,0,1,0,1,0,0,1,0,0,1,1,0,0,1,0,1,1,1,0,0,1,0,0,0,1]
100011111001000010010110010100100010110101001001100101110010001

which has a Hamming weigth of 28.

We proceed as follows:

1. We start by computing

   $m^x,(m^x{)}^x,(m^{x^2}{)}^x$.

2. Then, compute

   $m^p,m^{p^2},m^{p^3},(m^x{)}^p,(m^{x^2}{)}^p,(m^{x^3}{)}^p,(m^{x^2}{)}^{p^2}$ (use the Frobenius operator).

3. Next, we group the previous elements as follows:
   

   $$\begin{array}{rl}{y}_0& =m^p\cdot m^{p^2}\cdot m^{p^3},\\ y_1& =\overline{m},\\ y_2& =m^{x^2{p}^2},\\ y_3& =\overline{m^{xp}},\\ y_4& =\overline{m^x\cdot m^{x^{2}p}},\\ y_5& =\overline{m^{x^2}},\\ y_6& =\overline{m^{x^3}\cdot m^{x^{3}p}}.\end{array}$$
   (recall that at this point computing the conjugate is the same as computing the inverse)

4. Finally, compute the product:
   

   $$y_0\cdot y_1^2\cdot y_2^6\cdot y_3^{12}\cdot y_4^{18}\cdot y_5^{30}\cdot y_6^{36}$$
   Expand for correctness

   $$\begin{array}{rl}& y_0\cdot y_1^2\cdot y_2^6\cdot y_3^{12}\cdot y_4^{18}\cdot y_5^{30}\cdot y_6^{36}=\\ & (m^p\cdot m^{p^2}\cdot m^{p^3})\cdot {\left(\frac{1}{m}\right)}^2\cdot (m^{x^2{p}^2}{)}^6\cdot {\left(\frac{1}{m^{xp}}\right)}^{12}\cdot {\left(\frac{1}{m^x{m}^{x^{2}p}}\right)}^{18}\cdot {\left(\frac{1}{m^{x^2}}\right)}^{30}\cdot {\left(\frac{1}{m^{x^3}{m}^{x^{3}p}}\right)}^{36}\\ & =\underset{m^{{\lambda }_0}}{\underbrace{{\left(\frac{1}{m}\right)}^2\cdot {\left(\frac{1}{m^x}\right)}^{18}{\left(\frac{1}{m^{x^2}}\right)}^{30}{\left(\frac{1}{m^{x^3}}\right)}^{36}}}\cdot \underset{m^{{\lambda }_{1}p}}{\underbrace{m^p\cdot {\left(\frac{1}{m^{xp}}\right)}^{12}\cdot {\left(\frac{1}{m^{x^{2}p}}\right)}^{18}\cdot {\left(\frac{1}{m^{x^{3}p}}\right)}^{36}}}\cdot \underset{m^{{\lambda }_2{p}^2}}{\underbrace{m^{p^2}\cdot (m^{x^2{p}^2}{)}^6}}\cdot \underset{m^{{\lambda }_3{p}^3}}{\underbrace{m^{p^3}}}\\ & =m^{{\lambda }_0}\cdot m^{{\lambda }_{1}p}\cdot m^{{\lambda }_2{p}^2}\cdot m^{{\lambda }_3{p}^3}=m^{{\lambda }_0+{\lambda }_{1}p+{\lambda }_2{p}^2+{\lambda }_3{p}^3}=m^{\frac{p^4-p^2+1}{r}}\end{array}$$

   To compute this product as efficient as possible, we use the vectorial addition chain technique as explained in the last part of Section 5 of 2008/490:
   

   $$\begin{array}{rl}{T}_{0,1}& =y_6^2\cdot y_4\cdot y_5,\\ T_{1,1}& =T_{0,1}\cdot y_3\cdot y_5,\\ T_{0,2}& =T_{0,1}\cdot y_2,\\ T_{1,2}& =T_{1,1}^2\cdot T_{0,2},\\ T_{1,3}& =T_{1,2}^2,\\ T_{1,4}& =T_{1,3}\cdot y_0,\\ T_{0,3}& =T_{1,3}\cdot y_1,\\ T_{0,4}& =T_{0,3}^2\cdot T_{1,4}.\end{array}$$
   which requires only
   $9$ multiplications and
   $4$ squarings.

Costs:

$$\underset{\text{Step 1}}{\underbrace{(63·3)S+(x·3)M}}+\underset{\text{Step 2}}{\underbrace{5·(6\tilde{c}+5\tilde{m})+2·(5\tilde{m})}}+\underset{\text{Step 3}}{\underbrace{4M+5C}}+\underset{\text{Step 4}}{\underbrace{9M+4S}}$$

Frobenius Operator

As we have shown, in the final exponentation we must compute the first, second and third Frobenius operator

${\pi }_p,{\pi }_p^2,{\pi }_p^3$ over an element in

${\mathbb{F}}_{p^{12}}$. Recalling that

${\pi }_p^i(x)=x^{p^i}$, a naive implementation would take

$O(\log p)=O(254)$

${\mathbb{F}}_{p^{12}}$-operations, but we will do it much better.

Recall that we can write any element

$f\in {\mathbb{F}}_{p^{12}}$ as:

$$f=g+hw=g_0+h_{0}w+g_1{w}^2+h_1{w}^3+g_2{w}^4+h_2{w}^5$$
where

$g=g_0+g_{1}v+g_2{v}^2$ and

$h=h_0+h_{1}v+h_2{v}^2$ with

$g,h\in {\mathbb{F}}_{p^6}$ and

$g_i,h_i\in {\mathbb{F}}_{p^2}$.

Expand to see the proof

We will now show that raising an

${\mathbb{F}}_{p^2}$-element over any power of

$p$ is almost cost free and therefore so is

$f$.

• Take
  $b=b_0+b_{1}u\in {\mathbb{F}}_{p^2}$ with
  $b_0,b_1\in {\mathbb{F}}_p$. Then,
  $b^{p^{2i}}=b_0+b_1(u^{(p^2{)}^i})=b_0+b_{1}u=b$ for Lagrange theorem. Similarly,
  $b^{p^{2i-1}}=\frac{b^{p^{2i}}}{b^{p^2}\cdot b^p}=\frac{b}{b\cdot b^{-p}}=\frac{1}{b^{-p}}=b^p=\overline{b}$.
• Notice
  $w^p=(w^6{)}^{\frac{p-1}{6}}w=(9+u{)}^{\frac{p-1}{6}}w$ and therefore
  $(w^i{)}^p={\gamma }_{1,i}{w}^i$, where
  ${\gamma }_{1,i}=(9+u{)}^{i\frac{p-1}{6}}$.

Using all the previous observations, we have:

$$\begin{array}{rl}{f}^p& =(g_0+h_{0}w+g_1{w}^2+h_1{w}^3+g_2{w}^4+h_2{w}^5{)}^p\\ & ={\overline{g}}_0+{\overline{h}}_0{w}^p+{\overline{g}}_1{w}^{2p}+{\overline{h}}_1{w}^{3p}+{\overline{g}}_2{w}^{4p}+{\overline{h}}_2{w}^5\\ & ={\overline{g}}_0+{\overline{h}}_0{\gamma }_{1,1}w+{\overline{g}}_1{\gamma }_{1,2}{w}^2+{\overline{h}}_1{\gamma }_{1,3}{w}^3+{\overline{g}}_2{\gamma }_{1,4}{w}^4+{\overline{h}}_2{\gamma }_{1,5}{w}^5,\end{array}$$
which can be computed using a few multiplications and a few conjugations, assuming each

${\gamma }_{1,i}$ is pre-computed for

$i\in [5]$.

Similarly, notice

$w^{p^2}=(9+u{)}^{\frac{p^2-1}{6}}w=({\gamma }_{1,1}{)}^{1+p}w={\gamma }_{1,1}\cdot {\overline{\gamma }}_{1,1}w$ and therefore

$(w^i{)}^{p^2}={\gamma }_{2,i}{w}^i$, where

${\gamma }_{2,i}={\gamma }_{1,i}\cdot {\overline{\gamma }}_{1,i}$.

Lastly, we can use the identity

$p^3-1=p^2(p-1)+(p+1)(p-1)$ to show that

$w^{p^3}=(9+u{)}^{\frac{p^3-1}{6}}w=({\gamma }_{1,1}{)}^{p^2}\cdot {\gamma }_{2,1}w={\gamma }_{1,1}\cdot {\gamma }_{2,1}w$ and therefore

$(w^i{)}^{p^3}={\gamma }_{3,i}{w}^i$, where

${\gamma }_{3,i}={\gamma }_{1,i}\cdot {\gamma }_{2,i}$.

Then, we have:

$$\begin{array}{r}{f}^p={\overline{g}}_0+{\overline{h}}_0{\gamma }_{1,1}w+{\overline{g}}_1{\gamma }_{1,2}{w}^2+{\overline{h}}_1{\gamma }_{1,3}{w}^3+{\overline{g}}_2{\gamma }_{1,4}{w}^4+{\overline{h}}_2{\gamma }_{1,5}{w}^5,\\ f^{p^2}=g_0+h_0{\gamma }_{2,1}w+g_1{\gamma }_{2,2}{w}^2+h_1{\gamma }_{2,3}{w}^3+g_2{\gamma }_{2,4}{w}^4+h_2{\gamma }_{2,5}{w}^5,\\ f^{p^3}={\overline{g}}_0+{\overline{h}}_0{\gamma }_{3,1}w+{\overline{g}}_1{\gamma }_{3,2}{w}^2+{\overline{h}}_1{\gamma }_{3,3}{w}^3+{\overline{g}}_2{\gamma }_{3,4}{w}^4+{\overline{h}}_2{\gamma }_{3,5}{w}^5,\end{array}$$
with

${\gamma }_{1,i}=(9+u{)}^{i\frac{p-1}{6}}$,

${\gamma }_{2,i}=(9+u{)}^{i\frac{p^2-1}{6}}={\gamma }_{1,i}\cdot {\overline{\gamma }}_{1,i}$ and

${\gamma }_{3,i}=(9+u{)}^{i\frac{p^3-1}{6}}={\gamma }_{1,i}\cdot {\gamma }_{2,i}$. Hence, assuming the precomputation of

${\gamma }_{1,i},{\gamma }_{2,i},{\gamma }_{3,i}$ for

$i\in [5]$, we can compute

$f^p,f^{p^2}$ or

$f^{p^3}$ almost for free.

Costs:

$\underset{f^p}{\underbrace{6\tilde{c}+5\tilde{m}}},\underset{f^{p^2}}{\underbrace{5\tilde{m}}},\underset{f^{p^3}}{\underbrace{6\tilde{c}+5\tilde{m}}}$

Speeding up Computations in the Cyclotomic Subgroup

I follow 2010/542 closely for this section.

Recall that:

$$\begin{array}{rl}{\mathbb{F}}_{p^2}& ={\mathbb{F}}_p[u]/(u^2+1),\\ {\mathbb{F}}_{p^6}& ={\mathbb{F}}_{p^2}[v]/(v^3-(9+u)),\\ {\mathbb{F}}_{p^{12}}& ={\mathbb{F}}_{p^6}[w]/(w^2-v).\end{array}$$
and also:

$${\mathbb{F}}_{p^{12}}={\mathbb{F}}_{p^2}[w]/(w^6-(9+u))$$

For this optimitzation, we will further need the following

${\mathbb{F}}_{p^{12}}$ representation:

$$\begin{array}{rl}{\mathbb{F}}_{p^4}& ={\mathbb{F}}_{p^2}[w^3]/((w^3{)}^2-(9+u)),\\ {\mathbb{F}}_{p^{12}}& ={\mathbb{F}}_{p^4}[w]/(w^3-w^3),\end{array}$$
so that for

$f=(g_0+g_1{w}^3)+(g_2+g_3{w}^3)w+(g_4+g_5{w}^3)w^2\in C$, we can equivalently write

$f=g_0+g_{2}w+g_4{w}^2+g_1{w}^3+g_3{w}^4+g_5{w}^5\in B$.

This optimization is due to "compression" and "decompression" technique, which exploits the algebraic relationships arising from elements belonging to the cyclotomic subgroup

${\mathbb{G}}_{{\varphi }_6(p^2)}$.

Compression and Decompression

If

$f=g_0+g_{2}w+g_4{w}^2+g_1{w}^3+g_3{w}^4+g_5{w}^5$ then the compression function

$\mathcal{C}$ is:

$$\mathcal{C}(f)=[g_2,g_3,g_4,g_5].$$
and the decompression function

$\mathcal{D}$ is:

$$\mathcal{D}([{\tilde{g}}_2,{\tilde{g}}_3,{\tilde{g}}_4,{\tilde{g}}_5])={\tilde{g}}_0+{\tilde{g}}_{2}w+{\tilde{g}}_4{w}^2+{\tilde{g}}_1{w}^3+{\tilde{g}}_3{w}^4+{\tilde{g}}_5{w}^5,$$
where:

$$\{\begin{array}{ll}{\tilde{g}}_1=\frac{{\tilde{g}}_5^2(9+u)+3{\tilde{g}}_4^2-2{\tilde{g}}_3}{4{\tilde{g}}_2},{\tilde{g}}_0=(2{\tilde{g}}_1^2+{\tilde{g}}_2{\tilde{g}}_5-3{\tilde{g}}_3{\tilde{g}}_4)(9+u)+1,& \text{if }{\tilde{g}}_2\ne 0\\ {\tilde{g}}_1=\frac{2{\tilde{g}}_4{\tilde{g}}_5}{{\tilde{g}}_3},{\tilde{g}}_0=(2{\tilde{g}}_1^2-3{\tilde{g}}_3{\tilde{g}}_4)(9+u)+1,& \text{if }{\tilde{g}}_2=0\end{array}$$

Compressed Squaring

Now, on input

$f\in {\mathbb{G}}_{{\varphi }_6(p^2)}$ we will be able to compute its square on compressed form. That is,

$\mathcal{C}(f^2)=[h_2,h_3,h_4,h_5]$, where

$h=f^2=h_0+h_{2}w+h_4{w}^2+h_1{w}^3+h_3{w}^4+h_5{w}^5$:

$$\begin{array}{rl}{h}_2& =2(g_2+3(9+u)B_{4,5}),\\ h_3& =3(A_{4,5}-(10+u)B_{4,5})-2g_3,\\ h_4& =3(A_{2,3}-(10+u)B_{2,3})-2g_4,\\ h_5& =2(g_5+3B_{2,3}),\\ A_{i,j}& =(g_i+g_j)(g_i+(9+u)g_j),\\ B_{i,j}& =g_i\cdot g_j.\end{array}$$

Exponentiation in the Cyclotomic Subgroup

Given

$f\in {\mathbb{G}}_{{\varphi }_6(p^2)}\subset {\mathbb{F}}_{p^{12}}^{\ast }$ and

$e\ge 1$, the idea is to compute

$f^e$ using the squaring formula in compressed form.

Represent

$e$ in binary as

$e=e_{\ell -1}{e}_{\ell -2}\dots e_1{e}_0$ with

$e_{\ell -1}=1$ and define

$H_e=\{i:1\le i\le \ell -1\text{ and }{e}_i=1\}$. Then:

$$f^e=\prod _{i=0}^{\ell -1}{f}^{e_i\cdot 2^i}=f^{e_0}\prod _{i\in H_e}\mathcal{D}(\mathcal{C}(f{)}^{2^i}).$$

Therefore, we can use the square-and-multiply algorithm to compute such product. Since we are going to be squaring in compressed form, we start at

$\mathcal{C}(f)$ rather than the standard

$f$ in such algorithm.

def exp_cyclo(f,e):
    Cf = comp(f);
    result = 1;
    for i in range(len(e)):
        if (e[i] == 1) {
            result = decomp(Cf) * result;
        }

        Cf = square_comp(Cf);
    }

    return result;

Appendix A: Field Isomorphisms

Recall that:

$$\begin{array}{rl}{\mathbb{F}}_{p^2}& ={\mathbb{F}}_p[u]/(u^2+1),\\ {\mathbb{F}}_{p^6}& ={\mathbb{F}}_{p^2}[v]/(v^3-(9+u)),\\ A={\mathbb{F}}_{p^{12}}& ={\mathbb{F}}_{p^6}[w]/(w^2-v).\end{array}$$

We use two additional representations of

${\mathbb{F}}_{p^{12}}$.

The first one is used to efficiently compute the Frobenius operator and is as follows:

$$B={\mathbb{F}}_{p^{12}}={\mathbb{F}}_{p^2}[w]/(w^6-(9+u)).$$

The second one is used to speed up the squares and exponentations in the hard part of the final exponentiation:

$$\begin{array}{rl}{\mathbb{F}}_{p^4}& ={\mathbb{F}}_{p^2}[w^3]/((w^3{)}^2-(9+u)),\\ C={\mathbb{F}}_{p^{12}}& ={\mathbb{F}}_{p^4}[w]/(w^3-w^3).\end{array}$$

Isomorphism between A and B

We start with

$A={\mathbb{F}}_{p^6}[w]/(w^2-v)$ and

$B={\mathbb{F}}_{p^2}[w]/(w^6-(9+u))$, where we notice that if we take

$a=a_1+a_2\cdot w\in A$, we can use the identity

$v=w^2$ to rewrite it as:

$$\begin{array}{rl}a& =(a_{1,1}+a_{1,2}v+a_{1,3}{v}^2)+(a_{2,1}+a_{2,2}v+a_{2,3}{v}^2)\cdot w=\\ & =a_{1,1}+a_{2,1}\cdot w+a_{1,2}\cdot w^2+a_{2,2}\cdot w^3+a_{1,3}\cdot w^4+a_{2,3}\cdot w^5\in B\end{array}$$

This relationship induces the following field isomorphism

${\chi }_{AB}:A\to B$ and its inverse

${\chi }_{AB}^{-1}:B\to A$:

$$\begin{array}{rl}{\chi }_{AB}(a_{1,1},a_{1,2},a_{1,3},a_{2,1},a_{2,2},a_{2,3})& =(a_{1,1},a_{2,1},a_{1,2},a_{2,2},a_{1,3},a_{2,3}),\\ {\chi }_{AB}^{-1}(b_1,b_2,b_3,b_4,b_5,b_6)& =(b_1,b_3,b_5,b_2,b_4,b_6).\end{array}$$

Exercise: Check that

${\chi }_{AB}$ is a field isomorphism between

$A$ and

$B$ and that for all

$a\in A$ and

$b\in B$ it is true that

${\chi }_{AB}^{-1}({\chi }_{AB}(a))=a$ and

${\chi }_{AB}({\chi }_{AB}^{-1}(b))=b$.

Isomorphism between B and C

We continue with

$B={\mathbb{F}}_{p^2}[w]/(w^6-(9+u))$ and

$C={\mathbb{F}}_{p^4}[w]/(w^3-w^3)$, for which we notice that

$c=(c_1+c_2\cdot w^3)+(c_3+c_4\cdot w^3)w+(c_5+c_6\cdot w^3)\cdot w^2\in C$ can be equivalently written as

$c_1+c_3\cdot w+c_5\cdot w^2+c_2\cdot w^3+c_4\cdot w^4+c_6\cdot w^5\in B$.

This relationship induces the following field isomorphism

${\chi }_{BC}:B\to C$ and its inverse

${\chi }_{BC}^{-1}:C\to B$:

$$\begin{array}{rl}{\chi }_{BC}(b_1,b_2,b_3,b_4,b_5,b_6)& =(b_1,b_4,b_2,b_5,b_3,b_6),\\ {\chi }_{BC}^{-1}(c_1,c_2,c_3,c_4,c_5,c_6)& =(c_1,c_3,c_5,c_2,c_4,c_6).\end{array}$$

Exercise: Check that

${\chi }_{BC}$ is a field isomorphism between

$B$ and

$C$ and that for all

$b\in B$ and

$c\in C$ it is true that

${\chi }_{BC}^{-1}({\chi }_{BC}(b))=b$ and

${\chi }_{BC}({\chi }_{BC}^{-1}(c))=c$.

Appendix B: Finite Field Arithmetic

Fp2 Arithmetic

Addition

# INPUT: (a = a1 + a2·u), (b = b1 + b2·u) ∈ Fp2, where ai,bi ∈ Fp
# OUTPUT: (c1 + c2·u) = a + b ∈ Fp
c1 = a1 + b1
c2 = a2 + b2

Subtraction

# INPUT: (a = a1 + a2·u), (b = b1 + b2·u) ∈ Fp2, where ai,bi ∈ Fp
# OUTPUT: (c1 + c2·u) = a - b ∈ Fp
c1 = a1 - b1
c2 = a2 - b2

Multiplication

# INPUT: (a = a1 + a2·u), (b = b1 + b2·u) ∈ Fp2, where ai,bi ∈ Fp
# OUTPUT: (c1 + c2·u) = a·b ∈ Fp
c1 = a1·b1 - a2·b2
c2 = a1·b2 + a2·b1

Scalar Multiplication

# INPUT: b ∈ Fp, (a = a1 + a2·u) ∈ Fp2, where ai ∈ Fp
# OUTPUT: (c1 + c2·u) = a·b ∈ Fp
c1 = b·a1
c2 = b·a2

Square

# INPUT: (a = a1 + a2·u) ∈ Fp2, where ai ∈ Fp
# OUTPUT: (c1 + c2·u) = a² ∈ Fp
c1 = (a1 - a2)·(a1 + a2)
c2 = 2·a1·a2

Inverse

# INPUT: (a = a1 + a2·u) ∈ Fp2, where ai ∈ Fp
# OUTPUT: (c1 + c2·u) = a⁻¹ ∈ Fp
c1 = a1·(a1² + a2²)⁻¹
c2 = -a2·(a1² + a2²)⁻¹

Fp4 Arithmetic

Square

# INPUT: (a = a1 + a2·V) ∈ Fp4, where ai ∈ Fp2
# OUTPUT: (c1 + c2·V) = a² ∈ Fp4
c1 = a2²·(9 + u) + a1²
c2 = (a1 + a2)² - a1² - a2²

Fp6 Arithmetic

Addition

# INPUT: (a = a1 + a2·v + a3·v²),(b = b1 + b2·v + b3·v²) ∈ Fp6, where ai,bi ∈ Fp2
# OUTPUT: (c1 + c2·v + c3·v²) = a + b ∈ Fp6
c1 = a1 + b1
c2 = a2 + b2
c3 = a3 + b3

Subtraction

# INPUT: (a = a1 + a2·v + a3·v²),(b = b1 + b2·v + b3·v²) ∈ Fp6, where ai,bi ∈ Fp2
# OUTPUT: (c1 + c2·v + c3·v²) = a - b ∈ Fp6
c1 = a1 - b1
c2 = a2 - b2
c3 = a3 - b3

Multiplication

# INPUT: (a = a1 + a2·v + a3·v²),(b = b1 + b2·v + b3·v²) ∈ Fp6, where ai,bi ∈ Fp2
# OUTPUT: (c1 + c2·v + c3·v²) = a·b ∈ Fp6
c1 = [(a2+a3)·(b2+b3) - a2·b2 - a3·b3]·(9+u) + a1·b1
c2 = (a1+a2)·(b1+b2) - a1·b1 - a2·b2 + (9+u)·(a3·b3)
c3 = (a1+a3)·(b1+b3) - a1·b1 + a2·b2 - a3·b3

Scalar Multiplication

# INPUT: b ∈ Fp, (a = a1 + a2·v + a3·v²) ∈ Fp6, where ai ∈ Fp2
# OUTPUT: (c1 + c2·v + c3·v²) = a·b ∈ Fp6
c1 = b·a1
c1 = b·a2
c1 = b·a3

The following three algorithms help to compute the line equations. See the sparse multiplications algorithms of next section for a full understanding.

Sparse Multiplication A

# INPUT: (a = a1 + a2·v + a3·v²),b = b2·v ∈ Fp6, where ai,b2 ∈ Fp2
# OUTPUT: (c1 + c2·v + c3·v²) = a·b ∈ Fp6
c1 = b2·a3·(9+u)
c2 = b2·a1
c3 = b2·a2

Sparse Multiplication B

# INPUT: (a = a1 + a2·v + a3·v²),(b = b2·v + b3·v²) ∈ Fp6, where ai,bi ∈ Fp2
# OUTPUT: (c1 + c2·v + c3·v²) = a·b ∈ Fp6
c1 = [(a2 + a3)·(b2 + b3) - a2·b2 - a3·b3]·(9+u)
c2 = a1·b2 + a3·b3·(9+u)
c3 = (a1 + a3)·b3 - a3·b3 + a2·b2

Sparse Multiplication C

# INPUT: (a = a1 + a2·v + a3·v²),(b = b1 + b3·v²) ∈ Fp6, where ai,bi ∈ Fp2
# OUTPUT: (c1 + c2·v + c3·v²) = a·b ∈ Fp6
c1 = a1·b1 + a2·b3·(9+u)
c2 = a2·b1 + a3·b3·(9+u)
c3 = (a1 + a3)·(b1 + b3) - a1·b1 - a3·b3

Square

# INPUT: (a = a1 + a2·v + a3·v²) ∈ Fp6, where ai ∈ Fp2
# OUTPUT: (c1 + c2·v + c3·v²) = a² ∈ Fp6
c1 = 2·a2·a3·(9 + u) + a1²
c2 = a3²·(9 + u) + 2·a1·a2
c3 = 2·a1·a2 - a3² + (a1 - a2 + a3)² + 2·a2·a3 - a1²

Inverse

# INPUT: (a = a1 + a2·v + a3·v²) ∈ Fp6, where ai ∈ Fp2
# OUTPUT: (c1 + c2·v + c3·v²) = a⁻¹ ∈ Fp6
c1mid = a1² - (9 + u)·(a2·a3)
c2mid = (9 + u)·a3² - (a1·a2)
c3mid = a2² - (a1·a3)
c1 = (a1² - (9 + u)·(a2·a3))·(a1·c1mid + xi·(a3·c2mid + a2·c3mid))⁻¹
c2 = ((9 + u)·a3² - (a1·a2))·(a1·c1mid + xi·(a3·c2mid + a2·c3mid))⁻¹
c3 = (a2²-a1·a3)·(a1·c1mid + xi·(a3·c2mid + a2·c3mid))⁻¹

Fp12 Arithmetic

Addition

# INPUT: (a = a1 + a2·w),(b = b1 + b2·w) ∈ Fp12, where ai,bi ∈ Fp6
# OUTPUT: (c1 + c2·w) = a + b ∈ Fp12
c1 = a1 + b1
c2 = a2 + b2

Subtraction

# INPUT: (a = a1 + a2·w),(b = b1 + b2·w) ∈ Fp12, where ai,bi ∈ Fp6
# OUTPUT: (c1 + c2·w) = a - b ∈ Fp12
c1 = a1 - b1
c2 = a2 - b2

Multiplication

# INPUT: (a = a1 + a2·w),(b = b1 + b2·w) ∈ Fp12, where ai,bi ∈ Fp6
# OUTPUT: (c1 + c2·w) = a·b ∈ Fp12
c1 = a1·b1 + a2·b2·v
c2 = (a1+a2)·(b1+b2) - a1·b1 - a2·b2

The following two algorithms are used to compute the line equations avoiding all multiplications by zero that arise from the sparisity of such equations.

Sparse Multiplication A

# INPUT: (a = a1 + a2·w),(b = b1 + b2·w) ∈ Fp12, where ai ∈ Fp6, b1 = b12·v and b2 = b22·v + b23·v², with b12,b22,b23 ∈ Fp2
# OUTPUT: (c1 + c2·w) = a·b ∈ Fp12
c1 = a1·b1 + a2·b2·v
c2 = (a1+a2)·[(b12+b22)·v + b23·v²] - a1·b1 - a2·b2

Sparse Multiplication B

# INPUT: (a = a1 + a2·w),(b = b1 + b2·w) ∈ Fp12, where ai ∈ Fp6, b1 = b11 + b13·v² and b2 = b22·v, with b11,b13,b22 ∈ Fp2
# OUTPUT: (c1 + c2·w) = a·b ∈ Fp12
c1 = a1·b1 + a2·b2·v
c2 = (a1+a2)·(b11 + b22·v + b13·v²) - a1·b1 - a2·b2

Square

# INPUT: (a = a1 + a2·w) ∈ Fp12, where ai ∈ Fp6
# OUTPUT: (c1 + c2·w) = a² ∈ Fp12
c1 = (a1-a2)·(a1-a2·v) + a1·a2 + a1·a2·v
c2 = 2·a1·a2

Inverse

# INPUT: (a = a1 + a2·w) ∈ Fp12, where ai ∈ Fp6
# OUTPUT: (c1 + c2·w) = a⁻¹ ∈ Fp12
c1 = a1·(a1² - a2²·v)⁻¹
c2 = -a2·(a1² - a2²·v)⁻¹

First Frobenius Operator

# INPUT: (a = a1 + a2·w) = ((a11 + a12v + a13v²) + (a21 + a22v + a23v²)) ∈ Fp12, where ai ∈ Fp6 and aij ∈ Fp2
# OUTPUT: (c1 + c2·w) = a^p ∈ Fp12
c1 = a̅11     + a̅12·γ12·v + a̅13·γ14·v²
c2 = a̅21·γ11 + a̅22·γ13·v + a̅23·γ15·v²

Second Frobenius Operator

# INPUT: (a = a1 + a2·w) = ((a11 + a12v + a13v²) + (a21 + a22v + a23v²)) ∈ Fp12, where ai ∈ Fp6 and aij ∈ Fp2
# OUTPUT: (c1 + c2·w) = a^p² ∈ Fp12
c1 = a11     + a12·γ22·v + a13·γ24·v²
c2 = a21·γ21 + a22·γ23·v + a23·γ25·v²

Third Frobenius Operator

# INPUT: (a = a1 + a2·w) = ((a11 + a12v + a13v²) + (a21 + a22v + a23v²)) ∈ Fp12, where ai ∈ Fp6 and aij ∈ Fp2
# OUTPUT: (c1 + c2·w) = a^p³ ∈ Fp12
c1 = a̅11     + a̅12·γ32·v + a̅13·γ34·v²
c2 = a̅21·γ31 + a̅22·γ33·v + a̅23·γ35·v²

The following two algorithms are used for computing squarings and exponentiations over the cyclotomic subgroup

${\mathbb{G}}_{{\varphi }_6(p^2)}\subset {\mathbb{F}}_{p^{12}}$. Luckily for us, this is the case after performing the easy part of the final exponentiation.

Fast Squaring

# INPUT: (a = a1 + a2·w) ∈ GΦ6(p²), where ai ∈ Fp6
# OUTPUT: ((c11 + c12·v + c13·v²) + (c21 + c22·v + c23·v²)·w) = a² ∈ GΦ6(p²)
[t11,t22] = (a11 + a22·V)²
[t23,t12] = (a21 + a13·V)²
[t13,aux] = (a12 + a23·V)²
t21 = aux·(9+u)
c11 = -2·a11 + 3·t11
c12 = -2·a12 + 3·t23
c13 = -2·a13 + 3·t13
c21 = 2·a21 + 3·t21
c22 = 2·a22 + 3·t22
c23 = 2·a23 + 3·t12

Fast Exponentatiation

# INPUT: e, (a = a1 + a2·w) ∈ GΦ6(p²), where e ∈ [0,p¹²-2] ai ∈ Fp6
# OUTPUT: (c = c1 + c2·w) = a^e ∈ GΦ6(p²)
c = a
for i in range(L-2,-1,-1): # L is log2(e)
    c = c² # Use previous algorithm
    if e[i] == 1:
        c = c·a

return c

---

1. if you don't know why, make sure to review Galois theory. Start by noticing that Galois group

   $\text{Gal}({\mathbb{F}}_{p^{12}}/{\mathbb{F}}_{p^6})$ is a cyclic group of order
   $2$ generated by the Frobenius endomorphism
   $F(X)=X^{p^6}$. ↩︎

2. I recently discoverd that this is just a particular case of the more generic non-adjacent form (NAF) of a number. ↩︎