# BLS12-381 理论与实现

• 白菜
• 更新于 2024-05-07 16:34
• 阅读 1252

Pairings, KZG, SNARK

<br />

<br />

## 这里没有的

• group theory, field theory and homomorphism

相关基本概念在这里不会涵盖，详情请查阅任何abstract algebraic 相关的书籍

• divisors

相关基本概念在这里不会涵盖，对于了解Pairing 来说 Pairing for Beginners 已足够，如果你还想深入理解最好翻阅一下 algebraic geometry 相关的书籍

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

理论和算法部分这里不会涵盖，详情可以查阅 Guide to Elliptic Curve Cryptography

• hash to curve

bytes string 映射到$\mathbb{G}_1$或者$\mathbb{G}_2$ 上的点，简单说就是hash，是pairing 应用层面必备的一大模块，后续会详细补充这块内容

• non-affine coordinate

affine coordinate 其实只是椭圆曲线元素表达的需要，它的scalar multiplication 并不经济，所以实际计算上都会用non-affine coordinate 来替代，后续会补上这块内容

• advanced scalar multiplication algorithms GLV/GLS

特定的曲线上充分利用同态映射来加速scalar multiplication，同时还能(GPU)并行化处理也是当下硬件加速卖点，后续也会再补上

<br />

<br />

## 关于代码

• python implementation

主要集中在Pairing的计算逻辑上，包括Miller LoopFinal Exponentiation。目前已经完成验证。

Finite FieldElliptic Curves的算术运算并没有逐一实现，用的是Sagemath库自带的 Galois Field and Elliptic Curve.

• rust implementation

从零着手，从 Bigint 算术运算到 Finite Field 算术运算到 Elliptic Curve 算术运算，再到 Pairings Primitives。底层的逻辑已经验证完毕，目前在Pairings验证过程中 ...

<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$$ for: \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 />

## Pairing 的演进

### Weil Reciprocity

$g$ and $f$ 是两个定义在椭圆曲线上的divisor function, $f, g \in K(E)$，它们的divisor support 不存在交集, $supp((f)) \land supp((g)) = \emptyset$。然后我们就有:

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

<br />

Details of general definition of Weil Reciprocity, you can refer THEOREM 3.9 of Guide to Pairing-based Cryptography.

<br />

<br />

### Weil Pairing

<br />

<br />

<br />

<br />

<br />

<br />

<br />

#### 如何对divisor $D_P$ and $D_Q$ 进行evaluate

divisor $D_P = (P) - (\mathcal{O})$ 的evaluation 可以被进行一步简化:

$$f_{rD_Q}(DP) \equiv f{rD_Q}(P)$$

<br />

<br />

Miller Loop 使得divisor function的evaluation $f_{r((Q) - (\mathcal{O}))}(P)$ 变得更容易实现，是工程上的一大步。很明显 Weil Pairing 是几何上对称的, 它实际上需要运行两次 Miller Loop. 看起来并不太经济? 实际上单次就够了，这就是 Tate Pairing 要做的事情.

<br />

<br />

### Tate Pairing

<br />

#### 定义

<br />

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

<br />

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 />

<br />

<br />

### Ate Pairing

Ate Pairing中, 点 $Q$ 被严格约束在 $\mathbb{G}_2 = \pi[p]$中，同时点$P$ 也被约束在 $\mathbb{G}_1 = \pi[1]$中，即Frobenius Map： $$\pi(Q) = [p]Q, \pi(P) = [1]P$$ 充分利用Frobenius Map的特性，将大大降低pairing的计算成本。

#### Miller 算法的两个重要特性

\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} 关于这两个特性的proof，这里不再推演，熟悉divisor function 后很容易推导出来。

<br />

<br />

<br />

<br />

<br />

#### Trace of Frobenius

<br />

<span style="color:red">在Ate Pairing中double-add 作用在curve $\mathbb{G}_2$ 上的点$Q$ (line evaluation的点是$P$)，并没有跟Tate Pairing一样作用在$\mathbb{G}_1$ (定义在base prime field$F_p$上) 上的点$P$ (line evaluation的点是$Q$)，为什么？因为Ate Pairing 之所以能用缩短Miller Loop，根本原因是运用了Frobenius Map $\varphi$ 的特性，而$\mathbb{G}_1$ 上运用Frobenius Map 是trivial的，没有意义，$\varphi(P) = [1]P$。</span>

<br />

<span style="color:red">那么，计算上它会带来什么后果吗？ 有的，如果$Q$定义在$F_{p^{12}}$上，$Q \in \mathbb{G}2 = E(F{p^{12}})[r]$，那么它的Double-add 成本将会很大，所以这时候twist 就派上用场了，把$Q$映射 twist curve $\mathbb{G}2'$ (定义在域$F{p^2}$)上去，这样Double-add的成本就会最大程度地降低，但仍然要比Tate Pairing 中的Double-add成本要稍微大些（$F_{p^2}$ VS $F_p$），所幸的是Miller Loop 可以被大大缩短，这会在一定程度上覆盖掉Double-add 带来的损失。</span>

<br />

<span style="color:red">所以Ate Pairing最终的性能如何，还是要取决于$T$ 的bit length 与 $r$ 的bit length 之间的差距。但是大家通常所见到的都是Ate Pairing，为什么Tate Pairing没有出现过？原因其实很简单，通过polynomial function $r(x) = x^4 - x^2 + 1, t(x) = x + 1$ 很容易看出Ate Pairing 的Miller Loop 要大大小于Tate Pairing的Miller Loop，$\log(r) \approx 4 \log(x) \approx 4\log(T)$。因此上面提到的Double-add 带来的影响就基本上可以被忽略了。 </span>

<br />

p = 103
r = 7
k = 6

for i in range(1, k):
print('lambda[{0}] = {1}'.format(i, (p ** i) % r))

lambda[1] = 5
lambda[2] = 4
lambda[3] = 6
lambda[4] = 2
lambda[5] = 3

<br />

### Optimal Ate Pairing

<br />

<br />

Ate Pairing 中，我们有: $$f{\lambda, Q}^{x1} = f{r, Q}^m$$ 其中 $m = (\lambda^k - 1) // r$, and $x1 = \sum_{i = 0}^{k - 1} (k - 1 - i) \cdot p^i$.

<br />

<br />

<br />

Optimal Ate Pairing中把$\lambda$进行了更通用的定义：$\lambda = \sum_{i = 0}^l c_i \cdot p^i, l < k$, and $\lambda = m \cdot r$. 同样运用上面的 Miller 算法的特性, 我们把它的divisor function 进行展开:

\begin{aligned} f{\lambda, Q} &= f{\sum_{i = 0}^l ci \cdot p^i, Q} \ &= \prod{i = 0}^l f_{ci \cdot p^i, Q} \cdot \prod{j = 0}^{l - 1} \frac{l{[s{j + 1}]Q,[cj \cdot q^j]Q}}{v{[sj]Q}} \ &= \prod{i = 0}^l f_{p^i, Q}^{ci} \cdot (\prod{i = 0}^l f_{ci, Q}^{p^i} \cdot \prod{j = 0}^{l - 1} \frac{l{[s{j + 1}]Q,[cj \cdot q^j]Q}}{v{[sj]Q}}) \ \end{aligned} 其中: $$\prod{i = 0}^l f_{p^i, Q}^{ci} = \prod{i = 0}^l f_{p, Q}^{i \cdot ci \cdot p^{i - 1}} = f{p, Q}^{\sum_{i = 0}^l i \cdot c_i \cdot p^{i - 1}}$$

<br />

<br />

<br />

<br />

## 有限域上的算术运算

BLS12-381 曲线的定义是这样的: $$E(F{p^{12}}): y^2 = x^3 + 4$$ 其中 $F{p^{12}} = F_p[v] / X^{12} - 2 X^6 + 2$。但是这个extension field 是如何构建的呢？

<br />

### Pairing 中域的切换

<br />

• Miller Loop

line function $f_{r, P}$ 不会改变所在的域，$P$在哪个域，这个函数仍然在那个域，比如 $F_p$:

• Evaluation Line Function

比如，单个line function的evaluation： $$l_{r, P}(Q) = y_Q - y_T - \alpha \cdot (x_Q - x_T)$$

最终evaluation 的结果， $f{r, P}(Q)$ 不再定义在$Q$原本的域 $[r]F{p^{12}}^{\times}$ 上了.

• Final Exponentiation

• Easy Part

通过提指Mill Loop 的结果$f{r, P}(Q)$ 推进一个特殊的乘法group，这就是我们所说的 Cyclotomic Group， $F{\varPhi_{12}}$:

• Hard Part

再次通过提指Cyclotomic Group 拉到目标乘法group $F_{p^{12}}^{\times}[r]$:

<br />

### 域塔 Tower Fields

#### 定义

<br />

$$Fp[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 />

$F_{p^{12}}$模的由来

• $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 />

<br />

<br />

<br />

<br />

### Cyclotomic Group 上的算术运算

<br />

<br />

<br />

<br />

<br />

<br />

<br />

<br />

<br />

<br />

#### 更快的 Squaring 算子

Squaring 之后: \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 />

Tate/Ate PairingFinal Exponentiation中 $\alpha \in \mathbb{G}_{\varPhi_6(q)}$，根据上面刚刚推演出的结论: \begin{aligned} \alpha^{\varPhi_6(q)} = \alpha^{q^2 - q + 1} = 1 \ \Longrightarrow \alpha^{q^2} \cdot \alpha = \alpha^q \ \end{aligned}

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

根据Frobenius Map 的特性，我们很容易得到：$a^q = \bar{a}, b^q = \bar{b}, c^q = \bar{c}$ (先给出结论，在后面Frobenius Map 部分会进行推理)。因此上面的式子简化成: \begin{aligned} \alpha^q = \bar{a} + \bar{b} \cdot v^q + \bar{c} \cdot (v^q)^2 \end{aligned}

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

由于 $v^3 = u, u^2 = \xi, u \in F_{q^2}, \xi \in F_q$，因此： $$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$$ 其中 $\zeta_6$ 是域$F_q$上的 primitive 6-th root of unity，也就是说$\xi_6^6 \equiv 1 \mod q$ .

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

<br />

<br />

<br />

<br />

### Twist 的力量

<br />

#### Sextic Twist

<br />

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

<br />

<br />

<br />

#### Sextic Twist Map

• Twist Operation

把 $E_{12}$ 上的元素映射到 $E_2'$上: $$\varphi: (x, y) \longmapsto (x \cdot \sqrt[3]{\xi}, y \cdot \sqrt[2]{\xi})$$

• Untwist Operation

把$E2'$ 上的元素映射到 $E{12}$ 上: $$\varphi^{-1}: (x', y') \longmapsto (\frac{x'}{\sqrt[3]{\xi}}, \frac{y'}{\sqrt[2]{\xi}})$$

<br />

twist/untwist的过程是很便宜的，尤其是当我们把期间用到的常量 $\sqrt[2]{\xi}, \sqrt[3]{\xi}$ 预先算出来。最后明确一下，既然要用twist trick，那么将尽可能把运算都限制在低阶的域上，只是在必要的时候才通过untwist 把值转换到高阶域上。

<br />

### Frobenius Map 的力量

<br />

<br />

<br />

<br />

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

<br />

<br />

• 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 />

<br />

<br />

<br />

<br />

<br />

## Curve 上的算术运算

<br />

### $\mathbb{G}_1$ 上的算术运算

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

<br />

$\mathbb{G}_1$ 定义在 Base Prime Field $F_p$上，它只是 $E(F_p)$ 的一个$r$-torsion subgroup. 所以它的算术运算(Scalar Multiplication) 跟$E(F_p)$一样，定义在$F_p$。

<br />

### $\mathbb{G}_2$上的算术运算

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

<br />

<br />

## Python Implementation

### Instantiation of Curve BLS12-381

#### Trace/Unti-trace Map

Trace Map: $$Tr(Q) = \sum_{i = 0}^k (Q_x^{p^i}, Qy^{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)
## 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 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(Fp)$, 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, 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 &= \lambda1 \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))

## map P from curve E(Fp) (or E(Fp12)) into twisted curve E_t(Fp2)
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 maps fp2 into fp12 before line function evaluation
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 to be updated...

<br />

## References

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

<br />