
// Make newform 720.4.a.p in Magma, downloaded from the LMFDB on 18 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_720_a();"
// To make the coeffs of the qexp of the newform in the Hecke field type "qexpCoeffs();"
// To make the newform (type ModFrm), type "MakeNewformModFrm_720_4_a_p();".
// This may take a long time!  To see verbose output, uncomment the SetVerbose lines below.
// The precision argument determines an initial guess on how many Fourier coefficients to use.
// This guess is increased enough to uniquely determine the newform.
// To make the Hecke irreducible modular symbols subspace (type ModSym)
// containing the newform, type "MakeNewformModSym_720_4_a_p();".
// This may take a long time!  To see verbose output, uncomment the SetVerbose line below.
// The default sign is -1.  You can change this with the optional parameter "sign".
function ConvertToHeckeField(input: pass_field := false, Kf := [])
    if not pass_field then
        Kf := Rationals();
    end if;
    return [Kf!elt[1] : elt in input];
end function;


// To make the character of type GrpDrchElt, type "MakeCharacter_720_a();"
function MakeCharacter_720_a()
    N := 720;
    order := 1;
    char_gens := [271, 181, 641, 577];
    v := [1, 1, 1, 1];
    // chi(gens[i]) = zeta^v[i]
    assert UnitGenerators(DirichletGroup(N)) eq char_gens;
    F := CyclotomicField(order);
    chi := DirichletCharacterFromValuesOnUnitGenerators(DirichletGroup(N,F),[F|F.1^e:e in v]);
    return MinimalBaseRingCharacter(chi);
end function;

function MakeCharacter_720_a_Hecke(Kf)
    return MakeCharacter_720_a();
end function;


function ExtendMultiplicatively(weight, aps, character)
    prec := NextPrime(NthPrime(#aps)) - 1; // we will able to figure out a_0 ... a_prec
    primes := PrimesUpTo(prec);
    prime_powers := primes;
    assert #primes eq #aps;
    log_prec := Floor(Log(prec)/Log(2)); // prec < 2^(log_prec+1)
    F := Universe(aps);
    FXY<X, Y> := PolynomialRing(F, 2);
    // 1/(1 - a_p T + p^(weight - 1) * char(p) T^2) = 1 + a_p T + a_{p^2} T^2 + ...
    R<T> := PowerSeriesRing(FXY : Precision := log_prec + 1);
    recursion := Coefficients(1/(1 - X*T + Y*T^2));
    coeffs := [F!0: i in [1..(prec+1)]];
    coeffs[1] := 1; //a_1
    for i := 1 to #primes do
        p := primes[i];
        coeffs[p] := aps[i];
        b := p^(weight - 1) * F!character(p);
        r := 2;
        p_power := p * p;
        //deals with powers of p
        while p_power le prec do
            Append(~prime_powers, p_power);
            coeffs[p_power] := Evaluate(recursion[r + 1], [aps[i], b]);
            p_power *:= p;
            r +:= 1;
        end while;    
    end for;
    Sort(~prime_powers);
    for pp in prime_powers do
        for k := 1 to Floor(prec/pp) do
            if GCD(k, pp) eq 1 then
                coeffs[pp*k] := coeffs[pp]*coeffs[k];
            end if;
        end for;
    end for;
    return coeffs;
end function;


function qexpCoeffs()
    // To make the coeffs of the qexp of the newform in the Hecke field type "qexpCoeffs();"
    weight := 4;
    raw_aps := [[0], [0], [5], [-34], [18], [12], [106], [44], [56], [-270], [-204], [120], [-80], [-536], [-536], [-542], [-174], [186], [-332], [-132], [-602], [548], [-492], [1052], [482], [-1214], [-898], [1364], [218], [-1386], [814], [-1282], [-3066], [1332], [1470], [2592], [-3332], [748], [-2560], [-1206], [-1694], [3722], [2836], [-234], [3814], [-2352], [3660], [2646], [240], [-4698], [-3814], [-2148], [-3370], [6134], [4566], [-1920], [5802], [-1640], [-2792], [-1108], [-6028], [-7994], [-736], [-5380], [-1370], [-5770], [4172], [8206], [10848], [-1694], [6642], [-10376], [2198], [-12220], [10388], [10552], [8262], [2864], [12588], [10330], [-1250], [5670], [-12976], [-9050], [17528], [2568], [-12652], [-6230], [-5290], [-8110], [-2020], [-9684], [18426], [4558], [460], [-8568], [-16374], [-21620], [-16524], [-4990], [15224], [5698], [5976], [-16460], [18236], [20842], [-11772], [-4514], [25096], [16262], [-2262], [14216], [-2558], [17044], [-20980], [-7176], [2724], [-10392], [11958], [13366], [14698], [-7570], [-21378], [-15804], [22028], [1762], [-2474], [-32040], [12874], [28208], [-29068], [-28152], [29916], [32904], [21764], [-3570], [19486], [-19764], [-14390], [-28536], [-27732], [-8086], [-39854], [-17752], [23858], [13888], [-16568], [13034], [34356], [16016], [19780], [41036], [35108], [18648], [-21688], [-42064], [44420], [17124], [-11110], [-12962], [-25672], [-2082], [-5666], [28622], [24586], [40632], [-8768], [37212], [2342], [42306], [-10894], [-35026], [-9080], [20398], [-29180], [-41460], [-28692], [25734], [-2334], [39518], [64634], [-7706], [-13928], [-6810], [12672], [12590], [-34252], [4060], [27430], [27036], [-17138], [-552], [25580], [41766], [-32924], [15614], [-8990], [25708], [-32466], [24512], [52982], [-19916], [27652], [398], [-54830], [-41826], [-83996], [-12280], [-19956], [-70612], [20126], [69902], [-67354], [-36632], [-32036], [-2942], [48230], [38572], [-50632], [-92910], [22478], [-5624], [35160], [26578], [-39676], [25198], [-86754], [-44280], [44810], [-51990], [98396], [-46004], [-63224], [-956], [59820], [-15120], [11170], [-55458], [19442], [30428], [432], [68980], [-38718], [-57022], [-66782], [53820], [4070], [-83862], [95484], [38584], [29248], [-46464], [-31912], [-47710], [-102422], [-24234], [18970], [-56836], [22098], [-115402], [112622], [-34496], [-72098], [-117712], [1862], [39428], [105614], [-34664], [-57292], [3230], [-62530], [-92112], [79486], [-33604], [16398], [66582], [137616], [27498], [-125862], [15462], [81432], [-113800], [47240], [-32522], [-35968], [54396], [-65238], [82382], [-33184], [77280], [-25130], [130332], [-122766], [-66478], [30232], [68126], [-47452], [163794], [84954], [100144], [-76178], [29846], [-54760], [75660], [17436], [-50014], [97444], [-37410], [28412], [-31996], [-19952], [89922], [7580], [59288], [-17336], [64726], [131306], [-163164], [33438], [-51000], [-79652], [169686], [2874], [-123150], [53626], [106518], [48740], [46288], [-135696], [-1102], [-54362], [-80254], [-84268], [112444], [-99316], [-212972], [198582], [18186], [-202310], [126126], [-93620], [-71430], [-129206], [-3752], [47578], [86858], [11126], [167888], [-53308], [102258], [-36452], [-63926], [-96678], [66166], [-114022], [-47486], [-135408], [224098], [30430], [76928], [201260], [197472], [-29952], [-178094], [-194044], [-98062], [111594], [22150], [199478], [34650], [151636], [-107560], [111154], [40], [27232], [174666], [-197304], [12714], [32628], [-87406], [225702], [25822], [-1818], [-86502], [228196], [-166016], [18700], [185976], [220328], [-41552], [-184146], [-27926], [-100054], [138308], [-87856], [-186654], [106760], [37016], [-242108], [43606], [-121010], [-13218], [66246], [-226618], [179610], [81548], [185896], [39900], [-188212], [252462], [-172422], [-24834], [44876], [31892], [-164050], [-78110], [-42168], [184414], [-198378], [-240912], [-7186], [-16436], [30576], [233854], [73102], [239694], [29348], [-295632], [31220], [51216]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_720_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


// To make the newform (type ModFrm), type "MakeNewformModFrm_720_4_a_p();".
// This may take a long time!  To see verbose output, uncomment the SetVerbose lines below.
// The precision argument determines an initial guess on how many Fourier coefficients to use.
// This guess is increased enough to uniquely determine the newform.
function MakeNewformModFrm_720_4_a_p(:prec:=1)
    chi := MakeCharacter_720_a();
    f_vec := qexpCoeffs();
    Kf := Universe(f_vec);
    // SetVerbose("ModularForms", true);
    // SetVerbose("ModularSymbols", true);
    S := CuspidalSubspace(ModularForms(chi, 4));
    S := BaseChange(S, Kf);
    maxprec := NextPrime(2999) - 1;
    while true do
        trunc_vec := Vector(Kf, [0] cat [f_vec[i]: i in [1..prec]]);
        B := Basis(S, prec + 1);
        S_basismat := Matrix([AbsEltseq(g): g in B]);
        if Rank(S_basismat) eq Min(NumberOfRows(S_basismat), NumberOfColumns(S_basismat)) then
            S_basismat := ChangeRing(S_basismat,Kf);
            f_lincom := Solution(S_basismat,trunc_vec);
            f := &+[f_lincom[i]*Basis(S)[i] : i in [1..#Basis(S)]];
            return f;
        end if;
        error if prec eq maxprec, "Unable to distinguish newform within newspace";
        prec := Min(Ceiling(1.25 * prec), maxprec);
    end while;
end function;


// To make the Hecke irreducible modular symbols subspace (type ModSym)
// containing the newform, type "MakeNewformModSym_720_4_a_p();".
// This may take a long time!  To see verbose output, uncomment the SetVerbose line below.
// The default sign is -1.  You can change this with the optional parameter "sign".
function MakeNewformModSym_720_4_a_p( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_720_a();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,4,sign)));
    Vf := Kernel([<7,R![34, 1]>,<11,R![-18, 1]>],Snew);
    return Vf;
end function;