//****************************************************************************
//	Program to compute the tau-adic expansion of integers for an arbitrary 
//	curve. The characteristic, the degree of extension and the parameters
//	of the curve are defined by the user. The program chooses a random 
//	divisor D. Upon entry of an integer m the program computes the 
// 	tau-adic expansion of m and uses it to compute mD.
//****************************************************************************


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

print"This package demonstrates the computation of m-folds of a divisor class 
using the Frobenius endomorphism of the curve. The curve should be defined 
over a comparably small ground field (max. 5 elements). It is then considered 
over a large extension field. You can find appropriate curves together with 
their class numbers at
http://www.exp-math.uni-essen.de/~lange/KoblitzC.html.
In this program you can define yourself a curve over a PRIME field. If you'd 
better like to run a prepared example using the curve
	C:y^2 + (x^2 + x + 1)*y = x^5 + x^4 + 1 of genus 2 over \F_{2^{89}}) 
then break and load FrobExample.";

S<T>:=PolynomialRing(RingOfIntegers());
r<t>:=PolynomialRing(RealField());
//******************************************************************************
//	subroutines
//******************************************************************************

// **** subprogram to count the number of points for odd characteristic ****

NumOdd:=function(p,j,f)
	k:=0;
	if j eq 1 then return #{<x,y>:x,y in GF(p,j)|y^2 eq Evaluate(f,x)}; 
	else
	   for x in GF(p,j) do
	      if Evaluate(f,x)eq 0 then k:=k+1;
		elif Evaluate(f,x)^((p^j-1)div 2) eq GF(p,j)!1 then k:=k+2;
	      end if;
	   end for;	
	   return k;
	end if;
end function;

// **** subprogram to count the number of points for even characteristic ****

NumEven:=function(p,j,f,h)
	if j eq 1 then return #{<x,y>:x,y in GF(p,j)|y^2+Evaluate(h,x)*y eq 
	Evaluate(f,x)}; 
	else	
	   k:=0;
	   for x in GF(p,j) do
	       if Evaluate(h,x) eq 0 then k:=k+1;
		  elif Trace((Evaluate(f,x)/(Evaluate(h,x)^2))) eq 0 
		  then k:=k+2;
	       end if;
	   end for;
	   return k;
	end if;
end function;

// **** 
//	Program to count points on the hyperelliptic curve defined by 
//	C: y^2=f(x) for odd characteristic q and genus g
// 	returns the characteristic polynomial of the Frobenius endomorphism.
// ****

CountpointsOdd:=function(q,g,f)
a:=[IntegerRing()|];
M:=[IntegerRing()|];
R1:=[IntegerRing()|];
for i in [1..g] do
        M[i]:=NumOdd(q,i,f)-q^i;
	end for;
a[1]:=1;
a[2*g+1]:=q^g;
for i in [2..g] do
  	a[i]:=(&+[M[j]*a[i-j]: j in [1..i-1]]) div (i-1);
	a[2*g+2-i]:=q^(g+1-i)*a[i];
	end for;
a[g+1]:=(&+[M[j]*a[g+1-j]: j in [1..g]]) div g;
P:=S! Reverse(a);
if (IsIrreducible(P) eq false) then  
	print"The characteristic polynomial is not irreducible. Such curves 
are not recommended since the class number contains many prime factors.";
// **** Test on supersingularity
else
	R1:=Coefficients(P);
	if forall{l:l in [g+1..2*g]|(R1[l] mod q^(Ceiling((2*g-l+1)/2))) eq 0}
	  then print "The Jacobian variety of this curve is supersingular. 
We do not recommend this curve.";
		end if;
	end if;
return P;
end function;

// **** 
//	Program to count points on the hyperelliptic curve defined by 
//	C: y^2+h(x)y=f(x) for even characteristic q and genus g
// 	returns the characteristic polynomial of the Frobenius endomorphism.
// ****

CountpointsEven:=function(q,g,h,f);
a:=[IntegerRing()|];
M:=[IntegerRing()|];
R1:=[IntegerRing()|];
for i in [1..g] do
	M[i]:=NumEven(q,i,f,h)-q^i;
	end for;
a[1]:=1;
a[2*g+1]:=q^g;	
for i in [2..g] do
	a[i]:=(&+[M[j]*a[i-j]: j in [1..i-1]]) div (i-1);
	a[2*g+2-i]:=q^(g+1-i)*a[i];
	end for;
a[g+1]:=(&+[M[j]*a[g+1-j]: j in [1..g]]) div g;
P:=S! Reverse(a);
if (IsIrreducible(P) eq false) then 
	print"The characteristic polynomial is not irreducible. Such curves 
are not recommended since the class number contains many prime factors.";
// **** Test on supersingularity 
else
	R1:=Coefficients(P);
	if forall{l:l in [g+1..2*g]|(R1[l] mod q^(Ceiling((2*g-l+1)/2))) eq 0}
	   then print "The Jacobian variety of this curve is supersingular. 
We do not recommend this curve.";
		end if;
	end if;
return P;
end function;

// **** 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 for even characteristic


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

// **** subprogram to compute the tau-adic expansion for odd characteristic

ExpandOdd:=function(R2,P,G,Q)
	U:=[IntegerRing()|];
	j:=0;	
	while exists{i:i in [1..#R2]|R2[i] ne 0} do
		j:=j+1;
		if (R2[1] mod Q^G eq 0) then U[j]:=0;
			elif ((R2[1] mod Q^G) gt (Q^G-1)/2) then
				U[j]:=(R2[1] mod Q^G)-Q^G;
			else U[j]:=(R2[1] mod Q^G);
		end if;
		dummy:=(R2[1]-U[j]) div (Q^G);
		for i in [1..2*G-1] do
			R2[i]:=R2[i+1]-dummy*P[i+1];
		end for;
		R2[2*G]:=-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,p)
	a:=Coefficients(poly);
	for i in [1..#a] do
		a[i]:=a[i]^p;
	end for;
	return(a);
end function;

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


Take_odd:=function(J,f,K,KR,X)
	a:=Random(K);
	repeat 
		a:=Random(K);
	until IsSquare(Evaluate(f,a));
	b:=Sqrt(Evaluate(f,a));
	D:=J![X-KR!a,b];
	m:=Random(#K);
	D:=m*D;
	return(D);
end function;

Take_even:=function(J,h,f,K,KR,X,n)//we assume n to be odd
	repeat 
		a:=Random(K);
		A:=Evaluate(f,a)/(Evaluate(h,a)^2);
	until IsZero(Trace(A));
	b:=(&+[A^(2^(2*j)): j in [0..((n-1) div 2)]]);
	b:=b*Evaluate(h,a);
	D:=J![X-KR!a,b];
	m:=Random(#K);
	D:=m*D;
	return(D);
end function;

//**************************************************************************
//	MAIN
//**************************************************************************

// **** get characteristic and build finite field or BREAK
q:=12;
while (q ne 0) and (IsPrime(q) eq false) do
readi q, "Enter the cardinality of the finite field (enter 0 for break)";
end while;
if (q eq 0) then print"End";
else
R<x>:=PolynomialRing(GF(q));

// **** get equation of hyperelliptic curve, for even characteristic we 
//	 need two polynomials

print"Enter the defining equation of the hyperelliptic curve
	C: y^2+h(x)y=f(x),
where h and f are polynomials over GF(",q,"), deg h<= g, deg f=2g+1, where g 
is the genus of the curve. For security reasons we recommend choosing 
g<=4. h is only needed for even characteristic. The polynomials are 
entered coefficient wise.";
readi g,"Enter the genus of the curve.";
g;
if (q mod 2) eq 0 then
	print"Now enter the coefficients of h.";
	h:=[GF(q)|];
	for i in [0..g] do
		print "coefficient of x^",i;
		readi dummy;
		h[i+1]:=GF(q)!dummy;
	end for;
	h:=R!h;
end if;
print"Now enter the coefficients of f.";
	f:=[GF(q)|];
for i in [0..2*g+1] do
	print "coefficient of x^",i;
	readi dummy;
	f[i+1]:=GF(q)!dummy;
	end for;
	f:=R!f;

// **** definition of the hyperelliptic curve C, error and break if C is 
//	singular.

if (q mod 2) eq 0 then
	C:=HyperellipticCurve(h,f);
	else C:=HyperellipticCurve(f);
	end if;
print " The curve is given by C:",C;

// **** Computation of the characteristic polynomial P of the Frobenius 
//	endomorphism

if (q mod 2) eq 0 then
	P:=CountpointsEven(q,g,h,f);
	else P:=CountpointsOdd(q,g,f);
	end if;
print"The characteristic polynomial of the Frobenius endomorphism of your 
curve is ",P;
if Evaluate(P,1) le (q^g-1) then 
	print"Caution! For this curves there might occur periodic expansion. 
If one expansions takes unusually long please break.";
end if;


// **** The Jacobian variety of the curve is defined over an extension field,
//	a random point is chosen.

readi n, "Enter the degree of extension, for p=2 only odd n is allowed (for cryptographic applications n should be prime).";
// I don't want to spend too much time on how to find a point on J
// works without restriction in Magma 2.7 and Magma 2.10
K<xi>:=GF(q,n);
J:=BaseExtend(Jacobian(C),K);
KR<X>:=PolynomialRing(K);
//D:=Random(J); //does not work in Magma 2.8
	if (q mod 2) eq 0 then
		D:=Take_even(J,h,f,K,KR,X,n);
	else
		D:=Take_odd(J,f,K,KR,X);
	end if;
//while IsOrder(D,Evaluate(P,1))  do // does no longer work with 2.8
while IsZero(Evaluate(P,1)*D) do
//	D:=Random(J);
	if (q mod 2) eq 0 then
		D:=Take_even(J,h,f,K,KR);
	else
		D:=Take_odd(J,f,K,KR,X);
	end if;
end while;
D:=Evaluate(P,1)*D;
print"A random point is given by ";
print"D:=",D;


// **** when we reduce the length of the tau-adic representation we assume 
//	that D lies in a subgroup of large prime order. If you decide to 
//	take a curve with nontrivial cofactor, please change the above 
//	line D:=Evaluate(P,1)*D; to D:=D*cofactor;

print"If the Picard group is of almost primer order, i.e. factors as 
|Pic^0(C/F_q)|*large prime, then D is of this prime order. Otherwise 
there might occur some mistakes in the process of reducing the length - if 
so, please read the comments in the program.";

// **** from now on we only need the coefficients of P

P:=Coefficients(P);    
	print"Now we do some precomputations. For large genus and size of the 
ground field this might take a while. Please be patient or stop the program 
by typing 2 times Ctrl+C. (We do ~p^{g}/2, i.e.",Ceiling((q^(g)-1)/2)-1," 
relatively short precomputations.)";

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

t1:=Cputime();
Y:=[J|i*D: i in [1..Ceiling((q^g-1)/2)]];

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

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

for i in [1..n-1] do
	dalt:=d[2*g];
	for j:=2*g 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(r!e,r!P);
M:=Coefficients(polyM);
  for l in [Degree(polyM)+2..4] do M[l]:=0; end for;

print"The precomputations took ",Cputime(t1),".";

// **** loop to compute the expansion of m and mD

m:=1;
while m ne 0 do
readi m, "Please enter the multiplier m to expand (enter 0 for break)";
if m ne 0 then
	t1:=Cputime();
	M1:=[IntegerRing()|];	
	for i in [1..2*g] 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..2*g] do 
		R1[l]:=0; 
	end for;
	if (q mod 2) eq 0 then
		U:=ExpandEven(R1,P,g,q);	
		else U:=ExpandOdd(R1,P,g,q);
	end if;
	T1:=Cputime(t1);
	print"The tau-adic expansion of ",m, " is given by the following 
vector of length",#U,"(starting with the lowest power of tau):";
	U;

// **** initialization for computing mD using the tau-adic expansion.

	t2 := Cputime();
	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],q),KR!Frob(H[2],q)];
			if U[i] ne 0 then 
				H:=H+(Sign(U[i])*Y[Abs(U[i])]);
			end if;
		end for;
	end if;
	T2:=Cputime(t2);
	print" ",m,"D equals"; H;
	print" The computation of the tau-adic expansion took ",T1,"and the 
computation of mD was performed in ",T2," such that the total time was",T1+T2,".";
	print" To compare we now compute ",m,"D using Magma's version. ";
	time m*D;
end if;
end while;
end if;