
// Make newform 21.6.a.c 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_6_a_c();".
// 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_6_a_c();".
// 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_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 := 6;
    raw_aps := [[5], [9], [94], [-49], [52], [-770], [-2022], [1732], [-576], [5518], [6336], [-7338], [-3262], [5420], [864], [4182], [-11220], [-45602], [1396], [18720], [46362], [97424], [-81228], [-3182], [4914], [-166354], [157160], [-6764], [178398], [-45134], [-205056], [72964], [-94182], [-47796], [-124266], [-446296], [-159746], [247252], [-684488], [-610474], [662252], [154630], [486904], [620546], [-236570], [82104], [99892], [-186704], [336372], [-926314], [1257114], [-347016], [99170], [344428], [295130], [-1272464], [276774], [-1289936], [1716550], [-1472182], [1028812], [-1186066], [1518916], [212808], [-1894], [-1578978], [-3394708], [2027314], [3488852], [965566], [1153930], [1611104], [3677472], [649766], [320700], [-2361888], [-3533898], [4048110], [2076450], [2574314], [848148], [1436822], [2354376], [-3788078], [-3643224], [2483892], [-2631774], [-1161302], [2813846], [6842992], [3343140], [-4282480], [-8931752], [2753060], [4804084], [-6024648], [-8429866], [9250578], [5844940], [9225326], [-6443372], [3742126], [-14638364], [14180458], [-1251604], [5943778], [-6461924], [-2346054], [-13495920], [3878490], [-533488], [5146102], [-2378598], [16002316], [12345944], [-3437550], [-16219084], [11992856], [1580094], [6983580], [3696022], [1846882], [-7685010], [7121796], [-3237868], [-7391634], [-5333610], [-11456400], [-24953992], [-14339826], [922932], [-9389952], [-408032], [25960454], [18355410], [-747166], [20269166], [-4699820], [584710], [-16401270], [-304948], [34342838], [15668376], [-29688636], [-23070802], [23264248], [19151510], [-5346830], [39585836], [-25028376], [-5025890], [-26049030], [-6824620], [23383496], [39595916], [-4675704], [-4925944], [-32368774], [-33233718], [-26642650], [31466300], [-13451590], [-28496312], [18185772], [32094114], [15615384], [48449936], [-45433562]];
    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_6_a_c();".
// 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_6_a_c(:prec:=1)
    chi := MakeCharacter_21_a();
    f_vec := qexpCoeffs();
    Kf := Universe(f_vec);
    // SetVerbose("ModularForms", true);
    // SetVerbose("ModularSymbols", true);
    S := CuspidalSubspace(ModularForms(chi, 6));
    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_6_a_c();".
// 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_6_a_c( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_21_a();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,6,sign)));
    Vf := Kernel([<2,R![-5, 1]>],Snew);
    return Vf;
end function;