
// Make newform 14.8.a.b in Magma, downloaded from the LMFDB on 18 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_14_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_14_8_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_14_8_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
        Kf := Rationals();
    end if;
    return [Kf!elt[1] : elt in input];
end function;


// To make the character of type GrpDrchElt, type "MakeCharacter_14_a();"
function MakeCharacter_14_a()
    N := 14;
    order := 1;
    char_gens := [3];
    v := [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_14_a_Hecke(Kf)
    return MakeCharacter_14_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 := 8;
    raw_aps := [[8], [-66], [-400], [-343], [40], [-4452], [36502], [-46222], [-105200], [-126334], [-170964], [20954], [318486], [77744], [703716], [1603278], [-1171894], [-2068872], [-994268], [33280], [-2971454], [-2376168], [-2122358], [6920346], [4952710], [3230000], [-1799092], [-15642876], [-6318902], [-10228762], [6007248], [20639618], [4761986], [-5057234], [-27273642], [6489208], [-63081028], [83227112], [30691556], [-52733836], [84273852], [-103956300], [-124775256], [147588782], [155811574], [-13319260], [-204939964], [-68485816], [-19362654], [308156764], [35579722], [-230056104], [5034950], [-103282578], [-232282254], [416748960], [-314678668], [192136944], [-440393410], [359234794], [-81146726], [253416072], [872705554], [-671611384], [-192216410], [-133837218], [425298128], [1077025518], [-723763728], [-448132412], [-1499455086], [-256889896], [650424296], [-466127498], [285859576], [-1650749612], [1513041962], [-863794380], [114041742], [-1183278274], [227959494], [-39070050], [-2586200136], [-1788926482], [-458904888], [1384588684], [-2736106798], [2430532870], [3948844172], [2574529456], [298482218], [2620001836], [-4166618752], [2412999948], [-1040918876], [-1172730720], [-81381756], [-2774584074], [4992130226], [1630951038], [2009498976], [-4479586498], [-1507303010], [2330877174], [-2911013688], [-8648050526], [-6335131754], [-1701620838], [3019770024], [-5927081734], [-1456492416], [-6716071014], [7020270742], [5143517194], [-4415735056], [694175662], [9508087022], [-7732150044], [-5063212326], [8081130304], [6300892840], [-9626240902], [-9454287564], [248879244], [-3464118526], [5563217018], [8236965866], [5852120028], [15170608348], [-15596918520], [-6955729944], [-11780290968], [-6962001960], [-20711388486], [16539242782], [-23365045634], [7092628232], [-12302989254], [-3666498444], [-29660925226], [25127795602], [4577720206], [1931326360], [15809448468], [-2465355896], [25186129948], [10730560652], [-27933158466], [-19498261714], [-16355068552], [26887351186], [10891847946], [-39954189700], [13030553388], [40401450948], [19891905168], [-41099018216], [-19037392746], [39383049314], [-15918631144], [3572371568], [4051007578], [25579109264], [41032282794], [48727686138], [-6947617428], [-18356454496], [34495366616]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_14_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


// To make the newform (type ModFrm), type "MakeNewformModFrm_14_8_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_14_8_a_b(:prec:=1)
    chi := MakeCharacter_14_a();
    f_vec := qexpCoeffs();
    Kf := Universe(f_vec);
    // SetVerbose("ModularForms", true);
    // SetVerbose("ModularSymbols", true);
    S := CuspidalSubspace(ModularForms(chi, 8));
    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_14_8_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_14_8_a_b( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_14_a();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,8,sign)));
    Vf := Kernel([<3,R![66, 1]>],Snew);
    return Vf;
end function;