###

This is outdated!!

Please use the file eta_products.py instead

###

def EtaProduct(L,prec,silent=False):
    #L is a list of lists [[n1,exp1],...,[nk,expk]], prec is an integer. Produces the eta product (quotient) eta(n1 z)^exp1 * ... * exp(nk z)^expk to precision prec
    R.<q> = PowerSeriesRing(QQ, default_prec = prec)
    f = 1
    offset = 0
    weight = 0
    levels = []
    for k in range(len(L)):
        [n,exponent] = L[k]
        levels.append(n)
        offset = offset + exponent*n / 24
        weight = weight + exponent/2
        sign = sgn(exponent)
        absexponent = abs(exponent)
        exp_remainder = absexponent % 3
        jacobi_exponent = (absexponent - exp_remainder) / 3
        g = int(prec/n) + 2
        jacobi_bound = floor((sqrt(1 + 8*g) - 1) / 2) + 1
        jacobi_kfactor = 1
        for l in range(1,jacobi_bound):
            jacobi_kfactor = jacobi_kfactor + (-1)^l * (2*l+1)*q^ZZ(n*l*(l+1)/2)
        f = f * jacobi_kfactor^ZZ(sign*jacobi_exponent)
        if exp_remainder != 0:
            bound = floor(sqrt(1 + 24*g)/6) + 2
            kfactor = 1
            for l in range(1,bound):
                kfactor = kfactor + (-1)^l * q^ZZ(n*(3*l^2 + l)/2) + (-1)^l * q^ZZ(n*(3*l^2 - l)/2)
            f = f * kfactor^ZZ(sign*exp_remainder)
    f = f + O(q^prec)
    frac_offset = frac(offset)
    f = f * q^(ZZ(offset - frac_offset))
    if silent == False:
        print "Weight", weight
        print "Level", lcm(levels)
        if frac_offset != 0:
            print "Add", frac_offset, "to all exponents"
    return f

def IsItEtaProduct(f):
    #f should be a Laurent series. If exponents are offset by some fraction then remove the offset first. For example use IsItEtaProduct(1+q+q^3+q^6+q^10+O(q^15)) to find eta(2z)^2 / eta(z)
    prec = f.common_prec(f)
    fval = f.valuation()
    if prec == Infinity:
        prec = f.degree() + 1
    R.<q> = PowerSeriesRing(QQ,prec)
    g = f * q^(-fval)
    kludge_f = R(0) #command R(g) seems to be bugged if f has negative exponents
    for i in range(prec - f.valuation()):
        kludge_f = kludge_f + q^i * g[i]
    f = kludge_f
    prefactor = f[0]
    f = f / prefactor
    L = []
    logf = log(f)
    while not logf == 0:
        N = logf.valuation()
        if N >= prec - fval:
            logf = 0
            break
        if N > prec/2 and len(L) > prec/2: #arbitrary cutoff, change as needed
            print "Probably not."
            return
        exponent = -logf[N]
        if abs(exponent) > 20*prec: #arbitrary cutoff, change as needed
            print "Probably not."
            return
        if not exponent.is_integer():
            print "Definitely not."
            return
        L.append([N,exponent])
        g = EtaProduct([[N,exponent]],prec,silent=True)
        logf = logf - log(R(g/q^(g.valuation())))
    output = "It looks like "
    if prefactor != 1:
        output = output + str(prefactor)
        if len(L) > 0:
            output = output + "*"
    if len(L) == 0 and prefactor == 1:
        output = output + "1"
    for k in range(len(L)):
        [N,exponent] = L[k]
        output = output + "eta("
        if N != 1:
            output = output + str(N)
        output = output + "z)"
        if exponent < 0:
            output = output + "^("
        if exponent > 1:
            output = output + "^"
        if exponent != 1:
            output = output + str(exponent)
        if exponent < 0:
            output = output + ")"
        if k < len(L) - 1:
            output = output + "*"
    output = output + "."
    print output

#EXAMPLE: compute the eta quotient eta(tau)^3 / eta(4*tau) to precision q^(200).
EtaProduct([[1,3],[4,-1]],200)
#include silent=False to turn off weight,level,exponent offset info

#EXAMPLE: is 1 + 3q + 9q^2 + 21q^3 + 48q^4 + 99q^5 + 198q^6 + 375q^7 + 693q^8 + ... an eta product?
print "EXAMPLE 2"
R.<q> = PowerSeriesRing(QQ)
IsItEtaProduct(1 + 3*q + 9*q^2 + 21*q^3 + 48*q^4 + 99*q^5 + 198*q^6 + 375*q^7 + 693*q^8 + O(q^9))

#SANITY CHECK
print "EXAMPLE 3"
IsItEtaProduct(EtaProduct([[1,5],[2,4],[3,3],[4,-1]],1000,silent=True))