
// Make newform 33.4.e.a in Magma, downloaded from the LMFDB on 23 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_33_e();"
// To make the character of type GrpDrchElt with Codomain the HeckeField, type "MakeCharacter_33_e_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_33_4_e_a();".
// 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_33_4_e_a();".
// 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 := [1, -1, 1, -1, 1];
        Kf := NumberField(Polynomial([elt : elt in poly]));
        AssignNames(~Kf, ["nu"]);
    end if;
    Rfbasis := [Kf.1^i : i in [0..Degree(Kf)-1]];
    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_33_e();"
function MakeCharacter_33_e()
    N := 33;
    order := 5;
    char_gens := [23, 13];
    v := [5, 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;

// To make the character of type GrpDrchElt with Codomain the HeckeField, type "MakeCharacter_33_e_Hecke();"
function MakeCharacter_33_e_Hecke(Kf)
    N := 33;
    order := 5;
    char_gens := [23, 13];
    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 := [[4, -2, 4, 0], [0, 0, 0, -3], [1, -12, 12, -1], [0, 6, -25, 6], [-15, 35, -25, 41], [-22, -11, -22, 0], [78, -21, 21, -78], [-45, 45, 0, 83], [-15, 0, 44, -44], [0, 192, -30, 192], [-77, 33, -77, 0], [0, -22, -245, -22], [-34, 34, 0, -27], [-75, 0, -84, 84], [55, -55, 0, -200], [143, -241, 143, 0], [0, -451, 418, -451], [717, -438, 438, -717], [-243, 0, -561, 561], [-433, 708, -708, 433], [0, -318, 346, -318], [-702, 637, -702, 0], [-35, -431, 431, 35], [-1279, 0, 154, -154], [891, -282, 891, 0], [74, 843, 74, 0], [0, 44, -618, 44], [2, -2, 0, 1561], [132, 0, 1056, -1056], [-836, 836, 0, 1027], [806, -1592, 1592, -806], [432, 0, 193, -193], [1309, -671, 671, -1309], [0, 1335, 339, 1335], [-313, 699, -699, 313], [-316, 316, 0, 2050], [990, -990, 0, -1986], [-2409, 1141, -2409, 0], [269, -1262, 269, 0], [-1477, 1477, 0, 843], [1078, -1078, 0, -21], [-225, -1975, 1975, 225], [0, 616, 2035, 616], [4088, -4349, 4349, -4088], [1220, 0, 889, -889], [-297, 0, -792, 792], [1741, 807, 1741, 0], [-1408, 1408, 0, -381], [0, 1792, -1039, 1792], [-1628, -3874, -1628, 0], [-987, 1830, -987, 0], [299, -299, 0, -6172], [1566, 0, -2742, 2742], [957, -1443, 957, 0], [0, 451, -990, 451], [-645, 0, 63, -63], [-5010, 5934, -5934, 5010], [0, -1291, -2748, -1291], [5797, -4785, 5797, 0], [-4580, 4594, -4594, 4580], [-894, 894, 0, 7520], [0, 778, -1115, 778], [-1476, 0, 843, -843], [-4576, 4576, 0, -965], [117, -2801, 2801, -117], [2222, -1801, 2222, 0], [8329, 0, 5016, -5016], [0, 3310, -3756, 3310], [3808, -6356, 6356, -3808], [-1590, 1590, 0, -8115], [6564, 0, -1584, 1584], [0, -2178, -5082, -2178], [0, 2937, -6810, 2937], [-10429, 0, 548, -548], [-2551, -4863, 4863, 2551], [-2508, -2613, -2508, 0], [0, 473, 2581, 473], [5115, 0, 4950, -4950], [12324, -12357, 12357, -12324], [1580, -3188, 1580, 0], [873, 0, -6941, 6941], [-1397, 1397, 0, -4639], [6727, -10431, 6727, 0], [0, -10252, 6774, -10252], [-4197, 0, 5796, -5796], [-10978, 10978, 0, 5704], [2200, 5554, 2200, 0], [-962, -313, 313, 962], [-1545, 0, 243, -243], [10496, 0, -2409, 2409], [-7451, 3997, -3997, 7451], [-8472, 13055, -13055, 8472], [-4092, 4092, 0, 14397], [0, -3319, 10575, -3319], [0, -3861, 7906, -3861], [3734, -3734, 0, -13010], [6919, -6919, 0, -8956], [0, 7832, 692, 7832], [5952, -2001, 2001, -5952], [2717, -2574, 2717, 0], [10987, -10987, 0, -4410], [0, -2773, -2260, -2773], [-6700, -10785, -6700, 0], [6238, -6238, 0, 2766], [15515, 0, 7127, -7127], [-5050, -5444, 5444, 5050], [0, 0, -12561, 0], [8273, 0, -1945, 1945], [9794, 4880, -4880, -9794], [0, -3344, 2079, -3344], [7131, -5683, 7131, 0], [-2058, 2058, 0, 10632], [-13543, 0, 9526, -9526], [-9636, 9636, 0, 275], [0, -11, -8221, -11], [1375, -1375, 0, -6257], [-5614, 16273, -16273, 5614], [-9625, 23263, -9625, 0], [0, -10153, -13423, -10153], [29664, 0, 1050, -1050], [4246, 0, 20295, -20295], [2520, -6969, 2520, 0], [8276, -18784, 18784, -8276], [22275, 0, 12980, -12980], [-18986, 3302, -18986, 0], [-4031, 4031, 0, 19159], [-20254, 7857, -7857, 20254], [0, 13024, 3039, 13024], [-13868, 0, -21351, 21351], [0, 8846, -12266, 8846], [5046, 9489, 5046, 0], [2101, 4785, -4785, -2101], [6369, -6369, 0, -23229], [-9636, 28255, -9636, 0], [14442, -10191, 14442, 0], [-9235, 0, -30427, 30427], [-8822, 8822, 0, 23795], [-15468, 23802, -23802, 15468], [22494, -19656, 19656, -22494], [16783, -20496, 20496, -16783], [26751, -26751, 0, -26416], [0, 8374, -7707, 8374], [19338, -1871, 19338, 0], [-21238, -3890, -21238, 0], [0, 21813, -22703, 21813], [-20130, 20130, 0, 19035], [-7291, -8105, 8105, 7291], [-22279, 0, -32678, 32678], [-21105, 0, -6798, 6798], [-10558, 6290, -6290, 10558], [25628, -25628, 0, -13155], [-2064, 0, 18997, -18997], [3212, -3212, 0, 23650], [0, 31576, -24005, 31576], [23149, 644, -644, -23149], [-15004, 4477, -15004, 0], [-3405, 401, -401, 3405], [46809, -50703, 50703, -46809], [-2302, 33975, -2302, 0], [-26494, -14317, 14317, 26494], [-18021, 0, 15477, -15477], [0, 15238, -788, 15238], [14625, 0, -7787, 7787], [31020, -31020, 0, -4554], [9394, -623, 9394, 0], [0, 22132, 11694, 22132], [8433, 0, -19129, 19129], [0, -23075, 3920, -23075]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_33_e_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


// To make the newform (type ModFrm), type "MakeNewformModFrm_33_4_e_a();".
// 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_33_4_e_a(:prec:=4)
    chi := MakeCharacter_33_e();
    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_33_4_e_a();".
// 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_33_4_e_a( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_33_e();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,4,sign)));
    Vf := Kernel([<2,R![400, -200, 60, -10, 1]>],Snew);
    return Vf;
end function;