clear;

q:=47;
Fq:=GF(q);
a:=21; b:=15;
E:=EllipticCurve([Fq|a,b]);
r:=17;
k:=4;
Fq4<u>:=ExtensionField<Fq,x|x^4-4*x^2+5>;

P:=E![30,3]; //P:=Random(E); 
P:=(Order(P) div r)*P;
E:=BaseChange(E,Fq4);
P:=E!P;
Q:=E![5*u^3 + 37*u + 13,7*u^3 + 45*u^2 + 10*u + 7];//Q:=Random(E);
h:=#E div r^2;
Q:=h*Q;
pi:=FrobeniusMap(E);
O:=PointsAtInfinity(E)[1];

TraceMap:=function(Q,pi);
	return pi(pi(pi(Q)))+pi(pi(Q))+pi(Q)+Q;
end function;

Q:=k*Q-TraceMap(Q,pi); // Anti-trace map
pi(Q)-q*Q; //check

F<x,y>:=FunctionField(E);
fDBL:=function(P);
	lambda:=(3*P[1]^2+a)/(2*P[2]);
	c:=P[2]-lambda*P[1];
	l:=F!(y-(lambda*x+c));
	v:=F!(x-(lambda^2-2*P[1]));
	return l,v;
end function;

fADD:=function(P,Q);
	lambda:=(Q[2]-P[2])/(Q[1]-P[1]);
	c:=P[2]-lambda*P[1];
	l:=F!(y-(lambda*x+c));
	v:=F!(x-(lambda^2-P[1]-P[2]));
	return l,v;
end function;

Miller:=function(P,Q,r);
	T:=P; f:=1; 
	rbits:=IntegerToSequence(r,2);
	for i:=#rbits-1 to 2 by -1 do
		l,v:=fDBL(T);
		T:=2*T;
		f:=f^2*l(Q)/v(Q);
		"DBL";
		f;
		if rbits[i] eq 1 then
			l,v:=fADD(T,P);
			T:=T+P;
			f:=f*l(Q)/v(Q);
			"ADD";
			f;
		end if;
	end for;
	l,v:=fDBL(T);
	T:=2*T;
	f:=f^2*l(Q)/v(Q);	
	"DBL";
	f;
	l:=x-P[1];
	f:=f*l(Q);
	"ADD";
	l(Q);
	return f;
end function;

Tate:=function(P,Q,r);
	return Miller(P,Q,r)^((q^k-1) div r);
end function;

Tate(P,Q,r);

