
// Make newform 1521.4.a.i in Magma, downloaded from the LMFDB on 18 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_1521_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_1521_4_a_i();".
// 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_1521_4_a_i();".
// 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_1521_a();"
function MakeCharacter_1521_a()
    N := 1521;
    order := 1;
    char_gens := [677, 847];
    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_1521_a_Hecke(Kf)
    return MakeCharacter_1521_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 := [[3], [0], [-9], [-15], [-48], [0], [-45], [-6], [162], [144], [-264], [-303], [-192], [97], [111], [414], [522], [376], [36], [357], [1098], [-830], [-438], [-438], [852], [396], [-182], [612], [-1083], [-90], [-2086], [1467], [-414], [-2419], [-930], [1683], [1874], [1194], [-2388], [-1566], [-657], [1222], [-1260], [-342], [81], [-1996], [2833], [3507], [-228], [-5493], [3627], [6075], [-210], [7092], [-5805], [-792], [-5472], [-2331], [1384], [-4062], [3764], [4227], [-306], [-2106], [10051], [-2154], [-10770], [2171], [7047], [6873], [-9318], [4128], [-2536], [-92], [-10182], [-579], [2106], [1974], [11886], [-1254], [-5823], [7341], [-7485], [15203], [1762], [7317], [-5016], [-9870], [-14541], [2112], [3276], [-15453], [3660], [747], [15804], [12078], [16110], [-3915], [16184], [7923], [-14389], [-10383], [-16425], [12213], [6383], [-6426], [-21330], [12084], [-2394], [-21971], [-15406], [9630], [14748], [3672], [19875], [1710], [16452], [25902], [-18108], [32904], [-15318], [7729], [-19242], [22518], [-9168], [1170], [1662], [30960], [8372], [2739], [6756], [29643], [-18128], [-6410], [28290], [27960], [-5649], [-756], [31194], [-17055], [-35520], [1095], [2554], [21522], [13124], [-23424], [-31077], [-19422], [1744], [19179], [-29217], [-15633], [-30589], [25884], [12305], [-29772], [47644], [21924], [32398], [2097], [-20016], [24993], [40959], [48933], [47388], [16803], [-57526], [-25000], [21222], [50076], [-60264], [-13977], [6156], [12552], [49192], [33021], [30000], [-22314], [-16256], [-6878], [-51708], [54756], [-29612], [26394], [41760], [25110], [40066], [-67578], [-6504], [-63048], [-7099], [44232], [57215], [11418], [-63459], [-34938], [-69060], [-5038], [36342], [-8478], [75207], [35606], [16611], [-67039], [-17208], [-30510], [13971], [-43380], [-47928], [50303], [-27358], [71316], [-41042], [-16794], [33117], [88146], [-31952], [-2313], [19776], [-45519], [63920], [26955], [71556], [88563], [33741], [-42302], [45657], [-23958], [-23290], [32364], [-105878], [-20788], [85368], [28197], [14159], [22143], [98850], [5685], [-79956], [88596], [-7602], [97481], [5524], [-25821], [103980], [95706], [79779], [84660], [93516], [67554], [87525], [-12672], [-18687], [-24901], [80442], [43032], [75490], [-6126], [-14274], [109506], [95888], [-107145], [-91017], [-106432], [43968], [25765], [-90981], [-109950], [69756], [39852], [-62062], [69702], [-26004], [-142058], [-75206], [141840], [39474], [-90942], [-12144], [-93879], [55134], [-56697], [-125712], [-13302], [-60678], [-59670], [111373], [73233], [56457], [-78786], [-44784], [108135], [-124212], [-49632], [-43252], [-86580], [-25036], [15372], [-26415], [-6036], [-168617], [-80385], [-37316], [-45684], [-111229], [-88494], [-94653], [108344], [-134988], [-60424], [94446], [79797], [-119979], [-30980], [170757], [88945], [-6390], [-87876], [-182964], [37602], [-106652], [28128], [4590], [98211], [-125352], [-48998], [-161532], [-24732], [113418], [-47826], [-77562], [-174384], [-55550], [-148614], [-28512], [130056], [-5670], [-9096], [-176442], [31876], [48369], [-184599], [9159], [77578], [-44883], [93483], [107746], [143310], [195225], [-60300], [6834], [10812], [4995], [-197152], [153992], [93078], [-112383], [-43908], [-59319], [-71847], [84390], [36666], [-9684], [-19146], [-93476], [-147967], [101025], [106104], [-66287], [228645], [106567], [144099], [172170], [105676], [-50740], [-85848], [-132138], [2802], [-104922], [104578], [51630], [-93780], [96180], [16332], [146742], [119331], [-77031], [-79378], [-233562], [218934], [-214998], [193833], [-118410], [141361], [136359], [-93575], [115671], [122031], [-140060], [-255225], [-110571], [-70506], [-248577], [185010], [85269], [139138], [17391], [74340], [-119880], [216948], [109681], [-117558], [-142479], [-208631], [-161687], [79038], [149337], [192220], [-92556], [-215442], [-259110], [138810], [21336], [113085], [-101910], [149142], [18135], [-88344], [-226800], [-79884]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_1521_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


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