
// Make newform 1216.4.a.d in Magma, downloaded from the LMFDB on 25 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_1216_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_1216_4_a_d();".
// 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_1216_4_a_d();".
// 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_1216_a();"
function MakeCharacter_1216_a()
    N := 1216;
    order := 1;
    char_gens := [191, 837, 705];
    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_1216_a_Hecke(Kf)
    return MakeCharacter_1216_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], [1], [8], [-17], [70], [61], [83], [-19], [115], [-279], [-72], [34], [108], [192], [-392], [-131], [609], [-338], [461], [750], [1177], [-22], [810], [-476], [1426], [1230], [310], [1791], [779], [-1430], [466], [-2240], [-1145], [2336], [-3516], [-390], [-2306], [1272], [1768], [3726], [1048], [-2678], [1443], [890], [300], [3339], [2149], [3834], [329], [-430], [-438], [2099], [-2996], [-5518], [-4068], [4992], [1970], [-1861], [1268], [-8912], [-3302], [5435], [-1740], [-2837], [-1579], [-10065], [-2953], [-988], [4458], [3522], [8809], [-12611], [9412], [-6553], [-14111], [-9902], [-258], [-1224], [12104], [6076], [-8148], [13849], [-9322], [7090], [9740], [4658], [14314], [3641], [11540], [-13732], [-15978], [-4300], [-12326], [7236], [1148], [5739], [-8442], [744], [-10013], [-14582], [-17924], [-11960], [-1348], [-9820], [20904], [-15093], [16740], [3042], [27816], [19320], [17926], [-18294], [22538], [17558], [-12440], [-10842], [-322], [2995], [-14652], [-20613], [7875], [6216], [-13961], [7788], [11782], [22700], [-12162], [2199], [-18965], [5552], [-13136], [-13416], [-33348], [-10606], [-29061], [-15955], [36763], [-31655], [12945], [-20263], [-29477], [1756], [1721], [-3063], [25467], [2976], [-2126], [-4684], [-1370], [30630], [-19617], [1614], [-35258], [2166], [30419], [32864], [20225], [-29415], [3697], [55245], [-54296], [-25770], [20296], [-34476], [39952], [-56942], [45772], [24916], [23136], [-9460], [31071], [-571], [-31024], [19455], [-29820], [2726], [-10668], [-7510], [-26396], [-17500], [12821], [-10305], [-48694], [5768], [-67789], [-61638], [-8397], [-23632], [-16786], [-39245], [1402], [2592], [-29232], [-4845], [-35904], [7848], [-7015], [-58408], [-16525], [42477], [-64650], [-54112], [-24619], [33412], [12578], [17836], [-24473], [-79071], [42957], [-80379], [349], [4472], [34532], [63947], [79188], [25102], [33791], [-27830], [14022], [81114], [16254], [-70186], [-82614], [-54205], [26875], [-62634], [-35238], [102250], [39562], [-28662], [-49706], [14789], [98956], [-60354], [-40106], [11769], [39238], [-89232], [-15174], [732], [-38364], [37530], [-36767], [28477], [-2244], [60260], [-82304], [94940], [-29556], [30063], [-79474], [-86365], [-113209], [-4458], [-55148], [-40574], [56638], [46868], [70662], [-68186], [9972], [126084], [-3836], [23057], [39082], [-89279], [-12341], [47755], [11], [92058], [63877], [-68081], [-92599], [5685], [100106], [30365], [22098], [-4451], [31580], [111876], [82618], [-113905], [152427], [-4415], [-66864], [39800], [69057], [-57934], [-162731], [-48210], [-41250], [-48644], [-112706], [45532], [-11798], [156520], [-51002], [-76502], [-68429], [42054], [-130762], [-649], [25767], [-5987], [134252], [155614], [-23479], [41191], [54824], [-139548], [114984], [-48570], [74778], [96892], [-15720], [-150872], [-116806], [12018], [137427], [-67122], [33739], [169215], [69321], [116166], [-39585], [-50134], [114529], [50903], [-129948], [-23524], [30176], [-195777], [73775], [110898], [76854], [-113052], [-168697], [-34296], [199199], [-54259], [34103], [-176913], [169841], [82918], [32717], [104508], [-2262], [-112781], [60862], [-97716], [123424], [-214846], [-79256], [90144], [130578], [129876], [71999], [-42532], [205335], [133818], [-137704], [-119946], [-25466], [-224439], [-144362], [27712], [-95010], [31338], [-106610], [48342], [179243], [142983], [-56191], [224828], [-53576], [-51237], [-97248], [17902], [58269], [114188], [110898], [65258], [-91183], [-140762], [-15147], [-272718], [110897], [27676], [177358], [100312], [35226], [152615], [112218], [264390], [-117280], [-85538], [-38090], [215348], [-83348], [209732], [-181210], [31890], [92842], [-16320], [177822], [-77325], [-11110], [-197172], [-196713], [200268], [31884], [125686], [-86787], [114830], [21461], [-13164], [118966], [-70048], [-13309], [112292], [2450], [143323], [-211919], [32027], [-245066], [-220441], [280741], [-266469], [-204940], [-61941]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_1216_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


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