N:=50000;
p:=N;
repeat p:=NextPrime(p);
until ((p mod 4) eq 1) and IsPrime(ExactQuotient(p-1,4));
p;

function half_gcd(a,b)
	v0:=0; v1:=1;
	while b gt Abs(v1) do
		q,r:=Quotrem(a,b);
		a:=b;b:=r;
		v2:=v0-q*v1;
		v0:=v1;v1:=v2;
	end while;
	return b,v1;
end function;

F:=GF(p);
m:=SquareRoot(F!-1);
mm:=Integers()!m;
c,d:=half_gcd(p,mm);
m-F!(c/d);

A<i>:=GaussianIntegers();

L:=[];
B1:=17;B2:=17;S:=40;
//FB:=[2,3,5,7,11,13,1+i,2-i,2+i,3-2*i,3+2*i,4-i,4+i,i];
FB1:=[0];//pour d
FB2:=[i];
M1:=ZeroMatrix(Rationals(),0,#FB1);
M2:=ZeroMatrix(Rationals(),0,#FB2);

for a:=-S to S do
	for b:=0 to S do
		if GCD(a,b) eq 1 then
			F1:=Factorisation(d*a+c*b);
			F2:=Factorisation(a^2+b^2);
			if ((#F1 eq 0) or (F1[#F1][1] le B1)) and ((#F2 eq 0) or(F2[#F2][1] le B2)) then 
				Append(~L,[a,b]);
				//côté linéaire
				for fact in F1 do
					j:=Index(FB1,fact[1]);
					if j eq 0 then //on ajoute à la base
						Append(~FB1,fact[1]);
					end if;
				end for;
				if NumberOfColumns(M1) lt #FB1 then
					M1:=HorizontalJoin(M1,ZeroMatrix(Rationals(),NumberOfRows(M1),#FB1-NumberOfColumns(M1)));
				end if;				
				N1:=ZeroMatrix(Rationals(),1,#FB1);	
				N1[1,1]:=-1; //correspond à d
				for fact in F1 do
					N1[1,Index(FB1,fact[1])]:=fact[2];
				end for;
				M1:=VerticalJoin(M1,N1);
				//côté gaussien
				F3,ipow:=Factorisation(a+b*i);
				for fact in F3 do
					j:=Index(FB2,fact[1]);
					if j eq 0 then //on ajoute à la base
						Append(~FB2,fact[1]);
					end if;
				end for;
				if NumberOfColumns(M2) lt #FB2 then
					M2:=HorizontalJoin(M2,ZeroMatrix(Rationals(),NumberOfRows(M2),#FB2-NumberOfColumns(M2)));
				end if;				
				N2:=ZeroMatrix(Rationals(),1,#FB2);	
				N2[1,1]:=ipow[2]; //correspond à i
				if d*a+c*b lt 0 then N2[1,1]+:=2; end if; //-1 = i^2 
				for fact in F3 do
					N2[1,Index(FB2,fact[1])]:=fact[2];
				end for;
				M2:=VerticalJoin(M2,N2);
			end if;
		end if;
	end for;
end for;
#L;#FB1+#FB2;
FB1;FB2;
//algèbre linéaire
MM:=ChangeRing(HorizontalJoin(M1,M2),ResidueClassRing(ExactQuotient(p-1,4)));
K:=NullspaceOfTranspose(MM);
Dimension(K) eq 1;//sinon il faut continuer
Vlogs:=Basis(K)[1];

//descente

//challenge
/*repeat 
	g:=Random(F);
until g^4 ne 1 and g^(ExactQuotient(p-1,2)) ne 1;
s:=Random(p-1);
h:=g^s;*/
g:=F!3295;h:=F!5460;

//première relation
k1:=-1;
repeat 
k1+:=1;
n:=Integers()!(g*h^k1);
if n gt p/2 then n:=n-p; end if;
n1,n2:=half_gcd(p,n);
Fn1:=Factorisation(n1);
Fn2:=Factorisation(n2);
until ((#Fn1 eq 0) or (Fn1[#Fn1][1] le B1)) and ((#Fn2 eq 0) or(Fn2[#Fn2][1] le B1));
//1+k1*x = y*sum a_i Vlogs[i]=y*l1
l1:=0;
for fact in Fn1 do l1+:=fact[2]*Vlogs[Index(FB1,fact[1])]; end for;
for fact in Fn2 do l1-:=fact[2]*Vlogs[Index(FB1,fact[1])]; end for;

k2:=k1;
repeat 
k2+:=1;
n:=Integers()!(g*h^k2);
if n gt p/2 then n:=n-p; end if;
n1,n2:=half_gcd(p,n);
Fn1:=Factorisation(n1);
Fn2:=Factorisation(n2);
until ((#Fn1 eq 0) or (Fn1[#Fn1][1] le B1)) and ((#Fn2 eq 0) or(Fn2[#Fn2][1] le B1));
//1+k2*x = y*sum a'_i Vlogs[i]=y*l2
l2:=0;
for fact in Fn1 do l2+:=fact[2]*Vlogs[Index(FB1,fact[1])]; end for;
for fact in Fn2 do l2-:=fact[2]*Vlogs[Index(FB1,fact[1])]; end for;

k1*l2-k2*l1 ne 0;//sinon il faut continuer
x:=(l2-l1)/(k2*l1-k1*l2);

(g^4)^(Integers()!x) eq h^4; 





