Running the commands in this file (in MAGMA) computes the first 2 Hilbert modular polynomials in Chapter 2 of Chloe Martindale's PhD thesis over Q, F_p, or C for K_0 = Q(5^1/2).

load "hilbert_modular_polynomial_5.mgm";

//choose p such that FF:=F_p contains all of the l-th roots of unity (to speed up the code).

p:=2*3*5*7*11*13*17*19*23*29*31+1;
FF:=FiniteField(p); //can be any (reasonably large) finite field, or Q, or C.

//if computing both polynomials, must have (l,p) = 1

P<x>:=PolynomialRing(Rationals());
K0<a>:=NumberField(x^2-5);
u:=(a+1)/2;

//choose level of the polynomials

norm := 11; mu:=1/2*a + 7/2; //norm 11, must be a totally positive prime element of K0

//this choice of mu such that Norm(mu) = 11 minimises the absolute value (recommended)

//this function computes the repquired precision

if norm eq 5 then
    find:=[];
    arr:=Combi3(2,6,10,10+Norm(mu)*6);
    for i:= 1 to #arr do
        find:=Insert(find,i,arr[i][2]+2*arr[i][3]);
    end for;
    N:=Maximum(find) + 1;
else
    find:=[];
    arr:=Combi4(2,6,5,15,10+Norm(mu)*6);
    for i:= 1 to #arr do
	find:=Insert(find,i,arr[i][2]+arr[i][3]+2*arr[i][4]);
    end for;
    N:=Maximum(find) + 1;
end if;

//compute the coefficients of the modular polynomial (takes a long time)

time c_phi, c_psi1, c_psi2, Xphi, Xpsi, Gphi, Gpsi1, Gpsi2 :=compute_coefficients(N,K,FF,mu,u,a);

//create the first polynomial from the coefficients (c_phi)
ZZ:=IntegerRing();

if norm ne 5 then
    arrphi:=Combi4(2,6,5,15,(Norm(mu)+1)*6);
    Fun<g2,s6,s5,s15>:=FunctionField(FF,4);
    basis:=[];
    s5_lim:=ZZ!(Floor((Norm(mu)+1)*6/5));
    if IsEven(s5_lim) then
       g2_lim:=ZZ!(((Norm(mu)+1)*6-5*s5_lim)/2);
    else
       s5_lim:=s5_lim-1;
       g2_lim:=ZZ!(((Norm(mu)+1)*6-5*s5_lim)/2);
    end if;
    for i:= 1 to #arrphi do
        basis:=Insert(basis,i,(g2^arrphi[i][1]*s6^arrphi[i][2]*s5^arrphi[i][3]*s15^arrphi[i][4])/(s5^s5_lim*g2^g2_lim));
    end for;

    J1:=s6/g2^3; J2:=g2^5/s5^2; J3:=s15/s5^3;

    Pol<X1,X2,X3,Y,Z>:=PolynomialRing(FF,5);
    modpoly:=0;

    for j:= 0 to ZZ!(Norm(mu) + 1) do
        for i:= 1 to #arrphi do
            for aa:= 0 to Floor((Norm(mu)+1)*6/6) do
                for bb:= 0 to Floor((Norm(mu) + 1)*6/10) do
                    for cc:= 0 to 1 do
                        if basis[i] eq J1^aa*J2^bb*J3^cc then
                           modpoly := modpoly + c_phi[j+1][i]*X1^aa*X2^bb*X3^cc*Y^j;
                        end if;
                    end for;
                end for;
            end for;
        end for;
    end for;

//create the 2nd polynomial from the coefficients (c_psi1 and c_psi2)

s5_lim:=Floor((10+Norm(mu)*6)/5);
if IsEven(s5_lim) then
   g2_lim:=ZZ!((10+Norm(mu)*6-5*s5_lim)/2);
else
   s5_lim:=s5_lim-1;
   g2_lim:=ZZ!((10+Norm(mu)*6-5*s5_lim)/2);
end if;

    arrpsi:=Combi4(2,6,5,15,10+Norm(mu)*6);


Psi1:=0;
Psi2:=0;
basispsi:=[];
for i:= 1 to #arrpsi do
    basispsi:=Insert(basispsi,i,(g2^arrpsi[i][1]*s6^arrpsi[i][2]*s5^arrpsi[i][3]*s15^arrpsi[i][4])/(s5^s5_lim*g2^g2_lim));
end for;

for j:= 0 to ZZ!(Norm(mu)) do
    for i:= 1 to #arrpsi do
        for aa:= 0 to Floor((10+Norm(mu)*6)/6) do
            for bb:= 0 to Floor((10+Norm(mu)*6)/10) do
                for cc:= 0 to 1 do
                    if basispsi[i] eq J1^aa*J2^bb*J3^cc then
                       Psi1 := Psi1 + c_psi1[j+1][i]*X1^aa*X2^bb*X3^cc*Y^j;
                       Psi2 := Psi2 + c_psi2[j+1][i]*X1^aa*X2^bb*X3^cc*Y^j;
                    end if;
                end for;
            end for;
        end for;
    end for;
end for;

Psi:=Psi1*Z - Psi2;

end if;

if norm eq 5 then
    arrphi:=Combi3(2,6,10,(Norm(mu)+1)*6);
    Fun<g2,s6,s10>:=FunctionField(FF,3);
    basis:=[];
    for i:= 1 to #arrphi do
        basis:=Insert(basis,i,(g2^arrphi[i][1]*s6^arrphi[i][2]*s10^arrphi[i][3])/(s10^3*g2^3));
    end for;

    J1:=s6/g2^3; J2:=g2^5/s10;

    Pol<X1,X2,Y,Z>:=PolynomialRing(FF,4);
    modpoly:=0;

    for j:=0 to ZZ!(Norm(mu) + 1) do
        for i:= 1 to #arrphi do
            for aa:= 0 to Floor((Norm(mu)+1)*6/6) do
                for bb:= 0 to Floor((Norm(mu) + 1)*6/10) do
                    if basis[i] eq J1^aa*J2^bb then
                       modpoly:=modpoly + c_phi[j+1][i]*X1^aa*X2^bb*Y^j;
                    end if;
                end for;
            end for;
        end for;
    end for;

    arrpsi:=Combi3(2,6,10,10+Norm(mu)*6);
    s10_lim:=Floor((10+Norm(mu)*6)/10);
    if IsEven(s10_lim) then
       g2_lim:=ZZ!((10+Norm(mu)*6-10*s10_lim)/2);
    else
       s10_lim:=s10_lim-1;
       g2_lim:=ZZ!((10+Norm(mu)*6-10*s10_lim)/2);
    end if;

    Psi1:=0;
    Psi2:=0;
    basispsi:=[];
    for i:=1 to #arrpsi do
        basispsi:=Insert(basispsi,i,(g2^arrpsi[i][1]*s6^arrpsi[i][2]*s10^arrpsi[i][3])/(s10^s10_lim*g2^g2_lim));
    end for;
    
    for j:= 0 to ZZ!(Norm(mu)) do
        for i:=1 to #arrpsi do
            for aa:= 0 to Floor((10+Norm(mu)*6)/6) do
                for bb:= 0 to Floor((10+Norm(mu)*6)/10) do
                    if basispsi[i] eq J1^aa*J2^bb then
                        Psi1 := Psi1 + c_psi1[j+1][i]*X1^aa*X2^bb*Y^j;
                        Psi2 := Psi2 + c_psi2[j+1][i]*X1^aa*X2^bb*Y^j;
                    end if;
                end for;
            end for;
        end for;
    end for;

    Psi:= Psi1*Z - Psi2;

end if;

//save output

PrintFile("modpoly11_G",modpoly,"Default");
PrintFile("modpoly11_H",Psi,"Default");