# Open project: Let's compute the Fourier expansion of a non-solvable algebraic Maass form.

OK so let's see if I can use MO to explicitly compute an example of something, by getting other people to join in. Sort of "one level up"---often people answer questions here but I'm going to see if I can make people do a more substantial project. Before I start, note

(1) the computation might have been done already [I'd love to hear of a reference, if it has]

(2) the computation may or may not be worth publishing

(3) If it is worth publishing, there may or may not be a debate as to who the authors are.

I personally don't give a hang about (2) or (3) at this point, but others might. Let me get on to the mathematics. Oh---just a couple more things before I start---this project is related to the mathematics at this question, but perhaps pushes it a bit further (if we can get it to work). I had initially thought about these issues because I was going to give them to an undergraduate, but the undergraduate tells me today that he's decided to do his project on the holomorphic case, and it seemed a bit daft to let my initial investment in the problem go to waste, so I thought I'd tell anyone who was interested. If no-one takes the bait here, I'll probably just make this another UG project.

OK so here's the deal. Say $K/\mathbf{Q}$ is a finite Galois extension, and $\rho:Gal(K/\mathbf{Q})\to GL(2,\mathbf{C})$ is an irreducible 2-dimensional representation. General conjectures in the Langlands philosophy predict that $\rho$ comes from an automorphic form on $GL(2)$ over $\mathbf{Q}$. The idea is that we are going to "see" this form in an explicit example where general theory does not yet prove that it exists.

Now the determinant of $\rho$ is a 1-dimensional Galois representation, and it makes sense to ask whether $det(\rho(c))$ is $+1$ or $-1$, where $c$ is complex conjugation. It has to be one of these, because $c^2=1$. The nature of the automorphic form predicted to exist depends on the sign.

If the determinant is $-1$ then the form should be holomorphic, and a classical weight 1 cusp form. In this case the existence of the form is known, because it is implied by Serre's conjecture, which is now a theorem of Khare and Wintenberger.

If the determinant is $+1$ then the conjectural form should be a real-analytic function on the upper half plane, invariant under a congruence subgroup, and satisfying a certain differential equation (which is not the Cauchy-Riemann equations in this case). If the image of $\rho$ is a solvable group then the existence of this form is known by old work of Langlands and Tunnell.

So in summary then, the one case where the form is not known to exist is when the determinant of complex conjugation is $+1$ and the image of Galois is not solvable.

Here is an explicit example. The polynomial

g5=344 + 3106*x - 1795*x^2 - 780*x^3 - x^4 + x^5


has splitting field $L$, an $A_5$-extension of $\mathbf{Q}$ ramified only at the prime 1951. Now $A_5$ is isomorphic to the quotient of $SL(2,\mathbf{F}_5)$ by its centre $\pm 1$, and $L$ has a degree two extension $K_0$, also unramified outside 1951, with $Gal(K_0/\mathbf{Q})$ being $SL(2,\mathbf{F}_5)$. Turns out that $K_0$ can be taken to be the splitting field of the rather messier polynomial

g24 = 14488688572801 - 2922378139308818*x^2 + 134981448876235615*x^4 -   1381768039105642956*x^6 + 4291028045077743465*x^8 -   2050038614413776542*x^10 + 287094814384960835*x^12 - 9040633522810414*x^14 + 63787035668165*x^16 - 158664037068*x^18 +   152929135*x^20 - 50726*x^22 + x^24


I should perhaps say that David Roberts told me these polynomials in Jan 2008; they're in a paper by him and John Jones---but they learnt about them from a paper of Doud and Moore.

Now $SL(2,\mathbf{F}_5)$ has two faithful 2-dimensional complex representations; the traces of each representation take values in $\mathbf{Q}(\sqrt{5})$ and one is of course the conjugate of the other. The determinant of both representations is trivial so in fact they are $SL(2,\mathbf{C})$-valued. Oh---also, all the roots of g24 are real---and hence $K_0$ is totally real. So what we have here is a representation $$\rho_0:Gal(K_0/\mathbf{Q})\to GL(2,\mathbf{C})$$ which is conjectured to come from automorphic forms, but, as far as I know, the conjecture is not known in this case.

Unfortunately the conductor of $\rho_0$ is $1951^2$, which is a bit big. In fact let me say something more about what is going on at 1951. In the $A_5$ extension the decomposition and inertia groups at 1951 are both cyclic of order 5. If my understanding of what David Roberts told me is correct, in the $SL(2,\mathbf{F}_5)$ extension the decomposition and inertia groups are both cyclic of order 10 (in fact I just got magma to check this). But the upshot is that $\rho_0$ restricted to a decomposition group at 1951 is of the form $\psi+\psi^{-1}$ with $\psi$ of order 10. Which character of order 10? Well there are four characters $(\mathbf{Z}/1951\mathbf{Z})^\times\to\mathbf{C}^\times$ of order 10, and two of them will do, and two won't, and which ones will do depends on which 2-dimensional representation of $SL(2,\mathbf{F}_5)$ you chose.

The key point though is that if you get $\psi$ right, then the twist $\rho:=\rho_0\otimes\psi$ will have conductor 1951, which is tiny for these purposes.

Now, as Marty did in an $A_4$ example and as I did and Junkie did in a dihedral example in the Maass form question cited above, it is possible to figure out explicitly numbers $b_1$, $b_2$, $b_3$,..., with the property that

$$L(\rho,s)=\sum_{n\geq1}b_n/n^s.$$

If one had a computer program that could calculate $b_n$ for $n\geq1$, then there are not one but two ways that one could attempt to give computational evidence for the predictions given by the Langlands philosophy:

(A) one could use techniques that Fernando Rodriguez-Villegas explained to me a few months ago to try and get computational evidence that $L(\rho,s)$ had analytic continuation to the complex plane and satisfied the correct functional equation, and

(B) one could compute the corresponding real analytic function on the upper half plane, evaluate it at various places to 30 decimal places, and see if the function was invariant under the group $\Gamma_1(1951)$.

I don't know much about (A) but I once tried, and failed, to do (B), and my gut feeling is that my mistake is in the computer program I wrote to compute the $b_n$. But as junkie's response in the previous question indicates, there now seem to be several ways to compute the $b_n$ and one thing I am wondering is whether we can use the methods he/she indicated in this question.

Let me speak more about how I tried to compute the $b_n$. The character of the Maass form is $\psi^2$, the determinant of $\rho$. General theory tells us that $b_n$ is a multiplicative function of $n$ so we only need compute $b_n$ for $n$ a prime power. Again general theory (consider the local $L$-functions) says that for $p\not=1951$ one can compute $b_{p^n}$ from $b_p$. If $p=1951$ then $b_{p^n}=1$ for all $n$, because the decomposition and inertia groups coincide for 1951 in $K_0$ [EDIT: This part of the argument is wrong, and it explains why my programs didn't work. I finally discovered my mistake after comparing the output of mine and Junkie's programs and seeing where they differed. It's true that decomposition and inertia coincide in $K_0$ but when one twists by the order 10 character this stops being true. In fact $b_{1951}$ is a primitive 5th root of unity that I don't know how to work out using my method other than by trial and error.]. Finally, if the $L$-function of $\rho_0$ is $\sum_n a_n/n^s$ then $b_n=\psi(n)a_n$ for all $n$ prime to 1951, so it suffices to compute $a_p$ for $p\not=1951$.

To compute $a_p$ I am going to compute the trace of $\rho_0(Frob_p)$. I first compute the GCD of the degrees of the irreducible factors of g24 mod $p$. If $p$ doesn't divide the discriminant of g24 then this GCD is the order of $Frob_p$ in $SL(2,\mathbf{F}_5)$. If $p$ does divide the discriminant of g24 then rotten luck, I need to factorize $p$ in the ring of integers of the number field generated by a root of g24. I did this using magma. Here are the results:

                                                   prime   order
2       6
3      10
5       4
163       5
16061       1
889289      10
451400586583       2
1188493301983785760551727       2
120450513180827412314298160097013390669723824832697847       1


Now unfortunately just computing the order of the conj class of Frobenius is not enough to determine the trace of the Galois representation, because $SL(2,\mathbf{F}_5)$ contains two conj classes of elements of order 5, and two of order 10. However, in both cases, the conj classes remain distinct in $A_5$ and so it suffices to have an algorithm which can distinguish between the two conjugacy classes in $A_5$. More precisely, we have to solve the following problem: we label the conj classes of elements of order 5 in $A_5$ as C1 and C2, and we want an algorithm which, given a prime $p$ for which g5 is irreducible mod $p$, we want the algorithm to return "C1" or "C2" depending on which class $Frob_p$ is in. Here is a beautiful way of doing it, explained to me by Bjorn Poonen: if g5 is irreducible mod $p$ then its roots in an alg closure of $\mathbf{F}_p$ are $x,x^p,x^{p^2},...$. Set $x_i=x^{p^i}$ and compute $\prod_{i&lt;j}(x_i-x_j)$. This product is a square root of the discriminant of g5. Choose once and for all a square root of the discriminant of g5 in the integers; if the product is congruent to this mod $p$ return "C1", else return "C2".

That's it! I implemented this. I had a program which returned a bunch of $a_n$'s, and hence a bunch of $b_n$'s. I built the function on the upper half plane as the usual sum involving Bessel functions and so on, but it computationally did not come out to be invariant under $\Gamma_1(1951)$.

If anyone wants to take up the challenge of computing the $b_n$ that would be great. I have explained one way to do it above, but I am well aware that there might be other ways to compute the $b_n$ analogous to junkie's approach from the previous question.

• I should say: if anyone doesn't understand some stuff I wrote above, or wants to know how I did something, or wants me to justify a statement I made, then feel free to ask me via a comment (which I will notice because comments made here will be flagged to me) and I'll try to reply. – Kevin Buzzard May 5 '10 at 16:00
• Kevin, you might be interested in some papers of Andrew Booker, especially (maths.bris.ac.uk/~maarb/papers/modularity.pdf) where he computes the first 10^9 coefficients of an even A5 representation and tests its modularity. – David Hansen May 5 '10 at 16:03
• Hi James. Bump's big grey book is a good place to read about these things. What's going on is that local Langlands at the infinite place says that to rho:Gal(C/R)-->GL_2(C) there's a rep pi of GL_2(R). If rho is odd then pi is a limit of discrete series. If you consider pi as a rep of the Lie algebra of GL_2(R) then it's a rep of the univ enveloping algebra of this Lie algebra, and elts of the univ enveloping algebra are differential operators. In the odd case you can find an element annihilated by the diff op corresponding to the Cauchy-Riemann equations... – Kevin Buzzard May 6 '10 at 11:21
• Just a heads up: the David Roberts in the question is not me! I knew this would happen eventually...:) – David Roberts May 6 '10 at 14:12
• Of course it's not my area, but I think this is a fantastic use of MO. +1 (I wish it could be +10) – GS May 6 '10 at 14:26

Here is Magma code that gets you the answer in a few seconds. I made a special case for the bad primes, and did them by hand.

_<x> := PolynomialRing(Rationals());
f5 := 344 + 3106*x - 1795*x^2 - 780*x^3 - x^4 + x^5;
g24 := 14488688572801 - 2922378139308818*x^2 + 134981448876235615*x^4 -
1381768039105642956*x^6 + 4291028045077743465*x^8 -
2050038614413776542*x^10 + 287094814384960835*x^12 -
9040633522810414*x^14 + 63787035668165*x^16 - 158664037068*x^18 +
152929135*x^20 - 50726*x^22 + x^24;
K := NumberField(f5);
_,D := IsSquare(Integers()!Discriminant(f5));
prec := 30;
CHAR_TABLE := CharacterTable(GaloisGroup(g24));
chi := CHAR_TABLE;

[ <2,Polynomial([1,-1,1])>,
<3,Polynomial([1,-ComplexField(prec)!chi,1])>,
<5,Polynomial([1,0,1])>,
<7,Polynomial([1,0,1])>,
<71,Polynomial([1,0,1])>,
<137,Polynomial([1,1,1])>,
<163,Polynomial([1,-ComplexField(prec)!chi,1])>,
<1951,Polynomial()>,
<16061,Polynomial([1,-2,1])>,
<889289,Polynomial([1,-ComplexField(prec)!chi,1])> ];
FACTORS := [bf : bf in BAD_FACTORS];

function LOCAL(p,d : Precision:=prec)
R := Roots(ChangeRing(f5,GF(p)));
if #R eq 1 then return Polynomial([1,0,1]); end if;
if #R eq 2 then
ord := Lcm([Degree(f) : f in Factorization(Polynomial(GF(p),g24))]);
return Polynomial([1,ord eq 3 select 1 else -1,1]); end if;
if #R eq 5 then
ord := Lcm([Degree(f) : f in Factorization(Polynomial(GF(p),g24))]);
return Polynomial([1,ord eq 1 select -2 else 2,1]); end if;
r := Roots(ChangeRing(f5,GF(p^5)));
x := r;
prod := GF(p)!&*[x^(p^i)-x^(p^j) : j in [(i+1)..4], i in [0..4]];
wh := prod eq GF(p)!D;
ord := Lcm([Degree(f) : f in Factorization(Polynomial(GF(p),g24))]);
if ord eq 10 then class := wh select 8 else 9; // compatible with FACTORS
else class := wh select 6 else 5; end if;
return Polynomial([1,-ComplexField(prec)!chi[class],1]);
end function;

L := LSeries(1, [0,0], 1951^2, LOCAL : Precision:=prec);
// s->1-s, Gamma(s/2)^2
psi := DirichletGroup(1951, CyclotomicField(10)).1;
p1951 := Polynomial([1,-ComplexField(prec)!CyclotomicField(5).1]);
TP := TensorProduct(L, LSeries(psi : Precision:=prec), [<1951, 1, p1951>]);
CheckFunctionalEquation(TP);


Here is the special values:

ev := Evaluate(TP,0); // 2-1.453085056...
rel := PowerRelation(ev,4 : Al:="LLL");
NF := NumberField(rel);
Q5<zeta5> := CyclotomicField(5);
assert IsIsomorphic(NF,Q5);
Q5!NF.1;


So $L(\rho,0)=-4\zeta_5(1+\zeta_5)$ for Marty. I get $L(\rho_0,-1)=32(48723\sqrt{5} - 778741)$ as an algebraic. I get $L(\rho,-2)=8800\zeta_5^3 - 14444\zeta_5^2 + 35604\zeta_5 + 17412$ with more precision. I determined the TensorProduct factor at 1951 via trial and error, making the obvious guesses until one worked (the failure is at 100-110 digits). With this, I take it to 240 digits and I can even get $$L(\rho,-4)=-18475535360\zeta_5^3 - 11142861380\zeta_5^2 - 12091894020\zeta_5 - 7107607296$$ and $$L(\rho,-6)=25255057273186244\zeta_5^3 - 1015274469604000\zeta_5^2 - 15695788409197884\zeta_5 + 9459547822189412$$ The precision can go higher if you want more.

Finally, the Maass form:

function MaassEval(L,z)
x:=Real(z); y:=Imaginary(z);
printf "Using %o coefficients\n", Ceiling(11/y);
C := LGetCoefficients(L,Ceiling(11/y));
pi := Pi(RealField());
a := Sqrt(y)*&+[C[n]*KBessel(0,2*pi*n*y)*Sin(2*pi*n*x) : n in [1..#C]];
return a;
end function;
zz:=0.0001+0.0001*ComplexField().1;
MaassEval(TP,zz);
// Using 110000 coefficients
// -1.71477211817772949974178783985E-8 + 9.01673609747756708674470686948E-9*i
MaassEval(TP,zz/(1951*zz+1));
// Using 161297 coefficients
// -1.71477211817772949974179078240E-8 + 9.01673609747756708674496293450E-9*i

• Well done Junkie. I'll take a look at this on Monday or Tuesday when I have some time! Thanks a lot! – Kevin Buzzard May 9 '10 at 14:18
• You can also de-Magma-ize the code with little difficulty. The only hard work was computing the Euler factors at bad primes via the decomposition, now that it done once. The rest appeals to factoring polynomials over finite fields. The LSeries internals is not too deep, as you merely require the Mellin transforms for $\Gamma(s/2)^2$ and $\Gamma((s+1)/2)^2$, which are some variant of Bessel function I suppose, and the Magma code is nothing special for them but uses Dokchitser's general paper. Naming the special values at negative integers applies by LLL as with algdep in PARI. – Junkie May 10 '10 at 3:18
• Haven't looked at this in a while... this is an outstanding computation! Maybe it's well-known to experts that these L-values lie in the cyclotomic field generated by a 5th root of unity -- the coefficient field of the Artin motive, I suppose. Thanks for posting the code! – Marty May 16 '10 at 17:50
• Junkie: when I said above "I'll take a look at this on Monday or Tuesday" I obviously meant "I'll take a look at this in just under a year's time". Apologies for this. I took a look at it today. I absolutely believe you. What you didn't answer was where I had gone wrong! But I finally took the trouble to understand your code and compare output with mine. I computed the first 10,000 coefficients of the L-function according to you, and the first 10,000 coefficients according to me, and I compared. Three were different :-/ They were the coefficients corresponding to $q^{1951d}$ for $d=1,2,3$! – Kevin Buzzard Apr 28 '11 at 21:56
• [we got the same answer for $d=4,5$ because these coefficients vanish]. I finally realised my error. I gave an argument in the main question that $b_p=1$ for $p=1951$ and it's not right. If $\chi_1$ and $\chi_2$ are two finite order characters of the absolute Galois group of $\mathbf{Q}_p$ and they both cut out totally ramified extensions of $\mathbf{Q}_p$ then, of course, their product might not. Hence my argument that $b_{1951}=1$ is not correct. In fact it's a random 10th root of unity. Thanks so much for enabling me to find my error! I'll edit my original post and finally acceptyouranswer! – Kevin Buzzard Apr 28 '11 at 22:01

I have accepted Junkie's answer to this question. But on analysing his working program and comparing it with my program-that-didn't-work I could easily find the bug in my program: I erroneously thought that the coefficient of $q^{1951}$ was 1 and it's not, it's a primitive 5th root of unity. My code works now, and is written in PARI-GP, so I thought I would post it so that people who don't have access to magma can play with the miracle that is a cuspidal Maass eigenform that (conjecturally) corresponds to a non-solvable even Galois representation. I make my post community wiki because my post says nothing that Junkie's doesn't so I don't want any credit for it.

\\ nonsolvmaass.g: works out a maass form associated to an explicit
\\ nonsolvable finite image even irred 2-d Galois rep.
\\ Note: rho(c)=-1 so we use sines.

\\ N is number of coefficients of the L-function we will compute.

N=2000000;

\\ g5 has splitting field an A_5 field.

g5=344 + 3106*x - 1795*x^2 - 780*x^3 - x^4 + x^5

\\ g24 has splitting field a field with Galois group SL(2,5).

g24 = 14488688572801 - 2922378139308818*x^2 + 134981448876235615*x^4 -   1381768039105642956*x^6 + 4291028045077743465*x^8 -   2050038614413776542*x^10 + 287094814384960835*x^12 -   9040633522810414*x^14 + 63787035668165*x^16 - 158664037068*x^18 +   152929135*x^20 - 50726*x^22 + x^24

\\ The character table of SL(2,5) starts like this:
\\
\\
\\ -------------------------------------------
\\ Class |   1  2  3  4    5    6  7    8    9
\\ Size  |   1  1 20 30   12   12 20   12   12
\\ Order |   1  2  3  4    5    5  6   10   10
\\ -------------------------------------------
\\ X.1   +   1  1  1  1    1    1  1    1    1
\\ X.2   -   2 -2 -1  0   Z1 Z1#2  1  -Z1-Z1#2
\\ X.3   -   2 -2 -1  0 Z1#2   Z1  1-Z1#2  -Z1
\\
\\...
\\
\\ # denotes algebraic conjugation, that is,
\\ #k indicates replacing the root of unity w by w^k
\\
\\ Z1     = (CyclotomicField(5)) ! [ RationalField() | 0, 0, 1, 1 ]
\\
\\ If we let K be the splitting field of g24 then K contains a splitting
\\ field for g5. If we fix an isomorphism of Gal(K/Q) with SL(2,5) and
\\ use X.2 to define a 2-dimensional representation of Gal(K/Q), then
\\ the resulting representation has conductor 1951^2.
\\
\\ Here is a pari script to compute the order of Frob_p for the resulting
\\ representation (for p not 1951). First a function that only works for the primes
\\ not dividing disc(g24).

ordfrob0(p)=local(m);m=factormod(g24,p);lcm(vector(#m~,i,poldegree(m[i,1])))

\\ For the primes dividing disc(g24) other than 1951 I used magma to work
\\ out the answer by hand [computing a maximal order etc].

ordfrob(p)=if(p==1951,0,if(p==2,6,if(p==3,10,if(p==5,4,if(p==163,5,if(p==16061,1,if(p==889289,10,if(p==451400586583,2,if(p==1188493301983785760551727,2,if(p==120450513180827412314298160097013390669723824832697847,1,ordfrob0(p)))))))))))

\\ If the order of Frob_p in SL(2,5) is 5 or 10 then we have to decide which
\\ of the two conjugacy classes p is in! Classes 5 and 8 reduce to the same
\\ class in A_5, and classes 6 and 9 reduce to the other. So we have to
\\ be able to distinguish between the classes in A_5. Bjorn Poonen told
\\ me how to do this. The trick is to compute the product whose square
\\ is the discriminant mod p and to see which of the square roots of the
\\ discriminant that the product itself is isomorphic to.
\\ This will not work for any p dividing the discriminant of g5.
\\ Fortunately I checked explicitly that no prime p dividing disc(g5)
\\ has Frob_p of order 5 or 10. Hooray for good fortune!
\\ The function returns 0 or 1.

whichclass(p)=local(t,conjs,pro);t=Mod(x,Pol(Mod([1, -1, -780, -1795, 3106, 344],p)));conjs=[t];for(i=1,4,conjs=concat(conjs,conjs[i]^p));pro=1;for(i=1,5,for(j=i+1,5,pro=pro*(conjs[i]-conjs[j])));pro==518348075378;

\\ Now finally an algorithm to compute a_p.

zet10=exp(2*Pi*I/10);Z1=zet10^4+zet10^6;Z1h2=zet10^2+zet10^8;

ap(p)=if(p==1951,0,orp=ordfrob(p);if(orp==1,2,if(orp==2,-2,if(orp==3,-1,if(orp==4,0,if(orp==6,1,whp=whichclass(p);if(orp==5,if(whp==0,Z1,Z1h2),if(whp==0,-Z1,-Z1h2))))))))

\\ The problem with the Maass form associated to these eigenvalues is that, although it
\\ has trivial character, it has conductor 1951^2. At a prime above 1951, the
\\ inertia subgroup has order 10, as does the decomposition group. The representation
\\ sends a generator of this inertia subgroup to something with eigenvalues zet10
\\ and zet10^9. So we need to twist by a Dirichlet character of conductor 1951
\\ and order 10 and we will get a representation of conductor 1951.
\\ But which one?! Well, there are only four, and two will work, so we have
\\ a sporting chance! If we get it wrong we can just try another one.
\\ Note that 3 is a primitive root mod 1951.

twist=vector(1950,i,0);z0=zet10^3;z=z0;n=Mod(3,1951);for(i=1,1950,twist[lift(n)]=z;n=n*3;z=z*z0);

\\ The character of the form we are after is the square of "twist".

\\ Here's the trace of Frobenius on the conductor 1951 twist.
\\ Note that the value at p=1951 was found by trial and error.
\\ The error I made, that Junkie fixed, was that I thought that "pure thought" gave
\\ that bp(1951)=1.

bp(p)=if(p==1951,zet10^6,ap(p)*twist[p%1951])

\\ Here's the character of the form.

chi(p)=if(p==1951,0,twist[p%1951]^2);

\\ Now v will be the vector of Hecke eigenvalues.

v=vector(N,i,0);
v=1;
for(i=2,N,if(i%1000==0,print(N-i));fac=factor(i);k=matsize(fac);\
if(k>1,v[i]=prod(j=1,k,v[fac[j,1]^fac[j,2]]),\
if(fac[1,2]==1,v[i]=bp(i),\
p=fac[1,1];e=fac[1,2];v[i]=v[p]*v[p^(e-1)]-chi(p)*v[p^(e-2)]))\
);

\\ Now glue them together with the Maass form analogue of a q-expansion.
\\ The function q^n is replaced by sqrt(y)*K(2.Pi.n.y)*sin(2.Pi.n.x)
\\ where K=K_0 is a certain Bessel function.

F(z)=local(x,y,M);x=real(z);y=imag(z);M=ceil(11/y);if(M>N,error("y too small."));sqrt(y)*sum(n=1,M,if(v[n]==0,0,v[n]*besselk(1e-30*I,2*Pi*n*y)*sin(2*Pi*n*x)))


If you cut and paste that into a pari session then it will either run out of memory, or spend a long time counting down from 2,000,000 to 0 and then stop. If it runs out of memory then type something like "allocatemem(2000000000)" or some such thing. If you can't be bothered to wait for the countdown to finish then change the "$N=$" line to something smaller at the beginning. It's computing the first $N$ coefficients of the Maass form.

The form itself is the function $F$ defined on the last line. The larger you let $N$ be, the better precision $F$ will be defined to, and hence the nearer to the real line you'll be able to evaluate $F$. The point is that if $z$ is in the upper half plane but very close to the real line then you need a lot of terms of $F$ to evaluate $F(z)$ accurately.

As in Junkie's post, we may try a random verification, to see whether $F$ looks like it's invariant under $\Gamma_1(1951)$:

(22:47) gp > zz4=0.0001+0.0001*I
%7 = 0.0001000000000000000000000000000 + 0.0001000000000000000000000000000*I
(22:47) gp > F(zz4)
%8 = 0.000000008720947535106531893030584189 - 0.00000001730002410267070596763764527*I
(22:48) gp > F(zz4/(1951*zz4+1))
%9 = 0.000000008720947535106531893415409754 - 0.00000001730002410267070596828184932*I
(22:51) gp > %8-%9
%10 = -3.848255644 E-28 + 6.442040498 E-28*I


The point is that the difference between the two values of $F$ is a tiny number, so $F$ is, to within computational error, seemingly invariant under $\Gamma_1(1951)$. People can play with other examples, but be warned: the smaller $N$ is (i.e. the less time you could be bothered to wait at the beginning) the more careful one has to be (i.e. the more likely the program is to refuse to compute $F$ where you want it to compute it).

• I'll remark that the reason Junkie and I get completely different values for $F(0.0001+0.0001*I)$ is that in fact, by chance, we computed different forms: there are two conjugate eigenforms, I computed one and he computed the other. – Kevin Buzzard Apr 29 '11 at 9:15

Hmm, a problem. In 4 hours Magma can compute this:

> f:=x^24 - 50726*x^22 + 152929135*x^20 - 158664037068*x^18 + 63787035668165*x^16 -
9040633522810414*x^14 + 287094814384960835*x^12 - 2050038614413776542*x^10 +
4291028045077743465*x^8 - 1381768039105642956*x^6 + 134981448876235615*x^4 -
2922378139308818*x^2 + 14488688572801;
> L:=LSeries(NumberField(f) : Method:="Artin");


But the problem is that the decomposition of Lprod contains no 2-dimensional factors...

> [<Character(xparent),x> : x in Lprod | #xlpoles eq 0];
[
<( 3, 3, 0, -1, -zeta(5)_5^3 - zeta(5)_5^2, zeta(5)_5^3 + zeta(5)_5^2 + 1,
0, zeta(5)_5^3 + zeta(5)_5^2 + 1, -zeta(5)_5^3 - zeta(5)_5^2 ), 1>,
<( 3, 3, 0, -1, zeta(5)_5^3 + zeta(5)_5^2 + 1, -zeta(5)_5^3 - zeta(5)_5^2,
0, -zeta(5)_5^3 - zeta(5)_5^2, zeta(5)_5^3 + zeta(5)_5^2 + 1 ), 1>,
<( 5, 5, -1, 1, 0, 0, -1, 0, 0 ), 1>,
<( 6, -6, 0, 0, 1, 1, 0, -1, -1 ), 2>
]


So we achieved a decomposition $1+\rho_3+\bar\rho_3+\sigma_5+2\eta_6$ with this, but no degree 2 factor? Is this correct? Maybe the splitting field is necessary here (when all reps will appear). The guess is that the above computation already did this in essence, so I can re-run it after adding more intelligence.

• Even it this works, I don't know if Dokchitser's code will really be too fast in practice to get a billion terms. The splitting field (degree 120) was found in 1.5 hours. – Junkie May 6 '10 at 10:34
• Well let's have a think! All I know about the degree 24 poly is that its splitting field is the SL(2,5) extension. I don't know which subextension it corresponds to so let's try to work it out. We want a degree 24 subextension of a degree 120 extension so we want a subgroup fixed by a group of order 5, so it's a Sylow 5-subgroup. OK so what's going on is probably something like this: if L0 is the field obtained by adjoining one root of the deg 24 poly and K0 is the splitting field then K0/L0 has degree 5, and the L-function of L0 is then a product of five L-functions for L0, namely... – Kevin Buzzard May 6 '10 at 10:41
• ...the L-functions L(L0,alpha^n) where alpha is a non-trivial character Gal(K0/L0)-->C^* and 0<=n<=4. You only computed one of these factors and it will probably correspond to something like the representation Ind_P^SL(2,5)(1). So you'll see the 2-d rep here iff its restriction to P is trivial---which it isn't. So that's why you don't see the degree 2 rep. Can you compute L(deg 24,alpha) for alpha a character of the Galois group of order 5 with kernel K0? – Kevin Buzzard May 6 '10 at 10:43
• PS I don't know if this helps, but if you run the following in pari: g24 = (blah blah blah) X=polcompositum(g24,g24) then under 6 seconds later I get a vector consisting of eight polynomials, and the four of degree 120 are polynomials of degree 120 such that if you adjoin one root of them to Q you get the splitting field of g24. I mention this because you claimed that something similar took 1.5 hours for you. – Kevin Buzzard May 6 '10 at 10:45
• PPS of course I don't care if the calculation takes a week! ;-) – Kevin Buzzard May 6 '10 at 10:46

I will copy a comment of mine about Jehanne over here.

His paper is: http://dx.doi.org/10.1006/jnth.2001.2656

Jehanne has 4 totally real Examples 2-5 (page 353-6), including the one of Booker. He also gets the degree 24 field. The method to compute the $L$-series is listed in Section 6.

Jehanne gives Poonen's trick as being due to Serre (from Buhler's book). Booker says the same (page 332).

As Kevin Buzzard pointed out, the application to computational Stark's conjecture was only for complex $A_5$ though, so the explicit totally real computation might not have been done outside of Booker.

I don't think that Dokchitser does anything in his Magma code that is not in essence the same as you do here, as the splitting types are just realized in a alternative manner, usually more complicated due to generality. I think his code could probably do $A_4$ as a black box w/o much modification, but $A_5$ requires a degree 240 number field.

• I've never been able to discern philosophically the relation of Booker's smooth sums for modularity versus the alternative approach of testing the functional equation. For the latter, you have the completed L-function as a sum $\sum_n F(nxA) + (sign)\sum_n G(nx/A)$ for any $A$ as a 1-parameter family, though Rubinstein can widen the test functions even more I think. Then you just compute at various $A$ values, for instance $A=1.1, 1.2$ to see if the result is self-consistent. Booker draws graphs, but the statistical content of them is not directly transparent. – Junkie May 6 '10 at 2:53
• I guess I am wrong about the degree 240, as it is only 120. You want $\tilde A_5$ (2-extension) and not $\hat A_5$ (4-extension) as the latter gives quadratic determinant, to quote the terminology of Jehanne, with the former as trivial determinant. The 120 still looks too big. The octic fields for $163$ and $277$ are already in Bachoc and Kwon (3.2) matwbn.icm.edu.pl/ksiazki/aa/aa62/aa6211.pdf for $\tilde A_4$. The Magma Artin representation coding takes little time to compute this as a black box via the L-series, if it is what is correct that we are looking for in that case. – Junkie May 6 '10 at 3:43
• Yes---the Galois group is A_5-tilde, which is isomorphic to SL(2,5). I was hoping that you'd be able to get magma to do it for free! – Kevin Buzzard May 6 '10 at 6:39

Actually, I should read the help with Magma, for the direction to ArtinRepresentations is there.

  f:=14488688572801 - 2922378139308818*x^2 + 134981448876235615*x^4 -
1381768039105642956*x^6 + 4291028045077743465*x^8 -
2050038614413776542*x^10 + 287094814384960835*x^12 -
9040633522810414*x^14 + 63787035668165*x^16 - 158664037068*x^18 +
152929135*x^20 - 50726*x^22 + x^24;
K:=NumberField(f);
A:=ArtinRepresentations(K); // Four hours for computing
KartinrepdataFrob(7); // fails


Then A has degree 2, though the conductor is $1951^2$ as you say. Magma does tensor product L-functions, but maybe not automatically at bad primes. But the hassle is, it cannot compute the Frobenius for many primes, probably due to the lack of cycle type distinction. I don't think Magma is the right tool for such a specific problem. Maybe the internal Frob function can be overridden and it will work via the Serre-Poonen trick.

I got it to work!

First the preliminary data:

_<x> := PolynomialRing(Rationals());
f5 :=344 + 3106*x - 1795*x^2 - 780*x^3 - x^4 + x^5;
g24 :=14488688572801 - 2922378139308818*x^2 + 134981448876235615*x^4 -
1381768039105642956*x^6 + 4291028045077743465*x^8 -
2050038614413776542*x^10 + 287094814384960835*x^12 -
9040633522810414*x^14 + 63787035668165*x^16 - 158664037068*x^18 +
152929135*x^20 - 50726*x^22 + x^24;
K := NumberField(g24);
_,D := IsSquare(Discriminant(f5)); D:=Integers()!D;


Then the long computation of ArtinRepresentations. All the work is in the 7 or so bad primes you listed.

SetVerbose("ArtRep",3);
A:=ArtinRepresentations(K); // Four hours for computing
BAD:= Kartinrepdatabadprimes cat [p : p in Factorization(D)];
chi:=Character(A);
L:=LSeries(A);


Then the changed local function.

FUNC:=Lcffun;
function LOCAL(p,d : Precision:=30)
if p in BAD then return FUNC(p,d : Precision:=30); end if;
R:= Roots(ChangeRing(f5,GF(p)));
if #R eq 1 then return Polynomial([1,0,1]); end if;
if #R eq 2 then
ord := Lcm([Degree(f) : f in Factorization(Polynomial(GF(p),g24))]);
return Polynomial([1,ord eq 3 select 1 else -1,1]); end if;
if #R eq 5 then
ord := Lcm([Degree(f) : f in Factorization(Polynomial(GF(p),g24))]);
return Polynomial([1,ord eq 1 select -2 else 2,1]); end if;
r:= Roots(ChangeRing(f5,GF(p^5)));
x:=r;
prod:= GF(p)!&*[x^(p^i)-x^(p^j) : j in [(i+1)..4], i in [0..4]];
wh:=prod eq GF(p)!D;
ord:= Lcm([Degree(f) : f in Factorization(Polynomial(GF(p),g24))]);
if ord eq 10 then class:=wh select 8 else 9; // could be reversed, see BAD
else class:=wh select 6 else 5; end if;
return Polynomial([1,-ComplexField(30)!chi[class],1]);
end function;


This is not optimal with the GF(p) I think, but who cares. A difficulty is that the BAD primes force an ordering onto the conjugacy classes on 5-cycles, so there are two such functions I suspect.

Now change the coefficient function of L to be as desired.

Lcffun := LOCAL;
psi := DirichletGroup(1951, CyclotomicField(10)).1;
TP := TensorProduct(L,LSeries(psi), [<1951, 1>]); // tried psi^3 too
CheckFunctionalEquation(TP); // -4.7331654313E-30 - works!
LSetPrecision(L, 9); CheckFunctionalEquation(L); // L directly


I tried both $\psi$ and $\psi^3$ in the TensorProduct. The L-function is actually wrong at 1951 under this, but the computation only goes to $C\sqrt{1951}$ terms or about 700 so it does not matter. Also, you can check L directly, reducing the precision.

• Hey well done! I am going away for the weekend so can't play with this right now :-( If you put the L-series into the "F" formula in the other thread, possibly changing sin to cos or vice-versa (one corresponds to rho(c)=+1 and the other to rho(c)=-1) I wonder if you can find a Gamma-invariant function on the upper half plane? This is great. I think Marty will also be very interested in this L-function because I think he wants to look at its special values or something. – Kevin Buzzard May 7 '10 at 16:43

@Junkie: here are the first 100 coefficients according to my computer program that doesn't work:

[1, 0.3090 + 0.9511*I, -0.5000 + 1.539*I, 0, 0, -1.618, 0, 0.8090 + 0.5878*I, -1.309 - 0.9511*I, 0, 1.309 - 0.9511*I, 0, 0.8090 - 0.5878*I, 0, 0, -0.3090 + 0.9511*I, 0.5000 - 1.539*I, 0.5000 - 1.539*I, 0.5000 + 0.3633*I, 0, 0, 1.309 + 0.9511*I, 0, -1.309 + 0.9511*I, 0.8090 + 0.5878*I, 0.8090 + 0.5878*I, 0.8090 - 0.5878*I, 0, 0, 0, 0.6180, 0, 0.8090 + 2.490*I, 1.618, 0, 0, 1.309 + 0.9511*I, -0.1910 + 0.5878*I, 0.5000 + 1.539*I, 0, -0.5000 - 1.539*I, 0, -0.8090 + 0.5878*I, 0, 0, 0, 0, -1.309 - 0.9511*I, -0.3090 - 0.9511*I, -0.3090 + 0.9511*I, 2.118 + 1.539*I, 0, 0.5000 - 0.3633*I, 0.8090 + 0.5878*I, 0, 0, -0.8090 + 0.5878*I, 0, -0.3090 - 0.9511*I, 0, -1.618, 0.1910 + 0.5878*I, 0, 0.3090 + 0.9511*I, 0, -2.118 + 1.539*I, 0.8090 - 0.5878*I, 0, 0, 0, 0, -0.5000 - 1.539*I, 0, -0.5000 + 1.539*I, -1.309 + 0.9511*I, 0, 0, -1.309 + 0.9511*I, -0.1910 + 0.5878*I, 0, 0, 1.309 - 0.9511*I, -0.1910 + 0.5878*I, 0, 0, -0.8090 - 0.5878*I, 0, 1.618, 1.000, 0, 0, 0, -0.3090 + 0.9511*I, 0, 0, 0, 0.1910 + 0.5878*I, 0.8090 - 0.5878*I, -2.618, 0]


This is of course a comment but it's far too big for the comment box. Junkie: presumably the actual answers are different?

If the above vector is v then my script for computing the Maass form is this:

F(z)=local(x,y,M);x=real(z);y=imag(z);M=ceil(11/y);if(M>N,error("y too small."));sqrt(y)*sum(n=1,M,if(v[n]==0,0,v[n]*besselk(1e-30*I,2*Pi*n*y)*sin(2*Pi*n*x)))


One should get F(z)=F((az+b)/(cz+d)) for (a,b;c,d) in Gamma_1(1951) but it doesn't work out. I've made a slip somewhere. Hmm---I can't remember whether I knew for sure whether I'd twisted by the right order 5 character. Let me know how these figures compare with yours anyway!