
// Make newform 816.4.a.b in Magma, downloaded from the LMFDB on 13 September 2025.
// To make the character of type GrpDrchElt, type "MakeCharacter_816_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_816_4_a_b();".
// 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_816_4_a_b();".
// 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_816_a();"
function MakeCharacter_816_a()
    N := 816;
    order := 1;
    char_gens := [511, 613, 545, 241];
    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_816_a_Hecke( : Kf := Kf)
    return MakeCharacter_816_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]];
    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], [-3], [-7], [-4], [21], [-25], [17], [69], [-15], [58], [298], [72], [-369], [59], [138], [262], [-50], [-568], [-124], [-100], [-158], [-710], [-214], [-1016], [-1780], [32], [-233], [33], [-420], [653], [209], [-589], [1026], [-1958], [-1010], [-1872], [787], [-1166], [-2451], [2061], [-2242], [-254], [-4938], [294], [4465], [-2020], [3690], [-51], [2699], [-3890], [215], [160], [-448], [-124], [2124], [3648], [4555], [1985], [-4982], [-4276], [3106], [-60], [-6876], [6768], [-2028], [-778], [-3431], [-1894], [3372], [869], [-5502], [11508], [-2556], [3410], [13704], [8672], [-5944], [5976], [8651], [-9459], [4668], [3975], [-10312], [-11257], [-5516], [3774], [-8778], [-10179], [7114], [-8248], [3782], [9589], [6758], [1274], [-6778], [3511], [-19196], [-1033], [-15148], [4036], [-18144], [8446], [3942], [-6404], [8100], [5577], [10116], [-27470], [-4138], [17722], [12838], [9771], [-20418], [-10646], [-9719], [2773], [-29060], [-30482], [18225], [-13382], [-31943], [3234], [2479], [18147], [-30196], [-2664], [14598], [-8509], [-9944], [-21118], [9629], [-20408], [16554], [17093], [-9306], [-20663], [-39604], [25076], [-28464], [-42705], [-38870], [11505], [29488], [40699], [5470], [-16761], [-26246], [-42094], [12476], [-32184], [-21518], [27582], [37697], [-31017], [-190], [2701], [12201], [42525], [-9534], [6230], [-10420], [40356], [20291], [45514], [22278], [-34153], [46232], [9926], [240], [36882], [5406], [46253], [-12143], [32203], [45801], [28449], [-48304], [40755], [31807], [23591], [36623], [-61220], [-19718], [32904], [12158], [-57090], [21538], [11903], [8186], [12596], [55650], [42532], [-4171], [14752], [-13025], [218], [65778], [40108], [12038], [16046], [-41925], [71120], [-50919], [-35463], [-29268], [-18846], [-8608], [-8660], [43631], [-3596], [50266], [37122], [-55762], [72016], [-83353], [91628], [77503], [-18720], [61179], [-73212], [-30015], [57728], [-84972], [17558], [30210], [33470], [-61055], [-1223], [86539], [52776], [-33649], [-72354], [89908], [45854], [-56520], [59850], [-55736], [-38973], [-6439], [-45276], [6053], [-36932], [11549], [-70765], [-35463], [72777], [46818], [-26709], [-15313], [12492], [-79943], [76574], [15668], [38556], [-118910], [64604], [61708], [90964], [-59943], [-113306], [-86702], [-12372], [-110124], [-101328], [-78421], [-120287], [16470], [31692], [16640], [26958], [-9040], [-110357], [142807], [103408], [44766], [-42587], [145228], [-54475], [72229], [-97800], [44684], [-79422], [-5976], [-27282], [1346], [118260], [-14518], [-33938], [-39855], [114144], [-80799], [98597], [74326], [-111509], [101032], [-42843], [-69544], [-8412], [156809], [-156276], [8377], [96900], [-149080], [4073], [61144], [117486], [-3814], [110984], [-165340], [-125758], [26948], [88218], [39281], [-54367], [-134580], [5719], [-32230], [-9165], [-92176], [-24984], [100902], [-54986], [129582], [26657], [-70542], [-45565], [83106], [-18606], [-17800], [112335], [164752], [-80813], [131578], [-51122], [46748], [148773], [25857], [-114931], [-188014], [113463], [110885], [189796], [-187071], [-159251], [-47212], [-18059], [-128982], [150757], [-98799], [158130], [91975], [-23948], [207980], [29458], [-59977], [15912], [-113776], [-72108], [62401], [126462], [120836], [-32269], [152764], [192924], [239695], [-108439], [9391], [-109684], [68092], [-165392], [40320], [41759], [146256], [165440], [183328], [-48541], [-225516], [-174649], [44234], [-198551], [24157], [79252], [76400], [124573], [-106674], [244997], [152339], [-193826], [-13424], [-220684], [109778], [-86653], [-4992], [-165931], [-225486], [-32580], [13872], [70866], [-82988], [-10194], [-112168], [-169940], [242619], [117674], [53186], [117830], [150822], [-215988], [-238537], [98361], [180520], [268350], [-20708], [4628], [-111383], [228582], [132457], [2226], [17829], [-221794], [-11594], [-140986], [-48741], [-228870], [125708], [217940], [-273928], [146823], [49771], [24071]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_816_a_Hecke( : Kf := Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


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