
// Make newform 882.4.g.be in Magma, downloaded from the LMFDB on 25 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_882_g();"
// To make the character of type GrpDrchElt with Codomain the HeckeField, type "MakeCharacter_882_g_Hecke();"
// To make the coeffs of the qexp of the newform in the Hecke field type "qexpCoeffs();"
// To make the newform (type ModFrm), type "MakeNewformModFrm_882_4_g_be();".
// 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_882_4_g_be();".
// 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
        poly := [4, 0, 2, 0, 1];
        Kf := NumberField(Polynomial([elt : elt in poly]));
        AssignNames(~Kf, ["nu"]);
    end if;
    Rf_num := [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
    Rf_basisdens := [1, 1, 2, 2];
    Rf_basisnums := ChangeUniverse([[z : z in elt] : elt in Rf_num], Kf);
    Rfbasis := [Rf_basisnums[i]/Rf_basisdens[i] : i in [1..Degree(Kf)]];
    inp_vec := Vector(Rfbasis)*ChangeRing(Transpose(Matrix([[elt : elt in row] : row in input])),Kf);
    return Eltseq(inp_vec);
end function;


// To make the character of type GrpDrchElt, type "MakeCharacter_882_g();"
function MakeCharacter_882_g()
    N := 882;
    order := 3;
    char_gens := [785, 199];
    v := [3, 2];
    // 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;

// To make the character of type GrpDrchElt with Codomain the HeckeField, type "MakeCharacter_882_g_Hecke();"
function MakeCharacter_882_g_Hecke(Kf)
    N := 882;
    order := 3;
    char_gens := [785, 199];
    char_values := [[1, 0, 0, 0], [0, 0, 1, 0]];
    assert UnitGenerators(DirichletGroup(N)) eq char_gens;
    values := ConvertToHeckeField(char_values : pass_field := true, Kf := Kf); // the value of chi on the gens as elements in the Hecke field
    F := Universe(values);// the Hecke field
    chi := DirichletCharacterFromValuesOnUnitGenerators(DirichletGroup(N,F),values);
    return chi;
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, 0, -2, 0], [0, 0, 0, 0], [0, -1, 6, -1], [0, 0, 0, 0], [2, -6, 2, 0], [-24, 0, 0, 15], [-66, -11, -66, 0], [0, 46, -60, 46], [0, -102, 38, -102], [56, 0, 0, -150], [216, 54, 216, 0], [0, 180, -140, 180], [18, 0, 0, 127], [-64, 0, 0, 288], [0, 338, -132, 338], [134, 192, 134, 0], [-168, 298, -168, 0], [0, 353, 252, 353], [192, -144, 192, 0], [198, 0, 0, -342], [-156, 595, -156, 0], [0, 300, -424, 300], [-324, 0, 0, 80], [0, 175, -306, 175], [1092, 0, 0, -133], [810, 351, 810, 0], [0, 350, 264, 350], [0, -402, -694, -402], [680, 1008, 680, 0], [-246, 0, 0, 912], [-1028, 0, 0, -744], [0, 1008, 156, 1008], [1164, -1350, 1164, 0], [180, 0, 0, -844], [0, 216, -926, 216], [2092, -108, 2092, 0], [-2100, 1193, -2100, 0], [0, -636, -3008, -636], [108, 0, 0, 986], [0, -609, -390, -609], [2478, 810, 2478, 0], [408, 0, 0, 417], [0, -90, -1354, -90], [1918, -1104, 1918, 0], [-4042, 0, 0, 600], [2304, 876, 2304, 0], [-1116, 0, 0, 2244], [-3912, 0, 0, 1332], [-3264, -594, -3264, 0], [0, -237, 1632, -237], [0, 2106, -1500, 2106], [-1170, 0, 0, -3498], [1116, -3615, 1116, 0], [5664, 0, 0, 170], [0, 2253, 2778, 2253], [270, 3462, 270, 0], [-1506, -1681, -1506, 0], [0, 334, -6000, 334], [-3490, 864, -3490, 0], [3220, 0, 0, 2790], [5292, -746, 5292, 0], [2430, 0, 0, 4811], [-5460, 0, 0, -938], [-5100, 26, -5100, 0], [0, -3951, -1824, -3951], [0, 2016, 3886, 2016], [0, 840, 9988, 840], [1760, 0, 0, -5388], [6810, -798, 6810, 0], [5316, 0, 0, 3261], [7194, -2057, 7194, 0], [0, -3354, 82, -3354], [0, -4904, 0, 0], [0, -6072, 5506, -6072], [3148, 0, 0, 1560], [0, 2932, -4824, 2932], [-7784, -2886, -7784, 0], [0, 1081, 3612, 1081], [0, 618, 11508, 618], [-11196, -3309, -11196, 0], [-5016, 0, 0, 7846], [1530, 0, 0, -5064], [-3966, -162, -3966, 0], [5592, 0, 0, 117], [0, -2136, 6792, -2136], [0, 4158, 50, 4158], [4602, 0, 0, 2880], [0, 336, 9582, 336], [-6870, 0, 0, -1315], [-4460, 0, 0, 1884], [0, 2606, 3816, 2606], [-1836, 7078, -1836, 0], [8388, 2532, 8388, 0], [-4074, 0, 0, -1410], [0, 10020, -1336, 10020], [7656, 0, 0, -6472], [0, 277, -4662, 277], [-12510, -1961, -12510, 0], [0, -2580, 8148, -2580], [0, 2928, -170, 2928], [7156, 0, 0, 6972], [1266, -336, 1266, 0], [-9432, 1354, -9432, 0], [0, -11442, 2236, -11442], [-3940, 4284, -3940, 0], [10728, 8577, 10728, 0], [6552, 0, 0, -12306], [0, -11045, 4842, -11045], [8794, -846, 8794, 0], [-18696, 0, 0, 891], [0, 5988, -8232, 5988], [144, -18780, 144, 0], [-23132, 0, 0, -3294], [8892, 5060, 8892, 0], [-12028, 0, 0, 2148], [-22484, -6162, -22484, 0], [-7236, 0, 0, -13954], [-7236, 12138, -7236, 0], [0, 8526, -3912, 8526], [14286, 0, 0, -1902], [16620, 2841, 16620, 0], [-12, 0, 0, 1512], [0, -10749, -12582, -10749], [-5874, 17010, -5874, 0], [0, -6072, -15348, -6072], [10080, 0, 0, -7458], [0, -5448, 1120, -5448], [0, 7520, -10464, 7520], [7872, 0, 0, -2302], [0, 4203, 16632, 4203], [19152, 2676, 19152, 0], [-11738, 0, 0, 3702], [0, -468, -20536, -468], [-11560, 0, 0, -2544], [0, -2957, 31038, -2957], [7224, 0, 0, -3507], [13494, 2801, 13494, 0], [-1116, 22696, -1116, 0], [15774, 0, 0, 18591], [-21278, -14496, -21278, 0], [-1812, 0, 0, -10428], [0, 5208, 7134, 5208], [8600, 2388, 8600, 0], [9306, 0, 0, 26106], [-20136, 13759, -20136, 0], [-1764, 0, 0, -31094], [5796, 0, 0, 2121], [-26718, 2501, -26718, 0], [0, 20114, -16548, 20114], [0, 26334, -7426, 26334], [0, -5772, -2176, -5772], [31026, 0, 0, 4869], [-1648, 0, 0, -10356], [0, 12874, 27060, 12874], [5552, 15708, 5552, 0], [5270, 0, 0, 16830], [0, -9216, -5060, -9216], [0, -6901, 20382, -6901], [-8592, 0, 0, 5279], [-22410, -23153, -22410, 0], [0, 2130, -21798, 2130], [25774, 0, 0, -7584], [5220, 0, 0, -14736], [0, 16448, -1572, 16448], [20388, -12918, 20388, 0], [10920, 26576, 10920, 0], [-14636, 8904, -14636, 0], [-10788, -11279, -10788, 0]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_882_g_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


// To make the newform (type ModFrm), type "MakeNewformModFrm_882_4_g_be();".
// 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_882_4_g_be(:prec:=4)
    chi := MakeCharacter_882_g();
    f_vec := qexpCoeffs();
    Kf := Universe(f_vec);
    // SetVerbose("ModularForms", true);
    // SetVerbose("ModularSymbols", true);
    S := CuspidalSubspace(ModularForms(chi, 4));
    S := BaseChange(S, Kf);
    maxprec := NextPrime(997) - 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_882_4_g_be();".
// 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_882_4_g_be( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_882_g();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,4,sign)));
    Vf := Kernel([<5,R![1156, 408, 110, 12, 1]>,<11,R![4624, 272, 84, -4, 1]>,<13,R![126, 48, 1]>],Snew);
    return Vf;
end function;