
// Make newform 21.10.a.b in Magma, downloaded from the LMFDB on 19 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_21_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_21_10_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_21_10_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
        poly := [-588, -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_21_a();"
function MakeCharacter_21_a()
    N := 21;
    order := 1;
    char_gens := [8, 10];
    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_21_a_Hecke(Kf)
    return MakeCharacter_21_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 := 10;
    raw_aps := [[5, -1], [-81, 0], [550, 70], [-2401, 0], [-73148, 550], [43318, -108], [-103590, -22662], [-110468, -288], [-1020768, 5754], [4890142, -24032], [198912, -193824], [-6951450, -56916], [-21051550, -259450], [-1850020, -1063872], [-24491376, 704268], [-55575594, 2170836], [-93087756, -2201292], [7936198, 3849696], [29648404, 10977588], [-189205968, -3632250], [94807674, 2170548], [-35412976, -1766196], [86737404, 14519424], [32720674, -22510394], [875388930, -23923764], [-94812298, 65626322], [794953640, -2818296], [28974100, 11290150], [-152645922, -82463616], [78982162, -32872424], [-866624304, 100079244], [417746188, -95774888], [-3493207686, 122916084], [-942247164, 282576960], [-5513460906, -8915628], [1593138488, -191841984], [8700166294, -176236488], [4737841732, -166563684], [3173662024, -523909940], [-9460359826, 289085306], [1115595500, 477892250], [14256533758, -290215188], [-16908468152, -706614866], [-10235992654, 102579156], [24772747462, 659667604], [-12937973592, 338596200], [17154603988, -84355992], [-35315685680, 791185752], [-12310409028, 2837698548], [-5850316018, 1811003436], [23994569370, 351870564], [19150918104, -1188963654], [43836672914, 1059732108], [-71751337820, 845701684], [34634396234, -1648676218], [-22731799040, 3138932266], [-63972621954, 639639126], [-89292819872, -2613735432], [4276309510, 845313300], [84228289130, 2382956636], [11273011396, -2496997728], [-33519963946, -5239537822], [155634338044, 1298376072], [-99024073080, -4311449844], [257124670730, -317073024], [-39338423874, -1738820412], [-80460976324, 5100771816], [-95344720142, -8167147416], [-116766383356, 5218023386], [-247543903850, -9096315552], [-175757315270, -7484081750], [10079956832, 104064902], [343520968224, -3878044848], [-322557853834, -4000608000], [-333621005076, 2337271704], [214170663504, 5937207384], [5986617990, -32913393456], [416774382438, 12264165000], [-157611845598, 9055989084], [-147587790598, 12718280124], [346817662476, 7269454500], [17608446614, 24974982036], [-408948707256, -17890993698], [311927006818, -26415196776], [-624402770328, 16774690752], [-791246822028, -20309886582], [905482376034, -3077160816], [-549351298742, 34411337856], [549547508366, 31766984786], [-93442589168, 11320922916], [1114898188332, -35180356452], [803488385120, 22369012532], [1175503871272, -8907153624], [486641509892, -36170286550], [-1217718470828, 21932170200], [821298242328, 22804665936], [-1698258643474, 8578327958], [376387138482, -53679688878], [-1788067780844, 14146290000], [-2616731690578, -4692515868], [76580708740, 126367456932], [-3037172642834, -51209780804], [1857235255804, -12254541284], [-1651541606870, -10114841564], [-1402048766596, -100398468132], [-1769795353406, 93180510648], [-346339897484, 76560657724], [-468985506390, -19581024582], [126387152160, 28164724206], [1220843506938, 61784073144], [2116670604080, 34009126416], [399371318326, -97174274832], [-5825468929734, -4480293636], [652039852852, -744092712], [200552319704, 40639774860], [-2373828179406, -167055390924], [-775183080580, -167517407088], [1049794001864, -114472487500], [558519938094, -191045350536], [-3195130882932, 134825716350], [-1032177547874, 42768073656], [-430290439790, 49191291900], [-757819198842, 207448085946], [-3285303658524, -10580692878], [-5932650936868, 188212967784], [814695712350, -429452018256], [4962065142102, 279087801264], [-865132991424, 164044051776], [2515371397208, -365404250952], [5869005120630, -357426888300], [-3840078955548, 5099658444], [12408407907840, -21195901818], [944255577712, 352605924432], [1198870773734, -236594430000], [-1919710478190, -26704410762], [-231612347374, -7491474072], [6380545786454, -63661130950], [-14483944920836, 226247408856], [-10090148278562, 89021023738], [8864965789386, -170298175872], [-11493671758732, -42775376184], [-6970249003306, 394906173884], [-250178001096, -466669334556], [-16670750149212, -96459520938], [8394253788662, 114610013676], [1587204055144, 54387002356], [-6491788275634, -218088159264], [-2833185991790, 372325868518], [21136133554868, -60738684552], [15937046752152, -192687901998], [-2260152067298, -562581709200], [11745592888266, 425625280722], [-5199185413996, -131669180712], [7297492528184, 8611651916], [10282614339884, -1249948040256], [9725195516040, 114166335150], [6667266885512, 5314782888], [-8427238709158, -650567888674], [-9808536417078, -258021203112], [-1183810213234, 518793202550], [-8070535671796, -278240144794], [-9951701528806, 175576993232], [7712575886104, 1200445286796], [-6353188994316, 158147295216], [1800165039714, -309126447540], [-13725161296824, -396241999320], [-31127269791472, 107354513160], [9525651748654, 571789589688]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_21_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


// To make the newform (type ModFrm), type "MakeNewformModFrm_21_10_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_21_10_a_b(:prec:=2)
    chi := MakeCharacter_21_a();
    f_vec := qexpCoeffs();
    Kf := Universe(f_vec);
    // SetVerbose("ModularForms", true);
    // SetVerbose("ModularSymbols", true);
    S := CuspidalSubspace(ModularForms(chi, 10));
    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_21_10_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_21_10_a_b( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_21_a();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,10,sign)));
    Vf := Kernel([<2,R![-568, -9, 1]>],Snew);
    return Vf;
end function;