
// Make newform 2166.4.a.f in Magma, downloaded from the LMFDB on 20 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_2166_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_2166_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_2166_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_2166_a();"
function MakeCharacter_2166_a()
    N := 2166;
    order := 1;
    char_gens := [1445, 1807];
    v := [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_2166_a_Hecke(Kf)
    return MakeCharacter_2166_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 := [[2], [-3], [2], [-21], [-40], [-17], [36], [0], [74], [-100], [-103], [-187], [128], [121], [410], [-230], [744], [-277], [231], [-578], [609], [-1259], [-696], [612], [1550], [-418], [1701], [1690], [-42], [1828], [2124], [1206], [310], [1335], [-1116], [-1328], [-2537], [101], [-240], [1428], [2046], [-1510], [2290], [-627], [-748], [2171], [3177], [-4751], [4706], [-4807], [546], [3966], [5051], [4370], [670], [-7002], [-2484], [3736], [574], [9278], [5244], [666], [-2388], [9460], [10682], [4542], [-1313], [-6649], [-10706], [3799], [-10926], [204], [-9119], [13102], [-2841], [14024], [-9438], [8267], [-6664], [-2038], [8832], [9826], [8850], [-10643], [10387], [6652], [16976], [-6367], [-14008], [14161], [-12356], [-10094], [7580], [-19462], [14383], [1716], [-3954], [11356], [-14627], [-6355], [-8077], [-10346], [6126], [-1444], [20951], [8678], [-11700], [-6968], [10248], [1505], [-22187], [9682], [2844], [-27809], [-28477], [-10926], [-5465], [-25876], [-21594], [-6392], [-346], [-20603], [14342], [15530], [-24108], [-2854], [18033], [11878], [17323], [-5082], [1673], [-31172], [2941], [6061], [26004], [-20837], [-33162], [497], [-26466], [27090], [29376], [5916], [40784], [-5768], [30735], [19234], [-22293], [-29226], [-8315], [-45930], [-26815], [-37662], [14881], [14684], [51484], [38568], [-37591], [8138], [35543], [14696], [47286], [-20816], [-7839], [-38390], [41686], [-34692], [32375], [29131], [-25642], [-18480], [-15840], [37753], [28422], [16361], [-49035], [17598], [61361], [-45246], [-44947], [45899], [-41132], [-65128], [32075], [23858], [-21006], [-38970], [-43382], [52007], [-11642], [-8508], [28127], [-69676], [-17437], [10506], [20830], [-35472], [-1978], [73629], [-31480], [55050], [-32532], [17432], [27839], [8959], [-61834], [-37804], [36164], [-50412], [-16404], [51037], [89483], [49058], [66753], [34740], [62442], [49209], [46272], [41574], [22400], [59628], [-41253], [82159], [-8258], [-67021], [-47768], [-837], [51778], [22656], [75736], [-27180], [-27003], [22351], [35880], [104816], [55076], [-49404], [62409], [29576], [32600], [-15522], [-27946], [70696], [-8537], [88405], [15180], [-89040], [-21704], [85624], [36481], [-19800], [10009], [22370], [63416], [-28850], [55652], [22144], [-1047], [23305], [56880], [103555], [-57715], [-110364], [-44475], [-34155], [-88596], [79259], [-78768], [-17400], [45445], [74004], [-1899], [33368], [-134094], [-76960], [-64169], [-3093], [65634], [-122443], [82553], [-58386], [43962], [101752], [-38286], [11949], [67744], [12240], [-32595], [49842], [-16595], [59218], [-4014], [125610], [-89418], [-62342], [-15155], [77600], [29176], [-142084], [-59128], [2289], [-87866], [159598], [164881], [60250], [87963], [-70481], [106554], [123115], [4638], [-32209], [-4958], [173616], [106562], [-89919], [55134], [-100758], [140724], [-137494], [71225], [118120], [-7675], [-94509], [-31872], [101833], [132342], [7383], [105485], [103145], [87128], [-21340], [150709], [-125028], [139580], [139594], [-21012], [71438], [83649], [163914], [-29815], [101179], [2150], [109086], [40728], [74572], [-96236], [-38722], [185838], [76044], [-87806], [5682], [-141931], [-26921], [21032], [-88872], [147862], [8906], [31100], [-63474], [-212096], [-184548], [-150530], [-57850], [-149136], [39592], [-139636], [7510], [-79968], [114693], [-63679], [75246], [-180363], [-35128], [-64534], [-69019], [222221], [53710], [-90748], [213758], [176454], [98679], [-68250], [-66070], [-20096], [61326], [157239], [-116834], [201049], [-145322], [-70772], [-49342], [231475], [48646], [-245670], [130891], [12546], [20830], [8897], [36694], [-81319], [-89868], [-244931], [99934], [-130320], [125738], [-135750], [-196447], [-179075], [210454], [56420], [125540], [-57943], [32514], [199598], [-287983], [-251606], [-135480], [-144390], [-44204], [-305664], [192638], [93462], [128449], [274100], [4544], [-130123], [-232944], [-130724], [19966], [-72421], [-1114]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_2166_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


// To make the newform (type ModFrm), type "MakeNewformModFrm_2166_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_2166_4_a_f(:prec:=1)
    chi := MakeCharacter_2166_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_2166_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_2166_4_a_f( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_2166_a();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,4,sign)));
    Vf := Kernel([<5,R![-2, 1]>,<13,R![17, 1]>],Snew);
    return Vf;
end function;