
// Make newform 784.4.a.d in Magma, downloaded from the LMFDB on 20 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_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_784_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_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], [-4], [-12], [0], [-12], [76], [-8], [100], [56], [-166], [232], [-414], [72], [452], [-424], [-18], [-444], [-284], [-524], [1008], [896], [40], [-1388], [448], [-824], [-1916], [824], [-1796], [-1238], [626], [160], [724], [-1370], [-852], [750], [440], [1660], [-2532], [-1048], [-884], [-396], [-228], [-888], [-2942], [1086], [1544], [-4700], [-3440], [4244], [-3700], [3270], [1856], [936], [3852], [-6720], [-3072], [2556], [7616], [1934], [-6602], [4708], [3732], [-4812], [-8480], [6568], [4710], [7796], [-4434], [-1068], [-6708], [-9392], [13496], [-4976], [-5074], [-9852], [-2568], [10370], [11836], [7950], [-1816], [-11012], [2526], [-5136], [3656], [-5376], [-11236], [-1790], [-8966], [-13732], [9368], [10004], [10568], [9688], [-19044], [-1932], [-15824], [8348], [-11208], [-1372], [-3834], [-20516], [-18346], [6108], [2678], [17748], [-13840], [11964], [26032], [-25440], [-13280], [9504], [-20254], [-3386], [-13196], [8976], [-20802], [-24964], [6920], [9002], [-20260], [-2780], [-3106], [5788], [-12708], [-31084], [-31270], [-36766], [14216], [1656], [-18700], [-11012], [22024], [8224], [626], [-19368], [-26216], [-24156], [27988], [-3812], [11482], [3764], [-23762], [16640], [10172], [-11804], [-26568], [-3604], [34536], [24348], [17032], [-11126], [-1344], [29252], [9080], [27052], [-23424], [11376], [31864], [44064], [16260], [25756], [-20502], [2696], [-4060], [-28786], [-24808], [51656], [40484], [338], [16412], [-52932], [47180], [-6752], [21030], [35560], [-53040], [-3652], [-5022], [25552], [-8780], [-51888], [-1524], [46818], [-2928], [68944], [9052], [3526], [40548], [46214], [-32968], [7856], [30988], [-13516], [-45332], [24596], [-52848], [-32306], [62746], [-4320], [19968], [-23274], [27800], [20948], [74720], [-49020], [-19980], [24760], [36156], [67674], [-45124], [33102], [71748], [57752], [644], [63552], [83672], [-16472], [-86984], [21888], [-58214], [-7742], [47808], [-78270], [-73472], [29316], [59314], [-102848], [-49904], [16432], [9316], [-59510], [-2164], [-4664], [-10086], [35700], [75744], [-83624], [-23442], [38316], [15504], [-84068], [-20092], [37648], [-17930], [87840], [-18232], [-113344], [108036], [-101836], [51384], [-105786], [31952], [51440], [31360], [69348], [-28060], [71090], [93244], [-27996], [-55096], [-8608], [-58628], [-6524], [-79148], [38504], [-63460], [28150], [-41824], [88588], [-12606], [48300], [-71156], [-47056], [-85088], [23776], [-33200], [-103196], [-3770], [-48358], [6884], [-117840], [-26368], [46008], [66772], [-43468], [-59984], [-137678], [83982], [-128776], [-115280], [-56810], [-15260], [99478], [-146428], [-58698], [101244], [-12992], [-92900], [143100], [116236], [134792], [28682], [9176], [79948], [-54108], [-47906], [-37036], [96340], [74456], [106846], [-64368], [57618], [13598], [81292], [30288], [47488], [95836], [121112], [-16256], [-91634], [21964], [-11978], [-165868], [193448], [151578], [-133944], [23884], [-7324], [115568], [59330], [-81994], [77946], [-204424], [-39652], [102876], [-105788], [-96134], [24848], [-6256], [-79976], [-58898], [-31242], [85052], [-70944], [159834], [142220], [105228], [78580], [149808], [-71804], [187740], [-16102], [-227370], [106008], [-115826], [93456], [1960], [-100332], [-366], [133640], [-156002], [189600], [168192], [-125556], [-106012], [82182], [127596], [-100024], [-30614], [-172588], [44076], [107688], [190254], [28256], [-108154], [176228], [156352], [-43072], [-33464], [-156152], [55308], [154106], [-81920], [137570], [18508], [121920], [-2424], [220884], [181708], [-167976], [4126], [33348], [121372], [168012], [171176], [-2698], [-223896], [-123528], [-195700], [106162], [240308], [10562], [-115496], [-116280], [149356], [-105992], [-197194], [228306], [-192580], [-232356], [145760], [-132526], [36084], [-224332], [253446], [199852], [-12360], [-143624], [-119704], [-277328], [115674], [9796], [60152], [280268], [114280], [-245196], [-233044], [-170058], [260076], [-130104]];
    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_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_784_4_a_d(: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_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_784_4_a_d( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_784_a();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,4,sign)));
    Vf := Kernel([<3,R![4, 1]>,<5,R![12, 1]>],Snew);
    return Vf;
end function;