/*This supplementary file contains details for some of the calculations in "Six-dimensional sphere packing and linear programming" by Matthew de Courcy-Ireland, Maria Dostert, and Maryna Viazovska. Date: July 27, 2023.
The calculations are done in Pari/GP, available at https://pari.math.u-bordeaux.fr/
We used Version 2.11.2 of the GP/Pari Calculator. The estimated runtimes are taken from our machine: Intel(R) Core(TM) i7-4650U CPU @1.70 GHz processor base frequency, 8.00 GB installed RAM, operating system Windows 10 Pro version 1709.
Two other files "checklist-an-67250.txt" and "checklist-bn-4329.txt" are also available on arXiv. They list the 67250 indices n where we need to compute the cuspidal part in order to check the sign of a_n, or respectively the 4329 indices for b_n. For other n, positivity already follows from Deligne's bound, so there is no need to compute the cuspidal part. The reader can download these checklists and open them in Pari in order to reproduce the check that all terms are positive. On Windows, the command L=read(path,"file.txt") will copy the indices into a list L, where "path" describes where the file is saved. For instance, the path could start with "C:\Users\...". On Linux, a shorter command L=read("file.txt") does the same task, assuming "file" is saved in the current directory.
We recommend doing the calculation in stages with reference to this supplement. Copy-pasting the whole thing into Pari may lead to errors depending on where the checklists are saved.
Contents of this file:
First, we describe how to obtain x_1,...,x_44 as quadratic irrationals. These are converted to floating-point, keeping 38 digits, which are the values we used to obtain Proposition 5.2 (upper bound for the cuspidal part), as well as Proposition 7.1 and Table 7.1 (lower bound for the Eisenstein part).
Second, we substitute these values as described on page 22 to obtain Proposition 5.2.
Third, we derive Table 7.1. This involves several cases, which we unify by writing matrices that change the appropriate variables x_j to 0. This way we can use formula (7.17) throughout, applied to different vectors x. Taking the minimum of these bounds gives Proposition 7.1.
Fourth, using the values from Propositions 5.1 and 7.1, we solve equation (2.14) numerically. For n larger than this solution, we have a_n >=0 as required (or similarly for b_n >= 0). We call this initial bound "bound1".
Fifth, for n until bound1, we compute the number of divisors of n and the Eisenstein part (which is usually much larger than the worst case eps*n^2). If the Eisenstein part is larger than Deligne's bound, then positivity follows. Adding the remaining values n to a to-do list yields "checklist-an-67250.txt" and "checklist-bn-4329.txt". We call the largest value on the to-do list "bound2".
Sixth, we compute the cuspidal part using (5.4) for n in "checklist-an-67250.txt" and compare it to the Eisenstein part. The result is that all a_n are positive, except for the vanishing coefficients from (4.2) and when n=1 mod 4. A similar, but smaller, calculation shows that b_n >= 0 for all n with equality only for n from (4.3).
*/


/* Section 4: Here we set up the system of equations from page 18, equations (4.1),(4.2),(4.3) 
First, we write the eleven equations forcing a_n=0 for n=1 mod 4. A1mod4 is a matrix representing this system. We generate it by initializing eleven vectors at the zero vector of size 44, and then changing the corresponding entries. For instance, v10 will be one row of A1mod4, corresponding to the equation x_10=0. The vectors are labeled by the indices of the variables that enter, such as v2329 corresponding to x_23+x_29=0
*/

v10=vector(44,i,0);
v10[10]=1;
v41=vector(44,i,0);
v41[41]=1;
v43=vector(44,i,0);
v43[43]=1;
v111=vector(44,i,0);
v111[1]=1;
v111[11]=1;
v612=vector(44,i,0);
v612[6]=1;
v612[12]=1;
v2329=vector(44,i,0);
v2329[23]=1;
v2329[29]=1;
v2531=vector(44,i,0);
v2531[25]=1;
v2531[31]=-1;
v3538=vector(44,i,0);
v3538[35]=1;
v3538[38]=-1;
v132021=vector(44,i,0);
v132021[13]=1;
v132021[20]=1;
v132021[21]=3;
v162022=vector(44,i,0);
v162022[16]=1;
v162022[20]=1;
v162022[22]=2;
v182122=vector(44,i,0);
v182122[18]=1;
v182122[21]=-1;
v182122[22]=-2;
A1mod4=matconcat([v10;v41;v43;v111;v612;v2329;v2531;v3538;v132021;v162022;v182122]);

/* Next, we write a matrix corresponding to the equations a_n=0 or b_n=0 for the specific values n from (4.2) and (4.3)
"force" lists the values from (4.2) where we force a_n=0, except that the positions are increased by 1. This is because we start counting from q^0 in the q-series, but the vectors are indexed starting from 1. Thus "force" includes 1933 whereas (4.2) stops at n=1932. 
We also include 1, that is n=0, because (4.1) fixes the value a_0.
"transforce" lists the places where we force b_n=0, i.e. the transform is forced to vanish. Again, note the shift by 1 to [2,3,4,...,83] instead of [1,2,3,...,82].
"prec" is the number of coefficients we compute in the q-series. We set prec=2000 because we only need to go up to q^1932
M3,M4 are the subspaces of modular forms of level 48 and weight 3 for the two characters we used. B3 and B4 are the bases described in Section 3, with B3=(f1,...,f22) and B4=(f23,...,f44). The Atkin-Lehner transforms are computed using the data initialized in "atk3" and "atk4".
A3 is a matrix, initially empty. We then add a row for each of the basis functions from B3, listing the coefficients up to q^prec. Likewise A4 does the same thing for the other subspace, recording the coefficients of the basis functions from B4.
"dual3" and "dual4" multiply by the matrices from (3.12) to compute the transforms.
"dual4" involves a factor 1/sqrt(12). We represent the square root exactly as a polynomial Mod(t,t^2-12) in the quotient ring modulo t^2-12.
"sys3","sys4" are the submatrices listing only the coefficients q^n for n in (4.2), or equivalently the indices in "force".
"leh3" and "leh4" (for Lehner) do the same thing for the transforms, recording only the indices in "transforce".
"A" then concatenates these four matrices together to write the system of equations (4.1),(4.2),(4.3).
Computing dual4 takes more memory because it involves operations with polynomials instead of integers. It may be necessary to allocate memory, for instance allocatemem(16000000) was enough for us but the initial size of 8 Mbytes was not.
*/

allocatemem(16000000)
force=[1, 3, 4, 5, 7, 9, 11, 12, 13, 23, 27, 33, 39, 61, 65, 89, 91, 93, 107, 165, 1933];
transforce=[2, 3, 4, 5, 8, 9, 10, 11, 14, 15, 37, 83];
prec=2000;
M3=mfinit([48,3,-3]);
M4=mfinit([48,3,-4]);
B3=mfbasis(M3);
B4=mfbasis(M4);
atk3=mfatkininit(M3,48);
atk4=mfatkininit(M4,48);
A3=[];
for(j=1,length(B3),A3=matconcat([A3,mfcoefs(B3[j],prec)~]));
A4=[];
for(j=1,length(B4),A4=matconcat([A4,mfcoefs(B4[j],prec)~]));
dual3=A3*atk3[2]/4;
dual4=A4*atk4[2]/Mod(t,t^2-12);
sys3=matrix(length(force),length(B3),{i},{j},A3[force[i],j]);
sys4=matrix(length(force),length(B4),{i},{j},A4[force[i],j]);
leh3=matrix(length(transforce),length(B3),{i},{j},dual3[transforce[i],j]);
leh4=matrix(length(transforce),length(B4),{i},{j},dual4[transforce[i],j]);
A=matconcat([sys3,sys4;leh3,leh4]);

/* Lastly, we concatenate A1mod4 and A to get the complete system of 44 equations in 44 unknowns.
The right side of the equation has a 1 corresponding to a_0=1, and 0 in all other entries corresponding to a_n=0 or b_n=0. The equation a_0=1 corresponds to row 12 because we first write 11 equations from A1mod4. 
The transpose rhs~ is a column vector rather than a row.
*/

total=matconcat([A1mod4;A]);
rhs=vector(44,i,0);
rhs[12]=1;
x=matinverseimage(total,rhs~);

/* Output: typing "x" or "print(x)" should give the following output, which is a vector listing x_1,...,x_44 as exact elements of Q(sqrt(12))

[Mod(-5933985634895097866727312605651608472538244579636113805857/11160420408110884613435465825574808257902091756958813074304*t + 13305835801389714487894701329762261844642710724107326137255/5580210204055442306717732912787404128951045878479406537152, t^2 - 12), Mod(13171000968786893989488139647060543539330001150557422035099/11160420408110884613435465825574808257902091756958813074304*t - 29741032077010592910877726160603830062567867730204045307757/5580210204055442306717732912787404128951045878479406537152, t^2 - 12), Mod(10132495874189428928090966852619558774241896818329744462023/2790105102027721153358866456393702064475522939239703268576*t - 22292047118768347600029751974280125808098820591543738212525/1395052551013860576679433228196851032237761469619851634288, t^2 - 12), Mod(814026652012869083620458982244706792248044613750691468387/348763137753465144169858307049212758059440367404962908572*t - 1711011024264086482792146424552345491552564840439615696317/174381568876732572084929153524606379029720183702481454286, t^2 - 12), Mod(263567153417509190329320343543021790669072496873146686632/87190784438366286042464576762303189514860091851240727143*t - 983965156692193542319291976604278249460018244714403584608/87190784438366286042464576762303189514860091851240727143, t^2 - 12), Mod(-2745919821869583758552049373060736525798014416696636655621/100443783672997961520919192430173274321118825812629317668736*t - 23502822471985112282969606830752930319270029191986053175/16740630612166326920153198738362212386853137635438219611456, t^2 - 12), Mod(-3047672686969449855323754208360647216023554499866603403053/33481261224332653840306397476724424773706275270876439222912*t + 240558088070039545206788197900466642576699738522224083233/5580210204055442306717732912787404128951045878479406537152, t^2 - 12), Mod(-5517032999709669238020008899494318924515776276133032608197/25110945918249490380229798107543318580279706453157329417184*t + 5336488169994454549837231082175980522542493295768760295197/4185157653041581730038299684590553096713284408859554902864, t^2 - 12), Mod(-605625130304386369752208901277281388547982505328027071003/348763137753465144169858307049212758059440367404962908572*t + 1120940151926879699971180303832582709364168366395305489361/174381568876732572084929153524606379029720183702481454286, t^2 - 12), Mod(0, t^2 - 12), Mod(5933985634895097866727312605651608472538244579636113805857/11160420408110884613435465825574808257902091756958813074304*t - 13305835801389714487894701329762261844642710724107326137255/5580210204055442306717732912787404128951045878479406537152, t^2 - 12), Mod(2745919821869583758552049373060736525798014416696636655621/100443783672997961520919192430173274321118825812629317668736*t + 23502822471985112282969606830752930319270029191986053175/16740630612166326920153198738362212386853137635438219611456, t^2 - 12), Mod(-1426190309080687265007923744926769810663547620624569770905/16740630612166326920153198738362212386853137635438219611456*t - 896643704179693339832127950724646066150918706953523948679/2790105102027721153358866456393702064475522939239703268576, t^2 - 12), Mod(-2480219716894363087304574521184681706502131796569464525/174381568876732572084929153524606379029720183702481454286*t + 8782265378193457685235334588751937702810273511753210545/174381568876732572084929153524606379029720183702481454286, t^2 - 12), Mod(-153039165054267610462262996232307115971783403466816145233/87190784438366286042464576762303189514860091851240727143*t + 571335252022573469495287623233404129474409716726681663452/87190784438366286042464576762303189514860091851240727143, t^2 - 12), Mod(-1954017610707785801673793647269224544943335807657407657239/11160420408110884613435465825574808257902091756958813074304*t + 3604525000541974009440676832555762261487191891529196696899/5580210204055442306717732912787404128951045878479406537152, t^2 - 12), Mod(1601431679212579497405818054822214374081433728912523682473/2790105102027721153358866456393702064475522939239703268576*t - 4243200568646273955690618974257202679305650430012061446667/1395052551013860576679433228196851032237761469619851634288, t^2 - 12), Mod(-369294707713032534207807001086684956178418603791481541439/11160420408110884613435465825574808257902091756958813074304*t + 706615785968010167399492585560970242663027459995833740407/5580210204055442306717732912787404128951045878479406537152, t^2 - 12), Mod(-164097681559806670461665624532779756100994708803890523633/2790105102027721153358866456393702064475522939239703268576*t + 408751531961889884011640721366936988786571852868155632395/1395052551013860576679433228196851032237761469619851634288, t^2 - 12), Mod(2970273935493592444118781665632090641427860618036142791239/16740630612166326920153198738362212386853137635438219611456*t - 348129217224080182839257886028153293129650645020876356705/697526275506930288339716614098425516118880734809925817144, t^2 - 12), Mod(-257347271068817529851809653450886805127385499568595503389/8370315306083163460076599369181106193426568817719109805728*t + 763053524358671357063053164945753079556507095679009791833/2790105102027721153358866456393702064475522939239703268576, t^2 - 12), Mod(-78495038863827483216182389456507648025713813100062610761/66962522448665307680612794953448849547412550541752878445824*t - 819491262749332546726613744330535916449986731362185843259/11160420408110884613435465825574808257902091756958813074304, t^2 - 12), Mod(10807063769758740400626509384771581374585972529760214991307/43944155356936608165402146688200807515489486293025326480072*t - 1732584480024296855852316383390254238122707411632683503514/1831006473205692006891756112008366979812061928876055270003, t^2 - 12), Mod(-510318139695774051766380132532520671673889630134983764265/784717059945296574382181190860728705633740826661166544287*t + 4235883172963293580861386543573681472741270805825058948910/1831006473205692006891756112008366979812061928876055270003, t^2 - 12), Mod(553417968418361081537810462397597051264777747377480253153/2441341964274256009189008149344489306416082571834740360004*t - 710626614778322692764546272122106694551527589354719978080/610335491068564002297252037336122326604020642958685090001, t^2 - 12), Mod(-1224180889292205091406274783824358576928025678515027413960/610335491068564002297252037336122326604020642958685090001*t + 5124242633266170779011837591667174935442674116660421020864/610335491068564002297252037336122326604020642958685090001, t^2 - 12), Mod(-521895255451610272432474728484337891782930815999086020827/174381568876732572084929153524606379029720183702481454286*t + 8517869123108574299605122800807366484968545173716970683181/610335491068564002297252037336122326604020642958685090001, t^2 - 12), Mod(6258385265770172502845373064233750340189569032043090233244/610335491068564002297252037336122326604020642958685090001*t - 24611313159124722413159544840035552288092765140756074302872/610335491068564002297252037336122326604020642958685090001, t^2 - 12), Mod(-10807063769758740400626509384771581374585972529760214991307/43944155356936608165402146688200807515489486293025326480072*t + 1732584480024296855852316383390254238122707411632683503514/1831006473205692006891756112008366979812061928876055270003, t^2 - 12), Mod(1791682177206184586371179832799925814034617478461617398321/6277736479562372595057449526885829645069926613289332354296*t - 1971451158233066058695131614622533826900625188617562396528/1831006473205692006891756112008366979812061928876055270003, t^2 - 12), Mod(553417968418361081537810462397597051264777747377480253153/2441341964274256009189008149344489306416082571834740360004*t - 710626614778322692764546272122106694551527589354719978080/610335491068564002297252037336122326604020642958685090001, t^2 - 12), Mod(-15642286620520251780689356039142963588614947335085282759/1220670982137128004594504074672244653208041285917370180002*t + 18874324690421761821155154735665186885294175404325275706/610335491068564002297252037336122326604020642958685090001, t^2 - 12), Mod(-904512630300165417999544896003988506257382332006005360533/348763137753465144169858307049212758059440367404962908572*t + 7245626133429246403592717906335845311936640951261051571119/610335491068564002297252037336122326604020642958685090001, t^2 - 12), Mod(3982034714418952109205942493982417312298905762199881488677/1220670982137128004594504074672244653208041285917370180002*t - 11016474717997110868944983100532479461431320199233292645179/610335491068564002297252037336122326604020642958685090001, t^2 - 12), Mod(591104979895227299877935043039993733234955953830791661821/4882683928548512018378016298688978612832165143669480720008*t - 2168849440179559866933651832071419098735163496136076020655/2441341964274256009189008149344489306416082571834740360004, t^2 - 12), Mod(111650893441555058861488282013051158492007880464957194627/174381568876732572084929153524606379029720183702481454286*t - 2155370118507177384270842294348871264403260927385921515715/610335491068564002297252037336122326604020642958685090001, t^2 - 12), Mod(1844796611267723774926172431754667680566704615300279097900/610335491068564002297252037336122326604020642958685090001*t - 7334127124390146203257799821958616772197943389948584678216/610335491068564002297252037336122326604020642958685090001, t^2 - 12), Mod(591104979895227299877935043039993733234955953830791661821/4882683928548512018378016298688978612832165143669480720008*t - 2168849440179559866933651832071419098735163496136076020655/2441341964274256009189008149344489306416082571834740360004, t^2 - 12), Mod(529826174851738759618749660769841051533531007089862423299/697526275506930288339716614098425516118880734809925817144*t - 8435422424938455727598828092621655077316868661614608673655/2441341964274256009189008149344489306416082571834740360004, t^2 - 12), Mod(224678751381599950786492736474636280403894321816104205396/610335491068564002297252037336122326604020642958685090001*t - 204467917753821226838271594759441987634554724282777376024/87190784438366286042464576762303189514860091851240727143, t^2 - 12), Mod(0, t^2 - 12), Mod(1547878207279216966411867453673435218909977962703059346123/174381568876732572084929153524606379029720183702481454286*t - 4252834879046645182610368630252807407080103043428862903119/87190784438366286042464576762303189514860091851240727143, t^2 - 12), Mod(0, t^2 - 12), Mod(222468843388209209251464647056995626985567725610824359232/261572353315098858127393730286909568544580275553722181429*t - 1207204136558719277034054288103014283443918973665398147893/261572353315098858127393730286909568544580275553722181429, t^2 - 12)]~

*/

/* Substituting a 38-digit approximation to sqrt(12) for t gives the values from Table 4.1 (or rather, the values we truncated to make Table 4.1). See the vector "x" below with entries x1,...,x44.
[-5933985634895097866727312605651608472538244579636113805857/11160420408110884613435465825574808257902091756958813074304*sqrt(12) + 13305835801389714487894701329762261844642710724107326137255/5580210204055442306717732912787404128951045878479406537152,  13171000968786893989488139647060543539330001150557422035099/11160420408110884613435465825574808257902091756958813074304*sqrt(12) - 29741032077010592910877726160603830062567867730204045307757/5580210204055442306717732912787404128951045878479406537152,  10132495874189428928090966852619558774241896818329744462023/2790105102027721153358866456393702064475522939239703268576*sqrt(12) - 22292047118768347600029751974280125808098820591543738212525/1395052551013860576679433228196851032237761469619851634288,  814026652012869083620458982244706792248044613750691468387/348763137753465144169858307049212758059440367404962908572*sqrt(12) - 1711011024264086482792146424552345491552564840439615696317/174381568876732572084929153524606379029720183702481454286,  263567153417509190329320343543021790669072496873146686632/87190784438366286042464576762303189514860091851240727143*sqrt(12) - 983965156692193542319291976604278249460018244714403584608/87190784438366286042464576762303189514860091851240727143,  -2745919821869583758552049373060736525798014416696636655621/100443783672997961520919192430173274321118825812629317668736*sqrt(12) - 23502822471985112282969606830752930319270029191986053175/16740630612166326920153198738362212386853137635438219611456,  -3047672686969449855323754208360647216023554499866603403053/33481261224332653840306397476724424773706275270876439222912*sqrt(12) + 240558088070039545206788197900466642576699738522224083233/5580210204055442306717732912787404128951045878479406537152,  -5517032999709669238020008899494318924515776276133032608197/25110945918249490380229798107543318580279706453157329417184*sqrt(12) + 5336488169994454549837231082175980522542493295768760295197/4185157653041581730038299684590553096713284408859554902864,  -605625130304386369752208901277281388547982505328027071003/348763137753465144169858307049212758059440367404962908572*sqrt(12) + 1120940151926879699971180303832582709364168366395305489361/174381568876732572084929153524606379029720183702481454286,  0,  5933985634895097866727312605651608472538244579636113805857/11160420408110884613435465825574808257902091756958813074304*sqrt(12) - 13305835801389714487894701329762261844642710724107326137255/5580210204055442306717732912787404128951045878479406537152,  2745919821869583758552049373060736525798014416696636655621/100443783672997961520919192430173274321118825812629317668736*sqrt(12) + 23502822471985112282969606830752930319270029191986053175/16740630612166326920153198738362212386853137635438219611456,  -1426190309080687265007923744926769810663547620624569770905/16740630612166326920153198738362212386853137635438219611456*sqrt(12) - 896643704179693339832127950724646066150918706953523948679/2790105102027721153358866456393702064475522939239703268576,  -2480219716894363087304574521184681706502131796569464525/174381568876732572084929153524606379029720183702481454286*sqrt(12) + 8782265378193457685235334588751937702810273511753210545/174381568876732572084929153524606379029720183702481454286,  -153039165054267610462262996232307115971783403466816145233/87190784438366286042464576762303189514860091851240727143*sqrt(12) + 571335252022573469495287623233404129474409716726681663452/87190784438366286042464576762303189514860091851240727143,  -1954017610707785801673793647269224544943335807657407657239/11160420408110884613435465825574808257902091756958813074304*sqrt(12) + 3604525000541974009440676832555762261487191891529196696899/5580210204055442306717732912787404128951045878479406537152,  1601431679212579497405818054822214374081433728912523682473/2790105102027721153358866456393702064475522939239703268576*sqrt(12) - 4243200568646273955690618974257202679305650430012061446667/1395052551013860576679433228196851032237761469619851634288,  -369294707713032534207807001086684956178418603791481541439/11160420408110884613435465825574808257902091756958813074304*sqrt(12) + 706615785968010167399492585560970242663027459995833740407/5580210204055442306717732912787404128951045878479406537152,  -164097681559806670461665624532779756100994708803890523633/2790105102027721153358866456393702064475522939239703268576*sqrt(12) + 408751531961889884011640721366936988786571852868155632395/1395052551013860576679433228196851032237761469619851634288,  2970273935493592444118781665632090641427860618036142791239/16740630612166326920153198738362212386853137635438219611456*sqrt(12) - 348129217224080182839257886028153293129650645020876356705/697526275506930288339716614098425516118880734809925817144,  -257347271068817529851809653450886805127385499568595503389/8370315306083163460076599369181106193426568817719109805728*sqrt(12) + 763053524358671357063053164945753079556507095679009791833/2790105102027721153358866456393702064475522939239703268576,  -78495038863827483216182389456507648025713813100062610761/66962522448665307680612794953448849547412550541752878445824*sqrt(12) - 819491262749332546726613744330535916449986731362185843259/11160420408110884613435465825574808257902091756958813074304,  10807063769758740400626509384771581374585972529760214991307/43944155356936608165402146688200807515489486293025326480072*sqrt(12) - 1732584480024296855852316383390254238122707411632683503514/1831006473205692006891756112008366979812061928876055270003,  -510318139695774051766380132532520671673889630134983764265/784717059945296574382181190860728705633740826661166544287*sqrt(12) + 4235883172963293580861386543573681472741270805825058948910/1831006473205692006891756112008366979812061928876055270003,  553417968418361081537810462397597051264777747377480253153/2441341964274256009189008149344489306416082571834740360004*sqrt(12) - 710626614778322692764546272122106694551527589354719978080/610335491068564002297252037336122326604020642958685090001,  -1224180889292205091406274783824358576928025678515027413960/610335491068564002297252037336122326604020642958685090001*sqrt(12) + 5124242633266170779011837591667174935442674116660421020864/610335491068564002297252037336122326604020642958685090001,  -521895255451610272432474728484337891782930815999086020827/174381568876732572084929153524606379029720183702481454286*sqrt(12) + 8517869123108574299605122800807366484968545173716970683181/610335491068564002297252037336122326604020642958685090001,  6258385265770172502845373064233750340189569032043090233244/610335491068564002297252037336122326604020642958685090001*sqrt(12) - 24611313159124722413159544840035552288092765140756074302872/610335491068564002297252037336122326604020642958685090001,  -10807063769758740400626509384771581374585972529760214991307/43944155356936608165402146688200807515489486293025326480072*sqrt(12) + 1732584480024296855852316383390254238122707411632683503514/1831006473205692006891756112008366979812061928876055270003,  1791682177206184586371179832799925814034617478461617398321/6277736479562372595057449526885829645069926613289332354296*sqrt(12) - 1971451158233066058695131614622533826900625188617562396528/1831006473205692006891756112008366979812061928876055270003,  553417968418361081537810462397597051264777747377480253153/2441341964274256009189008149344489306416082571834740360004*sqrt(12) - 710626614778322692764546272122106694551527589354719978080/610335491068564002297252037336122326604020642958685090001,  -15642286620520251780689356039142963588614947335085282759/1220670982137128004594504074672244653208041285917370180002*sqrt(12) + 18874324690421761821155154735665186885294175404325275706/610335491068564002297252037336122326604020642958685090001,  -904512630300165417999544896003988506257382332006005360533/348763137753465144169858307049212758059440367404962908572*sqrt(12) + 7245626133429246403592717906335845311936640951261051571119/610335491068564002297252037336122326604020642958685090001,  3982034714418952109205942493982417312298905762199881488677/1220670982137128004594504074672244653208041285917370180002*sqrt(12) - 11016474717997110868944983100532479461431320199233292645179/610335491068564002297252037336122326604020642958685090001,  591104979895227299877935043039993733234955953830791661821/4882683928548512018378016298688978612832165143669480720008*sqrt(12) - 2168849440179559866933651832071419098735163496136076020655/2441341964274256009189008149344489306416082571834740360004,  111650893441555058861488282013051158492007880464957194627/174381568876732572084929153524606379029720183702481454286*sqrt(12) - 2155370118507177384270842294348871264403260927385921515715/610335491068564002297252037336122326604020642958685090001,  1844796611267723774926172431754667680566704615300279097900/610335491068564002297252037336122326604020642958685090001*sqrt(12) - 7334127124390146203257799821958616772197943389948584678216/610335491068564002297252037336122326604020642958685090001,  591104979895227299877935043039993733234955953830791661821/4882683928548512018378016298688978612832165143669480720008*sqrt(12) - 2168849440179559866933651832071419098735163496136076020655/2441341964274256009189008149344489306416082571834740360004,  529826174851738759618749660769841051533531007089862423299/697526275506930288339716614098425516118880734809925817144*sqrt(12) - 8435422424938455727598828092621655077316868661614608673655/2441341964274256009189008149344489306416082571834740360004,  224678751381599950786492736474636280403894321816104205396/610335491068564002297252037336122326604020642958685090001*sqrt(12) - 204467917753821226838271594759441987634554724282777376024/87190784438366286042464576762303189514860091851240727143,  0,  1547878207279216966411867453673435218909977962703059346123/174381568876732572084929153524606379029720183702481454286*sqrt(12) - 4252834879046645182610368630252807407080103043428862903119/87190784438366286042464576762303189514860091851240727143,  0,  222468843388209209251464647056995626985567725610824359232/261572353315098858127393730286909568544580275553722181429*sqrt(12) - 1207204136558719277034054288103014283443918973665398147893/261572353315098858127393730286909568544580275553722181429]~
*/

/*Section 5: bounding the cuspidal part
The following calculation yields the numbers 21.6161 and 24.0266 from Proposition 5.2.
W is the matrix from (3.12) carrying out the Atkin-Lehner transform; first we initialize M3, M4, atk3, atk4 as above.
y=W*x~ is the vector (column vector) of coefficients needed to write the transform as \sum_j y_j f_j
u,v are the vectors with entries |x_i| or |y_i| appearing from the triangle inequality.
c13,c16,...,c44 are the numbers from page 22 (Proof of Proposition 5.2). The basis function f_i has n-th Fourier coefficient of size at most c_i times Deligne's bound n*sigma_0(n)
The corresponding linear combinations are bounds depending on gcd(n,48), in fact only on gcd(n,12). If n is divisible by 12, we have a12. The others are labeled so that a6 applies when gcd(n,48)=6, and so on. The bounds for b_n are obtained the same way, but with the vector v in place of u.
The maximum values are a12 and b12, which we round up to obtain Proposition 5.2. This is computed using "proposition52a" and "proposition52b", where vecsort sorts the given values in increasing order. Each has 6 entries a1,a2,a3,a4,a6,a12 or b1,b2,b3,b4,b6,b12, so vecsort puts the maximum in the sixth position, hence the final command "[6]". 
*/

M3=mfinit([48,3,-3]);
M4=mfinit([48,3,-4]);
B3=mfbasis(M3);
B4=mfbasis(M4);
atk3=mfatkininit(M3,48);
atk4=mfatkininit(M4,48);
W=matconcat([atk3[2]/4,0;0,atk4[2]/sqrt(12)]);
x=[0.54260880498140096513867653943544603187, -1.2415641990507267991184361399538105655, -3.3991905564732487151430228250758675037, -1.7265328337275850241386961266251401995, -0.81363822218490473489145562398280189476, -0.096105123401520619760841407689890651869, -0.27221493503052813722741181942257235873, 0.51401353564796162790999644420796668884, 0.41269645847740238263636892690450016577, 0, -0.54260880498140096513867653943544603187, 0.096105123401520619760841407689890651869, -0.61648396750108306711304367224663420226, 0.0010926169103219513034252187992627946530, 0.47243563002773065942505045637470491711, 0.039437084243136701318315982461106041564, -1.0533219944523385506109234372658713817, 0.012002879245312781223893749659513914570, 0.089262594448894739650344612278568599360, 0.11554101925893365487160361897119358347, 0.16698098274738313741381335109181353960, -0.077489051751035178094959800716149812514, -0.094330188558529879564553410275673781724, 0.060639194906309794623804168766865979012, -0.37905807314629534412480287497367879036, 1.4476556752440492028213652494914163684, 3.5885598403833256663248165887542015140, -4.8033101379755878918290347229723305568, 0.094330188558529879564553410275673781724, -0.088040011215475046862305376778881572847, -0.37905807314629534412480287497367879036, -0.013466217519949019225100408189494107099, 2.8874442794684106808734151911975535527, -6.7493834708978057073626986721978412391, -0.46901482838125511797246129933033163853, -1.3134988693595305093541618640494079449, -1.5459762991088610300938595986576085849, -0.46901482838125511797246129933033163853, -0.82398184122267052703537196006578932030, -1.0698466798734160754078802044282581803, 0, -18.027492128156189368154902560159294763, 0, -1.6689434159500542284358558005305559402];
y=W*x~;
u=vector(44,i,abs(x[i]));
v=vector(44,i,abs(y[i]));
c13=1;
c16=2;
c18=6;
c20=24/7;
c21=24/7+2*12*(2+sqrt(2))/7;
c22=8*sqrt(2);
c35=2;
c38=4;
c41=1;
c43=2;
c44=2*sqrt(3);
a12=(u[13]+u[14]/2+u[15]/4)*c13+(u[16]+u[17]/2)*c16+(u[18]+u[19]/2)*c18+(u[35]+u[36]/2+u[37]/4)*c35+(u[38]+u[39]/2+u[40]/4)*c38+(u[41]+u[42]/3)*c41+u[43]*c43+u[44]*c44;
a6=(u[13]+u[14]/2)*c13+(u[16]+u[17]/2)*c16+(u[18]+u[19]/2)*c18+(u[35]+u[36]/2)*c35+(u[38]+u[39]/2)*c38+(u[41]+u[42]/3)*c41+u[43]*c43+u[44]*c44;
a4=(u[13]+u[14]/2+u[15]/4)*c13+(u[16]+u[17]/2)*c16+(u[18]+u[19]/2)*c18+(u[35]+u[36]/2+u[37]/4)*c35+(u[38]+u[39]/2+u[40]/4)*c38+(u[41])*c41+u[43]*c43+u[44]*c44;
a2=(u[13]+u[14]/2)*c13+(u[16]+u[17]/2)*c16+(u[18]+u[19]/2)*c18+(u[35]+u[36]/2)*c35+(u[38]+u[39]/2)*c38+u[41]*c41+u[43]*c43+u[44]*c44;
a3=(u[13])*c13+(u[16])*c16+(u[18])*c18+u[20]*c20+u[21]*c21+u[22]*c22+(u[35])*c35+(u[38])*c38+(u[41]+u[42]/3)*c41+u[43]*c43+u[44]*c44;
a1=(u[13])*c13+(u[16])*c16+(u[18])*c18+u[20]*c20+u[21]*c21+u[22]*c22+(u[35])*c35+(u[38])*c38+(u[41])*c41+u[43]*c43+u[44]*c44;
b12=(v[13]+v[14]/2+v[15]/4)*c13+(v[16]+v[17]/2)*c16+(v[18]+v[19]/2)*c18+(v[35]+v[36]/2+v[37]/4)*c35+(v[38]+v[39]/2+v[40]/4)*c38+(v[41]+v[42]/3)*c41+v[43]*c43+v[44]*c44;
b6=(v[13]+v[14]/2)*c13+(v[16]+v[17]/2)*c16+(v[18]+v[19]/2)*c18+(v[35]+v[36]/2)*c35+(v[38]+v[39]/2)*c38+(v[41]+v[42]/3)*c41+v[43]*c43+v[44]*c44;
b4=(v[13]+v[14]/2+v[15]/4)*c13+(v[16]+v[17]/2)*c16+(v[18]+v[19]/2)*c18+(v[35]+v[36]/2+v[37]/4)*c35+(v[38]+v[39]/2+v[40]/4)*c38+(v[41])*c41+v[43]*c43+v[44]*c44;
b2=(v[13]+v[14]/2)*c13+(v[16]+v[17]/2)*c16+(v[18]+v[19]/2)*c18+(v[35]+v[36]/2)*c35+(v[38]+v[39]/2)*c38+v[41]*c41+v[43]*c43+v[44]*c44;
b3=(v[13])*c13+(v[16])*c16+(v[18])*c18+v[20]*c20+v[21]*c21+v[22]*c22+(v[35])*c35+(v[38])*c38+(v[41]+v[42]/3)*c41+v[43]*c43+v[44]*c44;
b1=(v[13])*c13+(v[16])*c16+(v[18])*c18+v[20]*c20+v[21]*c21+v[22]*c22+(v[35])*c35+(v[38])*c38+(v[41])*c41+v[43]*c43+v[44]*c44;
proposition52a=vecsort([a1,a2,a3,a4,a6,a12])[6];
proposition52b=vecsort([b1,b2,b3,b4,b6,b12])[6];
proposition52=["For equation (5.1)",proposition52a;"For equation (5.2)",proposition52b]


/* Output: printing "proposition52" should give the following output. These numbers are rounded up to obtain Proposition 5.2, equations (5.1) and (5.2).
[21.616080382225432441406830453981038151]

[24.026534641158941800404376491923693398]
a12 and b12 give the numerical values from Proposition 5.2, equations (5.1) and (5.2). We computed the others hoping for a better ratio C/epsilon depending on gcd(n,48), but it turns out that the worst case for both the cuspidal bound and the Eisenstein bound is when n is a multiple of 12. 
*/

/* Section 7: lower bound for the Eisenstein part
The next stage in the calculation gives the numbers from Proposition 7.1, Table 7.1, and Proposition 7.2. 
low3,lowr,uppr are numerical values for the products over primes from Proposition 7.2. They were obtained using Pari's command "prodeuler".
We write matrices "Tab" indexed by the exponents a,b in n=2^a*3^b*n_0. Their effect is to change the appropriate variables x_j to 0 so that we can use the same formula (7.17) even when n is not sufficiently divisible by 2 or 3.
We check signs via g = if(t2 > 0, lowr, uppr) to determine whether to bound the ratio by lowr or uppr.
The bounds are written as a_{ij} or b_{ij}, where i,j are the exponents from n=2^i*3^j*n_0 (what we called a,b above), and the a or b now refers to a_n or b_n. Thus a10 and so on give the values "epsilon" in Table 7.1, whereas b10 etc. give "epsilon*" 
*/


/* We compute the products over primes up to 10^7, and then round up for lower bounds or round down for upper bounds. To see that primes up to this point are enough, the error in truncating the products can be estimated in the standard way by taking logarithms and comparing sums to integrals.
prodeuler(p=5,10000000,1-(Mod(p,3)==-1)/p^2)
0.94290850249008198537799490879916227120
prodeuler(p=7,10000000,1+(Mod(p,12)==7)*(-1+(p^2-1)/(p^2+1)))
0.94999004332694256771696237556537370703
prodeuler(p=5,10000000,1+(Mod(p,12)==5)*(-1+(p^2+1)/(p^2-1)))
1.0969540080545367239896571654940750770
*/

low3=0.9429085;
lowr=0.94999;
uppr=1.096954;

s1(v)=4*v[1]+v[2]+v[3]/4+v[4]/4^2+v[5]/4^3;
s2(v)=9*(v[23]+v[24]/4+v[26]/4^2) + v[25]+v[27]/4+v[28]/4^2;
s3(v)=4*v[6]-v[7]+v[8]/4-v[9]/4^2+v[10]/4^3;
s4(v)=v[23]+v[24]/4+v[26]/4^2-v[25]-v[27]/4-v[28]/4^2;
s5(v)=v[1]-v[2]+v[3]-v[4]+v[5];
s6(v)=9*(v[29]+v[30]+v[32])-v[31]-v[33]-v[34];
s7(v)=v[6]+v[7]+v[8]+v[9]+v[10];
s8(v)=v[29]+v[30]+v[32]+v[31]+v[33]+v[34];

X3(a,b,v,c) = 1/5*(4^a*9^b*s1(v)+4^a*(-1)^a*c*s3(v)+9^b*(-1)^a*s5(v)+c*s7(v));
X4(a,b,v,c) = 1/10*(4^a*9^b*s2(v)+4^a*(-1)^b*s4(v)+ 9^b*(-1)^b*c*s6(v)+c*s8(v));


/* a = 0 b = 0 */
T00= matrix(44, 44,{i},{j},0);
T00[1,1] = 1;
T00[6,6] = 1;
T00[23,23] = 1;
T00[29,29] = 1;

/* a = 0 b >= 0 */
T01 = T00;
T01[25,25] = 1;
T01[31,31] = 1;

/* a = 1 b = 0 */
T10 = T00;
T10[2,2] = 1;
T10[7,7] = 1;
T10[24,24] = 1;
T10[30,30] = 1;

/* a = 1 b >= 1 */
T11 = T10;
T11[25,25] = 1;
T11[31,31] = 1;
T11[27,27] = 1;
T11[33,33] = 1;

/* a= 1 b = 0 */
T10 = T00;
T10[2,2] = 1;
T10[7,7] = 1;
T10[24,24] = 1;
T10[30,30] = 1;

/* a = 2 b = 0 */
T20 = T10;
T20[3,3] = 1;
T20[8,8] = 1;
T20[26,26] = 1;
T20[32,32] = 1;

/* a = 3 b = 0 */
T30 = T20;
T30[4,4] = 1;
T30[9,9] = 1;

/* a = 4 b = 0 */
T40 = T30;
T40[5,5] = 1;
T40[10,10] = 1;

/* a = 3 b >= 1 */
T31 = matid(44);
T31[5,5] = 0;
T31[10,10] = 0;

/* a = 2 b >= 1 */
T21 = T31;
T21[4,4] = 0;
T21[9,9] = 0;

/* We load the variables x1,...,x44 and take the transform the same way as above for the cuspidal part.*/
M3=mfinit([48,3,-3]);
M4=mfinit([48,3,-4]);
B3=mfbasis(M3);
B4=mfbasis(M4);
atk3=mfatkininit(M3,48);
atk4=mfatkininit(M4,48);
W=matconcat([atk3[2]/4,0;0,atk4[2]/sqrt(12)]);
x=[0.54260880498140096513867653943544603187, -1.2415641990507267991184361399538105655, -3.3991905564732487151430228250758675037, -1.7265328337275850241386961266251401995, -0.81363822218490473489145562398280189476, -0.096105123401520619760841407689890651869, -0.27221493503052813722741181942257235873, 0.51401353564796162790999644420796668884, 0.41269645847740238263636892690450016577, 0, -0.54260880498140096513867653943544603187, 0.096105123401520619760841407689890651869, -0.61648396750108306711304367224663420226, 0.0010926169103219513034252187992627946530, 0.47243563002773065942505045637470491711, 0.039437084243136701318315982461106041564, -1.0533219944523385506109234372658713817, 0.012002879245312781223893749659513914570, 0.089262594448894739650344612278568599360, 0.11554101925893365487160361897119358347, 0.16698098274738313741381335109181353960, -0.077489051751035178094959800716149812514, -0.094330188558529879564553410275673781724, 0.060639194906309794623804168766865979012, -0.37905807314629534412480287497367879036, 1.4476556752440492028213652494914163684, 3.5885598403833256663248165887542015140, -4.8033101379755878918290347229723305568, 0.094330188558529879564553410275673781724, -0.088040011215475046862305376778881572847, -0.37905807314629534412480287497367879036, -0.013466217519949019225100408189494107099, 2.8874442794684106808734151911975535527, -6.7493834708978057073626986721978412391, -0.46901482838125511797246129933033163853, -1.3134988693595305093541618640494079449, -1.5459762991088610300938595986576085849, -0.46901482838125511797246129933033163853, -0.82398184122267052703537196006578932030, -1.0698466798734160754078802044282581803, 0, -18.027492128156189368154902560159294763, 0, -1.6689434159500542284358558005305559402];
y=W*x~;


/* computations of the bounds in Table 7.1
The cases where n is even (that is, a>0) can all be treated using equation (7.17). Depending on the factors of 2 and 3, we set certain variables to 0 if those terms should not be present. In each case, we let p be the vector x with those entries changed to 0. Likewise for the transform, q is the vector y with the same entries changed to 0. Multiplying by the matrix M changes the appropriate entries to 0.
The cases a=0, b=0 and a=0, b=1 involve extra terms from the quadratic twists, as explained in subsections 7.1 and 7.2. For b00, we take the minimum over four cases: n could be 1 mod 12, 5 mod 12, 11 mod 12, or 7 mod 12. For a00, we already know the result is 0 for n=1 mod 4 so we just take the minimum of two cases (n mod 3 could be 1 or 2). 
*/

/* Values n coprime to 48 */
a00=low3*min(x[1]-x[11]+x[6]-x[12]+uppr*(x[23]-x[29]),x[1]-x[11]-(x[6]-x[12])+uppr*(x[23]-x[29]));
b001=low3*(y[1]+y[11]+y[6]+y[12]+lowr*(y[23]+y[29]));
b005=low3*(y[1]+y[11]-y[6]-y[12]+lowr*(y[23]+y[29]));
b0011=low3*(y[1]-y[11]-y[6]+y[12]+lowr*(y[23]-y[29]));
b007=low3*(y[1]-y[11]+y[6]-y[12]+lowr*(y[23]-y[29]));
b00=vecsort([b001,b005,b0011,b007])[1];

/* Odd multiples of 3 */
a01=low3*(x[1]-x[11]-abs(x[6]-x[12])/9-uppr*0.1*(abs(9*x[23]-9*x[29]+x[25]+x[31])+abs(x[29]+x[25]+x[31])/9));
b01=low3*(-abs(y[1]+y[11])-abs(y[6]+y[12])/9+lowr*0.1*(9*y[23]+9*y[29]+y[25]-y[31]-abs(y[23]+y[29]-y[25]+y[31])/9));

/* Multiples of 2 */

a = 1;
b = 0;
M = T10;
p = M*x~;
q = M*y;
t1 = min(X3(a,b,p,1), X3(a,b,p,-1));
t2 = min(X4(a,b,p,1), X4(a,b,p,-1));
g = if(t2 > 0, lowr, uppr);
a10 = low3*(t1+g*t2)/(4^a*9^b);

t1 = min(X3(a,b,q,1), X3(a,b,q,-1));
t2 = min(X4(a,b,q,1), X4(a,b,q,-1));
g = if(t2 > 0, lowr, uppr);
b10 = low3*(t1+g*t2)/(4^a*9^b);

a = 2;
b = 0;
M = T20;
p = M*x~;
q = M*y;
t1 = min(X3(a,b,p,1), X3(a,b,p,-1));
t2 = min(X4(a,b,p,1), X4(a,b,p,-1));
g = if(t2 > 0, lowr, uppr);
a20 = low3*(t1+g*t2)/(4^a*9^b);

t1 = min(X3(a,b,q,1), X3(a,b,q,-1));
t2 = min(X4(a,b,q,1), X4(a,b,q,-1));
g = if(t2 > 0, lowr, uppr);
b20 = low3*(t1+g*t2)/(4^a*9^b);

a = 3;
b = 0;
M = T30;
p = M*x~;
q = M*y;
t1 = min(X3(a,b,p,1), X3(a,b,p,-1));
t2 = min(X4(a,b,p,1), X4(a,b,p,-1));
g = if(t2 > 0, lowr, uppr);
a30 = low3*(t1+g*t2)/(4^a*9^b);

t1 = min(X3(a,b,q,1), X3(a,b,q,-1));
t2 = min(X4(a,b,q,1), X4(a,b,q,-1));
g = if(t2 > 0, lowr, uppr);
b30 = low3*(t1+g*t2)/(4^a*9^b);

a = 2;
b = 1;
M = T21;
p = M*x~;
q = M*y;
t1 = min(X3(a,b,p,1), X3(a,b,p,-1));
t2 = min(X4(a,b,p,1), X4(a,b,p,-1));
g = if(t2 > 0, lowr, uppr);
a21 = low3*(t1+g*t2)/(4^a*9^b);

t1 = min(X3(a,b,q,1), X3(a,b,q,-1));
t2 = min(X4(a,b,q,1), X4(a,b,q,-1));
g = if(t2 > 0, lowr, uppr);
b21 = low3*(t1+g*t2)/(4^a*9^b);

a = 4;
b = 0;
p = x;
q = y~;
t1 = min(X3(a,b,p,1), X3(a,b,p,-1));
t2 = min(X4(a,b,p,1), X4(a,b,p,-1));
g = if(t2 > 0, lowr, uppr);
a40 = low3*(t1+g*t2)/(4^a*9^b);

/* a >= 4, b = 0 */
a = 4;
t1 = -1/5*(abs(s1(y)) + abs(s3(y)) + abs(s5(y))/4^a + abs(s7(y))/4^a);
t2 = 1/10*(s2(y) - abs(s4(y)) - abs(s6(y)+s8(y))/4^a);
g = if(t2 > 0, lowr, uppr);
b40 = low3*(t1+g*t2) ;
  
  
/* a >= 5, b = 0, worst case a = 5 */
a = 5;
t1 = -1/5*(abs(s1(x)) + abs(s3(x)) + abs(s5(x))/4^a + abs(s7(x))/4^a);
t2 = 1/10*(s2(x) - abs(s4(x)) - abs(s6(x)+s8(x))/4^a);
g = if(t2 > 0, lowr, uppr);
a50 = low3*(t1+g*t2) ;
 
 
/* a >= 4, b >= 1 here the worst case is when b = 1, a =4 */
a = 4;
b = 1;
t1 = -1/5*(abs(s1(x)) + abs(s3(x))/9^b + abs(s5(x))/4^a + abs(s7(x))/(9^b*4^a));
t2 = 1/10*(s2(x) - abs(s4(x))/9 - abs(s6(x))/4^a - abs(s8(x))/(9^b*4^a));
g = if(t2 > 0, lowr, uppr);
a41 = low3*(t1+g*t2) ;

t1 = -1/5*(abs(s1(y)) + abs(s3(y))/9^b + abs(s5(y))/4^a + abs(s7(y))/(9^b*4^a));
t2 = 1/10*(s2(y) - abs(s4(y))/9^b - abs(s6(y))/4^a - abs(s8(y))/(9^b*4^a));
g = if(t2 > 0, lowr, uppr);
b41 = low3*(t1+g*t2) ;

/* a = 3, b >= 1 here the worst case is when b = 1 */
a = 3;
b = 1;
M = T31;
p = M*x~;
q = M*y;
t1 = -1/5*(abs(s1(p)) + abs(s3(p))/9^b + abs(s5(p))/4^a + abs(s7(p))/(9^b*4^a));
tmp(c) = 1/10*(s2(p) + c*s6(p)/4^a - abs(s4(p) + c*s8(p)/4^a)/9^b);
t2 = min(tmp(-1),tmp(1));
g = if(t2 > 0, lowr, uppr);
a31 = low3*(t1+g*t2) ;

t1 = -1/5*(-s1(q) + abs(s3(q))/9^b + s5(q)/4^a + abs(s7(q))/(9^b*4^a));
tmp(c) = 1/10*(s2(q) + c*s6(q)/4^a - abs(s4(q) + c*s8(q)/4^a)/9^b);
t2 = min(tmp(-1),tmp(1));
g = if(t2 > 0, lowr, uppr);
b31 = low3*(t1+g*t2) ;

/* a = 1, b >= 1 here the worst case is when b = 1 */
a = 1;
b = 1;
M = T11;
p = M*x~;
q = M*y;
t1 = 1/5*(s1(p) - s5(p)/4^a  - abs(s3(p) - s7(p)/4^a)/9^b);
tmp(c) = 1/10*(s2(p) + c*s6(p)/4^a - abs(s4(p) + c*s8(p)/4^a)/9^b);
t2 = min(tmp(-1),tmp(1));
g = if(t2 > 0, lowr, uppr);
a11 = low3*(t1+g*t2) ;

p = q;
t1 = 1/5*(s1(p) - s5(p)/4^a  - abs(s3(p) - s7(p)/4^a)/9^b);
tmp(c) = 1/10*(s2(p) + c*s6(p)/4^a - abs(s4(p) + c*s8(p)/4^a)/9^b);
t2 = min(tmp(-1),tmp(1));
g = if(t2 > 0, lowr, uppr);
b11 = low3*(t1+g*t2) ;

/* a = 2, b >= 2 here the worst case is when b = 1 */
a = 2;
b = 2;
M = T21;
p = M*x~;
q = M*y;
t1 = 1/5*(s1(p) + s5(p)/4^a - abs(s3(p) + s7(p)/4^a)/9^b);
tmp(c) = 1/10*(s2(p) + c*s6(p)/4^a - abs(s4(p) + c*s8(p)/4^a)/9^b);
t2 = min(tmp(-1),tmp(1));
g = if(t2 > 0, lowr, uppr);
a22 = low3*(t1+g*t2) ;

p = q;
t1 = 1/5*(s1(p) + s5(p)/4^a - abs(s3(p) + s7(p)/4^a)/9^b);
tmp(c) = 1/10*(s2(p) + c*s6(p)/4^a - abs(s4(p) + c*s8(p)/4^a)/9^b);
t2 = min(tmp(-1),tmp(1));
g = if(t2 > 0, lowr, uppr);
b22 = low3*(t1+g*t2) ;

/* These numbers can be collected into Table 7.1 (except for the top row, which is treated separately in subsection 7.1 "Values n coprime to 48".)*/
table71=["b","a","epsilon","epsilon*";0,0,a00,b00;0,1,a10,b10;0,2,a20,b20;0,3,a30,b30;0,">= 4",min(a50,a40),b40;">= 1",0,a01,b01;">=1",1,a11,b11;">=1",2,a21,b21;">=1",3,a31,b31;">=1",">=4",a41,b41];
table71

/* Output: table71 should produce a higher-precision version of Table 7.1. These values are rounded down to the lower bounds stated in the paper.
[   "b"    "a"                                   "epsilon"                                 "epsilon*"]

[     0      0    0.64688766410471365715772609276335246519   0.16123767062435152951060983573993186125]

[     0      1  0.0037425552420029692389021159466196536326 0.0013580156001057262959747897746256134386]

[     0      2 0.00081346866993550601194991335944608850908  0.039028397743033393261930498049722824180]

[     0      3 0.00082643623125661282441846131330347052273   0.20381834858746950339881901391610645722]

[     0 ">= 4" 0.00026493542037053511895745232954803905552  0.011815972622709339403590807107693349445]

[">= 1"      0    0.74145799248145997864723275290709429132   0.27574721654059712010349701562523516988]

[ ">=1"      1  0.0067586286271784969513931046423677237831  0.054621211926853610953527760185745409313]

[ ">=1"      2 8.7536485524854605206408567012412760857 E-6  0.013633813680670959611910598565623054005]

[ ">=1"      3   0.015243808557392452178125135003082696944   0.13779452980908966998835824464357377433]

[ ">=1"  ">=4"   0.016353057226842559991824918631237739911   0.10085805690469353803871886699274511855]
*/

/* Reducing to a finite check
Let C be as in the upper bound for the cuspidal part (Proposition 5.2), and eps as in the lower bound for the Eisenstein part (Proposition 7.1).
Let R be the exponential rate from Nicolas-Robin, as described in Section 2.
"crossing" is the numerical solution of equation (2.14).
We round it to bound1, so that a_n >= 0 for n > bound1 by the estimates from Section 2.
*/

x=[0.54260880498140096513867653943544603187, -1.2415641990507267991184361399538105655, -3.3991905564732487151430228250758675037, -1.7265328337275850241386961266251401995, -0.81363822218490473489145562398280189476, -0.096105123401520619760841407689890651869, -0.27221493503052813722741181942257235873, 0.51401353564796162790999644420796668884, 0.41269645847740238263636892690450016577, 0, -0.54260880498140096513867653943544603187, 0.096105123401520619760841407689890651869, -0.61648396750108306711304367224663420226, 0.0010926169103219513034252187992627946530, 0.47243563002773065942505045637470491711, 0.039437084243136701318315982461106041564, -1.0533219944523385506109234372658713817, 0.012002879245312781223893749659513914570, 0.089262594448894739650344612278568599360, 0.11554101925893365487160361897119358347, 0.16698098274738313741381335109181353960, -0.077489051751035178094959800716149812514, -0.094330188558529879564553410275673781724, 0.060639194906309794623804168766865979012, -0.37905807314629534412480287497367879036, 1.4476556752440492028213652494914163684, 3.5885598403833256663248165887542015140, -4.8033101379755878918290347229723305568, 0.094330188558529879564553410275673781724, -0.088040011215475046862305376778881572847, -0.37905807314629534412480287497367879036, -0.013466217519949019225100408189494107099, 2.8874442794684106808734151911975535527, -6.7493834708978057073626986721978412391, -0.46901482838125511797246129933033163853, -1.3134988693595305093541618640494079449, -1.5459762991088610300938595986576085849, -0.46901482838125511797246129933033163853, -0.82398184122267052703537196006578932030, -1.0698466798734160754078802044282581803, 0, -18.027492128156189368154902560159294763, 0, -1.6689434159500542284358558005305559402];
u=vector(44,i,abs(x[i]));
c13=1;
c16=2;
c18=6;
c20=24/7;
c21=24/7+2*12*(2+sqrt(2))/7;
c22=8*sqrt(2);
c35=2;
c38=4;
c41=1;
c43=2;
c44=2*sqrt(3);
C=(u[13]+u[14]/2+u[15]/4)*c13+(u[16]+u[17]/2)*c16+(u[18]+u[19]/2)*c18+(u[35]+u[36]/2+u[37]/4)*c35+(u[38]+u[39]/2+u[40]/4)*c38+(u[41]+u[42]/3)*c41+u[43]*c43+u[44]*c44;
low3=0.9429085;
lowr=0.94999;
uppr=1.096954;
s1(v)=9*(13*v[1]+3*v[2]+1*v[3]);
s3(v)=13*v[6]-3*v[7]+v[8];
s2(v)=8*(16*v[23]+4*v[24]+1*v[26])+16*v[25]+4*v[27]+v[28];
s4(v)=-8*(v[29]+v[30]+v[32])+v[31]+v[33]+v[34];
eps=low3*(s1(x)-s3(x)+lowr*(s2(x)+s4(x)))/144;
R=log(sigma(6983776800,0))*log(log(6983776800))/log(6983776800);
f(n)=R/log(log(n))+log(C/eps)/log(n)-1;
crossing=solve(x=10^9,10^10,f(x));
bound1=floor(crossing);
["It is enough to check a_n until...", bound1]
M3=mfinit([48,3,-3]);
M4=mfinit([48,3,-4]);
B3=mfbasis(M3);
B4=mfbasis(M4);
atk3=mfatkininit(M3,48);
atk4=mfatkininit(M4,48);
W=matconcat([atk3[2]/4,0;0,atk4[2]/sqrt(12)]);
y=W*x~;

c13=1;

c16=2;

c18=6;

c20=24/7;

c21=24/7+2*12*(2+sqrt(2))/7;

c22=8*sqrt(2);

c35=2;

c38=4;

c41=1;

c43=2;

c44=2*sqrt(3);

v=vector(44,i,abs(y[i]));

C=(v[13]+v[14]/2+v[15]/4)*c13+(v[16]+v[17]/2)*c16+(v[18]+v[19]/2)*c18+(v[35]+v[36]/2+v[37]/4)*c35+(v[38]+v[39]/2+v[40]/4)*c38+(v[41]+v[42]/3)*c41+v[43]*c43+v[44]*c44;

low3=0.9429085;

lowr=0.94999;

uppr=1.096954;

eps=((4*y[23]+y[24]+y[29]+y[30])*lowr+3*y[1]+y[2]+3*y[6]-y[7])*low3/4;
R=log(sigma(6983776800,0))*log(log(6983776800))/log(6983776800);

f(n)=R/log(log(n))+log(C/eps)/log(n)-1;

crossing=solve(x=10^6,10^7,f(x));

bound1=floor(crossing);

["It is enough to check b_n until...",bound1]

/* Making a to-do list
We make an empty list "todo", and loop over n up to bound1. This may take several days of computer time. We ignore the progression n=1 mod 4 where we know a_n=0 already. For the other n, we compare the Eisenstein part to Deligne's bound. If the total a_{n,eis}-C*n*sigma(n,0) is negative, we add n to the to-do list. For these, we will have to compute the cuspidal part exactly instead of using Deligne's theorem.
In practice, instead of "listinsert", it may be more efficient to write the entries to a separate file. The result should be the same as "checklist-an-67250.txt".
x=[0.54260880498140096513867653943544603187, -1.2415641990507267991184361399538105655, -3.3991905564732487151430228250758675037, -1.7265328337275850241386961266251401995, -0.81363822218490473489145562398280189476, -0.096105123401520619760841407689890651869, -0.27221493503052813722741181942257235873, 0.51401353564796162790999644420796668884, 0.41269645847740238263636892690450016577, 0, -0.54260880498140096513867653943544603187, 0.096105123401520619760841407689890651869, -0.61648396750108306711304367224663420226, 0.0010926169103219513034252187992627946530, 0.47243563002773065942505045637470491711, 0.039437084243136701318315982461106041564, -1.0533219944523385506109234372658713817, 0.012002879245312781223893749659513914570, 0.089262594448894739650344612278568599360, 0.11554101925893365487160361897119358347, 0.16698098274738313741381335109181353960, -0.077489051751035178094959800716149812514, -0.094330188558529879564553410275673781724, 0.060639194906309794623804168766865979012, -0.37905807314629534412480287497367879036, 1.4476556752440492028213652494914163684, 3.5885598403833256663248165887542015140, -4.8033101379755878918290347229723305568, 0.094330188558529879564553410275673781724, -0.088040011215475046862305376778881572847, -0.37905807314629534412480287497367879036, -0.013466217519949019225100408189494107099, 2.8874442794684106808734151911975535527, -6.7493834708978057073626986721978412391, -0.46901482838125511797246129933033163853, -1.3134988693595305093541618640494079449, -1.5459762991088610300938595986576085849, -0.46901482838125511797246129933033163853, -0.82398184122267052703537196006578932030, -1.0698466798734160754078802044282581803, 0, -18.027492128156189368154902560159294763, 0, -1.6689434159500542284358558005305559402];
u=vector(44,i,abs(x[i]));
c13=1;
c16=2;
c18=6;
c20=24/7;
c21=24/7+2*12*(2+sqrt(2))/7;
c22=8*sqrt(2);
c35=2;
c38=4;
c41=1;
c43=2;
c44=2*sqrt(3);
C=(u[13]+u[14]/2+u[15]/4)*c13+(u[16]+u[17]/2)*c16+(u[18]+u[19]/2)*c18+(u[35]+u[36]/2+u[37]/4)*c35+(u[38]+u[39]/2+u[40]/4)*c38+(u[41]+u[42]/3)*c41+u[43]*c43+u[44]*c44;
low3=0.9429085;
lowr=0.94999;
uppr=1.096954;
s1(v)=9*(13*v[1]+3*v[2]+1*v[3]);
s3(v)=13*v[6]-3*v[7]+v[8];
s2(v)=8*(16*v[23]+4*v[24]+1*v[26])+16*v[25]+4*v[27]+v[28];
s4(v)=-8*(v[29]+v[30]+v[32])+v[31]+v[33]+v[34];
eps=low3*(s1(x)-s3(x)+lowr*(s2(x)+s4(x)))/144;
R=log(sigma(6983776800,0))*log(log(6983776800))/log(6983776800);
f(n)=R/log(log(n))+log(C/eps)/log(n)-1;
crossing=solve(x=10^9,10^10,f(x));
bound1=floor(crossing);
allocatemem(16384000000)
x3=vector(22,i,x[i]);
x4=vector(22,i,x[22+i]);
x3eis=vector(12,i,x[i]);
x4eis=vector(12,i,x[22+i]);
E3=mfinit([48,3,-3],{3});
E4=mfinit([48,3,-4],{3});
e3=mfbasis(E3);
e4=mfbasis(E4);
eis3=mflinear(e3,x3eis);
eis4=mflinear(e4,x4eis);
todo=List();
for(n=1,bound1,if(!(Mod(n,4)==1),if(mfcoef(eis3,n)+mfcoef(eis4,n)-C*n*sigma(n,0)<0, listinsert(todo,n,1))));
*/

/*
A smaller calculation is enough to verify b_n >= 0. The list can be stored in working memory instead of a separate file, and the calculation should take hours rather than days. We can apply Pari's built-in command "mfcoefs" directly to the basis functions, instead of a formula like (5.4), implemented as c(n) below.
We let y=W*x as above, and split the coefficients into four parts for cuspidal/Eisenstein subspaces of the two characters.
To study b_n instead of a_n, we take the values C* and epsilon* from the paper instead of C and epsilon (though we still call them C and eps below). We use them to compute a smaller value of "bbound1" than for a_n.
We check for n until "bbound1" whether positivity can already be concluded by comparing the Eisenstein part to Deligne's bound. If not, we add n to the to-do list "todo". In this case, the list is small enough that it can be stored in working memory. There is no need to write it to a separate file, although that would use less memory.
Finally, for n in the to-do list, we check whether b_n < 0. In case of numerical imprecision, we checked a less demanding condition b_n < 10^(-16) instead.
The output is written as pairs [n,b_n]. It shows that b_n > 0, except for the zeros we forced in (4.3) at n = 1, 2, 3, 4, 7, 8, 9, 10, 13, 14, 36, 82
*/

allocatemem(16384000000)
M3=mfinit([48,3,-3]);
M4=mfinit([48,3,-4]);
B3=mfbasis(M3);
B4=mfbasis(M4);
atk3=mfatkininit(M3,48);
atk4=mfatkininit(M4,48);
W=matconcat([atk3[2]/4,0;0,atk4[2]/sqrt(12)]);
x=[0.54260880498140096513867653943544603187, -1.2415641990507267991184361399538105655, -3.3991905564732487151430228250758675037, -1.7265328337275850241386961266251401995, -0.81363822218490473489145562398280189476, -0.096105123401520619760841407689890651869, -0.27221493503052813722741181942257235873, 0.51401353564796162790999644420796668884, 0.41269645847740238263636892690450016577, 0, -0.54260880498140096513867653943544603187, 0.096105123401520619760841407689890651869, -0.61648396750108306711304367224663420226, 0.0010926169103219513034252187992627946530, 0.47243563002773065942505045637470491711, 0.039437084243136701318315982461106041564, -1.0533219944523385506109234372658713817, 0.012002879245312781223893749659513914570, 0.089262594448894739650344612278568599360, 0.11554101925893365487160361897119358347, 0.16698098274738313741381335109181353960, -0.077489051751035178094959800716149812514, -0.094330188558529879564553410275673781724, 0.060639194906309794623804168766865979012, -0.37905807314629534412480287497367879036, 1.4476556752440492028213652494914163684, 3.5885598403833256663248165887542015140, -4.8033101379755878918290347229723305568, 0.094330188558529879564553410275673781724, -0.088040011215475046862305376778881572847, -0.37905807314629534412480287497367879036, -0.013466217519949019225100408189494107099, 2.8874442794684106808734151911975535527, -6.7493834708978057073626986721978412391, -0.46901482838125511797246129933033163853, -1.3134988693595305093541618640494079449, -1.5459762991088610300938595986576085849, -0.46901482838125511797246129933033163853, -0.82398184122267052703537196006578932030, -1.0698466798734160754078802044282581803, 0, -18.027492128156189368154902560159294763, 0, -1.6689434159500542284358558005305559402];
y=W*x~;

E3=mfinit([48,3,-3],{3});
E4=mfinit([48,3,-4],{3});
S3=mfinit([48,3,-3],{1});
S4=mfinit([48,3,-4],{1});
e3=mfbasis(E3);
e4=mfbasis(E4);
C3=mfbasis(S3);
C4=mfbasis(S4);
y3eis=vector(12,i,y[i]);

y3cusp=vector(10,i,y[12+i]);

y4eis=vector(12,i,y[22+i]);

y4cusp=vector(10,i,y[34+i]);

yeis3=mflinear(e3,y3eis);

yeis4=mflinear(e4,y4eis);

cusp3=mflinear(C3,y3cusp);

cusp4=mflinear(C4,y4cusp);

c13=1;

c16=2;

c18=6;

c20=24/7;

c21=24/7+2*12*(2+sqrt(2))/7;

c22=8*sqrt(2);

c35=2;

c38=4;

c41=1;

c43=2;

c44=2*sqrt(3);

v=vector(44,i,abs(y[i]));

C=(v[13]+v[14]/2+v[15]/4)*c13+(v[16]+v[17]/2)*c16+(v[18]+v[19]/2)*c18+(v[35]+v[36]/2+v[37]/4)*c35+(v[38]+v[39]/2+v[40]/4)*c38+(v[41]+v[42]/3)*c41+v[43]*c43+v[44]*c44;

low3=0.9429085;

lowr=0.94999;

uppr=1.096954;

eps=((4*y[23]+y[24]+y[29]+y[30])*lowr+3*y[1]+y[2]+3*y[6]-y[7])*low3/4;
R=log(sigma(6983776800,0))*log(log(6983776800))/log(6983776800);

f(n)=R/log(log(n))+log(C/eps)/log(n)-1;

crossing=solve(x=10^6,10^7,f(x));

bbound1=floor(crossing);

["It is enough to check until...",bbound1]
["Next compute the Eisenstein part for b_n until that point."]
["The cuspidal part only needs to be computed when the Eisenstein part is smaller than Deligne's bound. We put those values on a to-do list."]
["It may take around half an hour to compute the Eisenstein part and make the list. Directly computing the cuspidal part up to the end of the list may take an hour. A quick loop through the list then shows that b_n is positive ecept for the values forced to be 0 in equation (4.3)."]
["A faster alternative is to use the Hecke basis to compute b_{n,cusp} for n on the list."]
yeis3coefs=mfcoefs(yeis3,bbound1);

yeis4coefs=mfcoefs(yeis4,bbound1);

bneis=yeis3coefs+yeis4coefs;

todo=[];

for(n=1,bbound1,if(bneis[n+1]-C*n*sigma(n,0)<0,todo=concat(todo,[n])));



bbound2=todo[length(todo)]+1;
bncusp=mfcoefs(cusp3,bbound2)+mfcoefs(cusp4,bbound2);

for(i=1,length(todo),if(bncusp[todo[i]+1]+bneis[todo[i]+1]<1E-16,print([todo[i],bncusp[todo[i]+1]+bneis[todo[i]+1]])))
["This calculation shows that all b_n are non-negative. The calculation for a_n may take several days. See checklist-an-67250.txt for the list of values where the cuspidal part needs to be computed, and the rest of supplement.txt for explanations of how to generate and check the list."]

/* Output: we find the same list (4.3) containing the n for which b_n=0. Here some of the b_n appear to be small negative numbers because we used a 38-digit approximation to the exact coefficients x1,...,x44. We conclude that b_n >= 0 for all n.


[1, 1.6456920911512025112 E-37]
[2, -1.3812058622161878219 E-37]
[3, -9.403954806578300064 E-38]
[4, -2.586087571809032518 E-37]
[7, -1.1284745767893960076 E-36]
[8, -4.231779662960235029 E-37]
[9, -3.761581922631320025 E-37]
[10, -1.6985893369382054490 E-36]
[13, -7.523163845262640051 E-37]
[14, -3.667542374565537025 E-36]
[36, -2.3133728824182618157 E-35]
[82, -1.2761166672526753187 E-34]
*/


/* The final check.
First allocate memory. 16384 Mbytes was enough for us (and the largest allowed on a 64-bit laptop).
Take x as above. We separate into Eisenstein/cuspidal parts for both characters to get x3, x4, x3eis, x4eis, x3cusp, x4cusp.
E3, E4 are the spaces of Eisenstein series, with basis e3 and e4 respectively.
S3, S4 are the cuspidal subspaces, with bases C3 and C4.
eis3,eis4,cusp3,cusp4 are the components of the candidate g in each of the four subspaces.
h1,h2,h7,h9 are the Hecke eigenforms from Section 5, page 20.
L is the list of indices to check, read from a file "checklist-an-67250.txt"
c(n) is a_{n,cusp} computed according to (5.4). Note that this formula is only valid for even values of n. For odd n, there would be a further contribution from the quadratic twists f_{20},f_{21},f_{22}.
We check whether a_{n,eis}+a_{n,cusp}<0.001, in particular checking whether it is negative, for n not congruent to 1 mod 4 (those coefficients have been forced to be 0 already). The verification took 42 hours and 43 minutes on our device. The resulting values are just the ones from (4.1),(4.2),(4.3) where we forced zeros. It follows that a_n >= 0 for all n.
*/

allocatemem(16384000000)
x=[0.54260880498140096513867653943544603187, -1.2415641990507267991184361399538105655, -3.3991905564732487151430228250758675037, -1.7265328337275850241386961266251401995, -0.81363822218490473489145562398280189476, -0.096105123401520619760841407689890651869, -0.27221493503052813722741181942257235873, 0.51401353564796162790999644420796668884, 0.41269645847740238263636892690450016577, 0, -0.54260880498140096513867653943544603187, 0.096105123401520619760841407689890651869, -0.61648396750108306711304367224663420226, 0.0010926169103219513034252187992627946530, 0.47243563002773065942505045637470491711, 0.039437084243136701318315982461106041564, -1.0533219944523385506109234372658713817, 0.012002879245312781223893749659513914570, 0.089262594448894739650344612278568599360, 0.11554101925893365487160361897119358347, 0.16698098274738313741381335109181353960, -0.077489051751035178094959800716149812514, -0.094330188558529879564553410275673781724, 0.060639194906309794623804168766865979012, -0.37905807314629534412480287497367879036, 1.4476556752440492028213652494914163684, 3.5885598403833256663248165887542015140, -4.8033101379755878918290347229723305568, 0.094330188558529879564553410275673781724, -0.088040011215475046862305376778881572847, -0.37905807314629534412480287497367879036, -0.013466217519949019225100408189494107099, 2.8874442794684106808734151911975535527, -6.7493834708978057073626986721978412391, -0.46901482838125511797246129933033163853, -1.3134988693595305093541618640494079449, -1.5459762991088610300938595986576085849, -0.46901482838125511797246129933033163853, -0.82398184122267052703537196006578932030, -1.0698466798734160754078802044282581803, 0, -18.027492128156189368154902560159294763, 0, -1.6689434159500542284358558005305559402];
x3=vector(22,i,x[i]);
x4=vector(22,i,x[22+i]);
x3eis=vector(12,i,x[i]);
x4eis=vector(12,i,x[22+i]);
x3cusp=vector(10,i,x[12+i]);
x4cusp=vector(10,i,x[34+i]);
E3=mfinit([48,3,-3],{3});
E4=mfinit([48,3,-4],{3});
S3=mfinit([48,3,-3],{1});
S4=mfinit([48,3,-4],{1});
e3=mfbasis(E3);
e4=mfbasis(E4);
C3=mfbasis(S3);
C4=mfbasis(S4);
eis3=mflinear(e3,x3eis);
eis4=mflinear(e4,x4eis);
cusp3=mflinear(C3,x3cusp);
cusp4=mflinear(C4,x4cusp);
h1=C3[1];
h2=mflinear(C3,[0,0,0,(-1+sqrt(8)*I)/4/sqrt(2)/I,0,1/4/sqrt(2)/I,0,0,0,0]);
h7=mflinear(C4,[(1+sqrt(3)*I)/sqrt(12)/I,0,0,1/sqrt(12)/I,0,0,0,0,0,0]);
h9=mflinear(C4,[0,0,0,0,0,0,0,0,1/2,-1/sqrt(12)/I]);
L=read("checklist-an-67250.txt"); /* for Windows operating system, change this to the path specifying where "checklist-an-67250.txt" is saved in your directory, as in L=read("C:/Users/.../checklist-an-67250.txt")*/
c(n)=
{
	a=valuation(n,2);
	b=valuation(n,3);
	n0=n/2^a/3^b;
	f=factor(n0);
	l=matsize(f)[1];
	H1=prod(j=1,l,mfcoef(h1,f[j,1]^f[j,2]));
	H2=prod(j=1,l,mfcoef(h2,f[j,1]^f[j,2]));
	H7=prod(j=1,l,mfcoef(h7,f[j,1]^f[j,2]));
	H9=prod(j=1,l,mfcoef(h9,f[j,1]^f[j,2]));
	part1=(x[13]*mfcoef(h1,2^a)+(a>=1)*x[14]*mfcoef(h1,(a>=1)*2^(a-1))+(a>=2)*x[15]*mfcoef(h1,(a>=2)*2^(a-2)))*mfcoef(h1,3^b)*H1;
	part2=2*real( ( x[16]*mfcoef(h2,2^a)+x[18]*(1+sqrt(8)*I)*mfcoef(h2,2^a)+(a>=1)*(x[17]+x[19]*(1+sqrt(8)*I))*mfcoef(h2,(a>=1)*2^(a-1)))*mfcoef(h2,3^b)*H2);
	part7=2*real((x[35]*mfcoef(h7,2^a)+x[36]*(a>=1)*mfcoef(h7,(a>=1)*2^(a-1))+x[37]*(a>=2)*mfcoef(h7,(a>=2)*2^(a-2))+(-1+sqrt(3)*I)*(x[38]*mfcoef(h7,2^a)+x[39]*(a>=1)*mfcoef(h7,(a>=1)*2^(a-1))+x[40]*(a>=2)*mfcoef(h7,(a>=2)*2^(a-2))))*mfcoef(h7,3^b)*H7);
	part1+part2+part7;
}
for(i=1,length(L),if(mfcoef(eis3,L[i])+mfcoef(eis4,L[i])+c(L[i])<0.001,print([L[i],mfcoef(eis3,L[i])+mfcoef(eis4,L[i])+c(L[i])])))

/*
Output: one should find the following approximate zeros, corresponding to (4.2). The indices n=3 and n=11 from (4.2) do not appear because this calculation only checks even n.
[2, -9.110081218872728187 E-38]
[4, 2.938735877055718770 E-38]
[6, -3.056285312137947521 E-37]
[8, -6.964804028622053485 E-37]
[10, -9.521504241660528815 E-37]
[12, -1.2225141248551790083 E-36]
[22, -5.689392657979871539 E-36]
[26, -7.617203393328423052 E-36]
[32, -8.698658196084927559 E-36]
[38, -1.6927118651840940115 E-35]
[60, -2.0077443512044670636 E-35]
[64, -4.137740114894452028 E-35]
[88, -7.104687856369905698 E-35]
[90, -8.275480229788904056 E-35]
[92, -7.222237291452134449 E-35]
[106, -1.1134282490988707275 E-34]
[164, -2.535306215853509697 E-34]
[1932, -2.723911933454244083 E-32]
*/

/* The last odd index on the list is n=315. We already know from setting up the linear program that the coefficients have the correct sign up to that point, but we can check again as follows. We only need to check n congruent to 3 mod 4, hence n=4*i+3 with i at most 78.
*/

x=[0.54260880498140096513867653943544603187, -1.2415641990507267991184361399538105655, -3.3991905564732487151430228250758675037, -1.7265328337275850241386961266251401995, -0.81363822218490473489145562398280189476, -0.096105123401520619760841407689890651869, -0.27221493503052813722741181942257235873, 0.51401353564796162790999644420796668884, 0.41269645847740238263636892690450016577, 0, -0.54260880498140096513867653943544603187, 0.096105123401520619760841407689890651869, -0.61648396750108306711304367224663420226, 0.0010926169103219513034252187992627946530, 0.47243563002773065942505045637470491711, 0.039437084243136701318315982461106041564, -1.0533219944523385506109234372658713817, 0.012002879245312781223893749659513914570, 0.089262594448894739650344612278568599360, 0.11554101925893365487160361897119358347, 0.16698098274738313741381335109181353960, -0.077489051751035178094959800716149812514, -0.094330188558529879564553410275673781724, 0.060639194906309794623804168766865979012, -0.37905807314629534412480287497367879036, 1.4476556752440492028213652494914163684, 3.5885598403833256663248165887542015140, -4.8033101379755878918290347229723305568, 0.094330188558529879564553410275673781724, -0.088040011215475046862305376778881572847, -0.37905807314629534412480287497367879036, -0.013466217519949019225100408189494107099, 2.8874442794684106808734151911975535527, -6.7493834708978057073626986721978412391, -0.46901482838125511797246129933033163853, -1.3134988693595305093541618640494079449, -1.5459762991088610300938595986576085849, -0.46901482838125511797246129933033163853, -0.82398184122267052703537196006578932030, -1.0698466798734160754078802044282581803, 0, -18.027492128156189368154902560159294763, 0, -1.6689434159500542284358558005305559402];
x3=vector(22,i,x[i]);
x4=vector(22,i,x[22+i]);
x3eis=vector(12,i,x[i]);
x4eis=vector(12,i,x[22+i]);
x3cusp=vector(10,i,x[12+i]);
x4cusp=vector(10,i,x[34+i]);
E3=mfinit([48,3,-3],{3});
E4=mfinit([48,3,-4],{3});
S3=mfinit([48,3,-3],{1});
S4=mfinit([48,3,-4],{1});
e3=mfbasis(E3);
e4=mfbasis(E4);
C3=mfbasis(S3);
C4=mfbasis(S4);
eis3=mflinear(e3,x3eis);
eis4=mflinear(e4,x4eis);
cusp3=mflinear(C3,x3cusp);
cusp4=mflinear(C4,x4cusp);
for(i=0,78,if(mfcoef(eis3,4*i+3)+mfcoef(eis4,4*i+3)+mfcoef(cusp3,4*i+3)+mfcoef(cusp4,4*i+3)<0.001,print([4*i+3,mfcoef(eis3,4*i+3)+mfcoef(eis4,4*i+3)+mfcoef(cusp3,4*i+3)+mfcoef(cusp4,4*i+3)])))

*/
Output: the result should be the remaining values n=3 and n=11 from (4.2).
[3, 6.582768364604810045 E-37]
[11, -3.761581922631320025 E-36]
*/


/* Formula (5.4) is only valid for even n. The following formula works for all n, but part9 and part11 are 0 in case n is even.
x=[0.54260880498140096513867653943544603187, -1.2415641990507267991184361399538105655, -3.3991905564732487151430228250758675037, -1.7265328337275850241386961266251401995, -0.81363822218490473489145562398280189476, -0.096105123401520619760841407689890651869, -0.27221493503052813722741181942257235873, 0.51401353564796162790999644420796668884, 0.41269645847740238263636892690450016577, 0, -0.54260880498140096513867653943544603187, 0.096105123401520619760841407689890651869, -0.61648396750108306711304367224663420226, 0.0010926169103219513034252187992627946530, 0.47243563002773065942505045637470491711, 0.039437084243136701318315982461106041564, -1.0533219944523385506109234372658713817, 0.012002879245312781223893749659513914570, 0.089262594448894739650344612278568599360, 0.11554101925893365487160361897119358347, 0.16698098274738313741381335109181353960, -0.077489051751035178094959800716149812514, -0.094330188558529879564553410275673781724, 0.060639194906309794623804168766865979012, -0.37905807314629534412480287497367879036, 1.4476556752440492028213652494914163684, 3.5885598403833256663248165887542015140, -4.8033101379755878918290347229723305568, 0.094330188558529879564553410275673781724, -0.088040011215475046862305376778881572847, -0.37905807314629534412480287497367879036, -0.013466217519949019225100408189494107099, 2.8874442794684106808734151911975535527, -6.7493834708978057073626986721978412391, -0.46901482838125511797246129933033163853, -1.3134988693595305093541618640494079449, -1.5459762991088610300938595986576085849, -0.46901482838125511797246129933033163853, -0.82398184122267052703537196006578932030, -1.0698466798734160754078802044282581803, 0, -18.027492128156189368154902560159294763, 0, -1.6689434159500542284358558005305559402];
S3=mfinit([48,3,-3],{1});
S4=mfinit([48,3,-4],{1});
C3=mfbasis(S3);
C4=mfbasis(S4);
h1=C3[1];
h2=mflinear(C3,[0,0,0,(-1+sqrt(8)*I)/4/sqrt(2)/I,0,1/4/sqrt(2)/I,0,0,0,0]);
h7=mflinear(C4,[(1+sqrt(3)*I)/sqrt(12)/I,0,0,1/sqrt(12)/I,0,0,0,0,0,0]);
h9=mflinear(C4,[0,0,0,0,0,0,0,0,1/2,-1/sqrt(12)/I]);
f41=C4[7];
ac(n)=
{
	a=valuation(n,2);
	b=valuation(n,3);
	k=kronecker(-4,n);
	n0=n/2^a/3^b;
	f=factor(n0);
	l=matsize(f)[1];
	H1=prod(j=1,l,mfcoef(h1,f[j,1]^f[j,2]));
	H2=prod(j=1,l,mfcoef(h2,f[j,1]^f[j,2]));
	H7=prod(j=1,l,mfcoef(h7,f[j,1]^f[j,2]));
	H9=prod(j=1,l,mfcoef(h9,f[j,1]^f[j,2]));
	H11=prod(j=1,l,mfcoef(f41,f[j,1]^f[j,2]));
	part1=((x[13]+k*(x[20]+3*x[21]))*mfcoef(h1,2^a)+(a>=1)*x[14]*mfcoef(h1,(a>=1)*2^(a-1))+(a>=2)*x[15]*mfcoef(h1,(a>=2)*2^(a-2)))*mfcoef(h1,3^b)*H1;
	part2=2*real( ( (x[16]+k*(x[20]+2*x[22]))*mfcoef(h2,2^a)+(x[18]+k*(-x[21]-2*x[22]))*(1+sqrt(8)*I)*mfcoef(h2,2^a)+(a>=1)*(x[17]+x[19]*(1+sqrt(8)*I))*mfcoef(h2,(a>=1)*2^(a-1)))*mfcoef(h2,3^b)*H2);
	part7=2*real((x[35]*mfcoef(h7,2^a)+x[36]*(a>=1)*mfcoef(h7,(a>=1)*2^(a-1))+x[37]*(a>=2)*mfcoef(h7,(a>=2)*2^(a-2))+(-1+sqrt(3)*I)*(x[38]*mfcoef(h7,2^a)+x[39]*(a>=1)*mfcoef(h7,(a>=1)*2^(a-1))+x[40]*(a>=2)*mfcoef(h7,(a>=2)*2^(a-2))))*mfcoef(h7,3^b)*H7);
	part9=2*real(x[43]*mfcoef(h9,2^a)*mfcoef(h9,3^b)*H9)+2*imag(x[44]*sqrt(3)*mfcoef(h9,2^a)*mfcoef(h9,3^b)*H9);
	part11=(x[41]*mfcoef(f41,3^b)+x[42]*(b>=1)*mfcoef(f41,(b>=1)*3^(b-1)))*mfcoef(f41,2^a)*H11;
	part1+part2+part7+part9+part11;
}
*/

/* Changing x to y gives a similar formula for the cuspidal part of the transform.
M3=mfinit([48,3,-3]);
M4=mfinit([48,3,-4]);
B3=mfbasis(M3);
B4=mfbasis(M4);
atk3=mfatkininit(M3,48);
atk4=mfatkininit(M4,48);
W=matconcat([atk3[2]/4,0;0,atk4[2]/sqrt(12)]);
x=[0.54260880498140096513867653943544603187, -1.2415641990507267991184361399538105655, -3.3991905564732487151430228250758675037, -1.7265328337275850241386961266251401995, -0.81363822218490473489145562398280189476, -0.096105123401520619760841407689890651869, -0.27221493503052813722741181942257235873, 0.51401353564796162790999644420796668884, 0.41269645847740238263636892690450016577, 0, -0.54260880498140096513867653943544603187, 0.096105123401520619760841407689890651869, -0.61648396750108306711304367224663420226, 0.0010926169103219513034252187992627946530, 0.47243563002773065942505045637470491711, 0.039437084243136701318315982461106041564, -1.0533219944523385506109234372658713817, 0.012002879245312781223893749659513914570, 0.089262594448894739650344612278568599360, 0.11554101925893365487160361897119358347, 0.16698098274738313741381335109181353960, -0.077489051751035178094959800716149812514, -0.094330188558529879564553410275673781724, 0.060639194906309794623804168766865979012, -0.37905807314629534412480287497367879036, 1.4476556752440492028213652494914163684, 3.5885598403833256663248165887542015140, -4.8033101379755878918290347229723305568, 0.094330188558529879564553410275673781724, -0.088040011215475046862305376778881572847, -0.37905807314629534412480287497367879036, -0.013466217519949019225100408189494107099, 2.8874442794684106808734151911975535527, -6.7493834708978057073626986721978412391, -0.46901482838125511797246129933033163853, -1.3134988693595305093541618640494079449, -1.5459762991088610300938595986576085849, -0.46901482838125511797246129933033163853, -0.82398184122267052703537196006578932030, -1.0698466798734160754078802044282581803, 0, -18.027492128156189368154902560159294763, 0, -1.6689434159500542284358558005305559402];
y=W*x~;
S3=mfinit([48,3,-3],{1});
S4=mfinit([48,3,-4],{1});
C3=mfbasis(S3);
C4=mfbasis(S4);
h1=C3[1];
h2=mflinear(C3,[0,0,0,(-1+sqrt(8)*I)/4/sqrt(2)/I,0,1/4/sqrt(2)/I,0,0,0,0]);
h7=mflinear(C4,[(1+sqrt(3)*I)/sqrt(12)/I,0,0,1/sqrt(12)/I,0,0,0,0,0,0]);
h9=mflinear(C4,[0,0,0,0,0,0,0,0,1/2,-1/sqrt(12)/I]);
f41=C4[7];
bc(n)=
{
	a=valuation(n,2);
	b=valuation(n,3);
	k=kronecker(-4,n);
	n0=n/2^a/3^b;
	f=factor(n0);
	l=matsize(f)[1];
	H1=prod(j=1,l,mfcoef(h1,f[j,1]^f[j,2]));
	H2=prod(j=1,l,mfcoef(h2,f[j,1]^f[j,2]));
	H7=prod(j=1,l,mfcoef(h7,f[j,1]^f[j,2]));
	H9=prod(j=1,l,mfcoef(h9,f[j,1]^f[j,2]));
	H11=prod(j=1,l,mfcoef(f41,f[j,1]^f[j,2]));
	part1=((y[13]+k*(y[20]+3*y[21]))*mfcoef(h1,2^a)+(a>=1)*y[14]*mfcoef(h1,(a>=1)*2^(a-1))+(a>=2)*y[15]*mfcoef(h1,(a>=2)*2^(a-2)))*mfcoef(h1,3^b)*H1;
	part2=2*real( ( (y[16]+k*(y[20]+2*y[22]))*mfcoef(h2,2^a)+(y[18]+k*(-y[21]-2*y[22]))*(1+sqrt(8)*I)*mfcoef(h2,2^a)+(a>=1)*(y[17]+y[19]*(1+sqrt(8)*I))*mfcoef(h2,(a>=1)*2^(a-1)))*mfcoef(h2,3^b)*H2);
	part7=2*real((y[35]*mfcoef(h7,2^a)+y[36]*(a>=1)*mfcoef(h7,(a>=1)*2^(a-1))+y[37]*(a>=2)*mfcoef(h7,(a>=2)*2^(a-2))+(-1+sqrt(3)*I)*(y[38]*mfcoef(h7,2^a)+y[39]*(a>=1)*mfcoef(h7,(a>=1)*2^(a-1))+y[40]*(a>=2)*mfcoef(h7,(a>=2)*2^(a-2))))*mfcoef(h7,3^b)*H7);
	part9=2*real(y[43]*mfcoef(h9,2^a)*mfcoef(h9,3^b)*H9)+2*imag(y[44]*sqrt(3)*mfcoef(h9,2^a)*mfcoef(h9,3^b)*H9);
	part11=(y[41]*mfcoef(f41,3^b)+y[42]*(b>=1)*mfcoef(f41,(b>=1)*3^(b-1)))*mfcoef(f41,2^a)*H11;
	part1+part2+part7+part9+part11;
}
*/