
// Make newform 256.4.a.f in Magma, downloaded from the LMFDB on 14 June 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_256_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_256_4_a_f();".
// 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_256_4_a_f();".
// 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_256_a();"
function MakeCharacter_256_a()
    N := 256;
    order := 1;
    char_gens := [255, 5];
    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_256_a_Hecke(Kf)
    return MakeCharacter_256_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]];
    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 := [[0], [8], [-12], [-32], [-8], [20], [-98], [-88], [32], [-172], [256], [-92], [102], [-296], [320], [-76], [408], [-636], [552], [-416], [138], [64], [392], [-582], [238], [-1468], [992], [-584], [740], [-302], [1664], [-2328], [1734], [3032], [-1788], [-480], [-2300], [-1592], [-2208], [-3948], [-2104], [1412], [-3712], [1614], [-684], [4064], [2120], [-2816], [3848], [-652], [3050], [6336], [-4610], [792], [-5374], [-3488], [-4764], [-1344], [8596], [2874], [-2888], [-6540], [-10584], [-6368], [8758], [-716], [4408], [1202], [-5160], [-4876], [4834], [-4128], [4416], [4180], [-13736], [-512], [1732], [10436], [-12130], [5014], [-8024], [-2348], [1728], [62], [14112], [-4488], [-2482], [5894], [-7068], [-7616], [-13080], [13568], [-1696], [3096], [19208], [-16224], [-11292], [-5178], [-6856], [13732], [-10968], [-25612], [9768], [22838], [9208], [-10878], [18040], [26994], [18336], [-9286], [17536], [-5868], [-19286], [5240], [15520], [654], [8232], [24672], [22052], [-12024], [-19100], [9902], [19684], [19864], [-3256], [-2876], [7300], [2880], [-8800], [21076], [19336], [13664], [-19520], [20004], [31478], [7054], [9604], [3144], [22084], [14950], [-23432], [1044], [18208], [-12488], [-30172], [-32544], [-26156], [18646], [-5800], [25088], [-3004], [43282], [27880], [-7392], [29080], [-26688], [-19680], [-48466], [13610], [30692], [-4824], [-22986], [-17184], [-2920], [-27042], [44192], [-29824], [-11612], [22014], [-16572], [-20744], [-58476], [-40608], [-19958], [-44096], [20794], [25880], [45380], [3680], [21844], [-27392], [23048], [5364], [-16182], [27072], [-61244], [1860], [-16952], [-6106], [-58368], [-41662], [32600], [32744], [660], [-12216], [27498], [4754], [49892], [-12990], [-77728], [11028], [-12992], [7028], [7842], [-71496], [27156], [6656], [67528], [-37430], [-52360], [-60514], [18084], [-46944], [3576], [-59808], [-69786], [-67904], [30558], [53664], [31828], [-70988], [-15456], [-15218], [97600], [-60632], [11252], [45690], [39040], [63200], [-91048], [64516], [75880], [70912], [44874], [-29672], [48704], [-39970], [104228], [-88136], [43104], [101512], [87480], [-62112], [68276], [85490], [-96], [-41728], [-102936], [-20200], [22208], [38676], [52354], [79200], [75402], [47748], [-60728], [-62460], [3992], [-118812], [-67594], [-80896], [-83064], [48612], [-98828], [-68562], [-39128], [-98652], [-104358], [-104936], [-12572], [63956], [39272], [86650], [128640], [-136558], [-7776], [11448], [1636], [31690], [-136152], [-89088], [-36896], [1440], [70068], [39384], [52672], [119006], [-54636], [-2016], [-21534], [-17020], [10664], [44954], [-15624], [108708], [55204], [-156928], [-132444], [39480], [106472], [51334], [-132012], [33728], [123528], [48824], [-22546], [121592], [-97900], [159264], [-106492], [-38208], [-4060], [-164562], [64744], [-166304], [-122198], [-25080], [-156160], [-142142], [-60046], [-15864], [-137770], [-64908], [40064], [38054], [-141058], [77064], [49272], [183808], [-20364], [7668], [-96780], [-77824], [-78680], [-34920], [-135144], [-97004], [40290], [-200662], [68416], [-5692], [-139530], [-220988], [-76448], [60516], [-55992], [29076], [-149000], [64192], [150308], [84680], [95494], [148052], [-73920], [130884], [99834], [-11648], [25560], [80510], [176928], [58660], [6986], [27840], [-139176], [-231416], [57706], [40116], [84832], [146518], [-83928], [9208], [-38336], [-103020], [-164960], [-166156], [13640], [-257280], [4450], [-54850], [-176970], [-236620], [-119162], [41504], [-110046], [129416], [105696], [127296], [31732], [-76136], [119040], [214594], [178052], [203448], [-130936], [-23136], [168342], [-59776], [71782], [179864], [186068], [-118124], [52302], [58688], [-222762], [-25004], [138208], [-98572], [-217730], [-89208], [272872], [-5390], [-57228], [139288], [225672], [204266], [-4492], [52992], [290912], [-95394], [164768], [92836], [-56380], [-228032], [300760], [17222], [-171468], [-18584], [-160134], [-222536], [26912]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_256_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


// To make the newform (type ModFrm), type "MakeNewformModFrm_256_4_a_f();".
// 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_256_4_a_f(:prec:=1)
    chi := MakeCharacter_256_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(2999) - 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_256_4_a_f();".
// 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_256_4_a_f( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_256_a();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,4,sign)));
    Vf := Kernel([<3,R![-8, 1]>,<5,R![12, 1]>,<7,R![32, 1]>],Snew);
    return Vf;
end function;