
// Make newform 1521.4.a.h in Magma, downloaded from the LMFDB on 19 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_h();".
// 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_h();".
// 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 := [[1], [0], [-7], [-10], [22], [0], [-37], [30], [162], [113], [196], [13], [-285], [-246], [462], [537], [-576], [-635], [202], [1086], [-805], [884], [-518], [-194], [-1202], [429], [-1302], [1338], [-1034], [-1077], [-988], [-560], [519], [-348], [645], [2914], [-2079], [1700], [-3680], [-4146], [-3674], [-3283], [596], [-393], [3522], [2018], [160], [4072], [5794], [6482], [-6890], [-2466], [-3617], [-4860], [-565], [498], [-5546], [-2256], [2309], [-5833], [1650], [-2991], [-2422], [3402], [2310], [257], [1028], [2487], [2850], [-2018], [5287], [7278], [-4202], [-1583], [-2052], [6872], [11653], [6134], [10795], [-8489], [-1496], [-11695], [-10590], [-13949], [-10726], [-16228], [-7538], [15539], [-4811], [562], [-4914], [3600], [-17130], [11838], [8976], [-1682], [-15167], [6783], [-13918], [-1335], [-3806], [1905], [4800], [-14678], [-586], [8939], [13792], [-9569], [5192], [-3677], [-10960], [-26027], [-17681], [3192], [7580], [27707], [11216], [2536], [-17730], [-18920], [5241], [20467], [70], [-6432], [-6666], [14054], [-71], [-3936], [34202], [-27363], [21776], [-2484], [32906], [-3914], [33038], [-17586], [-18314], [-42068], [4282], [-40221], [-7084], [-17338], [35496], [14992], [-20659], [-28716], [13377], [27419], [2422], [34522], [13733], [22759], [-2168], [-15888], [-11628], [12584], [17184], [-12777], [9191], [-50498], [-1560], [21498], [-418], [-18132], [-12501], [43708], [-39614], [-36503], [-60450], [9261], [-63114], [-3827], [-34684], [49923], [-7056], [-58893], [50322], [26409], [-58592], [4254], [1514], [-43008], [2417], [-3654], [-65638], [34137], [-19514], [28364], [-61741], [10890], [-15629], [15394], [55750], [-58966], [-26214], [-20802], [-67373], [-63407], [64379], [2100], [-30687], [-12808], [-11615], [-64141], [34302], [1913], [30394], [31396], [36695], [74440], [-30493], [34398], [-85786], [-7646], [50608], [-15337], [-33304], [-20701], [59172], [-91506], [40529], [70252], [-79781], [-29828], [-19374], [-56915], [71031], [-19060], [90316], [-55876], [-13959], [94820], [-105820], [51147], [-80192], [-13020], [-50146], [14197], [-86344], [4762], [-60962], [-112058], [-49222], [100701], [5367], [5964], [-76844], [-81660], [52646], [-112818], [-71674], [32930], [22776], [99307], [44058], [-11150], [-85282], [-102034], [52878], [-99377], [-11772], [41980], [43253], [-84779], [-70401], [-12542], [-56230], [20055], [40382], [10249], [49133], [42324], [-15185], [20348], [66726], [140698], [-39164], [-80015], [71623], [73386], [51102], [30008], [32956], [80057], [48566], [-84002], [15451], [33801], [72746], [83711], [54678], [33234], [-22329], [96528], [-87914], [128905], [102312], [163921], [-123436], [8342], [154079], [-36647], [-65570], [52652], [115730], [89151], [139550], [-371], [144516], [-71307], [107294], [-129786], [-61049], [10030], [-91342], [-5773], [109926], [43174], [-87554], [-71973], [157682], [96826], [135105], [138438], [-105618], [98691], [-190426], [20302], [5314], [173013], [161909], [116877], [65098], [-226], [21454], [-27230], [182125], [85766], [154079], [-980], [-123954], [198871], [-68467], [-46726], [-135946], [-42282], [-112595], [-110384], [192698], [-140954], [-202566], [98454], [159405], [8738], [110501], [-119245], [-63310], [-123972], [14523], [110878], [90957], [-83445], [118224], [17844], [-79424], [167446], [237606], [-6728], [-99573], [61188], [-132270], [-133982], [221625], [-54210], [94485], [190302], [79458], [-36398], [-256769], [138903], [139734], [86283], [-236458], [-162738], [-84404], [-103010], [109290], [171150], [42632], [-67636], [-69133], [49961], [202046], [95702], [66704], [-225842], [-105646], [20874], [-153788], [-80331], [-63670], [88278], [110032], [106239], [157841], [-234526], [9610], [71135], [45486], [-49152], [-10213], [167105], [132882], [-139532], [94107], [-46851], [-206380], [-42742], [-266670], [-30004], [170689], [67809], [142262], [-218806], [38582], [163861], [-26868], [-250161], [307022], [6900]];
    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_h();".
// 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_h(: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_h();".
// 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_h( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_1521_a();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,4,sign)));
    Vf := Kernel([<2,R![-1, 1]>,<5,R![7, 1]>,<7,R![10, 1]>],Snew);
    return Vf;
end function;