from sympy import *

print("This computes the matrix for multiplication with c_1")
print("acting on the four-dimensional eigenspace in")
print("the rational elliptic surface. It also determines")
print("the characteristic polynomial, and checks that")
print("against the predictions from the ODE")

q = symbols('q')
exppsi = 1 - 4*q**3 + 2*q**6 + 8*q**9 - 5*q**12 - 4*q**15 - 10*q**18 + 8*q**21 + 9*q**24 + 14*q**30 - 16*q**33 - 10*q**36 + O(q**39)
derivpsi = -(12*q**2 + 36*q**5 + 48*q**8 + 84*q**11 + 72*q**14 + 144*q**17 + 96*q**20 + 180*q**23 + 156*q**26 + 216*q**29 + 144*q**32 + 336*q**35 + O(q**38))
z = 9*q + 126*q**4 + 936*q**7 + 5220*q**10 + 23940*q**13 + 96012*q**16 + 346752*q**19 + 1153296*q**22 + 3582810*q**25 + 10511730*q**28 + 29360448*q**31 + O(q**34)

"""
the basis is 1, M, \delta E, and 9point
the a13,a23,a33,a43 are all multiplied by q, to avoid 1/q issues 
"""
a11 = 0
a12 = 1
a13 = 0
a14 = 0

a21 = series(4*z,q,n=34)
a22 = series(-derivpsi/exppsi,q,n=38)
a23 = series(1/exppsi,q,n=39)
a24 = 0

a31 = series(2*q*diff(z,q),q,n=34)
a32 = series((q*derivpsi**2-q*diff(derivpsi,q))/exppsi,q,n=35)
a33 = series((-q*derivpsi-1)/exppsi,q,n=36)
a34 = 1

a41 = series((3*diff(z,q)+4*q*diff(derivpsi,q)*z - 4*q*derivpsi**2*z +q*diff(diff(z,q),q))/exppsi,q,n=33)
prea42 = (2*derivpsi**2 + 2*q*derivpsi*diff(derivpsi,q)-2*diff(derivpsi,q)-q*diff(diff(derivpsi,q),q))/(exppsi**2)
a42 = series(prea42 - 2*q*diff(z,q),q,n=33)
a43 = series(4*z*q,q,n=35)
a44 = 0

t = symbols('t')
m11 = t-a11
m12 = -a12
m13 = -a13
m14 = -a14

m21 = -a21
m22 = t-a22
m23 = -a23
m24 = -a24

m31 = -a31
m32 = -a32
m33 = t*q-a33
m34 = -a34

m41 = -a41
m42 = -a42
m43 = -a43
m44 = t-a44

def det2(m):
  d = m[0,0]*m[1,1]-m[1,0]*m[0,1]
  return(d)

def det3(m):
  d1 = m[1,1]*det2(Matrix([[m[2,2],m[2,0]],[m[0,2],m[0,0]]]))
  d2 = m[1,2]*det2(Matrix([[m[2,1],m[2,0]],[m[0,1],m[0,0]]]))
  d3 = m[1,0]*det2(Matrix([[m[2,1],m[2,2]],[m[0,1],m[0,2]]]))
  return(d1-d2+d3)

def det4(m):
  d1 = m[1,1]*det3(Matrix([[m[2,2],m[2,3],m[2,0]],[m[3,2],m[3,3],m[3,0]],[m[0,2],m[0,3],m[0,0]]]))
  d2 = m[1,2]*det3(Matrix([[m[2,1],m[2,3],m[2,0]],[m[3,1],m[3,3],m[3,0]],[m[0,1],m[0,3],m[0,0]]]))
  d3 = m[1,3]*det3(Matrix([[m[2,1],m[2,2],m[2,0]],[m[3,1],m[3,2],m[3,0]],[m[0,1],m[0,2],m[0,0]]]))
  d4 = m[1,0]*det3(Matrix([[m[2,1],m[2,2],m[2,3]],[m[3,1],m[3,2],m[3,3]],[m[0,1],m[0,2],m[0,3]]]))
  return(d1-d2+d3-d4)
                      
charmatrix = Matrix([[m11,m12,m13,m14],[m21,m22,m23,m24],[m31,m32,m33,m34],[m41,m42,m43,m44]])
print(charmatrix)
print()
charpoly = expand(det4(charmatrix)/q,t)
print("Characteristic polynomial:")
print(charpoly)
print()

theta11 = 1 + 6*q**3 + 6*q**9 + 6*q**12 + 12*q**21 + 6*q**27 + O(q**36)
theta21 = -18*q**2 - 72*q**5 - 306*q**8 - 1008*q**11 - 2934*q**14 - 7704*q**17 - 19134*q**20 - 44496*q**23 - 99270*q**26 - 212256*q**29 - 439272*q**32 + O(q**35)
theta12 = -q - q**4 - 2*q**7 - 2*q**13 - q**16 - 2*q**19 - q**25 - 2*q**28 - 2*q**31 + O(q**37)
theta22 = 1 + 8*q**3 + 44*q**6 + 152*q**9 + 487*q**12 + 1352*q**15 + 3518*q**18 + 8480*q**21 + 19503*q**24 + 42768*q**27 + 90530*q**30 + 185192*q**33 + O(q**36)

invproduct = series(1/(theta11**3+27*theta12**3),q,n=36)
c2 = expand(3*(theta11**2*theta21+27*theta12**2*theta22)*invproduct)
c1 = expand(3*(theta11*theta21**2+27*theta12*theta22**2)*invproduct)
c0 = expand((theta21**3+27*theta22**3)*invproduct)
cubicp = t**3 - c2*t**2 + c1*t - c0
lambdainfty = series(theta22/(theta12/q),q,n=36)/q
quarticp = expand((t-lambdainfty)*cubicp)

print("The prediction for the characteristic polynomial is:")
print(quarticp)
print()

print("The difference between the two computations is:")
print(charpoly-quarticp)