Issue with Miller's Algorithm in section 5.4.4

[Basem815]

Hi there appears to be an error with algorithm 8 in chapter 5 section 5.4.4. If the algorithm is implemented as shown in the manual using sage, then it will return an incorrect result (or in most cases, a division by zero error) as demonstrated below:

    def MillersAlgorithm(P, Q, E, r):
        ''' Miller’s algorithm for Short Weierstrass curves y^2 = x^3 + ax + b
            Input: P, Q are elements of the full r-torison elliptic curve E[r] 
                   r = b_0*2^0 + .... + b_t*2^t  with b_t = 1  '''
        INF = E(0)
        a = E.a4()
        if P == INF or Q == INF or P == Q:
            return (-1)**r
        b = r.binary()[::-1]    # Reverses the string s.t: b[0] = b_0 
        (x_p, y_p) = P.xy()
        (x_q, y_q) = Q.xy()
        (x_t, y_t) = P.xy()
        f_1, f_2 = 1, 1
        t = len(b) - 1
        for j in range(t-1, -1, -1):
            m = (3*x_t^2+a)/(2*y_t)
            f_1 = f_1^2 * (y_q-y_t-m*(x_q-x_t))
            f_2 = f_2^2 * (x_q + (2*x_t) - m^2)
            x_2t = m^2 - 2*x_t
            y_2t = -y_t - m*(x_2t - x_t)
            (x_t, y_t) = (x_2t, y_2t)
            if b[j] == '1':
                m = (y_t - y_p)/(x_t - x_p)
                f_1 = f_1 * (y_q-y_t-m*(x_q-x_t))
                f_2 = f_2^2 * (x_q + (x_p + x_t) - m^2)
                x_t_plus_p = m^2 - x_t - x_p
                y_t_plus_p = -y_t - m * (x_t_plus_p - x_t)
                (x_t, y_t) = (x_t_plus_p, y_t_plus_p)
        f_1 = f_1 * (x_q - x_t)
        return f_1/f_2         

    
    def WeilPairing(P,Q,curve,r):
        return (MillersAlgorithm(P,Q,curve,r)/ MillersAlgorithm(Q,P,curve,r))*(-1)**r

The algorithm should return the same value for example 98 as was demonstrated in the book using the built-in sage command for computing the Weil Pairing, but it returns an error instead:

    print("Example 98 (MMM)")
    F13 = GF(13)
    F13t.<t> = F13[]
    P_MOD_4 = F13t(t^4+2)
    F13_4.<t> = GF(13^4, name = 't', modulus = P_MOD_4)
    TJJF13_4 = EllipticCurve(F13_4, [8,8])
    P = TJJF13_4([7,2])
    Q = TJJF13_4([9*t^2+7, 12*t^3+2*t])
    r = 5 
    print("Weil-pairing using sage: ", P.weil_pairing(Q,r)) 
        # Output: 7*t^3 + 7*t^2 + 6*t + 3
    print("Weil-pairing of P and Q:", WeilPairing(P,Q,TJJF13_4,r))
        # Output: ZeroDivisionError: division by zero in finite field (In particular when computing the gradient m = (y_t - y_p)/(x_t - x_p) )

Issue

The issue is evidently when we are dealing with the case where the slope is not defined (vertical line). While we can add extra conditions to make checks for it, it'll make the algorithm more convoluted than it already appears so. In particular, because we have separated the computations of the x and y coordinates of point T in the algorithm above, it is unclear what to assign these coordinates when our gradient is infinity (E(0).xy() returns an error as the point at infinity does not have finite x and y coordinates).

Solution

Since the manual did not intend to elaborate on the concept of divisors, I won't attempt to myself. However, I think it helps to get an intuitive idea of what Miller's algorithm is intended to do.

Referring to the resources mentioned in section 5.4.4 (in particular, to Hoffstein et al. [2008]), Miller's algorithm is described as an double-and-add method that can be used to efficiently compute the Weil pairing. It makes use of the line function between two points that is described by the following equation:

$$g_{P,Q} = \Bigg\{ \begin{array}{ll} \frac{y - y_{P} - \lambda(x-x_{P})}{x + x_{P} + x_{Q} - \lambda^2} \quad \quad if \quad \quad \lambda \neq \infty \\\ x - x_{P} \quad \quad \quad \quad if \quad \quad \lambda = \infty \end{array} \Bigg.$$

where $\lambda$ is the slope of the line connecting P and Q.

In my opinion, using this alternative definition of Miller's algorithm is more concise and comprehensible:

$$\begin{array}{llllllllll} \text{[1] \quad Set } T = P \text{ and } f = 1 \\\ \text{[2] \quad Loop } i = n-2 \text{ down to } 0 \\\ \text{[3]}\quad \quad \quad \text{ Set } f = f^2 . g_{T,T} \\\ \text{[4]}\quad \quad \quad \text{ Set } T = 2T \\\ \text{[5]}\quad \quad \quad \text{ If } m_i = 1 \\\ \text{[6]}\quad \quad \quad \quad \quad \text{ Set } f = f . g_{T,P} \\\ \text{[7]}\quad \quad \quad \quad \quad \text{ Set } T = T + P \\\ \text{[8]}\quad \quad \quad \text{ End If } \\\ \text{[9]}\quad \text{ End *i* loop } \\\ \text{[10]}\quad \text{Return the value f} \end{array}$$

This version addresses the issue of when the two points are inverses of each other (or when the line between the two points intersects the curve at the point at infinity in projective coordinates).
Although the book shows a different approach to computing the Weil-pairing using this definition of Miller's algorithm, we can use the original definition as provided in the manual (5.49) to obtain the correct result as demonstrated below:

def g(P, Q, X, E):
    """ First computes the slope (lambda) of the line between P and Q
        Then computes the line function (g_{P,Q} evaluated at point X """
    
    x_p, y_p = P.xy()
    x_q, y_q = Q.xy()
    (x,y) = X.xy() 
    if Q == -P: # Case: lambda = infinity (Handles the case when y_p = 0)
        # g_{P,Q}(X) = x - x_p
        return x - x_p
    if P == Q and P.xy()[1] != 0: # Case: Point doubling
        # Case: lambda = (3*x_p^2 + a)/(2*y_p)
        a = E.a4()
        m = (3*x_p^2 + a)/(2*y_p)
    else:
        # Case: lambda = Point addition
        # lambda = (y_q - y_p)/(x_q - x_p)
        m = (y_q - y_p)/(x_q - x_p)

    return ( y - y_p - m*(x - x_p)) / (x+x_p + x_q - m^2)

def millerAtX(P, X, E, r):
    a = E.a4()
    Fp = E.base_field()
    b = r.binary()[::-1]    # Reverses the string s.t: b[0] = b_0 
    n = len(b) 
    T = P
    f = Fp(1)
    for i in range (n-2, -1, -1):
        f = f^2 * g(T, T, X, E)
        T = 2*T
        if b[i] == "1":
            f = f * g(T, P, X, E)
            T = T + P
    return f


def weilPairing(P,Q,curve,r):
    return (millerAtX(P,Q,curve,r)/ millerAtX(Q,P,curve,r))*(-1)**r

Optional

For completion purposes, I will also demonstrate how to compute the Weil-Pairing using the approach outlined in section 5.8.3 Hoffstein et al. [2008].

The Weil-pairing, denoted as $e_m(P,Q)$ as defined in the book is given as:

$$e_m(P,Q) = \left. \frac{f_{P}(Q+S)}{f_P(S)} \middle/ \frac{f_{Q}(P-S)}{f_Q(-S)} \right.$$

Here, $S \in E$ is a point that is not in the subgroup spanned by P and Q (i.e. $S \notin { \infty, P, -Q, P-Q }$ ).
Using Sage we can define the following function:

def weilPairingAlt(P,Q,curve,r):
    return (millerAtX(P, Q + S, E, r) / millerAtX(P, S, E, r)) / (millerAtX(Q, P - S, E, r) / millerAtX(Q, -S, E, r)) 

The same example can then be computed with this alternative definition using:

print("Example 98 (MMM)")
# Test (Example 98)
F13 = GF(13)
F13t.<t> = F13[]
P_MOD_4 = F13t(t^4+2)
F13_4.<t> = GF(13^4, name = 't', modulus = P_MOD_4)
TJJF13_4 = EllipticCurve(F13_4, [8,8])
P = TJJF13_4([7,2])
Q = TJJF13_4([9*t^2+7, 12*t^3+2*t])
S = TJJF13_4.random_element()
INF = TJJF13_4(0) # 0 is the identity element
while S == INF or S == P or S == -Q or S == P - Q:
    S = TJJF13_4.random_element()
r = 5 
print("Weil-pairing of P and Q:", weilPairing(P,Q, S, TJJF13_4,r))
# Weil-pairing of P and Q: 7*t^3 + 7*t^2 + 6*t + 3
print("Weil-pairing using sage: ", P.weil_pairing(Q,r)) 
# Weil-pairing using sage:  7*t^3 + 7*t^2 + 6*t + 3

I have not been able to figure out how to derive one formula from the other yet, but any help with that would be appreciated!