
// Make newform 1856.4.a.f in Magma, downloaded from the LMFDB on 20 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_1856_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_1856_4_a_f();".
// 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_1856_4_a_f();".
// 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_1856_a();"
function MakeCharacter_1856_a()
    N := 1856;
    order := 1;
    char_gens := [639, 581, 321];
    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_1856_a_Hecke(Kf)
    return MakeCharacter_1856_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], [7], [15], [-18], [-27], [57], [-44], [-152], [-152], [29], [-173], [120], [-314], [-339], [-357], [59], [572], [420], [-660], [726], [1004], [361], [168], [58], [-1206], [-1440], [1858], [1914], [-989], [-278], [-600], [1652], [-1212], [-1328], [1777], [-1934], [-734], [3337], [-198], [2598], [1510], [2865], [-1936], [3850], [-3034], [246], [1493], [1402], [-1166], [1466], [847], [-444], [5297], [-4061], [-6843], [3191], [-7436], [-2755], [2410], [2235], [1648], [-2878], [-4903], [-32], [-5455], [9552], [2433], [-7460], [-4002], [2167], [-2354], [-3471], [-10880], [-8827], [876], [5446], [2216], [-6823], [-5589], [-8258], [-3946], [-1380], [-5820], [-14504], [14912], [-7180], [1398], [-5658], [-11410], [-2560], [-10395], [-5627], [8738], [-3177], [5816], [-2267], [-6481], [-21357], [6192], [19800], [-10526], [11834], [-21771], [20778], [5596], [-14908], [-8136], [-3521], [7749], [5054], [27389], [2375], [3412], [-17659], [12742], [-10988], [13190], [22924], [10278], [10607], [24698], [16801], [1258], [1256], [15996], [189], [-2431], [-3630], [-1712], [-26584], [-8125], [-31804], [-37504], [76], [-23374], [-5448], [16062], [-19182], [-5556], [2616], [9626], [-16007], [23216], [46765], [37820], [6747], [13274], [39069], [-35329], [-36378], [219], [-5598], [-5886], [-29673], [-13668], [-4963], [9954], [30930], [-38770], [30115], [-31319], [7623], [-6269], [34020], [-51291], [38687], [-3158], [21944], [21121], [-13492], [6690], [-32125], [61998], [-3386], [-21172], [42486], [24528], [1504], [-30639], [-43978], [-4627], [67261], [-22058], [-56359], [-7640], [52141], [-52224], [27716], [5348], [-61930], [-44195], [61219], [-26675], [-66056], [-19812], [52322], [15992], [-42218], [23337], [-13788], [-11934], [57848], [-10834], [66], [-16788], [-14171], [-27739], [70906], [-22741], [-37361], [83666], [87150], [9737], [43749], [59096], [-51811], [17648], [-3164], [12552], [-15786], [-44252], [-57346], [-1866], [-54827], [4584], [66094], [30000], [-14088], [-73037], [44546], [-2488], [-27056], [-89188], [39024], [-86938], [69120], [19780], [-68342], [-102414], [66939], [48691], [-16742], [-41558], [99852], [30882], [-5406], [94096], [-7168], [-35150], [67857], [47360], [-79297], [29217], [-73600], [117830], [-80036], [30518], [113855], [72826], [4529], [-11840], [-27939], [96470], [-104536], [-22812], [35785], [-43304], [-90048], [8012], [-21309], [-105601], [-45214], [5757], [49048], [-112536], [-103793], [-18516], [129993], [60220], [-125896], [-68334], [59966], [8216], [86942], [28236], [21124], [126248], [-57402], [85442], [101401], [-66203], [-57864], [-115851], [113211], [-71098], [-64983], [-157665], [-113189], [-140898], [-102340], [-52280], [123746], [-160456], [177437], [-138004], [104321], [49509], [100061], [94730], [-109235], [-166852], [47950], [-67326], [141472], [234], [126323], [-13780], [11042], [62601], [-98428], [111839], [56833], [83130], [-169051], [127963], [-176542], [174466], [15908], [41805], [-80093], [138268], [-164607], [40340], [51620], [-172249], [-112696], [-188838], [-68598], [16450], [-43780], [15840], [64002], [21684], [-7024], [-147642], [-10375], [41556], [4044], [-29811], [-49650], [61716], [-201634], [-152622], [136432], [124948], [-21492], [-119049], [20438], [-47388], [-99782], [73201], [220313], [33141], [-2670], [-45271], [-125710], [204336], [-206866], [160594], [37337], [-226178], [145278], [177290], [32236], [-146309], [71452], [139301], [-166160], [-21961], [-152107], [257512], [-48127], [-51329], [153036], [-224426], [232814], [140824], [-154211], [-210739], [21604], [718], [-13749], [190263], [233975], [-17289], [-226051], [-189884], [218020], [237962], [-92790], [-207606], [29730], [35172], [-81891], [214937], [-121964], [-97397], [109608], [-58627], [137060], [234942], [54171], [-119730], [107552], [-262630], [-92478], [53805], [-114782], [36124], [-81391], [-96282], [104670], [-73159], [-297165], [283543], [-23977], [41820], [178080], [-232980], [-139095]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_1856_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


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