
// Make newform 1386.4.a.b in Magma, downloaded from the LMFDB on 19 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_1386_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_1386_4_a_b();".
// 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_1386_4_a_b();".
// 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_1386_a();"
function MakeCharacter_1386_a()
    N := 1386;
    order := 1;
    char_gens := [155, 199, 1135];
    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_1386_a_Hecke(Kf)
    return MakeCharacter_1386_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], [0], [-3], [7], [11], [-16], [-6], [14], [51], [-54], [95], [-193], [-102], [284], [72], [102], [63], [-790], [-433], [-135], [-238], [770], [1008], [639], [11], [-1692], [-532], [-1350], [1616], [-597], [302], [462], [-3021], [-1594], [-2814], [2450], [3899], [-124], [30], [-2268], [-1731], [785], [-51], [1226], [-2754], [-1744], [-4210], [-5497], [3732], [-3883], [222], [2304], [-1276], [255], [-2370], [-2466], [-6774], [-4876], [-3268], [-5928], [2672], [-5322], [2816], [-8760], [1337], [-4269], [-4843], [-3790], [-1722], [-3166], [4167], [-1380], [755], [-8596], [-11287], [1185], [11643], [-3382], [7566], [86], [12672], [12398], [-12372], [-14929], [9146], [11253], [-6711], [16952], [5850], [-4957], [-3753], [16272], [10613], [-11562], [-9604], [15852], [-15669], [2817], [-12364], [2372], [-12832], [14166], [18540], [21894], [23204], [18479], [1116], [3654], [-10248], [18848], [-13198], [1712], [-20994], [25691], [25445], [-19095], [-24391], [5691], [25749], [9030], [28895], [16856], [19470], [14820], [-18451], [-6288], [-17629], [-2067], [11315], [-23644], [-862], [-9072], [9407], [-15946], [8220], [14372], [-14082], [-12742], [-7131], [-21966], [-10978], [27534], [-19591], [-12744], [23567], [-30219], [8654], [43284], [-21931], [22404], [4826], [-40203], [-41716], [44358], [29384], [-5160], [10436], [-40626], [-2608], [57150], [657], [-5508], [-1690], [-7263], [-1569], [9261], [19400], [-20740], [8960], [-29313], [-18636], [-41254], [-31482], [49316], [-33415], [45105], [-55438], [-58443], [-6010], [40466], [51884], [55224], [-33847], [-19404], [7353], [9618], [33212], [-13060], [-36376], [-5058], [37403], [29022], [-54385], [22806], [-31194], [-8277], [-4336], [25787], [-45486], [-29826], [-33930], [49868], [24446], [32330], [-73815], [68709], [-52423], [60174], [-17844], [8228], [76964], [-18753], [-22111], [-36924], [15360], [61817], [-84544], [26910], [-66684], [-33723], [-75460], [20666], [74211], [80147], [-38724], [48890], [38058], [-44856], [41108], [37788], [-29839], [86114], [-63592], [15540], [-97279], [92946], [89915], [-35958], [-62829], [-92181], [-81357], [50066], [12476], [-84175], [83622], [4830], [80840], [-30492], [98606], [35358], [-37984], [-68856], [60969], [-18154], [52716], [124290], [11246], [71738], [110973], [-49444], [812], [3072], [-64486], [74768], [130722], [107156], [-27435], [-18435], [-98782], [-74160], [-91381], [-21139], [103358], [-21388], [-78766], [39107], [-137001], [4202], [73154], [-131040], [114636], [50339], [36126], [85658], [-32002], [-142431], [7658], [46938], [2069], [-70056], [69771], [89436], [-24348], [76260], [12638], [-15336], [-8995], [111039], [68010], [145976], [83360], [122772], [-136474], [101844], [-156907], [56183], [5643], [8615], [-79944], [-52564], [47922], [69534], [-76554], [-120841], [-33882], [184742], [108105], [52854], [58598], [47154], [-79996], [35117], [-163818], [9119], [-92106], [-127915], [43733], [-7933], [-3438], [-86400], [35000], [147666], [21062], [-108486], [192854], [137823], [-120679], [-39630], [-50455], [13604], [31421], [-144831], [-75588], [88805], [-65691], [-208824], [-63343], [71489], [-75384], [203955], [-30472], [123455], [110070], [4322], [103004], [102252], [-23445], [25542], [-148248], [53889], [88454], [164106], [-223497], [-174336], [-48889], [-47563], [30744], [-68308], [90956], [-202311], [-171220], [68412], [-189408], [-171808], [-144427], [-140748], [85638], [6014], [-38256], [114752], [62526], [-111687], [-54430], [15966], [-53260], [222609], [103520], [119711], [-178810], [261381], [83153], [87519], [58233], [22955], [262191], [-51946], [140888], [174366], [-171364], [-62880], [-171574], [-135309], [246440], [-115866], [-259560], [-79510], [34922], [-258696], [-267556], [-255459], [85160], [31848], [-113220], [-138328], [231254], [-23535], [-43872], [-93721], [-16650], [-61698], [-39375], [114776], [4236], [110910], [-124810], [164175], [-10113], [-112464], [-173149], [117126]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_1386_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


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