//Compute a-numbers, bounds, and statistics for covers of an elliptic curve branched at two points

// E: y^2=x^3-x
anumber_E:=function(p)
    if ((p mod 4) eq 1) then
        return 0;
    else return 1;
    end if;
end function;

// z^p - z = f over function field F of the curve
anumber_Y:=function(F,f,p)
        S<z> := PolynomialRing(F);
        FF<beta> := FunctionField(z^p-z - f);
        M:=CartierRepresentation(FF);
        return Dimension(Kernel(M));	
end function;

//n_{Q,i}
compute_n := function(i,d,p);
	return Ceiling( (p-1-i)*d /p);
end function;

//is a + d*k = 1 mod p for k<=i
check_mults := function(a,d,i,p);
    k :=( (1-a)* Modinv(d,p)) mod p ;
    return not (k gt i);
end function;


//lower bounds over E: compute L(\pi,P^1) exactly, two poles
lower_bound := function(d1,d2,p);
    lowerbound:= 0;
    optimum:=0;
    
    for j in [0..p-1] do
        bnd:=0;
                
        aj := compute_n(j,d1,p);
        bj:= compute_n(j,d2,p);
        for i in [0..j] do
            n:=compute_n(i,d1,p);
            m:=compute_n(i,d2,p);
 	       bnd:=bnd+ (n - Ceiling((aj+d1*(j-i))/p))+(m - Ceiling((bj+d2*(j-i))/p));
	end for;
        
        if bnd gt lowerbound then
            lowerbound := bnd;
            optimum := j;   
	end if;
        
    end for;
    return lowerbound;
end function;

// upper bound over E: bound U(E,\pi) via computing #T, one pole, 
upper_bound := function(d1,d2,p);

    count1 := 0; //not in kernel  
    
    for i in [0..p-1] do
        n:= compute_n(i,d1,p);
        for j in [1..n] do 
            if (j mod p ne 1) and check_mults(j,d1,i,p) then 
                count1 := count1+1;
	    end if;	
	end for;
    end for;

     count2 := 0; //not in kernel  
    
    for i in [0..p-1] do
        m:= compute_n(i,d2,p);
        for j in [1..m] do 
            if (j mod p ne 1) and check_mults(j,d2,i,p) then 
                count2 := count2+1;
	    end if;	
	end for;
    end for;	



    dimension := 0;
    for i in [0..p-2] do
        n:=compute_n(i,d1,p);
        m:=compute_n(i,d2,p);
        dimension := dimension + (n - Ceiling(n/p))+(m - Ceiling(m/p));
    end for;
    
    

    upperbound :=  p*anumber_E(p) + dimension - count1-count2;
    return upperbound;
end function;

termlen:=function(f)
	S:=Coefficients(f);
	len:=0;
	for q in S do 
		len:=len + #Terms(Numerator(q));
	end for;
	return len;
end function;

//compute nice representative
minrep:=function(L)
	v:=RationalFunction(Random(L));
	m:=termlen(v);
	for w in L do
		n:=termlen(RationalFunction(w));
		if n lt m then
			v:=w;
			m:=n;
		end if;
	end for;
	return(v); 
end function;



RRsample:=function(d1,d2,p,sample)
	K:=GF(p);
	R<x>:=FunctionField(K);
	P<y>:=PolynomialRing(R);
	F<alpha>:=ext<R|y^2-(x^3-x)>;
	I:=InfinitePlaces(F);
	Q:=Divisor(I[1]);
	J:=Zeroes((F!x))[1];
	P:=Divisor(J);
	
	V,f:=RiemannRochSpace(d1*Q+d2*P);
	T:=[];
	count:=0;
	while (count le sample) do
		v:=Random(V);
		if not(v eq 0) then
	    		D:=Divisor(f(v));
	    		n:=Valuation(D,I[1]);
			m:=Valuation(D,J);
	    		if ( (n eq -d1) and (m eq -d2) ) then
	        			T:=T cat [f(v)];
				count:=count+1;
			end if;	
	    	end if;
	end while;
	
	return T;
end function;

//sample
statisticsRR:=function(d1,d2,p,sample);
    K:=GF(p);
    R<x>:=FunctionField(K);
    P<y> := PolynomialRing(R);
    F<alpha> := ext< R | y^2 - (x^3-x) >;
    L:=lower_bound(d1,d2,p);
    U:=upper_bound(d1,d2,p);
    anumbers:= AssociativeArray();
    fns:=AssociativeArray();

    if (d1*d2 mod p eq 0) then
	    return "Invalid";
    end if;

    T:=RRsample(d1,d2,p,sample);

    max:=0;
    pol:=0;

    // count := 0;
    for f in T do
	a:=anumber_Y(F,f,p); 
        if IsDefined(anumbers,a) then
	        anumbers[a] := anumbers[a]+1;
	        	fns[a]:=Append(fns[a],f);
        else
	        anumbers[a] := 1;
	        fns[a]:=[f];
        end if;
        if (a gt max) then
	        max:=a;
	        pol:=f;
        end if;	
     //   count := count+1;
      //  if ((count mod 100) eq 0) then
       //     count, max,pol;
        // end if;
    end for;


    "Statistics",d1,d2,p;
    "Lower:",L,"Found:",max,"Upper:",U;
    for n in [0..max] do
        if IsDefined(anumbers,n) then
		print n, anumbers[n],minrep(fns[n]);
	end if;
    end for;
    return anumbers,fns;
end function;



sampleloop:=function(S,lo,hi,sample);
	for p in S do
		for d in [lo..hi] do
			if not (d mod p eq 0) then
				statisticsRR(d,p,sample);
			end if;
		end for;
	end for;
return 0;
end function;

//statisticsRR(4,6,5,1000);

