//****************************************************************************
//	Program to compute the tau-adic expansion of integers for a fixed 
//	curve and to get statistical results on the timings. 
//****************************************************************************

//****************************************************************************
//	 Introduction
//****************************************************************************


// **** Display text

Help:=procedure();
print"This package demonstrates the computation of m-folds of a divisor class 
using the Frobenius endomorphism of the curve. We have chosen the binary curve
 of genus 2 
	C:y^2 + (x^2 + x + 1)*y = x^5 + x^4 + 1
and consider it over \F_{2^{89}}. The characteristic polynomial of the 
Frobenius endomorphism is 
	T^4-2T^3+3T^2-4T+4.
The class number over this field is
	2*191561942608242456073498418252108663615312031512914969.
The divisor class D chosen has full order.

to display D enter
	> D;

to get the tau-adic expansion of m enter
	> Tau(m);

to compute mD using this expansion enter
	> FrobMult(m);

to compute X times an m-fold of D using the tau-adic expansion of m to get 
the timing enter
	> Stat(X);

for comparison you can enter
	> Compare(X);
to compute X times an m-fold of D using the built-in routine of Magma.

To redisplay this introduction enter
	> Help();


If you want to investigate your own chosen curve load
	> load FrobSelf;
to load a more general program.
";
end procedure;

//******************************************************************************
//	Subroutines
//******************************************************************************


// **** Subprogram to compute the nearest integer, in case of ambiguity the
//	integer of least absolute value is chosen.


round:=function(wert);
	if wert ge 0 then wert:=Ceiling(wert-0.5);
	else wert:=Floor(wert+0.5);
	end if;
	return wert;
end function;

// **** subprogram to compute the tau-adic expansion 



Expand:=function(R2)
	P:=[4,-4,3,-2,1];
	U:=[IntegerRing()|];
	j:=0;	
	while exists(i){i:i in [1..#R2]|R2[i] ne 0} do
		j:=j+1;
		if (R2[1] mod 4 eq 0) then U[j]:=0;
			elif (R2[1] eq -2) then U[j]:=R2[1];
			elif ((R2[1] mod 4) gt 2) then
				U[j]:=(R2[1] mod 4)-4;
			else U[j]:=(R2[1] mod 4);
		end if;
		dummy:=(R2[1]-U[j]) div 4;
		for i in [1..4-1] do
			R2[i]:=R2[i+1]-dummy*P[i+1];
		end for;
		R2[4]:=-dummy;
	end while; 
	return U;
end function;


// **** subprogram to compute sigma(F), where sigma denotes the  Frobenius-
//	endomorphism. If the field is represented via a normal basis, this 
//	only means cyclic shifting. Unfortunately Magma does not support
//	such a representation. Hence, we compute the q-th power of each
//	coefficient.


Frob:=function(poly)
	a:=Coefficients(poly);
	for i in [1..#a] do
		a[i]:=a[i]^2;
	end for;
	return(a);
end function;

// **** subprogram to find random elements on Jacobian - no internal 
//	procedure in Magma2.8

Take_even:=function(J,h,f,K,KR,X)
	repeat 
		a:=Random(K);
		A:=Evaluate(f,a)/(Evaluate(h,a)^2);
	until IsZero(Trace(A));
	b:=(&+[A^(2^(2*j)): j in [0..44]]);
	b:=b*Evaluate(h,a);
	D:=J![X-a,b];
	m:=Random(2^89);
	D:=m*D;
	return(D);
end function;

//**************************************************************************
//	MAIN
//**************************************************************************
Help();

// **** General definitions

S<T>:=PolynomialRing(RingOfIntegers());
r<t>:=PolynomialRing(RealField());
R<x>:=PolynomialRing(GF(2));
h:=x^2+x+1;
f:=x^5+x^4+1;
C:=HyperellipticCurve(h,f);
K<xi>:=GF(2,89);
J:=BaseExtend(Jacobian(C),K);
KR<X>:=PolynomialRing(K);
//D:=Random(J); //does not work in Magma 2.8
D:=Take_even(J,h,f,K,KR,X);
//while IsOrder(D,2)  do //does not work in Magma 2.8
while IsZero(2*D) do
//	D:=Random(J);
D:=Take_even(J,h,f,K,X);
end while;
D:=2*D;
// print"We use D:=",D,".";
E:=J!0;
// **** Vector, containing the coefficients of the characteristic polynomial 
//	of the Frobenius endomorphism.

P:=[4,-4,3,-2,1];    


// **** precomputations to be used in tau-adic method

Y:=[J|i*D: i in [1..2]];

// **** precomputations to be used in tau-adic expansion

d:=[IntegerRing()|];e:=[IntegerRing()|];
d[1]:=1;e[1]:=1;
for s in [2..4] do d[s]:=0; e[s]:=0; end for;

for i in [1..88] do
	dalt:=d[4];
	for j:=4 to 2 by (-1) do
		d[j]:=d[j-1]-P[j]*dalt;
		e[j]:=e[j]+d[j];
	end for;
	d[1]:=-P[1]*dalt;
	e[1]:=e[1]+d[1];
end for;
polyr:=r!P;
polyd:=r!e; 
ggT,polyM,dummy:=XGCD(polyd,polyr);
M:=Coefficients(polyM);
  for l in [Degree(polyM)+2..4] do M[l]:=0; end for;

// ************************************************************************
// 	definition of procedures
// ************************************************************************


// ************************************************************************
//	Program to compute the tau-adic expansion 
// ************************************************************************


Tau:=function(m)
	M1:=[IntegerRing()|];	
	for i in [1..4] do M1[i]:=round(m*M[i]); end for;
	polyR1:=S!(m-((S!e*S!M1) mod S!P));
	R1:=[IntegerRing()|];
	for l in [1..Degree(polyR1)+1] do 
		R1[l]:=(Coefficients(polyR1))[l]; 
	end for;
	for l in [Degree(polyR1)+2..4] do 
		R1[l]:=0; 
	end for;
	U:=Expand(R1);	
	print"The tau-adic expansion of ",m, " is given by the following 
vector of length",#U,"(starting with the lowest power of tau):";
	return (U);
end function;

// ************************************************************************
//	Program to compute mD using the tau-adic expansion
// ************************************************************************

FrobMult:=procedure(m)

// **** Computation of the expansion

	M1:=[IntegerRing()|];	
	for i in [1..4] do M1[i]:=round(m*M[i]); end for;
	polyR1:=S!(m-((S!e*S!M1) mod S!P));
	R1:=[IntegerRing()|];
	for l in [1..Degree(polyR1)+1] do 
		R1[l]:=(Coefficients(polyR1))[l]; 
	end for;
	for l in [Degree(polyR1)+2..4] do 
		R1[l]:=0; 
	end for;
	U:=Expand(R1);			
	print"The tau-adic expansion of ",m, " is given by the following vector of length",#U," (starting with lowest power):";	
U;

// **** Initialization

	if #U eq 0 then H:=E; else

	H:=Sign(U[#U])*Y[Abs(U[#U])];

// **** Computation of mD

	if #U gt 1 then 
		for i:=(#U-1) to 1 by (-1) do
			H:=J![KR!Frob(H[1]),KR!Frob(H[2])];
			if U[i] ne 0 then 
				H:=H+(Sign(U[i])*Y[Abs(U[i])]);
			end if;
		end for;
	end if;
	end if;
	print" ";
	print" ",m,"D equals  "; 
	H;
end procedure;

// ************************************************************************
//	Program to compute X times mD using the tau-adic expansion of m
//	where m is chosen randomly leq 
//	191561942608242456073498418252108663615312031512914969
// ************************************************************************


Stat:=procedure(anz);

// **** Initialization 

expantime:=0;
multtime:=0;
length:=0;

// **** Loop to do the wanted number of multiplications.



for i in [1..anz] do	
	m:=Random(191561942608242456073498418252108663615312031512914969);
	
// **** First we compute the tau-adic expansion

	t1:=Cputime();
	M1:=[IntegerRing()|];	
	for i in [1..4] do M1[i]:=round(m*M[i]); end for;
	polyR1:=S!(m-((S!e*S!M1) mod S!P));
	R1:=[IntegerRing()|];
	for l in [1..Degree(polyR1)+1] do 
		R1[l]:=(Coefficients(polyR1))[l]; 
	end for;
	for l in [Degree(polyR1)+2..4] do 
		R1[l]:=0; 
	end for;
	U:=Expand(R1);

	expantime:=expantime+Cputime(t1);
	length:=length+#U;	

// **** Now we use the expansion to compute m*D.

	t2:=Cputime();
// **** Initialization
	if #U eq 0 then H:=E; else
	H:=Sign(U[#U])*Y[Abs(U[#U])];

// **** Computation of mD.

	if #U gt 1 then 
		for i:=(#U-1) to 1 by (-1) do
			H:=J![KR!Frob(H[1]),KR!Frob(H[2])];
			if U[i] ne 0 then 
				H:=H+(Sign(U[i])*Y[Abs(U[i])]);
			end if;
		end for;
	end if;
	end if;
	multtime:=multtime+Cputime(t2);
end for;

// **** Display results

print"The average length of the expansions is",length/anz,".";
print"To compute the expansion we needed on average ",expantime/anz,".";
print"To perform the multiplication we needed additionally",multtime/anz,
"on average.";
print"Hence, on average we needed ",(expantime+multtime)/anz," to compute 
an m-fold.";

end procedure;

//***************************************************************************
// 	Program to compute X time an m-fold using the built-in routine of 
//	Magma where again m is chosen randomly leq 
//	191561942608242456073498418252108663615312031512914969.
// ************************************************************************



Compare:=procedure(anz)


// **** Initialization

totaltime:=0;


for i in [1..anz] do	
	m:=Random(191561942608242456073498418252108663615312031512914969);
	t1:=Cputime();
	H:=m*D;
	totaltime:=totaltime+Cputime(t1);
end for;

// **** Display results.

print"To perform one multiplication we needed on average ",totaltime/anz,".";

end procedure;