
// Make newform 16.4.e.a in Magma, downloaded from the LMFDB on 19 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_16_e();"
// To make the character of type GrpDrchElt with Codomain the HeckeField, type "MakeCharacter_16_e_Hecke();"
// To make the coeffs of the qexp of the newform in the Hecke field type "qexpCoeffs();"
// To make the newform (type ModFrm), type "MakeNewformModFrm_16_4_e_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_16_4_e_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
        poly := [1024, -512, -64, 96, 56, -80, 14, 6, -1, -2, 1];
        Kf := NumberField(Polynomial([elt : elt in poly]));
        AssignNames(~Kf, ["nu"]);
    end if;
    Rf_num := [[1, 0, 0, 0, 0, 0, 0, 0, 0, 0], [-1, 4, 0, 0, 0, 0, 0, 0, 0, 0], [3072, -2752, 384, 536, -120, -170, 82, 7, -14, 1], [-7424, 3264, -3008, 8, -120, 210, -114, 101, -2, 3], [7424, -1984, -1472, 872, 120, -230, 14, 49, -18, 7], [-4864, -1856, 2752, -792, 520, 90, 86, -39, 38, -17], [2304, -2624, -512, 1432, -680, -90, 74, 39, -38, 17], [-1280, 1088, -1216, 136, 184, 18, -58, 29, 6, -5], [5376, -3776, 3712, -4472, 1800, -350, 46, -179, 158, -37], [-7168, 17088, -11136, 616, 2040, 10, -818, 297, 46, -49]];
    Rf_basisdens := [1, 1, 1280, 1280, 1280, 1280, 1280, 256, 1280, 1280];
    Rf_basisnums := ChangeUniverse([[z : z in elt] : elt in Rf_num], Kf);
    Rfbasis := [Rf_basisnums[i]/Rf_basisdens[i] : i in [1..Degree(Kf)]];
    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_16_e();"
function MakeCharacter_16_e()
    N := 16;
    order := 4;
    char_gens := [15, 5];
    v := [4, 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;

// To make the character of type GrpDrchElt with Codomain the HeckeField, type "MakeCharacter_16_e_Hecke();"
function MakeCharacter_16_e_Hecke(Kf)
    N := 16;
    order := 4;
    char_gens := [15, 5];
    char_values := [[1, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0, 0, 0]];
    assert UnitGenerators(DirichletGroup(N)) eq char_gens;
    values := ConvertToHeckeField(char_values : pass_field := true, Kf := Kf); // the value of chi on the gens as elements in the Hecke field
    F := Universe(values);// the Hecke field
    chi := DirichletCharacterFromValuesOnUnitGenerators(DirichletGroup(N,F),values);
    return chi;
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, 0, 0, 0, -1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, -1, 0, 0, 0, 0], [-1, 0, 0, 0, 1, 0, -1, -1, -1, 1], [-1, -1, 3, -1, 4, -1, -2, 2, 0, -1], [2, 2, 0, -1, -2, 0, -6, 2, 2, 0], [2, -1, 1, 0, -9, 4, 7, -1, 1, 4], [-3, 3, -5, -2, -6, 2, -16, 0, -2, -5], [0, 8, 10, 0, 10, 3, 18, 2, -2, 2], [-7, -7, -27, 3, 12, 3, -6, -10, 0, 1], [-16, -3, 13, 0, -11, -24, -11, 5, -5, 0], [34, 2, 6, -4, 4, 4, 8, 0, -12, 6], [-5, 4, -10, -28, -9, 0, -7, 1, 1, -1], [2, 2, 6, 24, -16, 24, -12, 4, 0, 6], [-72, -8, -80, 3, -8, 0, 40, 8, 8, -16], [-102, 10, -2, -4, 36, 4, -24, 0, 20, -2], [-5, 12, -38, 68, -85, 0, 37, -19, -19, -21], [164, -36, 192, -29, 84, 0, 60, 28, 28, -8], [94, 19, -79, 0, 67, -12, 51, -21, 21, -4], [176, 0, -202, 0, 30, -13, -74, 22, -22, -26], [1, 1, 333, -21, -20, -21, 10, -26, 0, 9], [-29, -29, 37, 14, 64, 14, -62, 10, 0, -3], [-452, 12, 36, -8, 40, 8, 48, 0, -56, 36], [-240, -24, 256, 0, -56, 55, 104, -72, 72, 40], [53, 53, 35, -14, 0, -14, 110, 102, 0, -53], [25, -49, -73, -10, -158, 10, -48, 0, 86, -73], [-173, -20, -126, 44, 99, 0, -19, 5, 5, 27], [-75, -75, 689, 51, 60, 51, -190, -50, 0, 45], [-474, 22, -416, -57, 106, 0, -194, -10, -10, 80], [40, 13, 53, 0, -99, -32, 221, -35, 35, 80], [-92, 98, -30, 84, -4, -84, -256, 0, -140, -30], [886, -74, 2, 180, -116, -180, 152, 0, 28, 2], [608, 16, -594, 0, 6, -61, -2, 30, -30, -2], [72, 72, -208, -20, -272, -20, 136, -120, 0, 64], [-724, 28, -864, 11, -220, 0, 108, 60, 60, -112], [-385, -24, -404, -144, -83, 0, 179, -45, -45, -43], [-3, -3, -935, 155, -228, 155, -302, 62, 0, 117], [-558, 103, 601, 0, 311, 204, 71, 15, -15, -60], [912, -128, -1136, 0, 80, 47, -304, -144, 144, -96], [183, 183, 1243, 223, -156, 223, 526, 50, 0, -105], [-382, -17, 193, 0, 199, -332, -489, 111, -111, -172], [-736, -56, 744, 0, -272, 87, -16, 32, -32, 64], [751, 168, 548, -272, -299, 0, -373, 107, 107, -35], [1306, 154, 78, 364, 212, -364, -152, 0, -252, 78], [-539, -37, 3, 110, -310, -110, 80, 0, -242, 3], [1087, -192, 1144, -48, 297, 0, 471, 183, 183, -135], [43, 43, 927, -153, -396, -153, -122, -102, 0, 155], [-816, 200, 780, 0, 708, 31, -236, 164, -164, -236], [810, -182, 62, -468, 84, 468, 488, 0, 324, 62], [-64, 16, 136, 0, -328, -257, -104, 248, -248, 56], [-1401, 72, -1364, 24, -19, 0, -269, -93, -93, 109], [139, 139, -227, 302, -256, 302, 466, -166, 0, -11], [-1420, -28, -212, 232, -648, -232, -368, 0, -168, -212], [161, -217, -97, -282, -590, 282, 240, 0, 38, -97], [-168, 336, -688, 443, -1040, 0, -304, 0, 0, -184], [1004, 58, -38, 68, 460, -68, -192, 0, 420, -38], [-211, -211, 681, -301, 732, -301, -238, 126, 0, -155], [424, -259, -587, 0, -467, 272, -83, -243, 243, 96], [118, -58, 338, -92, 188, 92, 792, 0, -372, 338], [-1485, -188, -1162, 156, 663, 0, 89, 17, 17, 135], [-325, -325, 589, -618, 1024, -618, -222, -54, 0, -187], [322, -102, 432, -457, 390, 0, 18, 170, 170, 8], [-2385, -64, -2360, -224, 217, 0, 39, 167, 167, -39], [880, 8, -570, 0, -490, 539, 718, -98, 98, 302], [187, 187, -49, -33, -1100, -33, -250, -102, 0, 363], [114, 114, 422, -592, 144, -592, -252, 852, 0, -186], [3172, -271, -3611, 0, 25, -192, -647, -231, 231, -168], [914, -78, 624, 223, -578, 0, 890, 2, 2, -368], [-414, 368, 272, -352, 1120, 352, -192, 0, -160, 272], [-2660, -252, -2480, 259, 604, 0, 404, 244, 244, -72], [-4544, -275, 3965, 0, 85, -120, -1131, -27, 27, -304], [-136, 6, 550, 252, 628, -252, 1088, 0, -484, 550], [253, 253, -551, -157, -164, -157, 242, 606, 0, -171], [3434, 666, -306, 156, 516, -156, -1944, 0, -204, -306], [1835, 92, 1706, 164, -981, 0, 613, -723, -723, -37], [-1964, -212, -1328, 323, 1172, 0, -324, -100, -100, 424], [-1324, 340, 412, -360, 1512, 360, 144, 0, 8, 412], [3995, 20, 4158, 580, -161, 0, 81, -487, -487, 183], [316, -735, -1491, 0, -383, -712, -2143, -207, 207, -440], [441, -481, 279, -234, -702, 234, 1520, 0, -298, 279], [114, 114, 614, 320, -1904, 320, -412, -1036, 0, 838], [-4192, -96, 4186, 0, -110, 411, 250, -262, 262, 90], [-4057, -152, -3964, 616, 789, 0, -181, 603, 603, -59], [1468, -916, 228, -104, -696, 104, 2288, 0, 680, 228], [-871, 287, -265, 150, 450, -150, -1104, 0, 406, -265], [255, 255, 3459, 1143, -60, 1143, 2014, -1054, 0, -225], [-4062, -86, -3840, -1081, 486, 0, -142, 42, 42, 136], [-2871, 223, -953, 278, -1598, -278, -2352, 0, -138, -953], [658, 658, 2486, -536, -1296, -536, 884, 452, 0, -10], [1028, 745, 141, 0, -71, -1784, 1625, 713, -713, 424], [5712, -352, -640, 272, -2160, -272, -576, 0, -176, -640], [6784, 264, -5738, 0, -346, 1483, 2782, -690, 690, 782], [122, -102, 302, 268, 884, -268, 808, 0, 484, 302], [-217, -217, -6101, -417, 2532, -417, 1518, 146, 0, -1049], [5972, 596, 5104, -1109, -2020, 0, -364, -284, -284, -272], [576, -280, -500, 0, -1692, -1809, -268, 420, -420, 356], [-671, -671, 1709, 803, 1836, 803, -854, 6, 0, -247], [5950, 403, -6047, 0, 659, 52, -1341, 1147, -1147, -500], [-480, -480, -3592, 636, 1072, 636, -488, -360, 0, -56], [-474, 134, 32, 439, 282, 0, -818, -730, -730, 640], [-7666, 595, 8193, 0, 1267, 1988, 995, 59, -59, -68], [-624, 240, 1910, 0, -1906, 707, 2278, 294, -294, 1046], [-450, -165, 937, 0, -1469, 556, 1139, -165, 165, 652], [1376, -1128, -3470, 0, -606, 227, -4470, 282, -282, -966], [-96, -96, -2936, 420, -1968, 420, -1112, -1240, 0, 1080], [826, -318, 512, -2089, 846, 0, 426, 1474, 1474, -632], [-2922, 756, 1620, -440, 4184, 440, 1728, 0, -568, 1620], [4838, -538, 5872, -265, 2058, 0, 94, -10, -10, 496], [-3226, -364, -140, -184, -552, 184, 448, 0, 456, -140], [503, 503, -3973, -1385, -668, -1385, 814, 530, 0, -169], [-303, -303, -2457, -1446, 1024, -1446, -1226, 1038, 0, -209], [-1520, 752, -464, 1024, 704, -1024, -2432, 0, 128, -464], [-7273, -744, -6548, -904, 1597, 0, 1379, 147, 147, -19], [125, 125, 2907, -6, -320, -6, -978, 1158, 0, 35], [-4054, -1286, -1952, -1793, 3334, 0, 1810, -870, -870, 816], [-91, -91, 737, 1723, -452, 1723, 130, -946, 0, 317], [3021, -189, -245, 2174, -2054, -2174, -112, 0, -1186, -245], [3920, 192, -3474, 0, -842, -2013, 174, 718, -718, 254], [441, 441, 8709, 795, -1268, 795, 1050, -554, 0, 193], [960, -331, -2299, 0, 2245, 640, -1787, -891, 891, -1008], [4912, 40, -4296, 0, -416, 103, 1888, -656, 656, 576], [-6081, 1160, -7740, 864, -4267, 0, -373, -949, -949, -499], [653, -605, -1109, 958, -838, -958, -1008, 0, 2590, -1109], [-8517, -124, -9098, -1756, -265, 0, 761, 897, 897, -705], [1512, 1320, 288, 3203, -2696, 0, -2584, -248, -248, 96], [-1904, -232, 1836, 0, 884, 3823, 1540, -1676, 1676, 164], [8330, -593, -7607, 0, -3681, -404, 1583, -137, 137, 1316], [7459, -212, 6922, -916, -517, 0, 1365, 557, 557, -749], [-14916, -628, -700, 376, -3160, -376, -144, 0, -504, -700], [-1413, -1413, -3473, -1233, 2164, -1233, -1530, -1958, 0, 331], [-7492, -487, 8037, 0, -1263, 896, 2865, -1775, 1775, 1032], [6704, 1120, -5620, 0, 1916, -1257, 1772, 396, -396, -36], [-133, -133, 11567, -361, 1332, -361, -474, 1274, 0, -533], [-3284, -1508, -844, -264, -920, 264, 1328, 0, 3784, -844], [3807, 120, 5124, 3056, 1669, 0, -2149, -965, -965, 1437], [-890, -890, 2962, -1464, -48, -1464, -2084, -1524, 0, 914], [6565, 1627, 35, 494, -182, -494, -3184, 0, -3506, 35], [8835, -724, 10234, -1428, 2251, 0, 645, -547, -547, 675], [-11712, -968, 9336, 0, 752, -3129, -4880, 128, -128, -1408], [-7802, -341, 7521, 0, -677, -940, -437, -125, 125, 60], [742, 742, -3294, 112, 2352, 112, 5228, 92, 0, -1918], [-6106, 630, -6736, 2119, -1318, 0, -1202, -58, -58, 0], [4571, 12, 3746, -1964, -1277, 0, 1229, 373, 373, -813], [-2037, -2037, 2911, -449, 2548, -449, -6042, 442, 0, 763], [-3686, 2402, -4688, -457, -3810, 0, -5798, -1806, -1806, 1400], [7646, 147, -8895, 0, 2995, -1516, -2589, 91, -91, -1396], [203, 203, 18783, 423, 1204, 423, 2182, -166, 0, -805], [-2229, 716, -2662, -1964, 747, 0, -3611, 1613, 1613, 283], [172, 172, -5988, 732, -2224, 732, -3536, 2000, 0, 940], [11908, -852, 10896, 1843, -652, 0, 4060, 1372, 1372, -1864], [13482, 1226, -546, 684, -1836, -684, -3544, 0, -3196, -546], [-868, 1361, 1613, 0, 3889, 1352, 1425, 65, -65, -616], [7926, 1220, -700, -24, 824, 24, -3840, 0, -216, -700], [-10400, -552, 10414, 0, -146, -629, 2118, -2090, 2090, 566], [2547, 2547, -25657, 1423, -6316, 1423, 4502, -630, 0, 611], [-10680, 344, -10048, -557, -760, 0, -616, -2024, -2024, 976], [-25354, 838, 50, -860, 700, 860, -1576, 0, -1076, 50], [289, 289, -3443, 115, -6292, 115, -3158, -1978, 0, 2857], [-2511, 615, 1903, -2842, 4402, 2842, 2576, 0, -634, 1903], [-1805, -1805, -6219, -3138, 2176, -3138, -3262, -1782, 0, 717], [2648, 989, -4555, 0, 8077, -2464, -3507, -307, 307, -2896], [1712, -1272, -3962, 0, -938, 2779, -4850, 350, -350, -978], [236, 236, -6628, 1260, 2384, 1260, 4112, -784, 0, -1428], [-1265, -1265, 7827, -2881, 2500, -2881, -5666, 3106, 0, 15], [8708, -940, 11824, -2149, 7164, 0, -3404, 932, 932, 2176], [-10067, -141, 459, 3294, -870, -3294, 1200, 0, -1506, 459], [2819, 2819, -15433, 415, -4140, 415, 4342, 2794, 0, -749], [1536, -544, -96, -2144, -992, 2144, 896, 0, 288, -96], [1939, 1804, -310, -2372, -4133, 0, -3083, 365, 365, -445]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_16_e_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


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