
// Make newform 1521.4.a.s in Magma, downloaded from the LMFDB on 25 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_1521_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_1521_4_a_s();".
// 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_1521_4_a_s();".
// 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 := [-14, 0, 1];
        Kf := NumberField(Polynomial([elt : elt in poly]));
        AssignNames(~Kf, ["nu"]);
    end if;
    Rfbasis := [Kf.1^i : i in [0..Degree(Kf)-1]];
    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_1521_a();"
function MakeCharacter_1521_a()
    N := 1521;
    order := 1;
    char_gens := [677, 847];
    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_1521_a_Hecke(Kf)
    return MakeCharacter_1521_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 := [[1, 1], [0, 0], [12, -2], [0, 2], [-22, -12], [0, 0], [-82, -4], [-24, -2], [-4, -48], [-202, 24], [-20, 26], [50, -28], [100, 94], [-308, 52], [-162, 32], [82, 120], [70, 40], [314, 136], [236, 170], [214, -84], [450, -76], [-216, -88], [-694, 64], [480, -190], [266, 220], [182, -28], [-280, 436], [-984, -40], [30, -280], [-702, -232], [24, -340], [736, 144], [164, 190], [460, -40], [-1080, 290], [420, -114], [-322, -528], [1976, -238], [310, -112], [2862, -260], [-1572, -496], [1286, -288], [-1736, 208], [-954, 392], [-4252, 102], [2680, -100], [-3852, 48], [-4436, -174], [-306, -620], [686, -632], [-3118, 660], [1942, 876], [-2766, -892], [-5256, -312], [-910, 572], [-2628, 840], [3606, 408], [-856, 1378], [-1274, -144], [-2988, -290], [4236, 564], [-5016, -902], [-316, -986], [-2484, -896], [-1290, 816], [-4692, -378], [-6540, -254], [5182, 1144], [-5324, -1640], [-382, 1084], [5252, -498], [-2490, -1920], [264, -2400], [178, 1200], [-4240, -1770], [7450, 696], [-7990, 420], [9730, 612], [-2544, 1642], [-1254, 348], [1312, -544], [8350, 1736], [-6470, 2024], [5990, 776], [-4032, -1836], [-7264, 2952], [6984, -1838], [5114, -1036], [2692, 1630], [5680, -486], [640, 2232], [-2630, 1820], [12300, 1642], [-5968, -568], [11892, -254], [340, 2384], [-5020, -1330], [-11398, -8], [-5020, 56], [2078, -1200], [-9036, 1240], [-6480, -2942], [-5780, 3384], [-590, -3780], [9652, -1380], [-5522, 4296], [-1278, -96], [16920, -774], [2388, -736], [1126, -3008], [-10696, 3480], [12762, -3580], [10316, -3554], [-19252, 1094], [-10404, 1294], [6974, 4588], [3628, -502], [5904, -1976], [886, -3088], [11748, 584], [-30, 636], [-12246, -5392], [1850, -208], [4858, -1380], [-13084, 2914], [9466, 1524], [22518, -3616], [-2252, -1656], [-14536, 1580], [-8794, -5968], [916, 4478], [-11498, -2628], [17056, -1428], [3594, -5400], [-7836, -5754], [6150, 6000], [23512, -2678], [-21020, -1298], [8850, -8228], [9114, -3928], [-20072, 1910], [-26624, 4402], [14912, -3108], [-31294, 3636], [15522, 6056], [19754, 6452], [14642, 1420], [-10510, -4148], [-7692, -10136], [1162, -848], [1042, 9548], [-8566, -7172], [-13260, 356], [1240, 8472], [-3932, 1452], [-26104, 1704], [10144, 2632], [-3772, 494], [-21418, 2984], [9928, 298], [22766, 2496], [11434, -524], [-20140, 6802], [-12268, -3208], [-24496, -6710], [-5222, 14060], [-32864, 6392], [1854, -4480], [-8974, 5148], [-30770, -1436], [28670, -780], [-7758, 13012], [16428, -4040], [35274, -5924], [-42376, -548], [-2, 4304], [39404, 198], [36232, 6490], [10952, 5172], [26298, 7904], [22984, 4234], [13948, 7360], [42810, 3256], [2552, 3474], [-58042, -3520], [-13910, 3860], [41454, -3440], [256, -13098], [-24622, -5396], [-25114, 2752], [27538, -2480], [21474, -792], [-22028, 8560], [28408, -5342], [-23852, -8064], [-7894, 18932], [-27890, -5144], [-26262, -4400], [44640, 98], [-33008, 11120], [-4492, 1774], [3376, -13608], [-1066, -18816], [-6222, -7296], [28102, 7616], [19290, 10192], [26196, -12702], [-47160, 1104], [38952, 8066], [23972, -8780], [6726, -20824], [33474, -8672], [1432, 4032], [-8630, 1068], [17426, 11648], [-52306, -8592], [-13656, -1520], [-27610, 2616], [30398, -3748], [16952, 7242], [51886, 5792], [43432, -11214], [40616, -1142], [-84548, 490], [11584, 5664], [-4834, -8712], [40990, 1344], [57332, 5976], [60840, 3748], [-47786, 376], [55278, -7280], [-3476, 1872], [60712, 3862], [-26874, 7884], [1836, -3776], [11446, -25444], [-39902, 1964], [-50792, -6350], [-94868, -864], [464, -23912], [13138, 24960], [-26932, 1628], [-20032, -24520], [-24726, 2692], [-18912, -15110], [80620, 9760], [41064, -6486], [-6898, -4120], [-19920, -16530], [-53932, -13272], [-2370, 4912], [2980, 13134], [-31158, 5164], [-36234, 8488], [91758, -7472], [-10790, 1420], [22898, 8368], [-54960, -418], [-76918, 13588], [-4914, -20000], [39288, -21820], [-24060, -12008], [115286, -472], [-15610, 18256], [5784, -7950], [-108660, -1120], [-74628, 6246], [-4140, -4546], [54764, 1966], [36966, -5628], [93426, 6352], [47356, -9230], [-5030, -27836], [-45328, 4836], [15790, 15088], [5624, -33634], [29678, -4228], [86730, 16340], [-35666, 3848], [-10864, -1336], [21896, 11312], [129936, 1434], [-122320, -7488], [27082, 23132], [33492, 6998], [-64, 4344], [-14426, -34472], [42888, -15358], [-69604, 14318], [-77870, 9220], [-20830, 7936], [-52356, 8184], [-4608, 7882], [70366, 16872], [526, 12144], [11258, -7876], [-66384, 14592], [60522, 8820], [69724, 2952], [13316, 16310], [114122, -10288], [49324, -7650], [-48112, -11716], [55292, 20320], [-21236, -3344], [34122, -6732], [-7948, 40648], [-1322, 13744], [57850, -4100], [-68966, 19768], [57188, -28392], [-1772, 26214], [16922, 21912], [-23772, -24440], [-12162, -1848], [-5722, 704], [-75770, 28260], [-55470, 8984], [82786, 5572], [62742, -29124], [75788, 5028], [56038, -18584], [-19486, -19352], [-13972, 7998], [-63256, 33218], [57094, 24456], [-45160, -22510], [-108124, 23402], [-67664, 21160], [151110, 5424], [167878, -2880], [36366, 43440], [-9392, -21376], [40874, -7652], [-176816, 3238], [-133854, -2584], [27882, -33084], [183388, -4770], [-177686, 5268], [8400, -27372], [87458, -6852], [50490, 22192], [-41736, 29802], [-65856, 5036], [66188, 16830], [88380, -33640], [27298, 31488], [93152, -23222], [67454, -39664], [-76502, -884], [21868, 43050], [-105758, 19508], [-57504, -12278], [-73944, -30796], [43722, 19216], [26834, 17416], [-168218, 4016], [-198610, -5796], [-92726, -6076], [154034, 1608], [-59014, -15060], [15454, 53908], [-19756, -46816], [4222, -804], [159524, 10420], [78034, -18192], [3708, -29178], [-16076, -8474], [-45982, 6048], [132720, -13304], [9324, 26212], [47158, -6840], [-79510, -46600], [-180592, 9928], [17622, 12704], [38786, -28232], [25416, -2936], [178238, 19208], [-33962, -61576], [-108110, -23432], [-26976, -7454], [-31332, 34470], [17744, -46806], [-45176, 48930], [-67624, 1538], [129090, 34236], [-35680, -40314], [-59878, -6216], [-68920, 7254], [-173084, 19360], [-34942, 32844], [25160, 27186], [69330, 5860], [-89316, -7560], [9902, -55760], [-103010, 15896], [8192, 68646], [28910, -13740], [142420, -29736], [-161140, -4186], [-49142, -31396], [90054, 23996], [140080, 13682], [68624, -6790], [193932, 10726], [47560, -57256], [-80794, -12560], [95824, 48026], [-70948, -10250], [93234, 21532], [-111182, 13840], [48154, -39376], [-47500, 1224], [-53388, 62060], [-54718, -58936], [94578, 13760], [-116050, 13500], [114648, -41608], [37016, -33414], [-49420, 48584], [34118, 22484], [88330, -31988], [-24690, 44388], [-69600, 20880], [-64450, -24616], [-81200, 54858], [85944, -26024], [160560, 7850], [284384, -5942], [-60712, 5784]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_1521_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


// To make the newform (type ModFrm), type "MakeNewformModFrm_1521_4_a_s();".
// 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_1521_4_a_s(:prec:=2)
    chi := MakeCharacter_1521_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_1521_4_a_s();".
// 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_1521_4_a_s( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_1521_a();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,4,sign)));
    Vf := Kernel([<2,R![-13, -2, 1]>,<5,R![88, -24, 1]>,<7,R![-56, 0, 1]>],Snew);
    return Vf;
end function;