
// Make newform 40.6.a.d in Magma, downloaded from the LMFDB on 16 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_40_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_40_6_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_40_6_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 := [-32, -1, 1];
        Kf := NumberField(Polynomial([elt : elt in poly]));
        AssignNames(~Kf, ["nu"]);
    end if;
    Rf_num := [[1, 0], [-2, 4]];
    Rf_basisdens := [1, 1];
    Rf_basisnums := ChangeUniverse([[z : z in elt] : elt in Rf_num], Kf);
    Rfbasis := [Rf_basisnums[i]/Rf_basisdens[i] : i in [1..Degree(Kf)]];
    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_40_a();"
function MakeCharacter_40_a()
    N := 40;
    order := 1;
    char_gens := [31, 21, 17];
    v := [1, 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_40_a_Hecke(Kf)
    return MakeCharacter_40_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 := [[0, 0], [-6, -1], [25, 0], [26, -3], [280, 6], [694, 12], [74, -84], [-500, 36], [-1226, 123], [670, 312], [-1124, 54], [-2970, 216], [11538, -156], [8842, 339], [-1454, 885], [-2706, -540], [31292, -504], [7054, -1104], [-42706, -543], [23604, -546], [-33726, 1308], [-32952, 2508], [54362, 711], [-27510, -1464], [73834, -4620], [51366, -4440], [-65474, -1305], [79230, 405], [13934, -1080], [137538, 4296], [-167886, 2877], [52920, -102], [154666, -4704], [36588, -11304], [-74986, 10104], [-396652, -2586], [215038, 6456], [43186, 19707], [107730, -13383], [10574, 7608], [172868, -9312], [8966, -25104], [203924, 2706], [565210, 13116], [148030, 5892], [-11800, -7968], [-181480, 9954], [-504106, -17625], [-915282, -4911], [629014, -32160], [-468270, 6804], [-1086888, -19548], [812570, -13452], [109704, -27354], [-1011150, 33240], [1166886, 6939], [-1102786, -14400], [253764, -42270], [48086, 11736], [-377006, -2676], [-1505206, -25413], [176854, -3792], [792670, -60063], [132188, 108426], [-816502, 13464], [-2324922, -46020], [853008, 129810], [747490, 58056], [2295662, 2421], [901326, 33936], [1834034, -5304], [-2354136, -66192], [-2343910, -36855], [-943466, -169464], [2587060, -113580], [1578390, -100041], [2722518, -122208], [1311094, -17748], [1559442, -91296], [2117842, 154092], [1631364, -40632], [-2799258, -49224], [1048628, 148458], [-4041766, 996], [487384, -259968], [-1990014, -230025], [-1475958, 56004], [-4414838, 42192], [-896658, 321480], [-1912938, -250161], [726230, -233667], [-126056, 226188], [-2153246, 2049], [5953176, -22890], [9136364, 22692], [4722558, -148929], [-894194, 242208], [863530, -293760], [-980486, -125829], [4366558, 84600], [-3202402, 341001], [5625422, 281496], [-5540582, 404919], [-3069550, 11772], [4328272, -489222], [-2655262, 61680], [-5765130, 285777], [4761714, -168528], [-1938432, -475308], [-7688846, 35964], [-3513054, -15771], [409326, -123156], [-2196078, -38604], [2580980, 261564], [-3056492, 510582], [-932694, 208788], [-2045310, 708627], [717378, 10785], [4690214, -224700], [-14536284, -207072], [-9311994, -490992], [-9709926, 25932], [-8290130, -204948], [-9175398, 26811], [-2718400, -325674], [-10069186, -418392], [11308070, 542232], [9238072, 671844], [-9914086, -334827], [3190062, 867792], [18651684, -240072], [-13676658, 518703], [-14853060, -103266], [7460294, 563904], [-12292902, 252864], [-305774, -1078488], [9768126, 572004], [7263542, 369165], [-18748362, -170460], [-1946550, -906432], [2048408, -108234], [-3190538, 856920], [-14154746, 595443], [-13433538, -264435], [34042894, -5904], [3500912, 147252], [-26031538, -33180], [5365962, 174216], [18529676, 826560], [3257790, -371733], [-6742850, 48816], [-2831670, -129300], [20219498, 257247], [29830666, 117021], [-18213690, -570567], [31315564, 198534], [-22033168, -415476], [-15899798, 723396], [20110114, -16644], [-13538690, -337536], [13675294, 1424001], [31621050, 705744], [-27708822, 208677], [29794592, -333318], [-44622742, 341820], [-46413786, 48435], [-6065236, -1355586], [-14705298, -1198260]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_40_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


// To make the newform (type ModFrm), type "MakeNewformModFrm_40_6_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_40_6_a_d(:prec:=2)
    chi := MakeCharacter_40_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_40_6_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_40_6_a_d( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_40_a();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,6,sign)));
    Vf := Kernel([<3,R![-480, 12, 1]>],Snew);
    return Vf;
end function;