\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ mestre.gp \\ \\ \\ \\ This program "quickly" computes the action of the Hecke operators \\ \\ T_2, T_3, T_5, and T_7 on the space M_2(Gamma_0(N)) of modular forms \\ \\ of weight 2 and level a prime N > 7. In the process it has the \\ \\ "side effect" of computing all supersingular j-invariants in \\ \\ in characteristic N. \\ \\ \\ \\ AUTHOR: William Stein, was@math.berkeley.edu, \\ \\ University of California, Berkeley \\ \\ \\ \\ ALGORITHM: The algorithm is was developed by J.F. Mestre and \\ \\ Oesterle in connection with their work on a modular curve of \\ \\ conductor 5077. It is described in Mestre, "La methode des \\ \\ graphes. Exemples et applications." which is in "Proceedings \\ \\ of the international conference on class numbers and fundamental \\ \\ units, Katata, Japan, June 1986." The idea is that the Hecke \\ \\ operators act on the supersingular points of the modular curve \\ \\ X_0(1) in characteristic N and that the action on these points \\ \\ is enough to see the Hecke operators. To actually compute the \\ \\ action use modular polynomials. \\ \\ \\ \\ \\ \\ MODULAR POLYNOMIALS: In the algorithm it is necessary \\ \\ to quickly determine when two curves are p-isogeneous. \\ \\ To do this we use the modular polynomails PHIp. Mestre's paper \\ \\ gives PHI2 and PHI3. Unfortunately PHI3 contains a typo. I obtained\\ \\ (correct) versions of PHI2,PHI3,PHI5, and PHI7 from \\ \\ "Uber die Berechnung der Fourierkoeffizienten der Funktion j(tau)" \\ \\ by Von Oskar Herrmann, Journal Fur der Reine Agnew, 76???. \\ \\ [Has anyone computed PHI11 yet? It'd be nice to have it!] \\ \\ \\ \\ To compute all p-isogeneous curves gives this problem: \\ \\ PROBLEM: Let f be a polynomial with coefficients in F_{p^2}. \\ \\ Quickly find all roots of f assuming that \\ \\ all roots already lie in F_{p^2}. \\ \\ The solution I devised is: Let g = f*frob(f) which has \\ \\ coefficients in F_p. Factor g over F_p using a standard \\ \\ alogorithm (perhaps Berkelekamp's). Since the roots of f \\ \\ lie in F_{p^2}, the roots of g do also and hence g factors \\ \\ into a product of linear and quadratic factors. Now use the \\ \\ quadratic formula to find the roots of the quadratic factors. \\ \\ Finally check which roots of g are in fact roots of f. \\ \\ \\ \\ CHARACTERISTIC POLYNOMIAL: pari's built in function for computing \\ \\ characteristic polynomials of a matrix M is no good if M is large \\ \\ (with say > 100 entries). Instead I use a method which essential \\ \\ consists of computing traces of powers of M. \\ \\ \\ \\ CUSP FORMS: If f(x) is the characteristic polynomial of T_p acting \\ \\ on M_2(Gamma_0(N)), N prime, then f(x) / (x-(p+1)) is the \\ \\ characteristic polynomial of T_p acting on the cusp \\ \\ forms S_2(Gamma_0(N)). \\ \\ \\ \\ \\ \\ LAST MODIFIED: April 21, 1998. \\ \\ \\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ table_of_modp_roots: \\ In order to optimize speed we precompute a table of all \\ squares modulo p. \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ { table_of_modp_roots (p, locals) = v = vector (p-1, X, 0) ; \\ set v[i^2] = i (mod p). for (i = 1, p - 1, v[ (i * i) % p ] = Mod(i,p) ) ; v ; } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ my_sqrtmodp: compute the square root of b modulo p \\ and express it in terms of the quadratic extension of F_p \\ defined by a. \\ \\ * b is an integer modulo p \\ * a is assumed to be of the form w^2 + c, \\ e.g., Mod(1, 5)*w^2 + Mod(2, 5) \\ \\ +++ THE VARIABLE "table" must be set in the global name space \\ using the function table_of_modp_roots above. \\ Thus give a line like \\ table = table_of_modp_roots (251) \\ before using this function. \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ { my_sqrtmodp (b, p, a, c, r) = b = polcoeff(lift(lift(b)),0) ; if (b == 0, \\ easy case r = Mod(0,a) ) ; \\ else if b != 0 if (b != 0 & table[b] != 0, \\ so that b is a square mod p r = Mod(table[b],a) ) ; \\ else if b is not a square \\ return w times the square root of -b/c where w^2 =-c. if (b != 0 & table[b] == 0, c = polcoeff(a,0) ; r = Mod (table[lift(-b/c)] * variable(a), a) ) ; r ; } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ my_factormodp: factor quickly over F_{p^2} using quadratic formula, etc. \\ \\ f is a polynomial over F_p all of whose roots lie in \\ F_{p^2}. This function returns a list of the factors \\ of f in terms of a root of the irreducible quadratic \\ polynoomial a \in Fp[w]. \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ { my_factormodp (f, p, a, factors, b,c,d,i,k,n, dis) = \\ first factor into quadratic and linear factors print("f = ",f); modpfactors = factormod (f,p) ; n = matsize(modpfactors)[1] ; \\ next count the number of quadratic factors, with multiplicity k = 0 ; for (i = 1, n, if (poldegree(modpfactors[i,1]) == 2, k++ ) ) ; factors = matrix( (n-k) + 2*k, 2,X,Y,0 ) ; d = 1 ; for (i=1,n, if (poldegree(modpfactors[i,1]) == 1, factors[d,1] = modpfactors[i,1] ; factors[d,2] = modpfactors[i,2] ; d++, \\ else \\ x^2 + bx + c. b = polcoeff(modpfactors[i,1],1) ; c = polcoeff(modpfactors[i,1],0) ; dis = my_sqrtmodp (b^2-4*c, p, a) ; factors[d,1] = x - (-b + dis ) / 2 ; factors[d,2] = modpfactors[i,2] ; d++ ; factors[d,1] = x - (-b - dis ) / 2 ; factors[d,2] = modpfactors[i,2] ; d++ ) ) ; factors ; } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ rootsfp2: find all roots of a polynomial over Fp^2 which lie in Fp^2. \\ + all roots ARE ASSUMED to lie in Fp^2!!!! If this is not true \\ the function will bomb. \\ + a is an irreducible quadratic polynomial mod p in terms of \\ which the solutions will be expressed. \\ for example, a = w^2+1, if -1 is not in Fp. \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ { rootsfp2 (f, p, a, d, g, h, factors, c, s, i, n, r) = \\ compute g = Frob_p(f). g = 0 ; for (d = 0, poldegree(f), g += x^d * polcoeff(f,d)^p ) ; \\ let h = f * Frob_p(f) \in Fp[x]. h = f * g ; factors = my_factormodp (f * g, p, a) ; \\ n is the number of distinct linear factors of f*g n = matsize (factors) [1] ; \\ r is the vector of zeros of f*g: \\ this is where the assumption that they are all in Fp^2 is used! r = vector (n, d, -polcoeff (factors[d,1],0) / polcoeff (factors[d,1],1) ); \\ make a vector whose entries are precisely the roots \\ of f*g which are also roots of f, with the correct multiplicities. \\ return a list of roots, counted with multiplicity c = 1; s = vector (poldegree(f), X,0) ; for (i = 1, n, while ( poldegree(f) > 0 & changevar (f, [ r[i] ]) == 0, s[c] = r[i] ; c = c + 1; f = f / (x - r[i])) ) ; vector (c-1, i, s[i]) ; } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ twoisog: given a supersingular curve E of j-invariant j in \\ characteristic p this function returns a vector \\ of all j-invariants of curves 2-isogeneous to E. { twoisog (j, N, a, phi2) = phi2 = x^3 + j^3 - x^2*j^2 + 1488 * x * j * (x + j) - 162000 * (j^2 + x^2) + 40773375 * j * x + 8748000000 * (x + j) - 157464000000000 ; rootsfp2 (phi2, N, a) ; } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ I typed in the following three modular polynomials from \\ Herrmann's paper. Be glad you didn't have to type these in! \\ ALSO, Obviously multiplying things out once for all would speed \\ this function up somewhat. { threeisog (j, N, a, phi3) = phi3 = x*(x+2^15*3*5^3)^3+j*(j+2^15*3*5^3)^3-x^3*j^3 +2^3*3^2*31*x^2*j^2*(x+j)-2^2*3^3*9907*x*j*(x^2+j^2) +2*3^4*13*193*6367*x^2*j^2 +2^16*3^5*5^3*17*263*x*j*(x+j)-2^31*5^6*22973*x*j ; rootsfp2 (phi3, N, a) ; } { fiveisog (j, N, b, phi5) = a = matrix(6,6,X,Y,0) ; a[1,1] = 2^90 * 3^18*5^3*11^9 ; a[2,1] = a[1,2] = 2^77 * 3^16 * 5^3 * 11^6 * 31 * 1193 ; a[2,2] = -2^62 * 3^13 * 11^3 * 26984268714163 ; a[3,1] = a[1,3] = 2^60 * 3^13 * 5^2 * 11^3 * 13^2 * 3167 * 204437 ; a[3,2] = a[2,3] = 2^47*3^10*5^4*53359 * 131896604713 ; a[3,3] = 2^30*3^8*5^4*7*13*1861*6854302120759 ; a[4,1] = a[1,4] = 2^48*3^9*5^2*31*1193*24203*2260451 ; a[4,2] = a[2,4] = -2^31*3^7*5^3*327828841654280269 ; a[4,3] = a[3,4] = 2^17*3^4*5^3*2311*2579*3400725958453 ; a[4,4] = -2^2*5^2*11*17*131*1061*169751677267033 ; a[5,1] = a[1,5] = 2^30*3^7*5*13^2*3167*204437 ; a[5,2] = a[2,5] = 2^20*3^4*5^3*12107359229837 ; a[5,3] = a[3,5] = 3*5^3*167*6117103549378223 ; a[5,4] = a[4,5] = 2^5*3*5^2*197*227*421*2387543 ; a[5,5] = 2^3*5^2*257*32412439 ; a[6,1] = a[1,6] = 2^17*3^4*5*31*1193 ; a[6,2] = a[2,6] = -2*3*5^2*1644556073 ; a[6,3] = a[3,6] = 2^5*5^2*13*195053 ; a[6,4] = a[4,6] = -2^2*3^2*5*131*193 ; a[6,5] = a[5,6] = 2^3*3*5*31 ; a[6,6] = -1 ; phi5 = x^(5+1) + j^(5+1) + sum (v=1,6, sum (u=1,6, a[v,u]*x^(v-1)*j^(u-1))) ; rootsfp2 (phi5, N, b) ; } { sevenisog(j, N, b, phi7) = a = matrix(8,8,X,Y,0) ; a[1,1] = 0 ; a[1,2] = a[2,1] = 0 ; a[2,2] = 2^91*3^27*5^18*11*13*17^9 ; a[1,3] = a[3,1] = 2^90*3^27*5^18*7^3*17^9 ; a[2,3] = a[3,2] = -2^76*3^25*5^15*7^2*17^7*947*22541 ; a[3,3] = -2^61*3^20*5^12*7^2*17^3*1644361*60057292681 ; a[1,4] = a[4,1] = 2^76*3^25*5^15*7^3*17^6*31*26891 ; a[2,4] = a[4,2] = -2^61*3^19*5^12*7^2*17^3*19*487*88980809456419 ; a[3,4] = a[4,3] = 2^46*3^16*5^9*7^2*409*2633*231491957605001610911 ; a[4,4] = -2^31*3^10*5^6*7^2*23131*216217*11116146635649891498353 ; a[1,5] = a[5,1] = 2^60*3^19*5^12*7^2*17^3*397*1323331291097 ; a[2,5] = a[5,2] = 2^46*3^17*5^9*7^2*13*12983*3769379869638077087 ; a[3,5] = a[5,3] = 2^31*3^11*5^6*7^2*17*62349740297426529782049295279 ; a[4,5] = a[5,4] = 2^16*3^7*5^4*7^2*37*43*661*3893394856539704079067727101 ; a[5,5] = 2*5*7^2*197*912019631831096138476499139089037899 ; a[1,6] = a[6,1] = 2^47*3^16*5^9*7^2*13*31*26891*683503*9854261 ; a[2,6] = a[6,2] = -2^33*3^11*5^6*7^2*3697*9447061867111661492633 ; a[3,6] = a[6,3] = 2^19*3^9*5^3*7^2*178299075699438778621099394269 ; a[4,6] = a[6,4] = -2^3*7^2*1523*1510253179371566698174112073613709 ; a[5,6] = a[6,5] = 2^5*3^4*5^2*7^2*11*113*41659*85554840812719128607 ; a[6,6] = -2^2*3^2*7^2*10374612889856513538191507 ; a[1,7] = a[7,1] = 2^30*3^10*5^6*7*397*1323331291097 ; a[2,7] = a[7,2] = 2^17*3^7*5^3*7^2*59*10020909155496489683 ; a[3,7] = a[7,3] = 2^2*7^2*29*1291*6221*234068219582309095597 ; a[4,7] = a[7,4] = 2^4*3^5*5*7^2*197*227*2113*2087377*85827811 ; a[5,7] = a[7,5] = 2^3*3*7^2*302587*12536290128459761 ; a[6,7] = a[7,6] = 2^4*7^3*3259*9901340156731 ; a[7,7] = 3^2*7^2*13*67*97*8389943 ; a[1,8] = a[8,1] = 2^16*3^7*5^3*7*31*26891 ; a[2,8] = a[8,2] = -2^3*7^2*13*2129*5107*631559 ; a[3,8] = a[8,3] = 2^4*3^4*7^2*43*1801*146437 ; a[4,8] = a[8,4] = -2*3*7^2*13*1067425727 ; a[5,8] = a[8,5] = 2^5*5^2*7^2*11*43*509 ; a[6,8] = a[8,6] = -2^2*3^3*7*13553 ; a[7,8] = a[8,7] = 2^3*3*7*31 ; a[8,8] = -1 ; \\ done typing! phi7 = x^(7+1) + j^(7+1) + sum (v=1,8, sum (u=1,8, a[v,u]*x^(v-1)*j^(u-1))) ; rootsfp2 (phi7, N, b) ; } { pisog (j, p, N, a, r) = if (p == 2, r = twoisog (j, N, a)) ; if (p == 3, r = threeisog (j, N, a)) ; if (p == 5, r = fiveisog (j, N, a)) ; if (p == 7, r = sevenisog (j, N, a)) ; if (p != 2 & p != 3 & p!=5 & p!=7, print ("ERROR!!, ", p,"-isogenies not yet implemented.")) ; r ; } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ returns the number of supersingular elliptic curves in characteristic p >= 5: { count_supersingular_j (p, n) = n = floor(p/12) ; if (p%12 == 5 | p%12 == 7, n++ ) ; if (p%12 == 11, n += 2 ) ; n ; } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ find a in vector v: if a is already in v then returns the position \\ of a in v, otherwise returns (length of v) + 1. { find (a, v, n) = n = 1; while (v[n] != a & n < length(v), n++ ) ; if (v[n] == a, n, \\ else length(v) + 1 ) } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ given a prime level N, this function returns a matrix representing \\ the action of the hecke operator T_p on M_2(Gamma_0(N)). \\ \\ j = a supersingular j-invariant in characteristic N. { tpmatrix (p, N, j, a, i, r, M, pos, last, basis ) = j = supersingular_j (N) ; \\ Initialize the mod N arithmetic by setting \\ up the global table of roots mod N. table = table_of_modp_roots(N) ; \\ ok, mixing "p" and "N" is confusing! \\ Compute a monic irreducible polynomial of degree 2 over F_N. a = ffinit ( N, 2, w ) ; \\ Compute the dimension of M_2. dim = count_supersingular_j (N) ; \\ basis will eventually contain a basis of supersingular curves. basis = vector (dim, X, 0) ; \\ M will contain a matrix representing the action \\ of the hecke operator T_2 on M_2(Gamma_0(N)). M = matrix (dim, dim, X, Y, 0) ; \\ Initialize our basis (so far) with the initial \\ supersingular j-invariant. We'll find a basis \\ for all of M_2 by hitting what we've got so \\ far by T_p. That this works uses an unpublished \\ remark of Serre, the proof of which someone needs \\ to find! (See the discussionf on page 223 of Mestre.) basis[1] = j ; last = 1 ; \\ For each basis vector, compute the column of the matrix of T_p. if (verbose, print (dim," supersingular j's in char ",N,": "); print1 ("j = ") ) ; for (i = 1, dim, if (verbose, print1 (lift (lift (basis[i]))); if (i < dim, print1 (", "), print("")) ) ; \\ Compute all curves p-isogeneous to basis[i], with multiplicity. r = pisog (basis[i], p, N, a) ; \\ For each element of r: \\ 1. If we haven't seen it before add it to our basis, \\ otherwise find it in our basis. \\ 2. Add it to the column of the matrix M of T_2 \\ which we are constructing. for (j = 1, length(r), pos = find (r[j], basis) ; if (pos > dim | pos>last, last ++ ; basis [last] = r[j] ; M [last, i] ++, \\ else M [pos, i] ++ ) ) ; ) ; [M, basis, a] ; } tppoly (p,N) = CharPolyAdjoint ( tpmatrix (p,N) ) ; \\ Some useful shorthands. t2mat (N) = tpmatrix(2,N) ; t2poly (N) = tppoly(2,N) ; t3mat (N) = tpmatrix(3,N) ; t3poly(N) = tppoly(3,N) ; t5mat (N) = tpmatrix(5,N) ; t5poly(N) = tppoly(5,N) ; t7mat (N) = tpmatrix(7,N) ; t7poly(N) = tppoly(7,N) ; \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ CharPolyPS: \\ compute the characteristic polynomial of a LARGE matrix M \\ by first computing the traces of the powers of M \\ tr(M), tr(M^2), tr(M^3), ..., \\ and using the fact that one can recover the elementary symmetric \\ functions of n numbers from the n powersums of those numbers. \\ (see Janos' paper.) \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ { PowVecToChar (PowVec, f, g, n) = n = length (PowVec) ; f = Polrev (PowVec, x) + O(x^n) ; g = truncate (exp (intformal (-f,x) ) ) ; x^n * subst (g, x, 1/x) } { CharPolyPS (M, N, PowVec, k) = N = M ; PowVec = vector (length(M), k, 0) ; for (k = 1, length (M), PowVec [k] = trace (N) ; N *= M ) ; PowVecToChar (PowVec) ; } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ CharPolyAdjoint: \\ Characteristic polynomial using the "adjoint matrix method". \\ See [Cohen, Algorithms in Algebraic Number Theory, Algorithm 2.2.7] { CharPolyAdjoint (M, C, a, i, n, f) = n = length (M) ; C = matid (n) ; f = x^n ; if (verbose, print1("CharPoly computation (",n," steps): ") ) ; for (i = 1, n, if (verbose, print1(i) ; if (i