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


// To make the character of type GrpDrchElt, type "MakeCharacter_1_a();"
function MakeCharacter_1_a()
    N := 1;
    order := 1;
    char_gens := [];
    v := [];
    // 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_1_a_Hecke(Kf)
    return MakeCharacter_1_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 := 16;
    raw_aps := [[216], [-3348], [52110], [2822456], [20586852], [-190073338], [1646527986], [1563257180], [9451116072], [-36902568330], [71588483552], [-1033652081554], [1641974018202], [-492403109308], [-3410684952624], [6797151655902], [9858856815540], [4931842626902], [-28837826625364], [125050114914552], [-82171455513478], [-25413078694480], [-281736730890468], [715618564776810], [612786136081826], [-817641571654098], [741114547982552], [-2514301452571644], [1268353947457190], [-2054162866352238], [2990675947730816], [-1626226733523348], [10592201511845946], [-18670911522208540], [-12555957651134850], [28758788173002152], [-14527638158544394], [16774137626235212], [64199924334659736], [-75986044070753178], [93374877047641020], [74317669765796702], [-98622390566317248], [-8911776556935358], [35417574134917326], [-286460988828497800], [375833826551569052], [-25307787891567328], [303691994772520716], [107991586981028270], [-790506217682068518], [352956492128946960], [68568967462179602], [1588056499493163252], [-828562059034685694], [1404454851260256312], [1435820360137377030], [509159783589159152], [568946210839495166], [-4061836641519738198], [2785059360043480532], [-3638027115775808658], [-975296087300763844], [336691659808831752], [3655511591562842762], [-797379754904038314], [-10158512678893878748], [-4814653881814264654], [4502751666049603476], [22432314991369546550], [8025097307604388002], [1615072843193774040], [-9977994448626462064], [-23686588714485006178], [18685121156982396020], [-30252066823473936768], [-1007143932951970290], [35632419864211771526], [39432703089652313202], [-52782330312986685910], [8626299410782459980], [-42924857411864654098], [50476373204365891152], [50573397093112460882], [24794607909003802760], [-130654278967707957708], [-77827995175927650750], [-118450824223090809094], [138643031119690399302], [175645004594544548912], [-136630964339996286564], [64105810232391205920], [-24134317584467009704], [-28090793884898388348], [-171993603615123328900], [-183496515279116728248], [267203608600830147990], [-201467622292352438598], [358988638501746249572], [202328193434804061302], [74096261798900732876], [2096258165482303206], [687352768991581988412], [-90521802682779479670], [205773783706772924852], [570778361125879587266], [-930363435275990480604], [354224807936897233842], [-330045214241748274200], [-335883995287326350998], [-1334375729047804699744], [5688435372944444462], [398915477694945040986], [-540016641807600756220], [959110676308468369352], [-925925342118025231998], [-765927858754026197908], [1360754271840232136376], [-271808672598939928698], [674316250030836216540], [1267268866679850853502], [-1139445987246513182878], [1744312019604807273966], [-1437390792807621086268], [-1775475477541366791748], [1431003389167119462102], [-2418403187913507119410], [4744438233844295461680], [-3592646945141378296984], [-2768235432318591148618], [3558241342219996245260], [-1940919448954530628008], [4755648944590489774352], [3621365428229853791006], [-3863611493551703142198], [-5393270765226373961470], [6570369572295558011022], [-372073581785586295204], [2615111134656369060726], [5348992550333367700890], [8464923563697895536452], [7993973312037568404702], [1968405114515095582472], [-14353875238140627185484], [-8833271554558074901930], [12669646577325254931960], [6005322570637837971902], [14758947243889334451306], [964956153039498126740], [3443909030568133847712], [-10985023771436785030234], [-7984619570820038841198], [5452357163027018171132], [16712728456718529373896], [-15338437699335420028444], [-14913372374629398837648], [-5866674181682656972120], [-16794638211772494268830], [5044656161649011646746], [16542548428439130835302], [-9815830250691060979524], [-5979138890537025630198], [-14475696971583064201864], [1779211406030716089252], [10408769654278916425266], [32646136192413800125032], [-7473270757865514401248], [-31080942682464818856274]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_1_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


// To make the newform (type ModFrm), type "MakeNewformModFrm_1_16_a_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_1_16_a_a(:prec:=1)
    chi := MakeCharacter_1_a();
    f_vec := qexpCoeffs();
    Kf := Universe(f_vec);
    // SetVerbose("ModularForms", true);
    // SetVerbose("ModularSymbols", true);
    S := CuspidalSubspace(ModularForms(chi, 16));
    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_1_16_a_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_1_16_a_a( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_1_a();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,16,sign)));
    Vf := Kernel([],Snew);
    return Vf;
end function;