We give algorithms for computing the Weil and the Tate pairing on the elliptic curves E/(Subscript[𝔽, p]) :  (y^2) = 
        (x^3) + a x (p  =  3 (mod 4)) and E/(Subscript[𝔽, p]) :  (y^2) = 
        (x^3) + b (p = 2 (mod 3)). Since these curves are class (1) supersingular curves, they have embedding degree k = 2. Therefore, we consider the field extension (Subscript[𝔽, (p^2)]) of degree 2. We apply the algorithms to the curves (Subscript[E, 1])/(Subscript[𝔽, 11]) : (y^2) = 
        (x^3) + 3  x and (Subscript[E, 2])/(Subscript[𝔽, 131]) :  (y^2) = 
        (x^3) + 31.

<< Algebra`FiniteFields`;

fld = GF[p, {1, 0, 1}];

We start with defining some functions and modules that make our life easier, such as the function Pt[a,b,c,d] that gives us the point (a + b i, c + d i), and the module Neg[P]that returns the point on the curve given by -P.

Pt[a_, b_, c_, d_] := {fld[{a, b}], fld[{c, d}]}

Neg[P_] := 
      Module[

  

{nP}, 
        If[P == 𝒪, 
          nP = 𝒪, , 
If[
            P[[2]] == fld[{0, 0}], nP = P, , 
            nP = {P[[1]], -P[[2]]}]]; 
nP]

On input of two points P, 
      Q ϵ E((Subscript[𝔽, (p^2)])), the module PointAddition[P,Q] computes the sum of P and Q.

PointAddition[P_, Q_] := 
      Module[

 (* 
          Input :  
            two points P and Q on E ((Subscript[𝔽, (p^2)])); 
Output :  R,  
          the sum of P and Q on E ((Subscript[𝔽, (p^2)])). *) 

{λ, R1, R2, 
          R}, 
If[
          P == 𝒪, 
If[
            Q == 𝒪, R = 𝒪, , 
            R = Q], , 
If[Q == 𝒪, 
            R = P, , 
If[
              P == 
                Neg[Q], 
R = 𝒪, , 
If[
                P == 
                  Q, λ = ((3 P[[1]]^2 + a)/(2 P[[2]])), , λ = ((Q[[2]] - P[[2]])/(Q[[1]] - P[[1]]))]; 
If[R1 = (λ^2) - P[[1]] - Q[[1]]; 
                R1 == 0, R1 = fld[{0, 0}]]; 
If[
                R2 = λ (P[[1]] - R1) - P[[2]]; 
                R2 == 0, 
                R2 = fld[{0, 0}]]; 
R = {R1, 
                  R2}]]]; 
R]

We use the module above to define the module PointMultiplication[P,k], which enables us to multiply a point \(TraditionalForm\`Pϵ\ \(E(𝔽\_\(p\^2\))\)\) with a scalar k.

PointMultiplication[P_, k_] := 
      Module[

  

{binexp, l, 
          H}, 
binexp = IntegerDigits[k, 2]; 
        l = Length[binexp]; 
        H = P; 
For[j = 2, j ≤ l, 
          H = PointAddition[H, H]; 
          If[binexp[[j]] == 1, 
            H = PointAddition[H, 
                P]]; j++]; 
H]

We rewrite the module PointAddition[P,Q] to the module DivisorAddition[D1,D2,x,y], which adds two divisors in canonical form and expresses the sum of those divisors in canonical form. So the module takes as input a pair of divisors (Subscript[D, 1]) = {(Subscript[P, 1]), 
            (Subscript[f, 1])} = ((Subscript[P, 1])) - (𝒪) + div((Subscript[f, 1])) and (Subscript[D, 2]) = {(Subscript[P, 2]), 
            (Subscript[f, 2])} = ((Subscript[P, 2])) - (𝒪) + div((Subscript[f, 2])), and the variables x, y for the functions, and outputs (Subscript[D, 1]) + 
          (Subscript[D, 2]) = (Subscript[D, 3]) = {(Subscript[P, 3]), 
              (Subscript[f, 3])} = ((Subscript[P, 3])) - (𝒪) + 
              div((Subscript[f, 3])).

DivisorAddition[D1_, D2_, x_, y_] := 
      Module[

 (* 
          Input :  
              two divisors in canonical form D1  =  {P1, 
                  f1} and D2  =  {P2, 
                  f2} and the variables x and y for the function; 
Output :  the divisor D3  =  {P3, f3}, 
          where D3 is the sum of D1 and D2 in canonical form. *) 

{λ, R, R1, R2, l, 
          v}, 
If[
          D1[[1]] == 𝒪, 
If[
            D2[[1]] == 𝒪, R = 𝒪; 
            l = 1; v = 1, , R = D2[[1]]; l = x - D2[[1, 1]]; 
            v = x - D2[[1, 1]]], , 
If[
            D2[[1]] == 𝒪, R = D1[[1]]; 
            l = x - D1[[1, 1]]; 
            v = x - D1[[1, 1]], , 
If[
              D1[[1]] == 
                Neg[D2[[1]]], 
R = 𝒪; 
              l = x - D1[[1, 1]]; 
              v = 1, , 
If[
                D1[[1]] == 
                  D2[[1]], λ = ((3 D1[[1, 1]]^2 + a)/(2 D1[[1, 2]])), , λ = ((D2[[1, 2]] - D1[[1, 2]])/(D2[[1, 1]] - D1[[1, 1]]))]; 
If[
                R1 = (λ^2) - D1[[1, 1]] - D2[[1, 1]]; 
                R1 == 0, R1 = fld[{0, 0}]]; 
If[
                R2 = λ (D1[[1, 1]] - R1) - D1[[1, 2]]; 
                R2 == 0, R2 = fld[{0, 0}]]; R = {R1, R2}; 
              l = y - λ x + (λ D1[[1, 1]] - 
                      D1[[1, 2]]); 
              v = x - R1]]]; 
{R, ((D1[[2]] D2[[2]] l)/v)}]

Similarly, the module DivisorMultiplication[D1,k,x,y] multiplies a divisor (Subscript[D, 1]) in canonical form with an integer k, and expresses the resulting divisor in canonical form as well.

DivisorMultiplication[D1_, k_, x_, y_] := 
      Module[

  

{binexp, l, 
          D2}, 
binexp = IntegerDigits[k, 2]; 
        l = Length[binexp]; 
        D2 = D1; 
For[j = 2, j ≤ l, 
          D2 = DivisorAddition[D2, D2, x, y]; 
          If[binexp[[j]] == 1, 
            D2 = DivisorAddition[D2, D1, x, 
                y]]; j++]; 
D2]

The following module randomly picks a point P  ≠ 𝒪 on the curve over the ground field E((Subscript[𝔽, p])). This is done by choosing an x-coordinate at random, and solving the resulting quadratic equation for y, if possible. We assume that the prime p satisfies p = 3 (mod 4), because this greatly simplifies the extraction of square roots in the finite field (Subscript[𝔽, p]).

RandomPointGroundField := 
      Module[

  

{P1, P2}, 
        While[P1 = Random[Integer, {0, p - 1}]; 
          JacobiSymbol[P1^3 + a P1 + b, 
              p] == -1]; 
If[
          Random[Integer] == 1, 
          P2 = Mod[(((P1^3 + a P1 + b))^(((p + 1)/4))), p], 
          P2 = Mod[-(((P1^3 + a P1 + b))^(((p + 1)/4))), 
              p]]; 
Pt[P1, 0, P2, 0]]

Let P  and Q  be two m-torsion points on E((Subscript[𝔽, (11^2)])). The Weil pairing requires two points T ≠ -P and U ≠ -Q such that T ≠ U, Q + U and P + T ≠ U, Q + U. We use the module RandomPointGroundField to pick random points on the curve over the ground field, and check whether they satisfy the conditions. Note that taking points on the curve over the ground field rather than over the extension field does not affect the outcome of the pairing, but it does speed up the computations.

PickTUGroundField[P_, Q_] := 
      Module[

  

{T, U, 
          PT, QU},  
While[T = RandomPointGroundField; 
          PT = PointAddition[P, T]; U = RandomPointGroundField; 
          QU = PointAddition[Q, U]; 
P == 
              Neg[T] ∨ Q == Neg[U] ∨ PT == U ∨ 
            PT == QU ∨ T == U ∨ 
            T == QU]; 
{T, U, PT, QU}]

For the Weil pairing (Subscript[e, m])(P, Q) to be non-trivial, the two m-torsion points P and Q need to be independent. The module DistortionMap[P] computes a point Q of the same order as, but independent from P. Note that ϕ : (x, y) -> (-x, i y) is a distortion map for the curve (y^2) = (x^3) + a x , whereas a distortion map for (y^2) = (x^3) + b is given by ϕ : (x, y) -> (ω x, 
          y), with ω a primitive third root of unity.

DistortionMap[P_] := 
      Module[

 (* 
          Input :  point P; 
Outup :  
            point ϕ (P),  
          where ϕ is the distortion map. *) 

{Q, ω}, 
        If[P == 𝒪, Q = 𝒪, , 
          If[b == 0, 
            Q = {If[P[[1]] == fld[{0, 0}], 
                  fld[{0, 0}], , -P[[1]]], 
                If[P[[2]] == fld[{0, 0}], fld[{0, 0}], , 
                  fld[{0, 
                        1}] P[[2]]]}, 
While[ω = fld[{Random[Integer, {1, p - 1}], Random[Integer, {0, p - 1}]}]^(((p^2) - 1)/3); ω == fld[{1, 0}]]; 
            Q = {If[P[[1]] == fld[{0, 0}], 
                  fld[{0, 0}], , ω P[[1]]], 
                P[[2]]}]]; 
Q]

On input of an m-torsion point P, an arbitrary T, the value m, and variables x, y, the module ComputeFunctionWeil[P,T,m,x,y] computes the function f(x, y) such that div(f) = m(P + T) - m(T).

ComputeFunctionWeil[P_, T_, m_, x_, y_] := 
      Module[

  

{mD1, 
          mD2}, 
mD1 = 
          DivisorMultiplication[{PointAddition[P, T], fld[{1, 0}]}, m, x, y]; 
        mD2 = DivisorMultiplication[{T, 1}, m, x, 
            y]; 
mD1[[2]]/mD2[[2]]]

We are now ready to give the algorithm for the computation of the Weil pairing (Subscript[e, m])(P, Q) of two m-torsion points. The module WeilPairing[P,Q,m] first picks points T and U, such that the resulting functions (Subscript[f, A]) and (Subscript[f, B]) are defined at Q + U, U and P + T, T respectively, and then computes \(TraditionalForm\`\(e\_m\)(P, 
          S) = \(\(\(f\_A\)(Q + U)\) \(\(f\_B\)(T)\)\)\/\(\(\(f\_A\)(U)\) \(\(f\_B\)(P + T)\)\)\).

WeilPairing[P_, Q_, m_] := 
      Module[

  

{T, 
          U, PT, QU, fA, fB, x, y, DenfA, 
          DenfB}, 
While[{T, U, PT, QU} = 
            PickTUGroundField[P, Q]; 
fA[x_, y_] = 
            ComputeFunctionWeil[P, T, m, x, y]; 
fB[x_, 
              y_] = ComputeFunctionWeil[Q, U, m, x, y]; 
          DenfA[x_, y_] = Denominator[fA[x, y]]; 
          DenfB[x_, y_] = Denominator[fB[x, y]]; 
          DenfA[U[[1]], U[[2]]] == 0 ∨ 
            DenfA[QU[[1]], QU[[2]]] == 0 ∨ 
            DenfB[T[[1]], T[[2]]] == 0 ∨ 
            DenfB[PT[[1]], PT[[2]]] == 
              0]; 
(fA[QU[[1]], QU[[2]]] fB[
                T[[1]], T[[2]]])/(fA[U[[1]], U[[2]]] 
              fB[PT[[1]], PT[[2]]])]

We next deal with the Tate pairing. On input of an m-torsion point P, the value m, and variables x, y, the following module computes the function g(x, y) such that div(g) = m(P) - m(𝒪). Note that ComputeFunctionWeil[P,𝒪,m,x,y] = ComputeFunctionTate[P,m,x,y].

ComputeFunctionTate[P_, m_, x_, y_] := 
      Module[

 (* 
          Input :  an m - torsion point P,  the value m,  
          and the variables x and y for the function; 
Output :  function g such that div (g)  =  
            m (P) - 
              m (𝒪). *) 

{mD1}, 
mD1 = 
          DivisorMultiplication[{P, fld[{1, 0}]}, m, x, 
            y]; 
mD1[[2]]]

Now the module TatePairing[P,Q,m] proceeds as follows. After computing a function g(x, y) with div(g) = m(P) - m(𝒪), the algorithm picks a point Q' such that g is defined at Q + Q' and Q', and computes ⟨P, Q⟩_m = g(
            Q + Q')/g(Q').

TatePairing[P_, Q_, m_] := 
      Module[

  

{S, QP, g, Deng, Numg, x, 
          y}, 
g[x_, y_] = 
          ComputeFunctionTate[P, m, x, y]; 
Deng[x_, y_] = 
          Denominator[g[x, y]]; 
Numg[x_, y_] = 
          Numerator[g[x, y]]; 
While[
          QP = RandomPointGroundField; S = PointAddition[Q, QP]; 
          S == 𝒪 ∨ 
            Deng[S[[1]], S[[2]]] == 0 ∨ 
            Numg[QP[[1]], QP[[2]]] == 0 ∨ 
            Numg[S[[1]], S[[2]]] == 0 ∨ 
            Deng[QP[[1]], QP[[2]]] == 
              0]; 
(Numg[S[[1]], S[[2]]] Deng[
                QP[[1]], QP[[2]]])/(Deng[S[[1]], 
                S[[2]]] Numg[QP[[1]], QP[[2]]])]

The outcome of the Tate pairing is an equivalence class in \(TraditionalForm\`\(𝔽\_\(p\^2\)\^*\)/\((\(𝔽\_\(p\^2\)\^*\))\)\^m\). Whenever a unique value is required, one can raise the outcome to the power ((p^2) - 1)/m, to eliminate m-th powers.

TatePairingUnique[P_, Q_, m_] := 
      TatePairing[P, Q, m]^(((p^2) - 1)/m)

Time for some numerical examples. We apply the Weil and the Tate pairing to the elliptic curve (Subscript[E, 1])/(Subscript[𝔽, 11]) : (y^2) = 
        (x^3) + 3  x.

Clear[fld, p, a, b];
a = 3; b = 0; p = 11;
fld = GF[p, {1, 0, 1}];
PrimeQ[p]
Mod[p, 4] == 3
FactorInteger[p + 1]

The point P = (3, 5) has order 3:

P = Pt[3, 0, 5, 0]
P2 = PointMultiplication[P, 2]
P3 = PointMultiplication[P, 3]

When we apply the Weil pairing to two dependent m-torsion points, the outcome is always 1:

WeilPairing[P, P, 3]
WeilPairing[P, P2, 3]
WeilPairing[P, P3, 3]

We apply the distortion map to P to find an independent 3-torsion point:

Q = DistortionMap[P]

We now compute the Weil pairing (Subscript[e, 3])(P, Q), which is indeed a third root of unity:

e = WeilPairing[P, Q, 3] 
(e^3)

Similarly, we can compute the Tate pairing ⟨P, Q⟩_3. The outcome of this pairing is an equivalence class, so the evaluation does not always give the same result. However, if we eliminate the third powers, as is done in the module TatePairingUnique[P,Q,m], we get the unique representative for that particular equivalence class:

TatePairing[P, Q, 3]
TatePairing[P, Q, 3]

TatePairingUnique[P, Q, 3]
TatePairingUnique[P, Q, 3]

For the Tate pairing, the second point is not required to be an m-torsion point. In fact, the Tate pairing can be evaluated for any point on the curve. If \(TraditionalForm\`R\ \ ϵ\ \ m\ \(E(𝔽\_\(p\^2\))\)\), then the outcome of the pairing ⟨P, R⟩_m is the equivalence class representing 1.

arb = 
      PointAddition[RandomPointGroundField, 
        DistortionMap[RandomPointGroundField]]
R = PointMultiplication[arb, 3]
TatePairing[P, R, 3]
TatePairingUnique[P, R, 3]

Furthermore, for \(TraditionalForm\`R\ \ ϵ\ \ m\ \(E(𝔽\_\(p\^2\))\)\) the result of the pairing \(TraditionalForm\`⟨P, Q + R⟩\_m\) is the equivalence class of ⟨P, Q⟩_m.

TatePairing[P, PointAddition[Q, R], 3]
TatePairingUnique[P, PointAddition[Q, R], 3]

Our next example is the curve (Subscript[E, 2])/(Subscript[𝔽, 131]) : (y^2) = 
        (x^3) + 31.

Clear[fld, p, a, b, P, Q];
a = 0; b = 31; p = 131;
fld = GF[p, {1, 0, 1}];
PrimeQ[p]
Mod[p, 4] == 3
Mod[p, 3] == 2
FactorInteger[p + 1]

Consider the point P = (78, 25), which has order 11:

P = Pt[78, 0, 25, 0]
PointMultiplication[P, 11]

Again, the Weil pairing of two dependent m-torsion points is trivial.

P7 = PointMultiplication[P, 7]
WeilPairing[P, P7, 11]

So we apply the distortion map to find an independent 11-torsion point, compute the Weil pairing (Subscript[e, 11])(P, Q) and show that the result is an 11-th root of unity

Q = DistortionMap[P]
e = WeilPairing[P, Q, 11]
(e^11)

For the Tate pairing, we perform the same computations as in the previous example.

TatePairing[P, Q, 11]
TatePairing[P, Q, 11]

TatePairingUnique[P, Q, 11]
TatePairingUnique[P, Q, 11]

arb = 
      PointAddition[RandomPointGroundField, 
        DistortionMap[RandomPointGroundField]]
R = PointMultiplication[arb, 11]
TatePairing[P, R, 11]
TatePairingUnique[P, R, 11]

TatePairing[P, PointAddition[Q, R], 11]
TatePairingUnique[P, PointAddition[Q, R], 11]