##Mirror computations for Pencil of Debarre-Voisin fourfolds
##Georg Oberdieck, 2021

def GWNLRelation(d):
  #GW of the fiber of degree d = sum_{s} NL_{s,d} GW_{s,d}.
  #We output all s which appear in the sum
  #For non-zero NL number we need: Delta(d,s) = 22s - d^2 <= 0  <=> s <= d^2 / 22
  #For non-zero GW BPS numbers we need: s >= -5/2
  #Finally, 2s in ZZ with (2s ) % 0 mod 4.
  #We iterate over t such that -5 <= t <= d^2/11 satisfying t%0 (4). Then s = t/2.
  return [Rational((t,2)) for t in range(-5, floor(Rational((d**2,11))) + 1) if ( (t%4 in [0,3]) and is_even(11*t-d**2) )  ]
  #return [ s for s in L if NL(s,d) != 0]

def GWNLRelation_refined(d):
  #GW of the fiber of degree d = sum_{m,s} NL_{m,s,d} GW_{m,s,d}.
  #We output all pairs (m,s/m^2) which need to be checked.
  #For non-zero NL number we need: Delta(d,s) = 22s - d^2 <= 0  <=> s <= d^2 / 22
  #For non-zero GW number we need: s/m^2 >= -5/2
  #Finally, 2s in ZZ with (2s ) % 0 mod 4.
  #We iterate over t' such that -5 <= t' <= (d/m)^2/3 satisfying t%0 (4). Then s = m^2 t/2.
  return [(m,Rational((t,2))) for m in divisors(d) for t in range(-5, floor(Rational((d**2,11*m**2))) + 1) if ( (t%4 in [0,3]) and is_even(11*t-(d/m)**2) ) ]
  #return [ (m,a) for (m,a) in L if NLrefined(1,a,d/m) != 0]




##Expected GW answer for the BPS fiber invariants
@cached_function
def bps_f(s):
  #returns the coefficient given by <beta,beta>=s of F = 1/4*Theta^2/Delta
  k=int(2*s)
  prec=max(k+10,2)
  q=PowerSeriesRing(QQ,'q').gen()
  alpha=sum( sigma(n)*q**n for n in range(1,prec) if is_odd(n) ) + O(q**prec)
  theta=sum( q**(n**2) for n in range(-floor(sqrt(prec+5)),ceil(sqrt(prec+5))+1 ) ) + O(q**prec)
  Delta = prod( (1-q**(4*n))**24 for n in range(1,floor(Rational((prec,4))) )) + O(q**prec)
  f=Rational((-1,4))/(alpha*theta*Delta)
  if k+4 in f.exponents():
    return (-1 if k%4==3 else 1)*f.coefficients()[f.exponents().index(k+4)]
  else:
    return 0

@cached_function
def bps_g(s):
  #returns the coefficient given by <beta,beta>=s of F = 1/4*Theta^2/Delta
  k=int(2*s)
  prec=max(k+10,2)
  q=PowerSeriesRing(QQ,'q').gen()
  alpha=sum( sigma(n)*q**n for n in range(1,prec) if is_odd(n) ) + O(q**prec)
  theta=sum( q**(n**2) for n in range(-floor(sqrt(prec+5)),ceil(sqrt(prec+5))+1 ) ) + O(q**prec)
  Delta = prod( (1-q**(4*n))**24 for n in range(1,floor(Rational((prec,4)))+1 )) + O(q**prec)
  G2 = Rational((-1,24))+sum(sigma(n)*q**(4*n) for n in range(1,floor(Rational((prec,4)))+1)) +O(q**prec)
  g=(theta**4 + 4*alpha + 24*G2)/(12*alpha*theta*Delta)
  if k+4 in g.exponents():
    return (-1 if k%4==3 else 1)*g.coefficients()[g.exponents().index(k+4)]
  else:
    return 0

##GW BPS of HK
def GW_X_1(s,d):
  #Return <H^3>_{g=0,s,d}
  return 66*d*bps_g(s)

def GW_X_2(s,d):
  #Return <psi_1^3>_{g=0,s,d}
  return 24*bps_f(s)


##GW For fiber expected BPS:
def GWFamily_1(d):
  #<H^3>_{d}
  #for s in GWNLRelation(d):
  #  print((s,d),GW_X_1(s,d)*NL(s,d))
  return sum( GW_X_1(s,d)*NL(s,d) for s in GWNLRelation(d))

##Expected GW For fiber BPS:
def GWFamily_2(d):
  #<psi_1^3>_{d}
  #for s in GWNLRelation(d):
  #  print((s,d),GW_X_2(s,d)*NL(s,d))
  return sum( GW_X_2(s,d)*NL(s,d) for s in GWNLRelation(d))

NLvalues=[0,  2,  6,  8,  10,  18,  22,  24,  28,  30, 32,  40,  44,  46,  50,  52,  54,  66,  68,  72,  88]
NLvariables=var(['NL_%s' % i for i in NLvalues])

def NL(s,d):
  #returns unrefined NL numbers
  t=2*s
  if not (t in ZZ and int(t)%4 in [0,3]):
    return 0
  x=-(22*s-d**2)*Rational((1,2))
  if x not in ZZ:
    return 0
  return NLvariables[ NLvalues.index( int(x) ) ]




##Computing GW Invariants of the fiber using mirror symmetry

#We define the Chow ring of A*(Gr(6,10) x P1) as quotient of Q[y1,y2,y3,y4,hh], where y_i -> c_i(U*), U universal sub, hh hyplerplane class on P1
#We know that the Chern classes of quotient bundle generate the Chow ring, hence also the first four chern classes of universal sub.
Ry = PolynomialRing(QQ,'y1,y2,y3,y4,hh')
y1,y2,y3,y4,hh=Ry.gens()

Rx=PolynomialRing(QQ,'x1,x2,x3,x4,x5,x6')
Rxw=PolynomialRing(Rx,'w')

RPowx=PowerSeriesRing(QQ, 'x1,x2,x3,x4,x5,x6')
RPowxw=PowerSeriesRing(RPowx,'w')

x1,x2,x3,x4,x5,x6,w = var('x1,x2,x3,x4,x5,x6,w')  #w the lift of the generator of A*(P1)
rel_poly=Rx(taylor(1/(1-x1+x2-x3+x4), (x1,0), (x2,0), (x3,0), (x4,0), 25))

def rels(i): return sum( c*y1**e[0]*y2**e[1]*y3**e[2]*y4**e[3] for c,e in zip(rel_poly.coefficients(),rel_poly.exponents()) if e[0]+2*e[1]+3*e[2]+4*e[3]==i)

I=Ry.ideal([rels(i) for i in range(7,25)] + [hh**2])
I2=Ry.ideal(I.groebner_basis())
Rc=QuotientRing(Ry,I2, names=['c1','c2', 'c3', 'c4', 'h'])
c1,c2,c3,c4,h=Rc.gens()

Fundamental_Class_of_Family=105*c3**4*c4**2 - 90*c2*c3**2*c4**3 + 27*c2**2*c4**4 - 39*c1*c3*c4**4 + 6*c4**5 + 620*c3**5*c4*h - 230*c2*c3**3*c4**2*h + 234*c2**2*c3*c4**3*h - 1248*c1*c3**2*c4**3*h + 42*c1*c2*c4**4*h + 552*c3*c4**4*h

def degree_component(f,i):
  #return class of degree k in Rc
  f=f.reduce(I2)
  flift=f.lift()
  return sum( c*c1**e[0]*c2**e[1]*c3**e[2]*c4**e[3]*h**e[4] for c,e in zip(flift.coefficients(),flift.exponents()) if e[0]+2*e[1]+3*e[2]+4*e[3]+e[4]==i)



z=1
@cached_function
def I_ab(d):
  print("I_ab(d) for d = ", d)
  #return I function of Fano inside Gr(2,6) viewed inside A*(P5 x P5). x_1 = c_1(O(1,0)), x_2=c_1(O(0,1)).
  x1,x2,x3,x4,x5,x6=RPowx.gens()
  x=[0,x1,x2,x3,x4,x5,x6]
  w=RPowxw.gen()

  Res = sum( Contr( (d1,d2,d3,d4,d5,d-d1-d2-d3-d4-d5) ) for d1 in range(0,d+1) for d2 in range(0, d+1-d1) for d3 in range(0, d+1-d1-d2) for d4 in range(0, d+1-d1-d2-d3) for d5 in range(0, d+1-d1-d2-d3-d4))

  Res0, Res1 = Res.coefficients()
  print(len(Res0.exponents()), len(Res1.exponents()) )
  print(Res0.prec(), Res1.prec())
  for i in range(1,6):
    for j in range(i+1,7):
      print((i,j))
      Res0 = Res0/ (x[i]-x[j])
      Res1 = Res1/ (x[i]-x[j])

  return (Res0, Res1)

@cached_function
def Contr(d):
  #Returns contribution of this degree
  #d is of the form (d1, ..., d6)
  print(d)
  dd=[0]+list(d)
  x1,x2,x3,x4,x5,x6=RPowx.gens()
  x=[0,x1,x2,x3,x4,x5,x6]
  w=RPowxw.gen()

  if dd==sorted(dd):
    a=(-1)**sum(d) * prod( x[i]-x[j]+dd[i]-dd[j] for i in range(1,6) for j in range(i+1,7) ) + RPowx.O(18)
    b=prod((x[i1] + x[i2] + x[i3] + k + RPowx.O(18)) + w + O(w**2) for i1 in range(1,5) for i2 in range(i1+1,6) for i3 in range(i2+1,7) for k in range(1,dd[i1]+dd[i2]+dd[i3]+1))
    c=prod( (x[i]+k + RPowx.O(18))**(-10) for i in range(1,7) for k in range(1,dd[i]+1))
    return a*b*c

  else:
    #We sort to the right degree
    ddd=[(dd[i],i) for i in range(len(dd))]
    sorted_d,permutation = zip(*sorted(ddd))

    xperm=[0]+[ x[permutation[i]] for i in range(1,7) ]
    corr_sign=Permutation(list(permutation)[1:]).sign()  #We need to correct by the sign of the permutation (coming from the a term)

    a,b=Contr(sorted_d[1:]).coefficients()
    anew=corr_sign*a.substitute({ x[i]:xperm[i] for i in range(1,7)})
    bnew=corr_sign*b.substitute({ x[i]:xperm[i] for i in range(1,7)})
    return anew+w*bnew


@cached_function
def I_gr(d):
  w=RPowxw.gen()
  f0,f1=I_ab(d)
  print(f0.parent())
  #F=f0.polynomial() + w*f1.polynomial()
  res=0
  for f,e in zip([f0,f1],[0,1]):
    if e>=2:
      continue
    symf=SymmetricFunctions(QQ).from_polynomial(f.polynomial())
    symf_e = SymmetricFunctions(QQ).elementary()(symf)
    symf_e = symf_e.restrict_parts(4)
    res+= sum(c*prod(Rc.gens()[i-1] for i in P) for P,c in symf_e)*h**e
  return res



@cached_function
def I(k,prec):
  q=PowerSeriesRing(Rc, 'q').gen()
  return (1 if k==0 else 0) + sum(degree_component(I_gr(d),k)*q**d for d in range(1,prec)) + O(q**prec)

@cached_function
def f0_and_Qq_variable_change(prec):
  #Returns the inverse of Q=q*exp(f_1/f_0)
  qq=PowerSeriesRing(QQ,'q').gen()
  f0 = 1 + sum( QQ(I_gr(d).lift().constant_coefficient())*qq**d for d in range(1,prec)) + O(qq**prec)
  f1 = sum( QQ(I_gr(d).lift().coefficient({y1:1, y2:0, y3:0, y4:0, hh:0}))*qq**d for d in range(1,prec)) + O(qq**prec)

  Q = qq*(f1/f0).exp(prec)
  return (f0,Q.reverse())

@cached_function
def JGW(d):
  #Returns coefficient q^d of the Gromov-Witten J-function
  prec=d+1
  f0, Qrev = f0_and_Qq_variable_change(prec)
  Itotal = Fundamental_Class_of_Family*sum(I(k,prec) for k in range(0,3)) #Since the pencil is of dimension 5, we only need to know it up to degree 5.
  JJ = (-I(1,prec)/f0).exp(prec)/f0*Itotal
  Rqc=PowerSeriesRing(Rc, 'q')
  return JJ(Rqc(Qrev)).padded_list()[d]

@cached_function
def GW1(d):
  #Compute <H^3>_{d}
  return degree_component(JGW(d),22)*c1**3/(c4**6*h)

#@cached_function
#def GW2(d):
#  #compute <psi_1^3>_{d}
#  return degree_component(JGW(d),9)/(c2**4*h)

def GW1bps(d):
  return sum( (-1)**(d+int(d/k))*moebius(k)/k**2*GW1(d/k) for k in divisors(d))

#def GW2bps(d):
#  return sum( (-1)**(d+int(d/k))*moebius(k)/k**5*GW2(d/k) for k in divisors(d)) 




###Results

#Equations:
#deg 1: 0 = 264*NL_6
#deg 2: 130680 = 3960*NL_2 + 132*NL_24
#deg 3: 0 = 792*NL_10,
#deg 4: 3020160 = 264*NL_30 + 7920*NL_8,
#deg 5:  0 = 1320*NL_18,
#deg 6:  11880*NL_18 + 396*NL_40,
#deg 7:  1848*NL_30 + 55440*NL_8,
#deg 8:  266112*NL_10 + 15840*NL_32 + 528*NL_54,
#deg 9:  1176120*NL_2 + 71280*NL_24 + 2376*NL_46


##Calculation of Noether-Lefschetz modular form

#from weilrep import *
#E8=IntegralLattice(['E',8])
#U=IntegralLattice('U')
#F=IntegralLattice(matrix([[-2]]))
#E8m=E8.twist(-1)
#MDVpart=IntegralLattice(matrix([[-2,-1],[-1,-6]]))
#MDV=E8m.direct_sum(E8m).direct_sum(U).direct_sum(U).direct_sum(MDVpart)
#wDV=WeilRep(MDV.gram_matrix())
#wDV.modular_forms_dimension(11)

#sage: f=-10*w.modular_forms_basis(11,20)[0]
#sage: sum( (Rational((1,2)) if b != 0 else 1)*q**(int(11*b))*c(q**11) for a,b,c in f.fourier_expansion())    
#-10 + 640*q^11 + 990*q^12 + 5500*q^14 + 11440*q^15 + 21450*q^16 + 198770*q^20 + 510840*q^22 + 803440*q^23 + 1842720*q^25 + 2716340*q^26 + 3972320*q^27 + 15821520*q^31 + 29586880*q^33 + 39799980*q^34 + 70523970*q^36 + 92912160*q^37 + 121147620*q^38 + 329532060*q^42 + 524690520*q^44 + 657650400*q^45 + 1015808640*q^47 + 1252350000*q^48 + 1540719840*q^49 + 3377216480*q^53 + 4891095840*q^55 + 5851299080*q^56 + 8310836160*q^58 + 9870495360*q^59 + 11665672150*q^60 + 22242198990*q^64 + 30257991720*q^66 + 35200149600*q^67 + 47239297600*q^69 + 54495337380*q^70 + 62861641040*q^71 + 108747856800*q^75 + 141483221440*q^77 + 160817995140*q^78 + 207146877960*q^80 + 234776423200*q^81 + 265163160900*q^82 + 426936678520*q^86 + 537281886240*q^88 + 602145714720*q^89 + 751994036640*q^91 + 838028387790*q^92 + 934630852320*q^93 + 1424029516800*q^97 + 1746482128320*q^99 + 1929198954210*q^100 + 2351731419180*q^102 + 2595220874400*q^103 + 2855691669160*q^104 + 4165068142070*q^108 + 5003843069720*q^110 + 5483268655920*q^111 + 6555200245440*q^113 + 7152084916020*q^114 + 7812321496320*q^115 + 10996944441120*q^119 + 12991369488000*q^121 + 14092107537640*q^122 + 16580382116430*q^124 + 17984657560320*q^125 + 19457648548800*q^126 + 26595690836220*q^130 + 30983115705700*q^132 + 33443999334240*q^133 + 38828204353840*q^135 + 41760642196680*q^136 + 44978772352800*q^137 + 59979108777280*q^141 + 69049778684480*q^143 + 73962259801650*q^144 + 84899669397180*q^146 + 90987688424160*q^147 + 97273646861490*q^148 + 127002971109960*q^152 + 144738724802760*q^154 + 154567740935520*q^155 + 175710758567520*q^157 + 187045957138480*q^158 + 199425339071200*q^159 + 255668972349120*q^163 + 288834158689920*q^165 + 306519795394560*q^166 + 345523947723240*q^168 + 367002382575840*q^169 + 388925856887700*q^170 + 490769311119120*q^174 + 550179831771220*q^176 + 582829994438560*q^177 + 652127016928640*q^179 + 688827480198660*q^180 + 728765282868480*q^181 + 906815929201600*q^185 + 1009759780265280*q^187 + 1064036212525880*q^188 + 1182806293124940*q^190 + 1247772870392400*q^191 + 1313401762401300*q^192 + 1614130483144080*q^196 + 1786638715152720*q^198 + 1880756738762880*q^199 + 2078594354485440*q^201 + 2182183396054860*q^202 + 2294891937184000*q^203 + 2789412476046240*q^207 + 3070893194846400*q^209 + 3217939659108900*q^210 + 3537826960441840*q^212 + 3711976777907200*q^213 + 3886111529103540*q^214 + 4676731998669720*q^218 + O(q^220)