calculate weil-pairing from 0 to 1

为了更好理解 weil-pairing 以及 miller 算法，这里使用 python 从 0 开始实现了一遍。

用到的依赖只有 pow 函数，math.inf 和深拷贝函数 copy.deepcopy，且使用了本人并不顺手的 OOP 编程范式。

虽然写了之后对 weil-pairing 和 OOP 的理解同时加深了（

从 ${\mathbb{F}}_p$ 开始

我们在高中竞赛时学过 ${\mathbb{F}}_p$ 的运算。它的加减乘运算都很简单，只需在普通运算（这里指代码中的加减乘运算）的基础上 $\mathrm{mod}p$ 即可。而对于除法，我们学过至少两种方法（费马小定理和 exgcd）求一个元素在 ${\mathbb{F}}_p$ 上的逆元，从而可以化除为乘。

然而今天没有必要再讲一遍 exgcd，因为先前有一篇 post 已经详细介绍了其原理，这里我选择直接使用 python 的语法糖：

pow(x, -1, p)

来求 $x^{-1}\mathrm{mod}p$.

于是对于分数 $\frac{a}{b}$，我们也可以用 $a\times b^{-1}\mathrm{mod}p$ 来表达，${\mathbb{F}}_p$ 上的运算我们就处理完了。

引入椭圆曲线

这里引出本 post 的主角，一个定义在 ${\mathbb{F}}_{631}$ 上的椭圆曲线 $E$：

$$y^2=x^3+30x+34$$

这里我们使用一个类 EpCurve，来完成对它的抽象。目前，这个类已经可以放入求逆元和除法两个方法，为了体现这个曲线的特性，我们可以加一个判断点 $(x,y)$ 是否在这个曲线的方法，于是我们可以写出完整的实现：

class EpCurve:
    def __init__(self, a, b, p) -> None:
        self.A = a
        self.B = b
        self.mod = p
    
    def inverse(self, x:int) -> int:
        return pow(x, -1, self.mod) 

    def frac(self, x:int, y:int) -> int:
        return x * self.inverse(y) % self.mod
    
    def isOnCurve(self, x, y) -> bool:
        return y*y % self.mod == (x*x*x + self.A*x + self.B) % self.mod

以及对应的实例 curve = EpCurve(30, 34, 631)。

然后我们可以在 EpCurve 的基础上，抽象椭圆曲线上的点 EpPoint：

class EpPoint(EpCurve):
    def __init__(self, x, y, curve) -> None: 
        super().__init__(curve.A, curve.B, curve.mod)
        if not self.isOnCurve(x, y):
            raise ValueError(f"Point ({x}, {y}) is not on the curve")
        self.x = x % self.mod
        self.y = y % self.mod
        self.O = False
        self.curve = curve

    def __str__(self) -> str:
        if self.O:
            return 'O'
        return f"({self.x}, {self.y})"

这里 self.x self.y 就是点对应的坐标 $(x,y)$，而 self.O 对应了它是不是无穷远点 $\mathcal{O}$。__str__ 方法就是 python 内置的一个魔术方法，用于在 print 这个实例时调用这个函数输出对应的字符串。

点的运算

在大二下期的应用密码学课上，我们系统学习了椭圆曲线上点的运算，知道了如何进行点加运算与倍点运算：

对于两个点$P=(x_1,y_1),Q=(x_2,y_2)$，我们定义：

$$\lambda =\{\begin{array}{l}\frac{y_2-y_1}{x_2-x_1}\mathrm{mod}p(P\ne Q)\\ \frac{3x_1^2+a}{2y_1}\mathrm{mod}p\text{else}\end{array}$$

于是我们得到 $P+Q$ 的坐标 $(x_3,y_3)$ 满足：

$$\{\begin{array}{l}{x}_3={\lambda }^2-x_1-x_2\mathrm{mod}p\\ y_3=\lambda (x_1-x_3)-y_1\mathrm{mod}p\end{array}$$

在几何意义上 $P+Q$ 等同于连接 $P,Q$ 的直线与 $E$ 的第三个交点沿 $x$ 轴的对称点。而 $-P$，也就是 $P$ 沿 $x$ 轴的对称点。

对于无穷远点 $\mathcal{O}$，$\mathcal{O}$ 是 $E$ 上所有点构成的循环群中的单位元，因此 $\mathcal{O}+P=P+\mathcal{O}=P$，写代码的时候稍微特判一下即可：

    def calcLam(self, other) -> int:
        x1, x2, y1, y2 = self.x, other.x, self.y, other.y
        
        if x2 == x1 and y2 == y1:
            # P == Q
            if y1 == 0:
                return math.inf
            return self.frac(3*x1*x1 + self.A, 2*y1) # s = lambda
        else:
            # P != Q
            if x2 == x1:
                return math.inf
            return self.frac(y2 - y1, x2 - x1)

    def add(self, other):        
        if self.O:  # is the point at infinity
            # O + P = P
            return EpPoint(other.x, other.y, self.curve)
        if other.O:
            # P + O = P
            return EpPoint(self.x, self.y, self.curve)
        
        s = self.calcLam(other)
        if s == math.inf:
            res = copy.deepcopy(self)
            res.O = True 
            return res
        
        x1, x2, y1, y2 = self.x, other.x, self.y, other.y
        x3 = (s*s - x1 - x2) % self.mod
        y3 = (s * (x1 - x3) - y1) % self.mod
        
        return EpPoint(x3, y3, self.curve)

对于点的乘法，当然可以使用类似于快速幂的算法解决。但只不过之后的内容跟点乘运算关系不大，这里就直接使用普通的累加解决问题，用于调试的正确性检查：

    def mul(self, x): # check only, no performance, calc x*P
        res = self
        for i in range(x - 1):
            res = res.add(self)
        return res

阶的计算

原理也很简单，对一个点进行累加，如果变成无穷远点就代表这个点的阶（也就是累加次数）找到了。

    def calcOrder(self) -> int:
        back = self; cnt = 1

        while True:
            cnt += 1
            back = back.add(self)
            if back.O == 1:
                return cnt

应用密码学教材的内容就到此为止了。

miller 算法

我们首先简要介绍一下有理函数和除子 $\text{div}$ 的概念：

对于任意单变量有理函数：

$$f(X)=\frac{a_0+a_{1}X+a_2{X}^2+\cdots +a_n{X}^n}{b_0+b_{1}X+b_2{X}^2+\cdots +b_m{X}^m}$$

根据代数基本定理，我们总可以在复平面上将其因式分解为：

$$f(X)=\frac{a(X-{\alpha }_1{)}^{e_1}(X-{\alpha }_2{)}^{e_2}\cdots (X-{\alpha }_r{)}^{e_r}}{b(X-{\beta }_1{)}^{d_1}(X-{\beta }_2{)}^{d_2}\cdots (X-{\beta }_s{)}^{d_s}}$$

这里 ${\alpha }_i,{\beta }_j$ 共 $r+s$ 个数互不相同。

为了简便，我们可以使用除子将式子简化：

$$\text{div}(f)=e_1[{\alpha }_1]+e_2[{\alpha }_2]+\cdots +e_r[{\alpha }_r]-d_1[{\beta }_1]-d_2[{\beta }_2]-\cdots -d_r[{\beta }_r]$$

这里的${\alpha }_i$被称为 zeros（零点），${\beta }_i$ 被成为 poles（极点），而前面的系数$e_i,d_i$则被称为 multiplicity（重数）。

当然，这里的 $X$ 大写就意味着 $f$ 的对象不一定只能是一个数，$X$ 的本身也可以由多变量的函数构成，因此 $f$ 也变成了多变量函数。进一步地，如果$f=f(x,y)$ 中的 $(x,y)$ 由椭圆曲线 $E$ 所约束，那么 $E$ 中所对应的点有的可以让分子为零，有的可以让分母为零，因此这些点也可以被分别称为零点和极点。

e.g. 我们令$E:y^2=(x-{\alpha }_1)(x-{\alpha }_2)(x-{\alpha }_3)$, $P\in E$, $f(P)=f(x,y)=y$.

那么 $P_1,P_2,P_3=({\alpha }_1,0),({\alpha }_2,0),({\alpha }_3,0)$ 肯定在 $E$ 上。因此：

$$\text{div}(f)=[P_1]+[P_2]+[P_3]-3[\mathcal{O}]$$

这里$[P_1]+[P_2]+[P_3]$ 很好理解，因为此时 $f(x,y)=y=0$。虽然$f(x,y)$ 并没有分母，但极点$\mathcal{O}$的重数是可以通过射影变换与无穷小分析证明是$3$的，这里篇幅原因不做展开。

经过上述例子我们可以猜想 $E$ 上 $\text{div}(f)$ 的一些性质：

1. $\text{div}(f)=\text{div}(g)$，当且仅当存在常数$c$，$f=cg$.
2. 展开式系数（度数）之和为$0$.
3. 如果对这几个点做点加运算，结果为 $\mathcal{O}$.
4. 特别地，如果不存在零点和极点，那么 $f$ 为常数。

前者被称为除子的唯一性定理，后三者被称为主除子定理。结论与有理函数$f$的选取无关。

---

基础知识铺垫完成，现在引入 miller 算法中位于核心地位的 line function $g_{P,Q}$：

$$g_{P,Q}=\{\begin{array}{l}\frac{y-y_P-\lambda (x-x_P)}{x+x_P+x_Q-{\lambda }^2}(\lambda \ne \mathrm{\infty })\\ x-x_P\text{else}\end{array}$$

显然，$g_{P,Q}$ 也是一个有理函数，我们很容易发现$P,Q$是两个零点，代入一下点加公式可以知道 $P+Q$ 使得分母为零且分子不为零，是极点并且重数为$1$。因此根据主除子定理：

$$\text{div}(g_{P,Q})=[P]+[Q]-[P+Q]-[\mathcal{O}]$$

回到代码实现，由于这里$g_{P,Q}$可以抽象为 $g(x,y)$，并且之后计算$g_{P,Q}$要依赖一个具体的点$S$，因此这里可以直接将$P,Q,S$都作为函数的参数：

    def lineFunc(self, other, S) -> int: # P, Q, S. S is the f_p(S)
        xp, xq, yp, yq = self.x, other.x, self.y, other.y

        if self.calcLam(other) == math.inf:
            return S.x - xp
        lam = self.calcLam(other)
        return self.frac(S.y - yp - lam*(S.x - xp), S.x + xp + xq - lam*lam)

---

然后来到了 miller 算法本体：

设正整数 $m$ 的二进制表示为 $m=\overline{b_{n-1}{b}_{n-2}...b_0},b_i\in \{0,1\},b_{n-1}=1$，以下算法可以生成一个有理函数 $f_P$，使得 $\text{div}(f_P)=m[P]-[mP]-(m-1)[\mathcal{O}]$：

T = P, f = 1
for i from (n-2) to 0:
    f = f^2 * g_{T,T}
    T = 2T
    if b[i] == 1:
        f = f * g_{T,P}
        T = T + P
return f

使用归纳法证明如下：

首先 $m_1=b_{n-1}$ 时，算法返回 $1$ 且 $T=P$，$1$是常数，既没有零点也没有极点。将$m=1$ 代入 $\text{div}(f_P)$，所有的项都消掉了，因此 $m_1=b_{n-1}$ 成立。

设 $m_i=\overline{b_{n-1}{b}_{n-2}...b_{n-i}}$。有 $T=m_{i}P$，$\text{div}(f_P)=m_i[P]-[m_{i}P]-(m_i-1)[\mathcal{O}]$ 成立。我们需要证明 $m_{i+1}=\overline{b_{n-1}{b}_{n-2}...b_{n-i-1}}$ 时的情况成立：

假设 $b_{n-i-1}=0$，则不走 if 分支，我们实际的运算为 $f^{\prime }=f^2\cdot g_{m_{i}P,m_{i}P},m_{i}P=2m_{i}P$。则此时新的 $\text{div}(f^{\prime })$：

$$=2(m_i[P]-[m_{i}P]-(m_i-1)[\mathcal{O}])+([m_{i}P]+[m_{i}P]-[2m_{i}P]-[\mathcal{O}])$$$$=m_{i+1}[P]-2[m_{i}P]-2(m_i-1)[\mathcal{O}]+2[m_{i}P]-[m_{i+1}]P-[\mathcal{O}]$$$$=m_{i+1}[P]-[m_{i+1}]P-(m_{i+1}-1)[\mathcal{O}]$$

结果正确，并且此时 $T^{\prime }=2T=2m_{i}P=m_{i+1}P$，也符合归纳要求，因此归纳在 $b_{n-i-1}=0$ 分支下成立。同理也可证得 $b_{n-i-1}=1$ 分支下成立。

综上所述，该算法可以生成有理函数 $f_P$, 使得 $\text{div}(f_P)=m[P]-[mP]-(m-1)[\mathcal{O}]$。Q.E.D

特别的，如果点 $P\in E[m]$（即 $P$ 在由椭圆曲线上的点构成的阶为 $m$ 的群上），那么 $\text{div}(f_P)=m[P]-m[\mathcal{O}]$，符合之后要提到的 weil pairing 的条件。

回到代码实现，由于 miller 算法是基于 line function $g_{P,Q}$ 的，并且需要一个辅助点 $S$ 进行具体的数值计算，因此编码思路与 line function 一致：

    def miller(self, S) -> int: # P, S, calc f_P(S)
        if self.O:
            return 1
        T = self; f = 1
        n = bin(self.calcOrder())[3:] # 0b1 01001..101

        for bit in n:
            f = f * f * T.lineFunc(T, S) % self.mod
            T = T.add(T)
            if bit == '1':
                f = f * T.lineFunc(self, S) % self.mod
                T = T.add(self) 
        return f

weil pairing

我们设有理函数 $f_P,f_Q$ 满足 $\text{div}(f_P)=m[P]-m[\mathcal{O}],\text{div}(f_Q)=m[Q]-m[\mathcal{O}]$。假设$P,Q\in E[m]$，且 $S$ 不在$P,Q$ 张成的群里，则定义 weil pairing：

$$e_m(P,Q)=\frac{f_P(Q+S)}{f_P(S)}/\frac{f_Q(P-S)}{f_Q(-S)}$$

weil pairing 的值除了与有理函数 $f_P,f_Q$ 以及 $S$ 的选取无关外，还有以下几种性质：

• 单位根：$e_m(P,Q{)}^m=1$
• 双线性：$e_m(P_1+P_2,Q)=e_m(P_1,Q)e_m(P_2,Q)$，$e_m(P,Q_1+Q_2)=e_m(P,Q_1)e_m(P,Q_2)$
• 非退化性：$\mathrm{\forall }P\in E[m],s.t.e_m(P,Q)=1\Rightarrow Q=\mathcal{O}$
• 交错性：$e_m(P,Q)=e_m(Q,P{)}^{-1}$, $e_m(P,P)=1$.

由于我们已经使用 miller 算法计算出了 $f_P(S)$，因此我们可以直接通过 miller 类实现 weil 类：

    def weil(self, Q, S) -> int: # P, Q, S
        negS = EpPoint(S.x, -S.y, self.curve)    # -S

        res1 = self.frac(self.miller(Q.add(S)), self.miller(S))
        res2 = self.frac(Q.miller(self.add(negS)), Q.miller(negS))

        return self.frac(res1, res2)

例如取：

P = EpPoint(36, 60, curve)    # Order 5
Q = EpPoint(121, 387, curve)  # Order 5
S = EpPoint(0, 36, curve)     # Order 130

这里的 $P$ 与 $Q$ “线性无关”，而 $S$ 不在 $P$ 与 $Q$ 的“线性组合”内（即找不到 $m,n$ 使得 $S=mP+nQ$），因此可以用来进行 weil pairing 的计算。

print('P, Q\'s weil pairing:', P.weil(Q, S))
# P, Q's weil pairing: 242
print(242**5 % 631) # property 1 holds
# 1
print('Q, P\'s weil pairing:', Q.weil(P, S))
# Q, P's weil pairing: 279
print(242*279 % 631) # property 4 holds
# 1

然后我们再举另一个例子：

P3 = P.mul(3)
Q4 = Q.mul(4)
print('P3, Q4\'s weil pairing:', P3.weil(Q4, S))
# P3, Q4's weil pairing: 512
print(242**12 % 631) # property 2 holds
# 512

举例说明性质 3 也很简单：

O = P.mul(5)
print('P, O\'s weil pairing:', P.weil(O, S)) # property 3 holds
# P, O's weil pairing: 1

tate pairing

tate pairing 与 weil pairing 的区别是其椭圆曲线定义在有限域 ${\mathbb{F}}_q$ 上，而后者可以是任意域。我们定义：

$$\hat{\tau }(P,Q)=(\frac{f_P(Q+S)}{f_P(S)}{)}^{(q-1)/l}\in {\mathbb{F}}_q$$

其中 $l$ 是质数， $P\in E({\mathbb{F}}_q)[l]，Q\in E({\mathbb{F}}_q)$, $q\equiv 1(\mathrm{mod}l)$.

满足 weil pairing 中的单位根与双线性这两个性质。

可以发现，tate pairing 的计算量比 weil pairing 少一半，所以在密码学中更受青睐。例如 tate pairing 的变种 ate pairing 在以太坊中被广泛使用。

    def tate(self, Q, S) -> int:
        q = self.mod; l = self.calcOrder()
        if q % l != 1:
            raise ValueError('q and l don\'t suuport q == 1 mod l')

        tau = self.frac(self.miller(Q.add(S)), self.miller(S))
        return pow(tau, self.frac(q-1, l), q)

然后我们验证一下性质：

print('P, Q\'s tate pairing:', P.tate(Q, S))
# P, Q's tate pairing: 279
print('Q, P\'s tate pairing:', Q.tate(P, S))
# Q, P's tate pairing: 228

显然$279$ 与 $228$ 都是五次单位根，满足单位根性质；但 $279\times 228\ne 1$，因此交错性不满足。

然后看一下双线性：

P3 = P.mul(3)
Q2 = Q.mul(2)

print('P, Q\'s tate pairing:', P.tate(Q, S))
# P, Q's tate pairing: 279
print('3P, 2Q\'s tate pairing:', P3.tate(Q2, S))
# 3P, 2Q's tate pairing: 279

因为 ${279}^6={279}^{6\mathrm{\%}5}=279$，因此满足双线性性质。

完整代码

# y^2 = x^3 + 30x + 34
import copy, math

class EpCurve:
    def __init__(self, a, b, p) -> None:
        self.A = a
        self.B = b
        self.mod = p
    
    def inverse(self, x:int) -> int:
        return pow(x, -1, self.mod) 

    def frac(self, x:int, y:int) -> int:
        return x * self.inverse(y) % self.mod
    
    def isOnCurve(self, x, y) -> bool:
        return y*y % self.mod == (x*x*x + self.A*x + self.B) % self.mod

class EpPoint(EpCurve):
    def __init__(self, x, y, curve) -> None: 
        super().__init__(curve.A, curve.B, curve.mod)
        if not self.isOnCurve(x, y):
            raise ValueError(f"Point ({x}, {y}) is not on the curve")
        self.x = x % self.mod
        self.y = y % self.mod
        self.O = False
        self.curve = curve

    def __str__(self) -> str:
        if self.O:
            return 'O'
        return f"({self.x}, {self.y})"

    def calcLam(self, other) -> int:
        x1, x2, y1, y2 = self.x, other.x, self.y, other.y
        
        if x2 == x1 and y2 == y1:
            if y1 == 0:
                return math.inf
            return self.frac(3*x1*x1 + self.A, 2*y1) # s = lambda
        else:
            if x2 == x1:
                return math.inf
            return self.frac(y2 - y1, x2 - x1)

    def add(self, other):        
        if self.O:
            return EpPoint(other.x, other.y, self.curve)
        if other.O:
            return EpPoint(self.x, self.y, self.curve)
        
        s = self.calcLam(other)
        if s == math.inf:
            res = copy.deepcopy(self)
            res.O = True
            return res
        
        x1, x2, y1, y2 = self.x, other.x, self.y, other.y
        x3 = (s*s - x1 - x2) % self.mod
        y3 = (s * (x1 - x3) - y1) % self.mod
        
        return EpPoint(x3, y3, self.curve)

    def mul(self, x): # check only, no performance
        res = self
        for i in range(x - 1):
            res = res.add(self)
        return res

    def calcOrder(self) -> int:
        back = self; cnt = 1

        while True:
            cnt += 1
            back = back.add(self)
            if back.O == 1:
                return cnt

    def lineFunc(self, other, S) -> int: # S is the f_p(S)
        xp, xq, yp, yq = self.x, other.x, self.y, other.y

        if self.calcLam(other) == math.inf:
            return S.x - xp
        lam = self.calcLam(other)
        return self.frac(S.y - yp - lam*(S.x - xp), S.x + xp + xq - lam*lam)
    
    def miller(self, S) -> int: # P, S
        if self.O:
            return 1
        T = self; f = 1
        n = bin(self.calcOrder())[3:]

        for bit in n:
            f = f * f * T.lineFunc(T, S) % self.mod
            T = T.add(T)
            if bit == '1':
                f = f * T.lineFunc(self, S) % self.mod
                T = T.add(self) 
        return f
    
    def weil(self, Q, S) -> int: # P, Q, S
        negS = EpPoint(S.x, -S.y, self.curve)    # -S

        res1 = self.frac(self.miller(Q.add(S)), self.miller(S))
        res2 = self.frac(Q.miller(self.add(negS)), Q.miller(negS))

        return self.frac(res1, res2)
    
    def tate(self, Q, S) -> int:
        q = self.mod; l = self.calcOrder()
        if q % l != 1:
            raise ValueError('q and l don\'t suuport q == 1 mod l')

        tau = self.frac(self.miller(Q.add(S)), self.miller(S))
        return pow(tau, self.frac(q-1, l), q)


curve = EpCurve(30, 34, 631)  
P = EpPoint(36, 60, curve)  # Order 5
Q = EpPoint(121, 387, curve)
S = EpPoint(0, 36, curve)

# print('P + Q =', P.add(Q))
# print('P + Q =', P.add(P))
# print('Q + S =', Q.add(S))
# print('P\'s order:', P.calcOrder())
# print('P\'s miller function related to S:', P.miller(S))
# print('P + Q\'s miller function related to S:', P.miller(Q.add(S)))
# print('P, Q\'s weil pairing:', P.weil(Q, S))

P3 = P.mul(3)
Q2 = Q.mul(2)
O = P.mul(5)

print('P, Q\'s tate pairing:', P.tate(Q, S)) 
print('3P, 2Q\'s tate pairing:', P3.tate(Q2, S))

参考资料

• https://github.com/WTFAcademy/WTF-zk/tree/main
• https://crypto.stanford.edu/pbc/notes/ep/
• An Introduction to Mathematical Cryptography