######################################################################
# Author: morain@lix.polytechnique.fr
# Version: 1.0
# Date: 2023/07/24
######################################################################

# This script is devoted to the computation and verification of several
# identities related to CCR polynomials using the notations of the preprint.

# one ring to rule them all
# s is used for U_ell^e and is s.t. s*(ell+1) = 0 mod 24.
R.<ell,s,E2,E4,E6,sigma,E4t,E6t,d4,d6,dS,dS4,dS6,d46,f,df,df4,df6>=PolynomialRing(Rationals(),18)

########## The CCR case

# returns ell^-4 * dS^-1 * (-12*ell*E4^2*d6 + ...)
def check_E4t():
#    R.<ell,E2,E4,E6,sigma,d4,d6,s,dS,E4t>=PolynomialRing(Rationals(),10)
    E4p=(E2*E4 - E6)/3
    E6p=(E2*E6-E4^2)/2
    E2p=(E2^2-E4)/12
    E2t=(E2+2*sigma/ell)/ell
    sigp=ell/24*(4*sigma^2/ell^2+4*sigma/ell*E2-(ell^2*E4t-E4))
    print("sigma'=", sigp)
    tmp=sigp*dS+E4p*d4+E6p*d6
    tmp=tmp.numerator()
    print("degree(tmp, E2)=", tmp.degree(E2))
    # check that coeff of E2 is zero
    c1=tmp.coefficient({E2:1})
    # is a multiple of (2*E4*d4 + 3*E6*d6 + f*df), hence 0
    print("c1=", c1.factor())
    # find sigma as a root of constant coefficient
    e4t=tmp.coefficient({E2:0})
    e4t=-e4t.coefficient({E4t:0})/e4t.coefficient({E4t:1})
    # sig contains the value of sigma
    return e4t.factor()

# returns
# ell^-6 * dS^-3 * sigma^-1 * E6^-1 * E4^-1 * (-18*ell^3*E4^5*E6*d6^2*dS+...)
def check_E6t():
    e4t=check_E4t()
    E4p=(E2*E4 - E6)/3
    E6p=(E2*E6-E4^2)/2
    E2p=(E2^2-E4)/12
    E2t=(E2+2*sigma/ell)/ell
    sigp=ell*(4*sigma^2/ell^2+4*sigma/ell*E2-(ell^2*e4t-E4))/24
    print("sigma'=", sigp)
    # more derivatives
    E4pp=1/3*(E2p*E4+E2*E4p-E6p)
    E6pp=1/2*(E2p*E6+E2*E6p-2*E4*E4p)
    # crucial values
    E4tp=1/3*(E2t*e4t-E6t)
    E2tp=(E2t^2-e4t)/12
    E2pp=1/12*(2*E2*E2p-E4p)
    E2tpp=1/12*(2*E2t*E2tp-E4tp)
    sigpp=ell*(ell^3*E2tpp-E2pp)/2
    # inject diagonal derivatives
    dSS = (ell*dS    -2*E4*dS4 -3*E6*dS6)/sigma
    d44 = ((ell-1)*d4-sigma*dS4-3*E6*d46)/(2*E4)
    d66 = ((ell-2)*d6-sigma*dS6-2*E4*d46)/(3*E6)
    # starting point
    tmp=      sigpp*dS+sigp*(sigp*dSS+E4p*dS4+E6p*dS6)
    tmp=tmp + E4pp*d4+E4p*(sigp*dS4+E4p*d44+E6p*d46)
    tmp=tmp + E6pp*d6+E6p*(sigp*dS6+E4p*d46+E6p*d66)
    tmp=tmp.numerator()
    c2=tmp.coefficient({E2:2})
    print("E6t.c2=", c2.factor())
    c1=tmp.coefficient({E2:1})
    print("E6t.c1=", c1)
    c0=tmp.coefficient({E2:0})
    e6t=-c0.coefficient({E6t:0})/c0.coefficient({E6t:1})
    return e6t.factor()

########## The case of U_ell^e (Atkin's variant)
def checke_sigma():
#    R.<ell,E2,E4,E6,sigma,d4,d6,f,df>=PolynomialRing(Rationals(),9)
    t = 24/s
    tmp=4*E4*d4+6*E6*d6+s*f*df
    E4p=(E2*E4 - E6)/3
    E6p=(E2*E6-E4^2)/2
    E2p=(E2^2-E4)/12
    E2t=(E2+2*sigma/ell)/ell
    fp=f/t*(ell*E2t+E2)
    tmp=fp*df+E4p*d4+E6p*d6
    tmp=tmp.numerator()
    # check that coeff of E2 is zero
    c1=tmp.coefficient({E2:1})
    # is a multiple of (4*E4*d4 + 6*E6*d6 + s*f*df), hence 0
    print("c1=", c1.factor())
    # find sigma as a root of constant coefficient
    sig=tmp.coefficient({E2:0})
    sig=-sig.coefficient({sigma:0})/sig.coefficient({sigma:1})
    # sig contains the value of sigma
    return sig.factor()

# returns
# (-1) * df^-3 * f^-2 * s^-2 * ell^-2 * E6^-1 * E4^-1
#       * (-12*ell*s^2*E4^5*d6*f*df^2+...)
def checke_E4t():
    sig=checke_sigma()
    t = 24/s
    E4p=(E2*E4 - E6)/3
    E6p=(E2*E6-E4^2)/2
    E2p=(E2^2-E4)/12
    E2t=(E2+2*sig/ell)/ell
    fp=f/t*(ell*E2t+E2)
    fpp=f/t*(ell*E2t+E2)^2+f/12*(ell^2*(E2t^2-E4t)+(E2^2-E4))
    fpp=fpp/t
    E4pp=1/3*(E2p*E4+E2*E4p-E6p)
    E6pp=1/2*(E2p*E6+E2*E6p-2*E4*E4p)
    # inject diagonal derivatives
    dff = (        s*ell*df-4*E4*df4 -6*E6*df6)/f/s
    d44 = ((s*(ell+1)-4)*d4-s*f*df4-6*E6*d46)/(4*E4)
    d66 = ((s*(ell+1)-6)*d6-s*f*df6-4*E4*d46)/(6*E6)
    tmp=       fpp*df+ fp*(fp*dff+E4p*df4+E6p*df6)
    tmp=tmp + E4pp*d4+E4p*(fp*df4+E4p*d44+E6p*d46)
    tmp=tmp + E6pp*d6+E6p*(fp*df6+E4p*d46+E6p*d66)
    tmp=tmp.numerator()
    print("degree(tmp, E2)=", tmp.degree(E2))
    c2=tmp.coefficient({E2:2})
    print("E4t.c2=", c2.factor())
    c1=tmp.coefficient({E2:1})
    print("E4t.c1=", c1)
    c0=tmp.coefficient({E2:0})
    e4t=-c0.coefficient({E4t:0})/c0.coefficient({E4t:1})
    return e4t.factor()