//The function N(t) defined in Theorem 1.2.
function N(t)
    return t^4 + 5*t^3 + 15*t^2 + 25*t + 25;
end function;

//Given a rational number n, this function counts the number of primes p such that the
//padic valuation of the numerator of n is 0 mod 5.
function mu5(n)
    i := 0;
    fac := Factorization(Numerator(n));
    for j in [1..#fac] do
    if fac[j][2] mod 5 eq 0 then
    i := i+1;
    end if;
    end for;
    return i;
end function;

//Code to search the for x-coordinates t of n*P such that mu5(N(t)) is non-zero and n in a set A.
function SearchMu5(A)
    E := EllipticCurve([0,5,0,10,7]);
    P := E![-1,-1];
    v := [];
    for n in A do
    if n mod 6 in [1,5] then
        n;
        t := (n*P)[1];
        if mu5(N(t)) ne 0 then
             Append(~v,t);
             t;
             n;
        end if;
    end if;
    end for;
    return v;
end function;

//Given a vector v in F_2^l and an l cycle phi in S_l, computes the vector
//obtained by permuting the component of v according to phi.
function PermVector(v, phi)
    VS := Parent(v);
    l := Dimension(VS);
    G:=Sym(l);
    GS:=GSet(Sym(l));
    m:=Action(GS);
    return Parent(v)! [ v[m(i,phi)] : i in [1..l]];
end function;
 
 
 
//Given a vector v in F_2^l and an l cycle phi in S_l, computes the matrix
//M_{infty,u} where u has signature v.
function MInfty(v,phi)
    VS := Parent(v);
    l := Dimension(VS);
    Cols := [PermVector(v, phi^i): i in [0..l-2]];
    return Transpose(Matrix(BaseField(VS), l-1,l, Cols));
end function;
 
 
//Computes the dimension of the column space of MInfty(v,phi).
function dInfty(v,phi)
    return Rank(MInfty(v,phi));
end function;
 
 
//Given a vector v in F_2^l returns the sum of the components of v
function CompSum(v);
    VS := Parent(v);
    l := Dimension(VS);
    F := BaseField(VS);
    s := F!0;
    for i in [1..l] do
s := s+v[i];
end for;
return s;
end function;


//Given a Sophie Germain prime returns the signature of u.
function SignatureOfu(q)
    p := (q-1) div 2;
    V:=VectorSpace(GF(2),p);
    v:=V![ 0 :  t in [1..p]];
    if p mod 4 eq 1 then
        for t:=1 to Integers()!((p-1) div 2) do
        v[t]:=1;
        end for;
    else
        for t:=1 to Integers()!((p+1) div 2) do
        v[t]:=1;
        end for;
    end if;
    return v;
end function;

//Given a Sophie Germain prime q returns the permutation phi in Section 2.2.
function GalGen(q)
    p := (q-1) div 2;
    if p mod 4 eq 1 then
    g:=[ Min(2*i mod q, -2*i mod q) : i in [1..p]];
    else
    g:=[(p+1) - Min(2*(p+1-i) mod q, q-(2*(p+1-i) mod q)) : i in [1..p] ];
    end if;
    return Sym(p)!g;
end function;



//Does the computations mentioned in Theorem 5.3 for primes in [ n1, n2 ].
procedure VerifyThm53(n1,n2)
s:=Cputime();
for q in PrimesInInterval(n1, n2) do
    p := (q-1) div 2;
    if p gt 2 and IsPrime(p) then
        if (Order(GF(p)!2) mod 2) eq 1 then
            o:=Order(GF(p)!2);
            if o ne (p-1)/2 then
                v := SignatureOfu(q);
                phi := GalGen(q);
                assert dInfty(v,phi) eq p-1;
            else
                if (o mod 2) eq 0 then
                        v := SignatureOfu(q);
                        phi := GalGen(q);
                        assert dInfty(v,phi) eq p-1;
                end if;
            end if;
        end if;
        q,p;
    end if;
end for;
t:=Cputime();
printf "***  Verified Theorem 5.3 for all primes in [ %o, %o ] in %o seconds ***\n",n1,n2,t-s;
end procedure;



