#Program to perform the calculations in the Appendix of

#G. Oberdieck, D. Piyaratne, Y. Toda, Donaldson-Thomas invariants of abelian threefolds and Bridgeland stability conditions, arXiv:1808.02735.

#Usage:
#Run sage in same folder as the file. Then type: load("Exterior_Algebra.py").


import re
g=3
n=2*g


#Define a polynomial ring.
def st(L):
  return "".join([str(x) for x in L])

S=[s for s in Subsets(range(n)) if is_even(len(s))] #Even wedges
Sstr=[ "".join(re.findall(r'\d+', str(s))) for s in S]

P=PolynomialRing(QQ, ['a%s' % s for s in Sstr])
a = { s : P('a%s' % s) for s in Sstr }

#Matrices for the generator of so(2*n)
E=ExteriorAlgebra(QQ, n)
e=E.gens()
def Wedge(a,b):
  #Wedge with e_i e_j
  return matrix([[(e[a]*e[b]*E.basis()[sj]).coefficient(si) for sj in S] for si in S])

def AntiWedge(a,b):
  #apply \partial/partial x_i \partial/partial x_j
  return matrix([[ ((E.basis()[sj]).interior_product(e[b]).interior_product(e[a])).coefficient(si) for sj in S] for si in S])

def WedgeAntiWedge(a,b):
  #apply x_i \partial/partial x_j - 1/2 delta_{ij}
  M32=MatrixSpace(QQ,32,32)
  return matrix([[( e[a]*((E.basis()[sj]).interior_product(e[b])) ).coefficient(si) for sj in S] for si in S]) - (Rational((1,2))*M32.identity_matrix() if a==b else M32.zero_matrix())
  

print "Create list of monomials whose span contains Delta."
Delta_variable_list=[]
for p in xmrange([32]*4):
  I=[S[x] for x in p]
  if any( len([T for T in I if i in T])!=2 for i in range(6)):
    continue
  Delta_variable_list.append( prod(P('a%s' % st(T)) for T in I) )
Delta_variable_list = list(set(Delta_variable_list))
N=len(Delta_variable_list) #256



#COMPUTATION OF EQUATIONS FOR Delta BELOW
@cached_function
def Test_variable_list(r1,r2):
  print "Create list of monomials whose span contains <dDelta,Wegdge(%s,%s)*gamma>." % (r1,r2)
  TL=[]
  for p in xmrange([32]*4):
    I=[S[x] for x in p]
    if any( len([T for T in I if i in T])!=2 for i in range(6) if i not in [r1,r2]):
      continue
    if any( len([T for T in I if i in T])!=1 for i in [r1,r2]):
      continue
    TL.append( prod(P('a%s' % st(T)) for T in I) )
  return list(set(TL))

@cached_function
def Delta(L):
  #L list of |Flist| coefficients.
  assert(len(L)==len(Delta_variable_list))
  return sum(l*p for l,p in zip(L,Delta_variable_list))

@cached_function
def Check1(r1,r2,L):
  gamma=vector([a[s] for s in Sstr])
  dDelta=vector([Delta(L).derivative(a[s]) for s in Sstr])
  res=dDelta*Wedge(r1,r2)*gamma
  #print "res", res
  TL=Test_variable_list(r1,r2)
  #define column
  return [(QQ(res.coefficient(tl)) if tl in res.monomials() else 0) for tl in TL]

def TestMatrix(r1,r2):
  return matrix([ Check1(r1,r2, tuple([0]*(i) + [1] + [0]*(255-i))) for i in range(256)]).transpose()


TestMatrix=block_matrix([[TestMatrix(0,1)], [TestMatrix(0,2)], [TestMatrix(0,3)], [TestMatrix(0,4)], [TestMatrix(0,5)], [TestMatrix(1,2)]])
print "Rank of TestMatrix: ", TestMatrix.rank()

Coef=TestMatrix.right_kernel().basis()[0]
Delta_final = -2*sum(c*p for c,p in zip(Coef, Delta_variable_list))

#Coef=[-2, 2, -2, -4, 4, -2, -2, 2, 2, -2, 2, 1, 4, 4, 2, -2, -4, 2, -4, 2, -4, -2, 2, -2, 4, 4, 4, 4, 4, 1, -4, 4, 2, -2, -4, -4, 4, -2, 2, -4, -2, 1, -4, -4, -2, 4, -2, 1, -4, 4, -4, 1, -2, 4, 4, -2, -4, 4, 2, 4, 4, 1, 4, -2, -2, 4, -2, -2, 4, 2, -2, -4, 4, -2, -2, -2, -4, -4, -4, 4, -2, -2, -2, -2, -4, -4, 4, -4, -4, 2, -4, 2, 1, 4, -4, 2, 4, 1, 4, 4, 2, -4, 4, -4, -2, -2, -2, -4, -2, 4, 4, 4, 4, -4, 1, -2, -4, -4, 4, 1, -2, 2, -4, 2, -4, 2, -4, -2, -2, 2, 4, -2, 2, -2, -2, 1, -2, 2, 4, 4, 4, 4, 2, -2, -4, 1, -4, -4, -4, 4, 2, -4, -2, -2, -4, 4, -2, 2, 1, 2, -4, -2, 2, -2, 4, 2, 4, -4, -2, -4, -2, -4, 4, -2, 2, -2, 4, 4, -2, -4, -4, 4, -2, -4, -2, -2, -4, 4, -2, -4, 1, 1, 2, 4, -2, 4, 2, 2, 2, 2, -4, 2, -4, -4, -4, 2, 4, -4, -4, -2, -2, 4, 2, -2, -2, 2, -2, 4, -2, 4, -2, 4, 2, 4, -4, 4, -2, -2, -2, 2, 4, 4, 4, 2, -2, 4, 4, -2, 2, 4, -4, -2, -2, 2, 1, 2, 4, 4, -2, 2, -2, 2, 2, -2, 2, -4]
#Delta_final = [c*p for c,p in zip(Coef, Flist)

print "Found Delta: ", Delta_final

#Delta_final = a45^2*a0123^2 - 2*a35*a45*a0123*a0124 + a35^2*a0124^2 + 2*a34*a45*a0123*a0125 - 2*a34*a35*a0124*a0125 + a34^2*a0125^2 + 2*a25*a45*a0123*a0134 - 2*a25*a35*a0124*a0134 - 2*a25*a34*a0125*a0134 + 4*a24*a35*a0125*a0134 - 4*a23*a45*a0125*a0134 + a25^2*a0134^2 - 2*a24*a45*a0123*a0135 + 4*a25*a34*a0124*a0135 - 2*a24*a35*a0124*a0135 + 4*a23*a45*a0124*a0135 - 2*a24*a34*a0125*a0135 - 2*a24*a25*a0134*a0135 + a24^2*a0135^2 - 4*a25*a34*a0123*a0145 + 4*a24*a35*a0123*a0145 - 2*a23*a45*a0123*a0145 - 2*a23*a35*a0124*a0145 + 2*a23*a34*a0125*a0145 + 2*a23*a25*a0134*a0145 - 2*a23*a24*a0135*a0145 + a23^2*a0145^2 - 2*a15*a45*a0123*a0234 + 2*a15*a35*a0124*a0234 + 2*a15*a34*a0125*a0234 - 4*a14*a35*a0125*a0234 + 4*a13*a45*a0125*a0234 - 2*a15*a25*a0134*a0234 - 2*a15*a24*a0135*a0234 + 4*a14*a25*a0135*a0234 - 4*a12*a45*a0135*a0234 + 2*a15*a23*a0145*a0234 - 4*a13*a25*a0145*a0234 + 4*a12*a35*a0145*a0234 + a15^2*a0234^2 + 2*a14*a45*a0123*a0235 - 4*a15*a34*a0124*a0235 + 2*a14*a35*a0124*a0235 - 4*a13*a45*a0124*a0235 + 2*a14*a34*a0125*a0235 + 4*a15*a24*a0134*a0235 - 2*a14*a25*a0134*a0235 + 4*a12*a45*a0134*a0235 - 2*a14*a24*a0135*a0235 - 2*a14*a23*a0145*a0235 + 4*a13*a24*a0145*a0235 - 4*a12*a34*a0145*a0235 - 2*a14*a15*a0234*a0235 + a14^2*a0235^2 + 4*a15*a34*a0123*a0245 - 4*a14*a35*a0123*a0245 + 2*a13*a45*a0123*a0245 + 2*a13*a35*a0124*a0245 - 2*a13*a34*a0125*a0245 - 4*a15*a23*a0134*a0245 + 2*a13*a25*a0134*a0245 - 4*a12*a35*a0134*a0245 + 4*a14*a23*a0135*a0245 - 2*a13*a24*a0135*a0245 + 4*a12*a34*a0135*a0245 - 2*a13*a23*a0145*a0245 + 2*a13*a15*a0234*a0245 - 2*a13*a14*a0235*a0245 + a13^2*a0245^2 - 4*a15*a24*a0123*a0345 + 4*a14*a25*a0123*a0345 - 2*a12*a45*a0123*a0345 + 4*a15*a23*a0124*a0345 - 4*a13*a25*a0124*a0345 + 2*a12*a35*a0124*a0345 - 4*a14*a23*a0125*a0345 + 4*a13*a24*a0125*a0345 - 2*a12*a34*a0125*a0345 + 2*a12*a25*a0134*a0345 - 2*a12*a24*a0135*a0345 + 2*a12*a23*a0145*a0345 - 2*a12*a15*a0234*a0345 + 2*a12*a14*a0235*a0345 - 2*a12*a13*a0245*a0345 + a12^2*a0345^2 + 2*a05*a45*a0123*a1234 - 2*a05*a35*a0124*a1234 - 2*a05*a34*a0125*a1234 + 4*a04*a35*a0125*a1234 - 4*a03*a45*a0125*a1234 + 2*a05*a25*a0134*a1234 + 2*a05*a24*a0135*a1234 - 4*a04*a25*a0135*a1234 + 4*a02*a45*a0135*a1234 - 2*a05*a23*a0145*a1234 + 4*a03*a25*a0145*a1234 - 4*a02*a35*a0145*a1234 - 2*a05*a15*a0234*a1234 - 2*a05*a14*a0235*a1234 + 4*a04*a15*a0235*a1234 - 4*a01*a45*a0235*a1234 + 4*a*a0145*a0235*a1234 + 2*a05*a13*a0245*a1234 - 4*a03*a15*a0245*a1234 + 4*a01*a35*a0245*a1234 - 4*a*a0135*a0245*a1234 - 2*a05*a12*a0345*a1234 + 4*a02*a15*a0345*a1234 - 4*a01*a25*a0345*a1234 + 4*a*a0125*a0345*a1234 + a05^2*a1234^2 - 2*a04*a45*a0123*a1235 + 4*a05*a34*a0124*a1235 - 2*a04*a35*a0124*a1235 + 4*a03*a45*a0124*a1235 - 2*a04*a34*a0125*a1235 - 4*a05*a24*a0134*a1235 + 2*a04*a25*a0134*a1235 - 4*a02*a45*a0134*a1235 + 2*a04*a24*a0135*a1235 + 2*a04*a23*a0145*a1235 - 4*a03*a24*a0145*a1235 + 4*a02*a34*a0145*a1235 + 4*a05*a14*a0234*a1235 - 2*a04*a15*a0234*a1235 + 4*a01*a45*a0234*a1235 - 4*a*a0145*a0234*a1235 - 2*a04*a14*a0235*a1235 - 2*a04*a13*a0245*a1235 + 4*a03*a14*a0245*a1235 - 4*a01*a34*a0245*a1235 + 4*a*a0134*a0245*a1235 + 2*a04*a12*a0345*a1235 - 4*a02*a14*a0345*a1235 + 4*a01*a24*a0345*a1235 - 4*a*a0124*a0345*a1235 - 2*a04*a05*a1234*a1235 + a04^2*a1235^2 - 4*a05*a34*a0123*a1245 + 4*a04*a35*a0123*a1245 - 2*a03*a45*a0123*a1245 - 2*a03*a35*a0124*a1245 + 2*a03*a34*a0125*a1245 + 4*a05*a23*a0134*a1245 - 2*a03*a25*a0134*a1245 + 4*a02*a35*a0134*a1245 - 4*a04*a23*a0135*a1245 + 2*a03*a24*a0135*a1245 - 4*a02*a34*a0135*a1245 + 2*a03*a23*a0145*a1245 - 4*a05*a13*a0234*a1245 + 2*a03*a15*a0234*a1245 - 4*a01*a35*a0234*a1245 + 4*a*a0135*a0234*a1245 + 4*a04*a13*a0235*a1245 - 2*a03*a14*a0235*a1245 + 4*a01*a34*a0235*a1245 - 4*a*a0134*a0235*a1245 - 2*a03*a13*a0245*a1245 - 2*a03*a12*a0345*a1245 + 4*a02*a13*a0345*a1245 - 4*a01*a23*a0345*a1245 + 4*a*a0123*a0345*a1245 + 2*a03*a05*a1234*a1245 - 2*a03*a04*a1235*a1245 + a03^2*a1245^2 + 4*a05*a24*a0123*a1345 - 4*a04*a25*a0123*a1345 + 2*a02*a45*a0123*a1345 - 4*a05*a23*a0124*a1345 + 4*a03*a25*a0124*a1345 - 2*a02*a35*a0124*a1345 + 4*a04*a23*a0125*a1345 - 4*a03*a24*a0125*a1345 + 2*a02*a34*a0125*a1345 - 2*a02*a25*a0134*a1345 + 2*a02*a24*a0135*a1345 - 2*a02*a23*a0145*a1345 + 4*a05*a12*a0234*a1345 - 2*a02*a15*a0234*a1345 + 4*a01*a25*a0234*a1345 - 4*a*a0125*a0234*a1345 - 4*a04*a12*a0235*a1345 + 2*a02*a14*a0235*a1345 - 4*a01*a24*a0235*a1345 + 4*a*a0124*a0235*a1345 + 4*a03*a12*a0245*a1345 - 2*a02*a13*a0245*a1345 + 4*a01*a23*a0245*a1345 - 4*a*a0123*a0245*a1345 - 2*a02*a12*a0345*a1345 - 2*a02*a05*a1234*a1345 + 2*a02*a04*a1235*a1345 - 2*a02*a03*a1245*a1345 + a02^2*a1345^2 - 4*a05*a14*a0123*a2345 + 4*a04*a15*a0123*a2345 - 2*a01*a45*a0123*a2345 + 4*a05*a13*a0124*a2345 - 4*a03*a15*a0124*a2345 + 2*a01*a35*a0124*a2345 - 4*a04*a13*a0125*a2345 + 4*a03*a14*a0125*a2345 - 2*a01*a34*a0125*a2345 - 4*a05*a12*a0134*a2345 + 4*a02*a15*a0134*a2345 - 2*a01*a25*a0134*a2345 + 4*a*a0125*a0134*a2345 + 4*a04*a12*a0135*a2345 - 4*a02*a14*a0135*a2345 + 2*a01*a24*a0135*a2345 - 4*a*a0124*a0135*a2345 - 4*a03*a12*a0145*a2345 + 4*a02*a13*a0145*a2345 - 2*a01*a23*a0145*a2345 + 4*a*a0123*a0145*a2345 - 2*a01*a15*a0234*a2345 + 2*a01*a14*a0235*a2345 - 2*a01*a13*a0245*a2345 + 2*a01*a12*a0345*a2345 + 2*a01*a05*a1234*a2345 - 2*a01*a04*a1235*a2345 + 2*a01*a03*a1245*a2345 - 2*a01*a02*a1345*a2345 + a01^2*a2345^2 + 4*a05*a14*a23*a012345 - 4*a04*a15*a23*a012345 - 4*a05*a13*a24*a012345 + 4*a03*a15*a24*a012345 + 4*a04*a13*a25*a012345 - 4*a03*a14*a25*a012345 + 4*a05*a12*a34*a012345 - 4*a02*a15*a34*a012345 + 4*a01*a25*a34*a012345 - 4*a04*a12*a35*a012345 + 4*a02*a14*a35*a012345 - 4*a01*a24*a35*a012345 + 4*a03*a12*a45*a012345 - 4*a02*a13*a45*a012345 + 4*a01*a23*a45*a012345 - 2*a*a45*a0123*a012345 + 2*a*a35*a0124*a012345 - 2*a*a34*a0125*a012345 - 2*a*a25*a0134*a012345 + 2*a*a24*a0135*a012345 - 2*a*a23*a0145*a012345 + 2*a*a15*a0234*a012345 - 2*a*a14*a0235*a012345 + 2*a*a13*a0245*a012345 - 2*a*a12*a0345*a012345 - 2*a*a05*a1234*a012345 + 2*a*a04*a1235*a012345 - 2*a*a03*a1245*a012345 + 2*a*a02*a1345*a012345 - 2*a*a01*a2345*a012345 + a^2*a012345^2



#We check Delta_final is invariant under so(12)

dDelta = vector([Delta_final.derivative(a[s]) for s in Sstr])
gamma=vector([a[s] for s in Sstr])

def res1(r1,r2):
  return dDelta*Wedge(r1,r2)*gamma

def res2(r1,r2):
  return dDelta*AntiWedge(r1,r2)*gamma

def res3(r1,r2):
  return dDelta*WedgeAntiWedge(r1,r2)*gamma 

print "Invariance under so(12):", all( all([res1(r1,r2)==0, res2(r1,r2)==0, res3(r1,r2)==0]) for r1,r2 in xmrange([6]*2) if r1 < r2)




##We calculate Delta(r_1 e^{omega_1} + r_2 e^{omega_2})
s01,s02,s03,s04,s05,s12,s13,s14,s15,s23,s24,s25,s34,s35,s45=var('s01,s02,s03,s04,s05,s12,s13,s14,s15,s23,s24,s25,s34,s35,s45')
t01,t02,t03,t04,t05,t12,t13,t14,t15,t23,t24,t25,t34,t35,t45=var('t01,t02,t03,t04,t05,t12,t13,t14,t15,t23,t24,t25,t34,t35,t45')
r1,r2=var('r1,r2')

#vector of e^{\omega} with respect to basis S. omega = \sum_{i<j} s_{ij} e[i]e[j].
Omegavector=[1, s01, s02, s03,s04,s05, s12,s13,s14,s15,s23,s24,s25,s34,s35,s45, \
 (s03*s12 - s02*s13 + s01*s23), \
 (s04*s12 - s02*s14 + s01*s24), \
 (s05*s12 - s02*s15 + s01*s25), \
 (s04*s13 - s03*s14 + s01*s34), \
 (s05*s13 - s03*s15 + s01*s35), \
 (s05*s14 - s04*s15 + s01*s45), \
 (s04*s23 - s03*s24 + s02*s34), \
 (s05*s23 - s03*s25 + s02*s35), \
 (s05*s24 - s04*s25 + s02*s45), \
 (s05*s34 - s04*s35 + s03*s45), \
 (s14*s23 - s13*s24 + s12*s34), \
 (s15*s23 - s13*s25 + s12*s35), \
 (s15*s24 - s14*s25 + s12*s45), \
 (s15*s34 - s14*s35 + s13*s45), \
 (s25*s34 - s24*s35 + s23*s45), \
 s05*s14*s23 - s04*s15*s23 - s05*s13*s24 + s03*s15*s24 + s04*s13*s25 - s03*s14*s25 + s05*s12*s34 - s02*s15*s34 + s01*s25*s34 - s04*s12*s35 + s02*s14*s35 - s01*s24*s35 + s03*s12*s45 - s02*s13*s45 + s01*s23*s45]

Omegavector2=[1, t01, t02, t03,t04,t05, t12,t13,t14,t15,t23,t24,t25,t34,t35,t45, \
 (t03*t12 - t02*t13 + t01*t23), \
 (t04*t12 - t02*t14 + t01*t24), \
 (t05*t12 - t02*t15 + t01*t25), \
 (t04*t13 - t03*t14 + t01*t34), \
 (t05*t13 - t03*t15 + t01*t35), \
 (t05*t14 - t04*t15 + t01*t45), \
 (t04*t23 - t03*t24 + t02*t34), \
 (t05*t23 - t03*t25 + t02*t35), \
 (t05*t24 - t04*t25 + t02*t45), \
 (t05*t34 - t04*t35 + t03*t45), \
 (t14*t23 - t13*t24 + t12*t34), \
 (t15*t23 - t13*t25 + t12*t35), \
 (t15*t24 - t14*t25 + t12*t45), \
 (t15*t34 - t14*t35 + t13*t45), \
 (t25*t34 - t24*t35 + t23*t45), \
 t05*t14*t23 - t04*t15*t23 - t05*t13*t24 + t03*t15*t24 + t04*t13*t25 - t03*t14*t25 + t05*t12*t34 - t02*t15*t34 + t01*t25*t34 - t04*t12*t35 + t02*t14*t35 - t01*t24*t35 + t03*t12*t45 - t02*t13*t45 + t01*t23*t45]

v=[r1*Omegavector[i] + r2*Omegavector2[i] for i in range(32)]

print factor(expand(Delta_final(v)))