
// Make newform 1521.4.a.r in Magma, downloaded from the LMFDB on 17 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_r();".
// 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_r();".
// 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 := [-4, -1, 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 := [[0, 1], [0, 0], [-2, 1], [10, -11], [34, 12], [0, 0], [-18, 17], [26, 32], [-104, 12], [70, -96], [26, 34], [-102, -5], [-126, 22], [72, 143], [278, -121], [74, -30], [-246, 124], [-434, -190], [-150, 232], [50, -231], [-98, -260], [-524, 40], [1070, -182], [-166, -388], [718, -508], [-558, -142], [-152, -304], [-1216, 140], [-782, 43], [-1246, 236], [328, -344], [-604, -567], [-2326, 202], [-676, 847], [-2134, -244], [2106, -441], [1802, -924], [-2610, 430], [-1718, 902], [1642, 1114], [80, 29], [-358, 558], [-2916, 874], [-794, -206], [-2306, 1395], [680, 670], [1284, -537], [3214, -1737], [694, 1018], [-470, 101], [-1334, 621], [-1006, -993], [-1466, 2536], [3516, -1092], [-894, -151], [-3996, -480], [-834, -2190], [874, 1193], [2830, -3402], [2290, 484], [-5316, 1056], [-1682, -1487], [2830, 1748], [-584, -2858], [1434, 3345], [378, -2580], [4322, 122], [-1686, 1841], [-7332, 1753], [-6518, -409], [-6566, 1650], [-2670, 54], [4380, 858], [-2246, 2010], [2834, -378], [3866, 2721], [-1730, 4392], [8326, 696], [3442, -530], [-4746, 5646], [3988, -6781], [-6262, 6139], [8482, -1417], [-13098, 1909], [2616, -2190], [-10028, 111], [-3830, 1276], [-8082, -2216], [7734, 1669], [-9194, -1734], [-11404, 1008], [-9146, -2137], [7058, 1070], [-5376, 6329], [-2134, -178], [5812, 4394], [-726, -112], [11030, 1303], [-4488, -4580], [-2830, -5763], [7172, 1301], [1534, -797], [-1744, 2025], [-3262, -8085], [1648, -3555], [4966, -4152], [-2994, 1566], [-21510, 4416], [1772, -5578], [5634, 1559], [9572, -1134], [-194, 9088], [14858, -4610], [-11226, 1534], [20666, 1235], [-23590, 3976], [12290, -6848], [7288, -11282], [1918, 7682], [-1360, 7028], [-2298, -9636], [-4238, 6007], [-7254, 1538], [19418, 1986], [6778, 5846], [-7954, -10266], [10558, -1964], [-8656, 4854], [22248, -1238], [-4926, -1199], [-5986, -8978], [-17278, -729], [-21152, -1704], [-6426, -4218], [15762, -7068], [5886, 474], [12606, -131], [-18738, 5978], [-14234, 10156], [7958, -10489], [-16702, 1522], [8574, 4201], [-5584, 1350], [-5606, -540], [46, -12764], [-20550, 6446], [-6590, -595], [-17242, 15436], [-8920, 16580], [-10542, -1397], [25946, -9371], [-6190, 2725], [-2552, -13307], [30020, 4164], [-22952, 4875], [33844, 690], [13260, -4756], [49458, 2450], [-12246, 3532], [-24546, 13027], [20698, 2964], [2514, 3859], [33094, -2509], [-30852, -7249], [19058, -7136], [-21818, 10859], [-13828, 8902], [-19586, 18784], [21230, 1680], [-14022, -17852], [28678, -13632], [-41638, -6607], [4024, -14792], [20486, 13418], [-5216, 4466], [-24254, -10693], [-28354, 10350], [24018, -3692], [33392, 6348], [-21614, 18394], [-11922, 24218], [-300, 16456], [17294, -880], [9778, -25848], [27450, -6022], [13534, -3364], [-28534, -2884], [-16934, -7548], [34390, 7814], [-63430, 1390], [-23958, 25679], [-18258, 24318], [-16260, -12595], [-67182, -4316], [29648, -3909], [41622, 680], [5766, 4130], [43490, -7258], [6194, 17834], [23380, 18554], [-1810, -27119], [10108, -27960], [29642, 3045], [-5286, 10731], [21562, -11194], [490, -23540], [-48042, -45], [-22800, 3876], [-58366, -1096], [-18800, 4025], [45826, -25100], [2914, -17828], [19924, -9576], [-11638, 7656], [-30822, -14707], [31862, 2586], [27548, -19528], [11218, 15453], [-1818, 37406], [20182, -31293], [68054, 2470], [43018, -1713], [47498, -17150], [52062, -19981], [-55336, -10521], [10502, -25692], [-6190, -4413], [-35900, -20556], [64256, 6350], [48714, 7012], [-79414, 15932], [-59984, 4860], [1938, 36518], [48414, -12645], [24260, 5993], [-56054, -20401], [-62962, 14068], [-13290, -18641], [11180, 1860], [-15004, 49810], [-1258, 32364], [73824, -18257], [-44772, 35596], [-43870, 4949], [-12758, 12724], [-26472, -9368], [54834, -8013], [31866, -9358], [26562, 24036], [-9704, -10530], [9938, 8531], [-78286, 9816], [2230, -30251], [12242, -17791], [72430, -28634], [17034, 22672], [36154, -9382], [40426, -2336], [39542, -44078], [-19666, -21568], [58576, 28180], [3940, 31369], [15618, -25631], [-36938, -37240], [38874, 9750], [72292, 3367], [69678, -39075], [41378, -41350], [-22562, -17530], [-13038, 11898], [-48490, 26234], [-102342, -9832], [-24018, 33296], [-55792, 37980], [-56514, 31664], [-84414, 13054], [94278, 10736], [16294, -40640], [123974, 1918], [-36912, 33791], [-51660, -40324], [61578, -22557], [-41588, -21318], [80510, -10100], [10654, -1298], [20992, 41958], [-53970, 25631], [-23818, 53951], [-34470, 33157], [-83250, -7496], [100370, -1958], [-42108, -7011], [71026, 2446], [-6502, 54276], [45274, -54348], [16830, -7186], [-66204, 49458], [-41418, 37374], [36860, -38589], [46686, 244], [-48978, -43643], [-37450, -46899], [128904, 8170], [-79504, -7268], [32844, -25843], [-26370, -59982], [134504, -23063], [115854, -10954], [13146, 40042], [40650, 20642], [68812, 4164], [45782, -54633], [51118, -57], [75320, -26236], [134478, -30351], [122850, 12703], [-34402, -56616], [-20506, -44374], [92166, -11830], [3018, -22656], [-57436, 30000], [-19098, -51214], [128450, -19274], [-56782, 40035], [-18054, 32414], [-93278, 23136], [-133194, 11332], [21426, -3578], [65952, 34000], [-83466, 7110], [33286, -17352], [143310, -24102], [79244, -8456], [84122, -72602], [-48246, 80606], [24062, -29806], [115998, 16776], [-72478, -35526], [34126, -14214], [85464, -23370], [127226, -4929], [48790, 20059], [152574, -16863], [-22472, 48514], [-52130, 85911], [152696, -18139], [-38510, 11220], [-37298, -61510], [70038, -37807], [-11970, -12044], [-48086, -12600], [-60706, 16132], [40006, 57913], [22456, 43000], [-71578, 38468], [31610, 3580], [-18434, 7493], [-56378, 42120], [-91814, 90191], [52666, 46629], [99614, 13620], [112482, -19736], [33564, -33910], [-66442, -8850], [55564, 16424], [-78542, -13029], [124998, -36129], [88638, -19510], [-28822, 86643], [-3524, -1421], [36104, -39151], [-66466, 50835], [-50142, -10216], [136492, -12676], [-20810, 55024], [-172594, 12148], [-93776, -7514], [-17826, 92638], [115342, -22168], [82062, 28520], [74290, -77192], [-32214, 57036], [100718, -12114], [-14470, 86856], [-70410, 83386], [94878, -24657], [-107074, 49485], [-47602, -77238], [-59038, -24468], [-61812, 48208], [21926, -96336], [66790, -3777], [-53246, -15116], [123516, -84297], [2366, -31637], [210566, -33281], [46814, -57537], [-153626, 63729], [-7448, 90660], [36318, -71707], [60986, -92093], [-131830, -45196], [192338, 17755], [152706, 34478], [-192842, -15257], [-83548, -58778], [-126926, 87917], [52570, -77386], [26614, -70180], [-39710, 44518], [29742, -9169], [-234062, 14408], [-194684, 9951], [73432, 25183], [33654, -51953], [75794, 7994], [-109118, -29547], [69724, -85508], [156130, 6600], [-57888, 55994], [-41238, 83546], [-59178, 72596], [122262, -67566], [5436, 2283], [20838, 100318], [156986, 29292], [22984, -120731], [82106, 32150], [84646, -14806], [-8276, 30246]];
    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_r();".
// 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_r(: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_r();".
// 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_r( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_1521_a();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,4,sign)));
    Vf := Kernel([<2,R![-4, -1, 1]>,<5,R![-2, 3, 1]>,<7,R![-494, -9, 1]>],Snew);
    return Vf;
end function;