#! /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 <david.madore@ens.fr> - 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