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

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


// 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;

//dim H^0(E,\ker V(n Q))
dim_kerVE:=function(n,p);
    if ((anumber_E(p) eq 0) or (n ge p)) then
        return (n-Ceiling(n/p));
    else    

    K := GF(p);
    E := EllipticCurve([K | -1, 0]);
    I:=PointsAtInfinity(E);
    Q:=Divisor(I[1]);

    D:=n*Q;
    V,f:=DifferentialSpace(-D);
    g:=Inverse(f);
    B:=Basis(V);

    S:=[];
    for w in B do
        S:=Append(S,Vector(g(Cartier(f(w)))));
    end for;
    C:=Matrix(S);

    return Dimension(Kernel(C));

    end if;
end function;

lower_bound := function(d,p);
    lowerbound:= 0;
    optimum:=0;
    
    for j in [0..p-1] do
        sum :=0;
                
        for i in [j..p-1] do
            sum +:= Floor( i *d /p) - Floor(i * d/p - (1-1/p) * j *d /p );
	    end for;
        if sum gt lowerbound then
            lowerbound := sum; 
            optimum := j;   
	    end if;
    end for;

    return lowerbound;    
end function;

tau := function(p,d,i);
    count :=0;
    for n in [1..Floor(i *d/p)] do
        m := (( -n) * Modinv(d,p) ) mod p;
        if (0 lt m) and (m le p-1-i) then
            count:= count+1;
        end if;
    end for;
    return count;
end function;

upper_bound := function(d,p);
    sum := p* anumber_E(p);

    for i in [1..p-1] do
        sum +:= (Floor( i * d/p) - Floor(i * d/p^2) -tau(p,d,i));
    end for;
    return sum;
end function;

bounds := function(d,p);
    trivU:=p*anumber_E(p) + (p-1)*(d-1)/2;
    trivL:=dim_kerVE(compute_n(0,d,p),p);
    L:=lower_bound(d,p);    
    U:=upper_bound(d,p);
    return trivL,L,U,trivU;
end function;

//get list of functions with order of vanishing -d at Q
sample:=function(d,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]);
	V,f:=RiemannRochSpace(d*Q);
	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]);
	    		if n eq -d then
	        			T:=T cat [f(v)];
				count:=count+1;
			end if;	
	    	end if;
	end while;
	
	return T;
end function;

//sample a-numbers using random function f with pole of degree d in z^p- z = f over E
statistics:=function(d,p,num);
    K:=GF(p);
    R<x>:=FunctionField(K);
    P<y> := PolynomialRing(R);
    F<alpha> := ext< R | y^2 - (x^3-x) >;
    L:=lower_bound(d,p);
    U:=upper_bound(d,p);
    anumbers:= AssociativeArray();
    fns:=AssociativeArray();

    if (d mod p eq 0) then
	    return "Invalid";
    end if;

    T:=sample(d,p,num);

    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",d,p;
    "Lower:",L,"Found:",max,"Upper:",U;
    for n in [0..max] do
        if IsDefined(anumbers,n) then
		print n, anumbers[n],Random(fns[n]);
	end if;
    end for;
    return anumbers,fns;
end function;

RRspace:=function(d,p)
	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]);
	V,f:=RiemannRochSpace(d*Q);
	T:=[];
	for v in V do
		if not(v eq 0) then
	    		D:=Divisor(f(v));
	    		n:=Valuation(D,I[1]);
	    		if n eq -d then
	        			T:=T cat [f(v)];
			end if;
		end if;
	end for;
	
	return T;
end function;



//do an exhaustive search
searchRR:=function(d,p);
    K:=GF(p);
    R<x>:=FunctionField(K);
    P<y> := PolynomialRing(R);
    F<alpha> := ext< R | y^2 - (x^3-x) >;
    L:=lower_bound(d,p);
    U:=upper_bound(d,p);
    anumbers:= AssociativeArray();

    if (d mod p eq 0) then
	    return "Invalid";
    end if;

    T:=RRspace(d,p);
        max:=0;
    	pol:=0;
    for f in T do
        a:=anumber_Y(F,f,p); 
        if IsDefined(anumbers,a) then
	        anumbers[a] := anumbers[a]+1;
        else
	        anumbers[a] := 1;
        end if;
        if (a gt max) then
	        max:=a;
	        pol:=f;
        end if;	
    end for;

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


"Table 3";
for d in [2..16] do
    if d mod 5 ne 0 then
        bounds(d,5);
    end if;
end for;

"Table 4";
for d in [2..16] do
    if d mod 7 ne 0 then
        bounds(d,7);
    end if;
end for;

"Example Statistics";
statistics(3,5,100);

//searchRR(6,7);

