clear;

x:=-1;
q:=36*x^4+36*x^3+24*x^2+6*x+1;
r:=36*x^4+36*x^3+18*x^2+6*x+1;

while not IsPrime(q) or not IsPrime(r) do
	x:=x+1;
	q:=36*x^4+36*x^3+24*x^2+6*x+1;
	r:=36*x^4+36*x^3+18*x^2+6*x+1;
end while;

Fq:=GF(q);
b:=1;
E:=EllipticCurve([Fq|0,b]);
while #E mod r ne 0 do
	b+:=1;
	E:=EllipticCurve([Fq|0,b]);
end while;
	
#E mod r;
O:=PointsAtInfinity(E)[1];

p1:=AllRoots(Fq!1,3)[2];
p2:=AllRoots(Fq!1,3)[3];

phi1:=function(P);
	return E![P[1]*p1,P[2]];
end function;

phi2:=function(P);
	return E![P[1]*p2,P[2]];
end function;


for i:=1 to 1 do
	P:=Random(E);
	if P eq O then
		P:=Random(E);
	end if;
	Log(P,phi1(P)),Log(P,phi2(P));
	lambda1:=Log(P,phi1(P)); lambda2:=Log(P,phi2(P));
end for;

(lambda1^2+lambda1+1) mod r;
(lambda2^2+lambda2+1) mod r;

IntegerToSequence(lambda1,2);
IntegerToSequence(lambda2,2);
x;

