
// Make newform 960.4.a.k in Magma, downloaded from the LMFDB on 25 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_960_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_960_4_a_k();".
// 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_960_4_a_k();".
// 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_960_a();"
function MakeCharacter_960_a()
    N := 960;
    order := 1;
    char_gens := [511, 901, 641, 577];
    v := [1, 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_960_a_Hecke(Kf)
    return MakeCharacter_960_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 := [[0], [-3], [5], [-20], [-56], [86], [-106], [4], [-136], [206], [152], [-282], [-246], [412], [-40], [126], [56], [2], [-388], [672], [1170], [-408], [668], [66], [-926], [198], [1532], [-444], [-62], [414], [996], [-264], [2278], [1812], [1534], [3016], [1814], [-1844], [-3768], [-938], [3968], [3514], [1480], [-2774], [3806], [856], [3020], [1684], [2004], [5042], [-3090], [-2136], [98], [-5040], [-1986], [-1416], [6670], [-48], [-6938], [-1694], [-6364], [3134], [-236], [-3776], [-7918], [-4362], [7980], [-8294], [-964], [-8670], [-2314], [1896], [-1484], [-12370], [5620], [-5880], [-2082], [1742], [-3270], [-6134], [-10392], [12690], [7408], [-5062], [7160], [17100], [8634], [2986], [2406], [14316], [-292], [-14056], [11204], [4608], [2468], [-12192], [-1714], [-18014], [-16748], [14018], [-412], [-18218], [23524], [23330], [-13124], [11714], [-17628], [-2802], [2664], [23962], [-11940], [-16794], [-20706], [10724], [5744], [27906], [20556], [10224], [12982], [-1512], [-16710], [7962], [-12226], [-8748], [-7324], [21934], [10690], [-13792], [24004], [-8562], [-13836], [22224], [-11544], [3814], [-25662], [30658], [30894], [21596], [8646], [24954], [40004], [-16570], [4388], [14364], [21170], [-10664], [3190], [20814], [-18988], [-11664], [8246], [22890], [33548], [32264], [51228], [2144], [-33584], [-3590], [-21686], [5174], [35524], [-16122], [-19188], [-38464], [-43930], [-17328], [18160], [9102], [14306], [5710], [-56056], [12170], [4360], [8450], [-31400], [-31806], [-60428], [-21170], [-2668], [62474], [-8756], [-8952], [25558], [-5018], [-23864], [-33698], [-36170], [-2308], [64490], [63080], [18314], [-36276], [380], [49318], [46108], [76454], [51666], [-7578], [-10066], [-29552], [77238], [63992], [29374], [-59694], [23160], [37318], [-20104], [-29084], [-40398], [-2356], [-47414], [-22322], [-31604], [40044], [-16672], [-74262], [-58780], [-37702], [-12432], [62438], [13330], [84680], [-92318], [-25284], [-67380], [-28990], [6662], [4768], [101524], [105272], [-54442], [21852], [12400], [-68030], [-26940], [94608], [-41070], [-52530], [-94040], [-38160], [66772], [52252], [41588], [-42094], [56534], [-96800], [-94828], [80240], [23548], [106528], [-94386], [21250], [48024], [56810], [-8762], [10968], [35242], [-2412], [-82866], [117906], [110324], [-83020], [-31518], [-54722], [82622], [-93836], [-56602], [28026], [-45708], [99614], [29530], [-9164], [-2814], [30040], [31322], [17932], [134956], [-111758], [133370], [120216], [50400], [-31656], [101776], [63698], [57004], [127064], [-29230], [-22674], [67704], [-4062], [38726], [121884], [69726], [125928], [159134], [31926], [-66680], [-139602], [-27848], [-25212], [-52638], [-90954], [69280], [-139300], [-102996], [-40526], [-121788], [83786], [5376], [88414], [-107584], [-65610], [-105366], [14260], [154096], [-11798], [-55616], [-75048], [104986], [89698], [169220], [125714], [92958], [-22612], [6334], [66018], [-163652], [123212], [-143576], [510], [-29318], [40454], [86984], [-37748], [-28172], [135348], [-26342], [172038], [-129062], [-201564], [42814], [-109458], [-73402], [-181832], [127958], [221920], [60002], [-119636], [3480], [-97842], [153676], [-111414], [-110242], [-196028], [112090], [-482], [179664], [-198056], [-82322], [-64088], [16174], [156442], [81304], [79320], [-24516], [48346], [-87546], [135916], [166922], [-167232], [-10876], [-5480], [71718], [-35976], [-210866], [-79680], [71776], [-201838], [67914], [45962], [-101738], [5350], [-68436], [81990], [-33300], [-12952], [101336], [216006], [-105828], [-134224], [187522], [19614], [-157848], [-37460], [-194552], [-131686], [117904], [-85734], [-165028], [-41106], [168042], [-138066], [-203500], [-31434], [-244074], [-131904], [-263986], [-9414], [129100], [-230472], [96802], [250158], [38852], [110116], [-106286], [226158], [-159520], [97196], [248190], [153720], [24150], [8846], [5672], [102736], [209010], [27030], [-153636], [131946], [-191332], [-45752]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_960_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


// To make the newform (type ModFrm), type "MakeNewformModFrm_960_4_a_k();".
// 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_960_4_a_k(:prec:=1)
    chi := MakeCharacter_960_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_960_4_a_k();".
// 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_960_4_a_k( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_960_a();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,4,sign)));
    Vf := Kernel([<7,R![20, 1]>,<11,R![56, 1]>],Snew);
    return Vf;
end function;