
// Make newform 648.4.a.a in Magma, downloaded from the LMFDB on 20 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_648_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_648_4_a_a();".
// 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_648_4_a_a();".
// 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_648_a();"
function MakeCharacter_648_a()
    N := 648;
    order := 1;
    char_gens := [487, 325, 569];
    v := [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_648_a_Hecke(Kf)
    return MakeCharacter_648_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], [36], [64], [-65], [59], [-28], [160], [-57], [164], [-321], [-246], [-8], [84], [478], [-32], [415], [-220], [884], [-77], [-80], [1268], [123], [1346], [-834], [1060], [300], [-557], [-1617], [-1136], [116], [615], [1452], [839], [-2056], [1099], [-72], [2124], [-4005], [-1428], [-2226], [3756], [-3637], [795], [2500], [-3944], [-680], [3924], [4615], [-4493], [-3124], [-309], [1352], [3695], [-980], [-5293], [-4912], [-1394], [-7165], [-9340], [-6873], [204], [4056], [1095], [10419], [3596], [-2910], [-5076], [8118], [-9262], [1248], [9880], [-9778], [3260], [1620], [6910], [-1705], [-4005], [-5385], [5640], [4331], [-11596], [-2765], [7932], [-12268], [-10478], [-16129], [-1034], [9080], [18504], [17012], [-6140], [-10944], [-5208], [1440], [2310], [-9910], [6640], [-15969], [500], [-3845], [-3188], [17263], [-4656], [-23209], [-17624], [17695], [2748], [19579], [17804], [-22690], [1755], [-6920], [10744], [9931], [4068], [17540], [-30090], [-8820], [13319], [14035], [-9410], [1060], [8540], [17259], [-24145], [-17792], [10620], [-8650], [4652], [33868], [-12852], [-9730], [-33361], [-26557], [-18333], [-37376], [16075], [32183], [-1424], [-26961], [-21100], [35120], [21238], [-26088], [-11122], [-2985], [-49548], [21688], [10935], [-40126], [-42748], [-27496], [-21680], [-6824], [26512], [-22597], [38115], [-43485], [-45420], [-32505], [14276], [-35844], [-24270], [19152], [-9164], [36259], [-23917], [-8402], [7908], [-34041], [6644], [-26165], [-22600], [10682], [-54760], [57411], [-43520], [-5857], [17984], [-4372], [22443], [-23910], [38832], [-10677], [26395], [54408], [31402], [-42180], [-30425], [23448], [8440], [36019], [33396], [79322], [6271], [39238], [-1489], [-44400], [28943], [4976], [28814], [-22817], [-82620], [-26490], [-54652], [52920], [31610], [-26228], [74610], [24231], [-59300], [38480], [49352], [45866], [80960], [56731], [22380], [-3705], [-21829], [-64188], [-64197], [-1840], [87696], [467], [99690], [63180], [-71136], [55300], [-36237], [-56012], [-37720], [-43205], [61500], [-12520], [-41413], [-818], [-15640], [39724], [-8952], [-75640], [58312], [11110], [-58830], [-18328], [52156], [50640], [-27128], [-44128], [8391], [-97401], [56], [-62813], [-56553], [-17700], [-120841], [-110060], [92799], [-63494], [-7952], [-12336], [-92573], [73830], [-97406], [-76132], [-20221], [6491], [-93400], [-27565], [38779], [122960], [26890], [15748], [23815], [-2548], [-64560], [122363], [-7277], [-119600], [-80648], [4000], [105820], [134615], [120776], [27364], [102827], [56526], [-136132], [-22585], [145126], [76140], [134347], [94356], [-43497], [-41717], [26384], [24843], [124872], [61040], [1415], [115574], [-6412], [-32720], [149736], [109634], [116736], [118742], [-111292], [18575], [4740], [-120882], [141779], [-54212], [70036], [-151078], [20568], [120696], [123263], [-130302], [51760], [-76406], [92819], [158152], [35770], [-85989], [-173688], [69120], [106900], [36030], [-8714], [-43869], [-179352], [103152], [53196], [100656], [152887], [45495], [-175366], [-132840], [-76037], [-163590], [82803], [207300], [-72293], [-74300], [-77961], [-36776], [136404], [208535], [152496], [7535], [-78965], [-128040], [185867], [87155], [-155948], [-976], [138610], [-112860], [-178225], [-126141], [126624], [100688], [21376], [-160245], [63510], [-9908], [-63014], [147464], [206976], [118520], [-126737], [156304], [11651], [-70392], [-38856], [143503], [134802], [13319], [164295], [-52225], [-198276], [206415], [-179888], [35472], [13980], [-243025], [114148], [17640], [-27785], [-206285], [-61320], [204864], [-122064], [151567], [25812], [-20841], [-116992], [11071], [102563], [42690], [-42500], [85199], [25243], [-50960], [135199], [-194877], [156152], [-144012], [243527], [-163745], [264788], [30904], [13275], [77215], [-248412], [74840], [91851], [-97840], [-136810], [178526], [255680], [288392], [100087], [118031], [42148], [-162069], [157416], [66760]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_648_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


// To make the newform (type ModFrm), type "MakeNewformModFrm_648_4_a_a();".
// 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_648_4_a_a(:prec:=1)
    chi := MakeCharacter_648_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_648_4_a_a();".
// 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_648_4_a_a( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_648_a();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,4,sign)));
    Vf := Kernel([<5,R![5, 1]>],Snew);
    return Vf;
end function;