calculate weil-pairing from 0 to 1

To better understand the Weil pairing and the Miller algorithm, I implemented them from scratch in Python here.

The only dependencies used are the pow function, math.inf, and the deep copy function copy.deepcopy. I also adopted the OOP programming paradigm, which I'm not entirely comfortable with.

Although after writing this, my understanding of both the Weil pairing and OOP has improved (

Starting from ${\mathbb{F}}_p$

In high school competitions, we learned about operations in ${\mathbb{F}}_p$. Its addition, subtraction, and multiplication are straightforward—simply apply $\mathrm{mod}p$ to the ordinary operations (here referring to the addition, subtraction, and multiplication in code). For division, we learned at least two methods (Fermat's Little Theorem and the extended Euclidean algorithm) to compute the multiplicative inverse of an element in ${\mathbb{F}}_p$, thereby converting division into multiplication.

However, there's no need to revisit the extended Euclidean algorithm today, as a previous post has already explained its principles in detail. Here, I choose to directly use Python's syntactic sugar:

pow(x, -1, p)

to compute $x^{-1}\mathrm{mod}p$.

Thus, for a fraction $\frac{a}{b}$, we can also express it as $a\times b^{-1}\mathrm{mod}p$. With this, we have fully covered operations in ${\mathbb{F}}_p$.

Introducing Elliptic Curves

Here we introduce the protagonist of this post, an elliptic curve $E$ defined over ${\mathbb{F}}_{631}$:

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

Here, we use a class EpCurve to abstract it. Currently, this class already includes methods for modular inversion and division. To reflect the characteristics of this curve, we can add a method to check whether a point $(x,y)$ lies on the curve. Thus, we can write the complete implementation:

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

And the corresponding instance curve = EpCurve(30, 34, 631).

Then, based on EpCurve, we can abstract points on the elliptic curve with 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})"

Here, self.x and self.y are the coordinates $(x,y)$ of the point, while self.O indicates whether it is the point at infinity $\mathcal{O}$. The __str__ method is a built-in magic method in Python, used to return the corresponding string when the instance is printed.

Point Operations

In the Applied Cryptography course during the second semester of my sophomore year, we systematically studied point operations on elliptic curves and learned how to perform point addition and point doubling:

For two points $P=(x_1,y_1)$ and $Q=(x_2,y_2)$, we define:

$$\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}$$

Then we obtain the coordinates $(x_3,y_3)$ of $P+Q$ satisfying:

$$\{\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}$$

Geometrically, $P+Q$ is equivalent to the reflection across the x-axis of the third intersection point of the line connecting $P$ and $Q$ with the curve $E$. Meanwhile, $-P$ is the reflection of $P$ across the x-axis.

For the point at infinity $\mathcal{O}$, $\mathcal{O}$ is the identity element in the cyclic group formed by all points on $E$. Therefore, $\mathcal{O}+P=P+\mathcal{O}=P$. In code, we can handle this with a simple special case:

    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)

For point multiplication, we can of course use an algorithm similar to fast exponentiation. However, since the subsequent content is not heavily related to point multiplication, we simply use ordinary accumulation here for correctness checking during debugging:

    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

Calculation of Order

The principle is also quite simple: continuously add a point to itself, and if it becomes the point at infinity, the order (i.e., the number of additions) has been found.

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

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

The content from the applied cryptography textbook ends here.

Miller's Algorithm

We first briefly introduce the concepts of rational functions and the divisor $\text{div}$:

For any univariate rational function:

$$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}$$

According to the fundamental theorem of algebra, we can always factorize it in the complex plane as:

$$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}}$$

Here, the ${\alpha }_i$ and ${\beta }_j$ are $r+s$ distinct numbers.

For simplicity, we can use the divisor to simplify the expression:

$$\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]$$

Here, ${\alpha }_i$ are called zeros, ${\beta }_i$ are called poles, and the coefficients $e_i,d_i$ are called multiplicities.

Of course, the uppercase $X$ implies that the object of $f$ is not necessarily just a number; $X$ itself can also be composed of multivariate functions, so $f$ becomes a multivariate function. Furthermore, if the $(x,y)$ in $f=f(x,y)$ are constrained by an elliptic curve $E$, then some points in $E$ can make the numerator zero, and others can make the denominator zero, so these points can also be called zeros and poles, respectively.

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

Then $P_1,P_2,P_3=({\alpha }_1,0),({\alpha }_2,0),({\alpha }_3,0)$ must lie on $E$. Therefore:

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

Here, $[P_1]+[P_2]+[P_3]$ is easy to understand because $f(x,y)=y=0$. Although $f(x,y)$ has no denominator, the multiplicity of the pole $\mathcal{O}$ can be proven to be $3$ through projective transformation and infinitesimal analysis, which we omit here due to space constraints.

From the above example, we can conjecture some properties of $\text{div}(f)$ on $E$:

1. $\text{div}(f)=\text{div}(g)$ if and only if there exists a constant $c$ such that $f=cg$.
2. The sum of the expansion coefficients (degrees) is $0$.
3. If point addition is performed on these points, the result is $\mathcal{O}$.
4. In particular, if there are no zeros or poles, then $f$ is a constant.

The first is called the uniqueness theorem of divisors, and the latter three are called the principal divisor theorem. The conclusions are independent of the choice of the rational function $f$.

---

With the basics covered, we now introduce the line function $g_{P,Q}$, which is central to Miller's algorithm:

$$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}$$

Clearly, $g_{P,Q}$ is also a rational function. We can easily see that $P$ and $Q$ are two zeros. Substituting the point addition formula, we find that $P+Q$ makes the denominator zero while the numerator is nonzero, so it is a pole with multiplicity $1$. Therefore, by the principal divisor theorem:

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

Returning to code implementation, since $g_{P,Q}$ can be abstracted as $g(x,y)$ and the subsequent calculation of $g_{P,Q}$ depends on a specific point $S$, we can directly take $P,Q,S$ as function parameters:

    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)

---

Now we come to the main body of Miller's algorithm:

Let the binary representation of the positive integer $m$ be $m=\overline{b_{n-1}{b}_{n-2}...b_0},b_i\in \{0,1\},b_{n-1}=1$. The following algorithm generates a rational function $f_P$ such that $\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

Proof by induction:

First, when $m_1=b_{n-1}$, the algorithm returns $1$ and $T=P$. $1$ is a constant with no zeros or poles. Substituting $m=1$ into $\text{div}(f_P)$, all terms cancel out, so $m_1=b_{n-1}$ holds.

Assume $m_i=\overline{b_{n-1}{b}_{n-2}...b_{n-i}}$. Let $T=m_{i}P$ and $\text{div}(f_P)=m_i[P]-[m_{i}P]-(m_i-1)[\mathcal{O}]$ hold. We need to prove the case for $m_{i+1}=\overline{b_{n-1}{b}_{n-2}...b_{n-i-1}}$:

Suppose $b_{n-i-1}=0$. Then the if branch is skipped, and the actual operation is $f^{\prime }=f^2\cdot g_{m_{i}P,m_{i}P},m_{i}P=2m_{i}P$. The new $\text{div}(f^{\prime })$ is:

$$=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}]$$

The result is correct, and $T^{\prime }=2T=2m_{i}P=m_{i+1}P$ also meets the induction requirement. Therefore, the induction holds for the $b_{n-i-1}=0$ branch. Similarly, it can be proven for the $b_{n-i-1}=1$ branch.

In conclusion, this algorithm generates a rational function $f_P$ such that $\text{div}(f_P)=m[P]-[mP]-(m-1)[\mathcal{O}]$. Q.E.D.

In particular, if the point $P\in E[m]$ (i.e., $P$ is on the subgroup of order $m$ on the elliptic curve), then $\text{div}(f_P)=m[P]-m[\mathcal{O}]$, which satisfies the condition for the Weil pairing to be discussed later.

Returning to code implementation, since Miller's algorithm is based on the line function $g_{P,Q}$ and requires an auxiliary point $S$ for specific numerical calculations, the coding approach is consistent with that of the 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

Let rational functions $f_P$ and $f_Q$ satisfy $\text{div}(f_P)=m[P]-m[\mathcal{O}]$ and $\text{div}(f_Q)=m[Q]-m[\mathcal{O}]$. Assume $P,Q\in E[m]$, and $S$ is not in the group generated by $P$ and $Q$. Then the Weil pairing is defined as:

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

Apart from being independent of the choice of rational functions $f_P$, $f_Q$, and $S$, the Weil pairing also has the following properties:

• Root of unity: $e_m(P,Q{)}^m=1$
• Bilinearity: $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)$
• Non-degeneracy: $\mathrm{\forall }P\in E[m],s.t.e_m(P,Q)=1\Rightarrow Q=\mathcal{O}$
• Alternating: $e_m(P,Q)=e_m(Q,P{)}^{-1}$, $e_m(P,P)=1$.

Since we have already computed $f_P(S)$ using Miller's algorithm, we can directly implement the Weil class via the Miller class:

    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)

For example, take:

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

Here, $P$ and $Q$ are "linearly independent", and $S$ is not in the "linear span" of $P$ and $Q$ (i.e., no $m,n$ exist such that $S=mP+nQ$), so it can be used to compute the 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

Now consider another example:

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

Demonstrating property 3 is also straightforward:

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

The difference between the Tate pairing and the Weil pairing is that the elliptic curve for the Tate pairing is defined over a finite field ${\mathbb{F}}_q$, while the latter can be over any field. We define:

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

where $l$ is a prime, $P\in E({\mathbb{F}}_q)[l]$, $Q\in E({\mathbb{F}}_q)$, and $q\equiv 1(\mathrm{mod}l)$.

It satisfies the two properties of the Weil pairing: being a root of unity and bilinearity.

It can be observed that the computational cost of the Tate pairing is half that of the Weil pairing, making it more favored in cryptography. For example, the variant of the Tate pairing, the ate pairing, is widely used in Ethereum.

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

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

Then we verify the properties:

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

Clearly, $279$ and $228$ are both fifth roots of unity, satisfying the root of unity property; however, $279\times 228\ne 1$, so the alternating property is not satisfied.

Now, let's check the bilinearity:

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

Since ${279}^6={279}^{6\mathrm{\%}5}=279$, the bilinearity property is satisfied.

Complete Code

# 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))

References

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