
// Make newform 21.6.e.a in Magma, downloaded from the LMFDB on 25 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_21_e();"
// To make the character of type GrpDrchElt with Codomain the HeckeField, type "MakeCharacter_21_e_Hecke();"
// 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_e_a();".
// 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_e_a();".
// 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 := [1, -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_e();"
function MakeCharacter_21_e()
    N := 21;
    order := 3;
    char_gens := [8, 10];
    v := [3, 2];
    // 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;

// To make the character of type GrpDrchElt with Codomain the HeckeField, type "MakeCharacter_21_e_Hecke();"
function MakeCharacter_21_e_Hecke(Kf)
    N := 21;
    order := 3;
    char_gens := [8, 10];
    char_values := [[1, 0], [-1, 1]];
    assert UnitGenerators(DirichletGroup(N)) eq char_gens;
    values := ConvertToHeckeField(char_values : pass_field := true, Kf := Kf); // the value of chi on the gens as elements in the Hecke field
    F := Universe(values);// the Hecke field
    chi := DirichletCharacterFromValuesOnUnitGenerators(DirichletGroup(N,F),values);
    return chi;
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 := [[2, -2], [0, 9], [-11, 11], [126, 7], [0, -269], [-308, 0], [0, -1896], [164, -164], [3264, -3264], [2417, 0], [0, -2841], [11328, -11328], [-16856, 0], [-7894, 0], [-21102, 21102], [0, 29691], [0, 8163], [-15166, 15166], [0, 32078], [-38274, 0], [0, -34866], [-13529, 13529], [-68103, 0], [114922, -114922], [154959, 0], [0, 107570], [-8936, 8936], [-193667, 193667], [0, -205110], [46664, 0], [-304365, 0], [13303, -13303], [0, 398262], [230286, 0], [97134, -97134], [0, 29047], [0, 576500], [265232, -265232], [-363790, 0], [164846, -164846], [0, -30628], [-651392, 0], [757360, -757360], [0, 160339], [-61738, 0], [0, -370908], [217450, 0], [589771, 0], [0, 387045], [232732, -232732], [-42096, 42096], [-313416, 0], [0, -857807], [454517, 0], [-878182, 878182], [0, -1960934], [0, 1053795], [105059, -105059], [0, 427592], [638878, 0], [0, 2451422], [1716169, 0], [1808968, 0], [0, -1521456], [1348399, -1348399], [49695, -49695], [-1587836, 1587836], [214825, 0], [0, -2588596], [-24878, 0], [0, 1736008], [-862426, 862426], [0, -3115419], [1796944, -1796944], [3454656, 0], [-2155038, 2155038], [0, 462774], [4066212, -4066212], [5068632, -5068632], [0, -2877335], [-3413424, 0], [-1307374, 0], [0, -1935468], [516670, 0], [-2915301, 2915301], [1783791, -1783791], [4001580, 0], [1167665, -1167665], [3613582, 0], [-1801112, 0], [2369748, -2369748], [0, -518146], [0, -2826127], [9347473, 0], [-8171854, 8171854], [-7377132, 0], [326315, -326315], [0, 2167026], [723404, -723404], [-5999642, 5999642], [7015696, 0], [0, 8918725], [0, 13348231], [-1102148, 1102148], [0, -1893482], [0, -2829511], [-10679893, 0], [14699688, -14699688], [0, -8491626], [8629467, 0], [-10580735, 10580735], [0, -3847840], [-15100062, 0], [0, -9931018], [-9252241, 0], [0, -5004282], [12613706, 0], [0, 12538310], [8660661, -8660661], [7941768, 0], [0, 2114162], [-442307, 0], [-10756071, 10756071], [0, 11488593], [-10138832, 10138832], [-19683879, 0], [20171706, -20171706], [-4157346, 4157346], [-15443323, 0], [6204138, -6204138], [0, -21898446], [-2754834, 0], [-12912079, 12912079], [-26431522, 0], [-12221388, 12221388], [6106541, 0], [0, 3022202], [0, 20828522], [23232839, 0], [0, -10866768], [-22263178, 0], [12388091, -12388091], [0, -1694808], [-378495, 0], [0, 10428748], [30408224, 0], [28031486, 0], [0, -18802958], [-7863230, 7863230], [11285766, -11285766], [-13415008, 13415008], [31854690, 0], [-30592168, 0], [4637720, -4637720], [0, -604376], [-24405918, 0], [-36709496, 36709496], [22908854, -22908854], [-5996109, 0], [0, 11651585], [1109264, -1109264], [10574314, 0], [6326665, 0], [-39239463, 39239463], [0, 1550742], [0, 48748428], [0, 1925515], [0, -54256390]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_21_e_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


// To make the newform (type ModFrm), type "MakeNewformModFrm_21_6_e_a();".
// 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_e_a(:prec:=2)
    chi := MakeCharacter_21_e();
    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_e_a();".
// 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_e_a( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_21_e();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,6,sign)));
    Vf := Kernel([<2,R![4, -2, 1]>],Snew);
    return Vf;
end function;