######################################################################
# Author: morain@lix.polytechnique.fr
# Version: 2.1
# Date: 2024/02/20
######################################################################

# This script is devoted to the computation and verification of several
# identities related to modular polynomials using the notations of the 
# various preprints.

R.<ell,r,s,E2,E4,E6,sigma,E4t,E6t,d4,d6,dS,dS4,dS6,d46,X,dX,dX4,dX6>=PolynomialRing(Rationals(),19)

########## 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 + s*X*dX), 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 Fr_{ell, r, s}

# one ring to rule them all
# r+s*ell = r*ell+s = 0 mod 24; w = (r+s)/2 should be even.

# returns (2) * dX^-1 * X^-1 * s^-1 * ell * (3*E4^2*d6 + 2*E6*d4)
def double_eta_check_sigma():
#    R.<ell,r,s,E2,E4,E6,sigma,d4,d6,X,dX>=PolynomialRing(Rationals(),11)
#    tmp=w*X*dX+4*E4*d4+6*E6*d6 # quid?
    # w = (r+s)/2 is the weight
    E4p=(E2*E4 - E6)/3
    E6p=(E2*E6-E4^2)/2
    E2p=(E2^2-E4)/12
    E2t=(E2+2*sigma/ell)/ell
    Xp=(X/24)*(r*E2+s*ell*E2t)
    tmp=Xp*dX+E4p*d4+E6p*d6
    tmp=tmp.numerator()
    # check that coeff of E2 is zero
    c1=tmp.coefficient({E2:1})
    # is a multiple of (w*X*dX + 4*E4*d4 + 6*E6*d6), 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) * dX^-3 * X^-2 * s^-2 * ell^-2 * E6^-1 * E4^-1 * (r+s)^(-1) *
#       * (-6*ell*r^2*s*E4^5*d6*X*dX^2+...)
def double_eta_check_E4t():
    w = (r+s)/2
    sig=double_eta_check_sigma()
    E4p=(E2*E4 - E6)/3
    E6p=(E2*E6-E4^2)/2
    E2p=(E2^2-E4)/12
    E2t=(E2+2*sig/ell)/ell
    Xp=(X/24)*(r*E2+s*ell*E2t)
    Xpp=(X/24^2)*((s*ell*E2t+r*E2)^2+2*(s*ell^2*(E2t^2-E4t)+r*(E2^2-E4)))
    E4pp=1/3*(E2p*E4+E2*E4p-E6p)
    E6pp=1/2*(E2p*E6+E2*E6p-2*E4*E4p)
    # inject diagonal derivatives
    dXX = (        w*ell*dX-4*E4*dX4 -6*E6*dX6)/X/w
    d44 = ((w*(ell+1)-4)*d4-w*X*dX4-6*E6*d46)/(4*E4)
    d66 = ((w*(ell+1)-6)*d6-w*X*dX6-4*E4*d46)/(6*E6)
    tmp=       Xpp*dX+ Xp*(Xp*dXX+E4p*dX4+E6p*dX6)
    tmp=tmp + E4pp*d4+E4p*(Xp*dX4+E4p*d44+E6p*d46)
    tmp=tmp + E6pp*d6+E6p*(Xp*dX6+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()