\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ 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 <stange at math dot brown dot edu> 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]);

}