
// Make newform 23.4.a.b in Magma, downloaded from the LMFDB on 24 April 2024.
// To make the character of type GrpDrchElt, type "MakeCharacter_23_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_23_4_a_b();".
// 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_23_4_a_b();".
// 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, -5, -16, -2, 1];
        Kf := NumberField(Polynomial([elt : elt in poly]));
        AssignNames(~Kf, ["nu"]);
    end if;
    Rf_num := [[1, 0, 0, 0], [0, 1, 0, 0], [-10, -20, -1, 1], [-1, 28, 5, -2]];
    Rf_basisdens := [1, 1, 3, 3];
    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_23_a();"
function MakeCharacter_23_a()
    N := 23;
    order := 1;
    char_gens := [5];
    v := [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_23_a_Hecke(Kf)
    return MakeCharacter_23_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, 0, 0, 1], [1, 1, 1, 0], [2, 2, -2, -2], [4, -4, 4, -2], [10, -12, -4, 2], [31, -11, -11, -10], [32, -2, -6, 10], [14, 32, 24, 24], [-23, 0, 0, 0], [-31, 35, 23, -26], [-35, -19, 9, 12], [46, 30, -54, 10], [-49, 21, -7, -18], [-24, 70, 14, 28], [-119, -71, 13, 36], [-118, 14, -26, -120], [-304, 36, 60, 28], [206, 14, -26, 36], [110, 92, 84, 110], [33, -143, -111, 68], [281, 127, -1, -66], [-214, -98, -50, -132], [96, -146, 6, -22], [604, 82, 82, 160], [388, -38, 150, 282], [-618, -72, -156, -80], [776, -244, -36, 50], [706, -12, 164, 28], [-640, -208, -16, -368], [424, -174, -146, -222], [479, -249, -445, -152], [553, -219, -7, 172], [-420, 14, -210, -172], [631, 291, -281, 92], [-1274, 354, -142, 44], [-287, 121, 513, -52], [-224, 316, 372, -180], [-177, 711, -185, 304], [0, -792, -432, -436], [-1158, 208, 532, -92], [31, 175, 235, 176], [-74, 442, 86, -174], [1722, 74, 274, -290], [2001, 355, -81, -178], [-1135, 683, 95, 618], [494, -26, 206, -314], [-324, 32, -336, 348], [-916, 532, -104, -1276], [72, -254, -254, 886], [-1484, -644, -812, -592], [719, 377, 1137, 778], [-471, -119, 153, 472], [1450, -728, -140, 134], [2296, -506, -170, -742], [-2509, 1089, 649, -494], [-242, 902, 318, 20], [-471, 51, 835, 498], [136, -160, 80, 208], [2717, -153, 291, -302], [-584, -494, 182, -522], [750, 152, -632, -70], [-560, -588, -436, -868], [-3696, -208, -624, -796], [1253, 525, -279, -104], [3624, 114, -982, -1056], [1134, -1676, -520, -680], [795, -1193, 115, 1092], [6996, -1254, -326, 1184], [-3664, 1832, 700, -1028], [-2855, 275, -709, 1078], [5005, -2393, -177, 494], [4142, -110, -854, -1388], [-7560, 1432, 384, -94], [618, 1130, 1982, 1666], [56, 946, 26, 96], [6084, 908, 740, 968], [-1224, -1136, 968, -660], [6295, -3387, -267, -126], [-5850, -1636, -1808, -434], [-1799, 1155, -1565, -434], [-4072, 3418, 1578, 1728], [672, -1704, -572, 14], [-6960, -1660, 2052, 724], [2838, 740, -1920, 518], [2167, -537, 671, -1256], [-207, 1921, 1937, 2764], [-4114, -144, -1712, 1416], [-412, -1590, 670, -2182], [-1927, -117, -437, -1922], [708, 228, 2344, 1604], [5126, -152, 992, 576], [-264, -20, -444, 1738], [4939, 1643, 1963, 2024], [9103, 2175, 2507, 740], [-12227, 2173, 9, -240], [9866, -154, -682, -1108], [4725, -1585, -2261, 274], [-11764, -1370, 318, 960], [6344, -3658, -4138, -2606], [-2203, 1263, -593, -130], [-133, 1959, 1187, 292], [2836, 3276, 860, -3068], [6820, -4870, -1374, 2044], [-2, -1500, -2496, 3494], [9804, -4318, -1310, 950], [-695, -1097, 531, -1250], [2609, -5339, -2163, -3132], [-10190, 4932, 2456, 1916], [2568, -2152, 1448, 560], [2153, -4713, -1449, -2150], [1772, 4612, -1172, 2508], [-5328, 164, -2568, 1134], [-3350, 6176, 2016, 2332], [-4272, -802, -1322, -2176], [9236, -1868, -3244, -3996], [-516, 4782, 2078, 4128], [-8822, -2444, 3044, 2718], [-4539, -3723, -1683, 576], [-4749, -2863, -5187, -2146], [-776, 4326, -1754, -2862], [8062, 2222, 4054, 132], [8079, 1213, -739, 1918], [14474, 942, 1942, -2020], [-19497, -389, 1543, 1608], [-9588, 6968, 2336, 2672], [-11992, -2908, -1360, -2350], [3828, -4428, -3712, -4810], [-15424, -1888, 2696, -816], [-1078, 5962, -1782, -2862], [16520, -3160, 444, 718], [-19877, -205, 1259, -564], [7760, 1644, 292, -3010], [3800, -96, 1656, -1804], [-3638, -938, 3674, -4246], [-4531, -2881, -5029, -4510], [5706, 3032, -3476, 206], [-13860, 6844, 1108, 3804], [8108, -6678, -534, 3768], [-13906, -2874, 5422, 4896], [-1890, -316, 1408, -644], [6361, 6181, 1149, 1472], [3590, 3432, 1684, 1308], [16265, -4207, 305, 7684], [-19768, -3838, 2714, 168], [21142, -68, 1000, 6344], [-4272, 5764, 2948, -6518], [13050, 3940, 2816, 6592], [14247, 801, -3207, -2630], [-7455, -3743, -4703, -1516], [-30525, -1213, -437, -4264], [1802, 1472, -5076, 1736], [-9776, -3134, -1274, -930], [-1364, 2600, 2848, -3400], [11989, -1667, -1595, -5012], [-4854, 3688, 8408, 4], [18722, -10994, -5274, -3248], [-9104, -392, -3584, -1738], [671, -1483, 6593, 1962], [22110, -4432, -1500, 274], [-24028, 4468, 5328, 2270], [30117, 5477, 3873, 2176], [-21402, 624, -4112, -4116], [-24989, -2981, 703, 5552], [30402, -3328, 1744, -1908], [-13322, 2916, -3160, 2694], [7598, 3858, -830, 4768], [14552, 272, -4264, 2680], [16602, -11840, -4916, -2548]];
    aps := ConvertToHeckeField(raw_aps);
    chi := MakeCharacter_23_a_Hecke(Universe(aps));
    return ExtendMultiplicatively(weight, aps, chi);
end function;


// To make the newform (type ModFrm), type "MakeNewformModFrm_23_4_a_b();".
// 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_23_4_a_b(:prec:=4)
    chi := MakeCharacter_23_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(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_23_4_a_b();".
// 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_23_4_a_b( : sign := -1)
    R<x> := PolynomialRing(Rationals());
    chi := MakeCharacter_23_a();
    // SetVerbose("ModularSymbols", true);
    Snew := NewSubspace(CuspidalSubspace(ModularSymbols(chi,4,sign)));
    Vf := Kernel([<2,R![2, 61, -24, -2, 1]>],Snew);
    return Vf;
end function;