clear;

//Doubling and line function 
function Dbl(X1,Y1,Z1,xQ,yQ,u,v,w,z)

	xQ:=3*xQ;

	A:=X1^2; 
	B:=Y1^2; 
	C:=Z1^2;  
	D:=2*C;
	E:=2*D;
	D:=2*E;
	D:=D+E;		
	D:=D*v;        
	E:=X1+Y1;          		
	E:=E^2;
	E:=E-A;
	E:=E-B;   
	F:=Y1+Z1;        		
	F:=F^2;
	F:=F-B;
	F:=F-C; 

	L1:=D-B;	
	L2:=xQ*A;
	L3:=yQ*F;
	L3:=-L3;

	l := -L2*z^2+ L1 +L3*(u*w);
	 
	G:=2*D;
	G:=G+D;
	X3:=B-G;
	X3:=E*X3;
	E:=2*X3;
	X3:=E+X3;      
	E:=G^2;
	E:=2*E;
	E:=2*E;   
	Y3:=B+G;
	Y3:=Y3^2;
	C:=2*Y3;
	Y3:=C+Y3;     		
	Y3:=Y3-E;   
	Z3:=B*F;
	Z3:=2*Z3;
	E:=2*Z3;
	Z3:=2*E;
	Z3:=Z3+E;   		                		

	return X3, Y3, Z3, l;

end function;

function Add(X1, Y1, Z1, X2, Y2, xQ, yQ,u,v,w,z)	

	A:=Z1*X2; 
	A:=X1-A; 
	B:=Z1*Y2; 
	B:=Y1-B; 
	C:=X2*B;
	D:=A*Y2;

	L1:=C-D;
	L2:=xQ*B;
	L2:=-L2;
	L3:=yQ*A;	

	l:=L1+L3*(u*w)+L2*(-z^2); 
	
	C:=A^2;  
	X3:=C*X1; 
	C:=A*C; 
	D:=B^2; 
	D:=D*Z1; 
	D:=C+D; 
	D:=D-X3; 
	D:=D-X3; 
	X3:=X3-D; 
	B:=B*X3; 
	Y3:=C*Y1;
	Y3:=B-Y3; 	
	X3:=A*D; 
	Z3:=Z1*C; 
	
	return X3, Y3, Z3, l;
	
end function;

function MillerAte(P, Q, r,u,v,w,z)

	x1 := P[1];	y1 := P[2];	z1:=P[3];
	x2 := P[1];	y2 := P[2];	z2:=P[3];
	xQ := Q[1]; 	yQ := Q[2];	

	f1 := 1; 
	B := IntegerToSequence(r, 2);

	for i:=#B-1 to 1 by -1 do
		x1,y1,z1, line := Dbl(x1,y1,z1,xQ,yQ,u,v,w,z);
		f1 := f1^2*line;
		if B[i] eq 1 then
			x1,y1,z1, line := Add(x1, y1,z1, x2, y2, xQ, yQ,u,v,w,z);
			f1 := f1*line;
		end if;
	end for;

	return f1;

end function;

unTwist := function(XT, YT,z)
	Xu :=  XT/z^2;
	Yu := YT/z^3;
	return Xu, Yu;
end function;

Twist := function(XT, YT, z)
	Xt := z^2*XT;
	Yt := z^3*YT;
	return Xt,Yt;
end function;

k:=24;
x0:= 9223372036854782449;
q:=(x0-1)^2*(x0^8-x0^4+1) div 3 + x0;
r:=(x0^8-x0^4+1);

Fq:=GF(q);
Fq2<u>:=ExtensionField<Fq,x|x^2+5>;
Fq4<v>:=ExtensionField<Fq2,x|x^2-u>;
Fq12<w>:=ExtensionField<Fq4,x|x^2-v>;
Fq24<z>:=ExtensionField<Fq12,x|x^3-w>;

b:= 1;
bt:= v;

E:=EllipticCurve([Fq|0,b]);
t:=TraceOfFrobenius(E);
T:=t-1;
Et:=EllipticCurve([Fq4|0,bt]);

P:=(#E div r)*Random(E);
Qt:=(#Et div r)*Random(Et);

xT,yT:=Twist(P[1],P[2],z);
Pt:=Et(Fq24)![xT,yT];

x,y:=unTwist(Qt[1],Qt[2],z);
Q:=E(Fq24)![x,y];

fact1:=Random(1,r);
fact2:=Random(1,r);

function HardPart(m,x,Fq);

m1:=m^x;
m2:=m1^x;
m1:=m1^2;
m1:=1/m1;
mu7:=m2*m1;
mu7:=mu7*m;
mu6:=mu7^x;
mu5:=mu6^x;
mu4:=mu5^x;
mu7i:=1/mu7;
mu3:=mu4^x;
mu3:=mu3*mu7i;
mu2:=mu3^x;
mu1:=mu2^x;
mu0:=mu1^x;
msqr:=m^2;
mu0:=mu0*msqr;
mu0:=mu0*m;

mu1:=Frobenius(mu1,Fq,1);
mu2:=Frobenius(mu2,Fq,2);
mu3:=Frobenius(mu3,Fq,3);
mu4:=Frobenius(mu4,Fq,4);
mu5:=Frobenius(mu5,Fq,5);
mu6:=Frobenius(mu6,Fq,6);
mu7:=Frobenius(mu7,Fq,7);

f:=mu0*mu1;
f:=f*mu2;
f:=f*mu3;
f:=f*mu4;
f:=f*mu5;
f:=f*mu6;
f:=f*mu7;

	return f;

end function;

//All on twist

f1:=MillerAte(fact1*Qt,fact2*Pt,T,u,v,w,z);
f2:=MillerAte(fact2*Qt,fact1*Pt,T,u,v,w,z);

F1:=1/f1;
f1:=Frobenius(f1,Fq,12);
f1:=f1*F1;
F1:=Frobenius(f1,Fq,4);
f1:=f1*F1;

F2:=1/f2;
f2:=Frobenius(f2,Fq,12);
f2:=f2*F2;
F2:=Frobenius(f2,Fq,4);
f2:=f2*F2;

HardPart(f1,x0,Fq)-HardPart(f2,x0,Fq);


