
// Make newform 784.4.a.i in Magma, downloaded from the LMFDB on 19 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_784_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_784_4_a_i();".
// 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_784_4_a_i();".
// 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_784_a();"
function MakeCharacter_784_a()
    N := 784;
    order := 1;
    char_gens := [687, 197, 689];
    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_784_a_Hecke(Kf)
    return MakeCharacter_784_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], [-2], [16], [0], [-24], [68], [-54], [-46], [-176], [-174], [-116], [74], [10], [480], [-572], [-162], [-86], [904], [-660], [-1024], [-770], [904], [682], [102], [218], [80], [108], [-1492], [-1510], [-2138], [-1616], [-414], [-926], [494], [-474], [-88], [-2972], [1480], [2500], [-980], [2580], [-2388], [1912], [-1938], [-122], [-3900], [1164], [1480], [3906], [-252], [-2070], [5352], [6298], [-1810], [334], [3488], [-3796], [2416], [6734], [2250], [3642], [7288], [-3790], [2760], [6714], [-6386], [-7168], [2478], [2144], [-7652], [-7442], [-2616], [-6616], [-3770], [6232], [12212], [7706], [9260], [2638], [9122], [-5370], [8046], [5064], [7922], [-1080], [-3740], [5458], [-3674], [-6252], [-19472], [-2806], [588], [7136], [7140], [12908], [-12192], [-420], [-18582], [17106], [-8786], [7952], [-9042], [-22818], [726], [2504], [16158], [2070], [-7002], [13272], [25990], [448], [11018], [-3786], [-342], [-26992], [-16978], [22510], [-9404], [-4182], [-6096], [-6632], [-1990], [22124], [-6428], [-23838], [-886], [16602], [27260], [-3652], [13432], [-34360], [5624], [-10392], [2666], [3138], [-38878], [-3976], [23258], [2604], [39926], [-17422], [3870], [22120], [13724], [11320], [18876], [33748], [-3230], [-43922], [-26968], [-11486], [-22026], [-8428], [-18996], [-20948], [1392], [-41432], [25674], [46238], [-17432], [-26880], [30970], [19344], [50122], [25850], [32012], [-18240], [23144], [-14562], [-5628], [12340], [-11004], [18656], [19722], [-1076], [-27074], [51724], [-29526], [-23856], [13896], [-14312], [-13506], [47122], [-68762], [26512], [42524], [70], [-36318], [-43390], [43300], [18622], [-44364], [16052], [-76444], [13268], [63294], [-78482], [58482], [-81914], [-2392], [79334], [-72428], [-5256], [-53818], [56482], [36024], [-80220], [-20808], [-84370], [84422], [-81026], [-23840], [-75608], [-16646], [37384], [-54414], [20816], [-19454], [-87336], [-86438], [27410], [83672], [10066], [-35320], [4530], [-55838], [56854], [24776], [-89664], [12984], [-46326], [44710], [14352], [-1590], [-7070], [-15744], [43210], [80214], [-92476], [-85120], [34848], [30130], [38144], [-31258], [-58602], [-90780], [-53648], [-6030], [-52112], [68752], [-32626], [69438], [-53424], [14694], [-85672], [112356], [-71790], [-3466], [131216], [1234], [83192], [-61944], [48912], [46552], [99034], [-48010], [-44834], [-16602], [-94436], [-112286], [-17652], [88768], [-30226], [85248], [-75794], [62704], [33376], [-22210], [-44178], [-16742], [-25432], [108968], [-53508], [31560], [-72190], [-144496], [92770], [-22842], [105500], [-17010], [32566], [-11494], [11718], [-96334], [75334], [-17844], [-135520], [59996], [13662], [-50066], [32226], [-50806], [92192], [-138936], [7184], [-158054], [-57072], [-109904], [104016], [105502], [165624], [34698], [-143542], [-51144], [23168], [-28442], [148542], [-42720], [-21978], [-155774], [70870], [145762], [-48680], [-113128], [-21786], [-44686], [-95536], [21058], [9008], [58482], [1566], [-121878], [-98604], [148598], [43804], [-60830], [134074], [128558], [176118], [-10076], [181862], [-22362], [55788], [-1248], [-165734], [131328], [51052], [-182824], [154824], [60648], [-143090], [51210], [-5562], [133812], [184966], [192846], [-185844], [-200870], [111474], [48832], [84438], [-76962], [112464], [-187300], [-136814], [-82586], [-117940], [36992], [107766], [-152220], [-65314], [-81424], [-112234], [-115048], [-95698], [-168094], [81592], [71718], [-81926], [-177406], [236508], [34230], [-118648], [69330], [-247314], [-241256], [-147448], [122192], [-83036], [167940], [31102], [82176], [108540], [-148778], [112888], [-33710], [-101764], [-96294], [-132068], [-143190], [232760], [115698], [119152], [-174174], [102148], [14380], [-169530], [8158], [110806], [-17250], [39878], [-55182], [-84760], [-56452], [184090], [-273420], [-216312], [-278692], [135546], [152448], [-41590], [-178056], [180272], [-80666], [-59542], [-158736], [241824], [195162], [-52122], [71764]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_784_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


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