#! /usr/bin/env python
# This Python program computes the coefficients of the modular
# invariant function (j) by using the action of the Hecke operators.
# David A. Madore - 2007-07-31 - Public Domain
# The basic relation is this: if
# j(q) = 1/q + c1 q + c2 q^2 + c3 q^3 + ...
# (we use the normalization 0 for the constant term), there exists a
# unique, easily computed, monic polynomial of degree n (the n-th
# Faber polynomial for j), F_n, such that F_n(j) starts like 1/q^n +
# terms of order at least 1 in q. Now modular theory tells us that
# F_n(j) actually gives (n times) the n-th Hecke operator T_n acting
# on j, i.e., something like
# 2 T_2(j) = 1/q^2 + 2 c2 q + (2c4+c1) q^2 + 2 c6 q^3 + (2c8+c2) q^4 + ...
# 3 T_3(j) = 1/q^3 + 3 c3 q + 3 c6 q^2 + (3c9+c1) q^3 + 3c12 q^4 + ...
# 4 T_4(j) = 1/q^4 + 4 c4 q + (4c8+2c2) q^2 + 4c12 q^3 + (4c16+2c4+c1) q^4 + ...
# etc. This relation F_n(j) = n T_n(j) can be used both ways:
# computing F_n(j) allows us to compute higher (divisible)
# coefficients from lower ones, but in the other direction, computing
# F_2(j) with the highest unknown coefficient as an indeterminate
# allows us to compute the latter from others.
# Set this to True to print coefficients as soon as they are computed,
# False to print them consecutively (e.g., do not print c(18) until
# c(17) has been computed).
earlyprint = False
# We bootstrap with a few known coefficients:
bootcoefs = [0, 196884, 21493760, 864299970, 20245856256, 333202640600,
4252023300096, 44656994071935, 401490886656000, 3176440229784420,
22567393309593600, 146211911499519294, 874313719685775360,
4872010111798142520, 25497827389410525184, 126142916465781843075,
593121772421445058560]
# lencomplete is the *first unknown* coefficient.
lencomplete = len(bootcoefs)
# coefs is the main storage of j coefficients; it is a dictionary so
# we don't have to store things contiguously.
coefs = dict()
for i in range(lencomplete):
coefs[i] = bootcoefs[i]
if i>0:
print "%d\t%d"%(i, coefs[i])
# Now we need an object class for doing simple operations on power
# series (actually Laurent series: index is the valuation):
class PowerSeries:
def __init__(self, index, coefs):
self.index = index
self.coefs = coefs
def precis(self): # Return the (big-O) precision we have
return self.index + len(self.coefs)
def multiply(s1, s2): # Return a _new_ series by multiplying s1 and s2
index = s1.index + s2.index
length = min(len(s1.coefs),len(s2.coefs))
coefs = []
for i in range(length):
x = 0
for j in range(i+1):
x += s1.coefs[j]*s2.coefs[i-j]
coefs.append(x)
return PowerSeries(index, coefs)
def addMonomial(self, c, k): # _Change_ a series by adding a monomial
if k < self.index:
self.coefs[:0] = (self.index - k) * [0]
self.index = k
self.coefs[k-self.index] += c
def addMonomialTimes(self, c, k, f): # Add a monomial times f
if k + f.precis() < self.index:
raise Exception("PowerSeries defined nowhere")
elif k + f.precis() < self.precis():
del self.coefs[(f.precis() - self.index):]
if k + f.index < self.index:
self.coefs[:0] = (self.index - k - f.index) * [0]
self.index = k + f.index
for i in range(len(self.coefs)):
if self.index + i >= k + f.index:
self.coefs[i] += c*f.coefs[i+self.index-k-f.index]
def coef(self, k): # Return a given coefficient
if k < self.index:
return 0
elif k < self.precis():
return self.coefs[k-self.index]
else:
raise Exception("PowerSeries insufficiently precise")
def strictDivisors(n): # Return the list of divisors of n between 1 and n-1
l = []
for i in range(1,n):
if (n%i)==0:
l.append(i)
return l
# An exception thrown in various places to signify we can't compute
# something yet.
class MissingCoef:
def __init__(self, value):
self.value = value
def __str__(self):
return repr(self.value)
# The main part!
while True:
# Start by forming a series with whatever contiguous coefficients we have.
jcoefs = []
for k in range(lencomplete):
jcoefs.append(coefs[k])
j = PowerSeries(0, jcoefs)
j.addMonomial(1, -1)
# Now compute the first Faber polynomials of j:
jFaber = [None, j]
for n in range(2, min(lencomplete,7)):
# ((Here 7 is a heuristic, meaning we compute the first 6 Faber
# polynomials: any value at least equal to 5 should work, but
# higher values are interesting only if you wish to earlyprint
# high coefficients.))
jn = PowerSeries.multiply(jFaber[n-1], j)
for k in range(n-2, 0, -1):
jn.addMonomialTimes(-jn.coef(-k),0,jFaber[k])
jn.addMonomial(-jn.coef(0),0)
jFaber.append(jn)
# At this point, jn is the n-th Faber polynomial of j.
ld = strictDivisors(n)
for k in range(1, jn.precis()):
if not coefs.has_key(k*n):
# Compute a coefficient from the action of the n-th
# Hecke operator (with just n*coefs[k*n], the one we
# will deduce, missing from the sum):
try:
v = 0
for d in ld:
a = n/d
if (k%a)==0:
kk = (k/a)*d
if coefs.has_key(kk):
v += (n/a)*coefs[kk]
else:
raise MissingCoef(kk)
w = jn.coef(k)-v
if (w%n)!=0:
raise Exception("Divisibility check failed!")
w /= n
coefs[k*n] = w
if earlyprint:
print "%d\t%d"%(k*n, w)
except MissingCoef, kk:
pass # (Actually never happens...)
# Now try the other way around: deduce some lower coefficients
# from the higher ones (known through the Hecke operators). We
# only use T_2 here, so we only deal with F_2(j), which is
# essentially j^2 (up to a constant -2*196884 we don't care about
# since we're interested only in one, higher, coefficient).
try:
while True:
# See if lencomplete can be increased.
while coefs.has_key(lencomplete):
if not earlyprint:
print "%d\t%d"%(lencomplete, coefs[lencomplete])
lencomplete += 1
# Try to compute the first unknown coefficient
# (lencomplete) by computing the previous coefficient in
# 2 T_2(j) and equating.
k = lencomplete-1
if not coefs.has_key(k*2):
raise MissingCoef(k*2)
v = 2*coefs[k*2]
if (k%2)==0:
if not coefs.has_key(k/2):
raise MissingCoef(k/2)
v += coefs[k/2]
# At this point, v is coefficient k of j^2, computed from
# the Hecke operators. Now we can deduce coefficient k+1
# of j from this:
for i in range(1, k):
v -= coefs[i] * coefs[k-i]
if (v%2)!=0:
raise Exception("Evenness check failed!")
v /= 2
coefs[k+1] = v
if earlyprint:
print "%d\t%d"%(k+1, v)
except MissingCoef, kk:
pass