\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ PARI/GP Script for Computation of Tate Pairing v. 1.1 \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ This script uses the algorithm presented in \\
\\ "The Tate Pairing via Elliptic Nets" by Katherine E. Stange \\
\\ See http://www.math.brown.edu/~stange/tatepairing/ \\
\\ or contact for info \\
\\ \\
\\ Note that the optimisations mentioned in the paper are not \\
\\ all implemented here; rather, the algorithm is implemented in \\
\\ its simplest form for clarity. \\
\\ \\
\\ v. 1.1 fixes a bug that affects pairings on points whose order is \\
\\ a power of two. \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Set debug = 1 to report steps of algorithm
global(debug);
debug = 0;
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ double_or_add
\\ Given a block V centred at k, and the initial data relevant to
\\ the elliptic net, returns either a block centred at 2k or 2k+1
\\ depending on whether variable "add" is 0 or 1 respectively.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
{
double_or_add(V, initial_data, add) =
local(doubleV, m, i, j);
\\ Create the output block
doubleV = matrix(2,8);
\\ initial_data contains the precomputed inverses
inverse_20 = initial_data[1]; \\ inverse of W(2,0)
inverse_11 = initial_data[2]; \\ inverse of W(1,1)
inverse_n1 = initial_data[3]; \\ inverse of W(-1,1)
inverse_2n = initial_data[4]; \\ inverse of W(2,-1)
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Fill out first vector of output block
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
for(j=-1,2,
i = j;
m = 4; \\ index to middle of block
doubleV[1,m + 2*i - add] = ((V[1,m+i])*(V[1,m+i+2])
*(V[1,m+i-1])^2 - (V[1,m+i])
*(V[1,m+i-2])*(V[1,m+i+1])^2)
*inverse_20;
\\ when we hit j=-1, if add=1, calculate
\\ W(2k+5,0) instead of W(2k-3,0)
if( i == -1 && add == 1, i = 3; );
doubleV[1,m + 2*i - 1 - add] = (V[1,m+i+1])
*(V[1,m+i-1])^3 - (V[1,m+i-2])
*(V[1,m+i])^3;
);
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\Fill out second vector of output block
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
m2 = 2; \\ index to middle of second vector of block
m1 = 4; \\ index to middle of first vector of block
if( add == 0,
doubleV[2,1] = ( V[2,m2+1]*V[2,m2-1]*V[1,m1-1]^2
- V[1,m1]*V[1,m1-2]*V[2,m2]^2 )*inverse_11;
);
doubleV[2,2-add] = ( V[2,m2-1]*V[2,m2+1]*V[1,m1]^2
- V[1,m1-1]*V[1,m1+1]*V[2,m2]^2);
doubleV[2,3-add] = ( V[2,m2+1]*V[2,m2-1]*V[1,m1+1]^2
- V[1,m1]*V[1,m1+2]*V[2,m2]^2)*inverse_n1;
if( add == 1,
doubleV[2,3] = ( V[1,m1+1]*V[1,m1+3]*V[2,m2]^2
- V[2,m2-1]*V[2,m2+1]*V[1,m1+2]^2 )*inverse_2n;
);
return(doubleV);
}
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ net_loop
\\ Given a starting block V centred at 1, initial data relevant to
\\ the elliptic net, and an integer m, returns the block centred at m
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
{
net_loop(V,initial_data, m) =
local(currentV, m_size, add, i, j);
\\ determine the number of steps in the double-and-add loop
m_size = ceil(log(m+1)/log(2));
\\ the variable storing the current block
currentV = V;
\\ ignore the first "1" in the binary expansion of m
m = m - 2^(m_size-1);
\\ step through the digits in the binary expansion of m
for(j=1,m_size-1,
i = m_size - j; \\kludgy version of "down to"
\\ determine if this is a double step or a
\\ double-and-add step based on current digit
\\ of m; set "add" accordingly
if( m - 2^(i-1) >= 0,
add = 1;
m = m - 2^(i-1);
,
add = 0;
);
\\ call the double or double-and-add function to
\\ update the current block
currentV = double_or_add(currentV, initial_data, add);
\\ print information about the current step
if( debug == 1,
print("The digit of $m$ for this step is ",
add, ".");
print("The resulting block from this step is ",
currentV);
);
);
return(currentV);
}
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ tate_pairing_alg
\\ Arguments: elliptic curve, two points, and integer.
\\ Returns: tate pairing of the two points with respect to the integer.
\\ Note: this is currently not implemented for curves in characteristic
\\ 2 or 3.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
{
tate_pairing_alg(elliptic_curve, point_a, point_b, torsion) =
local(V, initial_data, x1, y1, x2, y2, a4, a6, resultV);
\\ Create a starting block for the net
V = matrix(2,8);
\\ Create a vector to store the pre-computed inverses
initial_data = vector(4);
\\ Make sure the points are on the curve
if( ellisoncurve(elliptic_curve, point_a) == 0,
print("The first point is not on the curve!");
);
if( ellisoncurve(elliptic_curve, point_b) == 0,
print("The second point is not on the curve!");
);
\\ If the curve is not in nice Weierstrass form
\\ y^2 = x^3 + Ax + B, do a change of coordinates
if( elliptic_curve[1] != 0 || elliptic_curve[2] != 0
|| elliptic_curve[3] != 0,
print("Curve not in two-coeff. Weierstrass form!");
a1 = elliptic_curve[1];
a2 = elliptic_curve[3];
a3 = elliptic_curve[2];
a4 = elliptic_curve[4];
a6 = elliptic_curve[5];
coordinate_change = [1/6,-a2/3,-a1/2,-a3/2+a1*a2/6];
elliptic_curve = ellchangecurve(elliptic_curve,
coordinate_change);
point_a = ellchangepoint(point_a, coordinate_change);
point_b = ellchangepoint(point_b, coordinate_change);
print("New curve: ", elliptic_curve);
print("New first point: ", point_a);
print("New second point: ", point_b);
);
\\ Set up the usual variable names for elliptic curves
x1 = point_a[1];
y1 = point_a[2];
x2 = point_b[1];
y2 = point_b[2];
a4 = elliptic_curve[4];
a6 = elliptic_curve[5];
\\ Fill out the first vector of the block
V[1,4] = 1;
V[1,5] = 2*y1;
V[1,6] = 3*x1^4 + 6*a4*x1^2 + 12*a6*x1-a4^2;
V[1,7] = 4*y1*(x1^6 + 5*a4*x1^4 + 20*a6*x1^3 - 5*a4^2*x1^2
- 4*a4*a6*x1 - 8*a6^2 - a4^3);
V[1,8] = -((V[1,6])^3 - (V[1,5])^3*(V[1,7]));
V[1,3] = 0;
V[1,2] = - V[1,4];
V[1,1] = - V[1,5];
\\ Fill out the second vector of the block
V[2,1] = 1;
V[2,2] = 1;
V[2,3] = 2*x1+x2 - ((y2-y1)/(x2-x1))^2;
\\ Print starting block
if( debug == 1,
print("The beginning block is: ", V);
);
\\ Pre-compute the inverses
initial_data[1] = 1/(V[1,5]);
initial_data[2] = 1;
initial_data[3] = 1/(x1 - x2);
initial_data[4] = 1/((y1+y2)^2 - (2*x1 + x2)*(x1-x2)^2);
\\ Print the pre-computed info
if( debug == 1,
print("Precomputed inverses: ", initial_data);
);
\\ Call the net algorithm to obtain the block centred
\\ at "torsion"
resultV = net_loop(V, initial_data, torsion);
\\ Apply the Tate pairing formula and return result
return(resultV[2,3]/resultV[1,5]);
}