
// Make newform 65.4.a.d in Magma, downloaded from the LMFDB on 25 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_65_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_65_4_a_d();".
// 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_65_4_a_d();".
// 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 := [-3, 0, 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_65_a();"
function MakeCharacter_65_a()
    N := 65;
    order := 1;
    char_gens := [27, 41];
    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_65_a_Hecke(Kf)
    return MakeCharacter_65_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 := [[2, 1], [-5, -3], [-5, 0], [-18, -4], [-5, 3], [13, 0], [-52, -2], [-7, 17], [17, 3], [58, -26], [-53, -121], [48, 180], [-60, -142], [-361, 33], [-230, 84], [-276, -50], [171, 211], [-134, -90], [-4, 378], [117, -139], [-168, 100], [-434, 330], [566, -92], [1090, -292], [-126, -672], [534, 588], [-1361, -239], [819, -495], [418, -792], [152, 414], [-511, 511], [-720, 1128], [930, -908], [-770, -602], [254, -776], [-2433, -1], [790, 56], [-380, 30], [802, 1072], [-1256, 966], [468, -1876], [-914, -1638], [-1796, 84], [838, -864], [186, -1140], [-2540, -320], [-832, -1088], [-882, 1292], [-1848, 1898], [-248, 606], [-2686, 904], [105, 193], [-1052, 1026], [3770, 158], [-3522, 96], [-2751, 2555], [4754, 2280], [4027, 803], [-364, -1182], [-448, -1634], [2083, 3301], [2328, -672], [4814, 1792], [-1122, -5358], [3896, -4110], [-3708, 1260], [3197, 1837], [8788, -1014], [4097, 2767], [-3880, -1414], [-9644, 1740], [6157, 321], [-9109, 2033], [-7034, -1640], [1695, -1333], [-3054, 3736], [6654, -1888], [-5116, 784], [-2598, 3340], [688, 7126], [7550, 2834], [-13690, -796], [6971, -2993], [-5874, 2324], [-5064, 3960], [3395, -8531], [11748, -946], [13070, -1204], [-9444, 406], [-496, -5146], [13315, 697], [-7027, -2863], [2248, 1350], [3334, -1206], [733, 6557], [8679, 2441], [-6912, 1122], [-12634, 3678], [-13457, -3255], [-11716, 2042], [-12045, -1479], [312, -9412], [-18329, 2385], [-11678, -5370], [16598, -570], [-3960, 1172], [-8004, -11274], [12702, 4636], [-8198, 1182], [4786, -4024], [-9715, 5571], [-270, 12032], [-11018, -4960], [-3221, 619], [4149, 4185], [4542, -4844], [14746, -7632], [-11943, 3767], [-3126, -4352], [-6054, -13922], [17846, -276], [1792, -2038], [4460, 12862], [25444, 546], [5529, -17599], [3842, -9056], [7424, 10126], [-25476, -1696], [-14289, 1477], [8862, -15936], [1805, -823], [-12570, 9116], [-3714, 1254], [2460, -15838], [2218, 9688], [17742, -1356], [9612, 14184], [1706, -7628], [-6086, 2044], [21834, -3874], [3281, 8377], [6562, 2844], [-29781, 2965], [-5154, 20180], [-2746, -22046], [-14057, 10599], [8580, -16340], [2926, -15696], [8930, -410], [25538, 1488], [16858, -13824], [15610, 11250], [-29337, -6407], [13467, 6401], [4307, 905], [-6004, 21040], [4382, -1654], [-15252, 6670], [46670, 1652], [-32060, -10982], [-17910, 15704], [7786, -964], [-13768, -4706], [28284, 9396], [-14520, 12884], [-14784, -15066], [-11472, 19768], [-15616, 4914]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_65_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


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