##Mirror computations for Pencil of Fano variety of cubic fourfolds.
##Georg Oberdieck, 2021

##NL invariants
@cached_function
def NLdisc(d):
  #Returns the coefficient of q^{d/6} in the NL modular form.
  if d not in ZZ: return 0
  if d % 6 not in [0,2]: return 0
  q=PowerSeriesRing(QQ,'q').gen()
  alpha = 1 + 6*sum( sum( legendre_symbol(d,3) for d in divisors(n))*q**n for n in range(1,d+1)) + O(q**(d+1))
  beta = sum( sum((n/d)**2*legendre_symbol(d,3) for d in divisors(n))*q**n for n in range(1,d+1)) + O(q**(d+1))
  NL_modform=-alpha(q**6)**11 + 162*alpha(q**6)**8*beta(q**6) + 91854*alpha(q**6)**5*beta(q**6)**2 + 2204496*alpha(q**6)**2*beta(q**6)**3 - alpha(q**2)**11 + 66*alpha(q**2)**8*beta(q**2) - 1386*alpha(q**2)**5*beta(q**2)**2 + 9072*alpha(q**2)**2*beta(q**2)**3 + O(q**(d+1))
  return NL_modform.padded_list()[d]

@cached_function
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=-(6*s-d**2)*Rational((1,2))
  if x not in ZZ:
    return 0
  return NLdisc(int(x))

def NLrefined(m,s,d):
  #returns refined NL numbers of the family counting classes beta of divisibility m, square (beta,beta)=s, and degree beta*H = d.
  if m==1:
    return sum( moebius(k)*NL( s*Rational((1, k**2)), d/k ) for k in divisors(d) )
  elif d % m ==0:
    return NLrefined(1, s*Rational((1,m**2)), d/m )
  else:
    return 0

@cached_function
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) = 6s - d^2 <= 0  <=> s <= d^2 / 6
  #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/3 satisfying t%0 (4). Then s = t/2.
  L = [Rational((t,2)) for t in range(-5, floor(Rational((d**2,3))) + 1) if ( ((t % 4) in [0,3]) and is_even(3*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) = 6s - d^2 <= 0  <=> s <= d^2 / 6
  #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.
  L = [(m,Rational((t,2))) for m in divisors(d) for t in range(-5, floor(Rational((d**2,3*m**2))) + 1) if ( ((t % 4) in [0,3]) and is_even(3*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 18*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 GWExp1(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 GWExp2(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))


##Computing GW Invariants of the fiber using mirror symmetry

#We define the Chow ring of A*(Gr(2,6) x P1) as quotient of Q[y1,y2,hh], where y_i -> c_i(U*), hh hyplerplane class on P1
Ry = PolynomialRing(QQ,'y1,y2,hh')
y1,y2,hh=Ry.gens()

#We determine relations of the Grassmannian (generated by c_i(Q) = 0 for i>4, where 0 -> U -> C^6 -> Q -> 0 is canonical sequence).
Rx=PolynomialRing(QQ,'x1,x2') 
Rxw=PolynomialRing(Rx,'w')
x1,x2,w = var('x1,x2,w')  #
rel_poly=Rx(taylor(1/(1-x1+x2), (x1,0), (x2,0), 10)) 

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

I=Ry.ideal([rels(5),rels(6),rels(7), rels(8), hh**2])
I2=Ry.ideal(I.groebner_basis())
Rc=QuotientRing(Ry,I2, names=['c1','c2', 'h'])
c1,c2,h=Rc.gens()

#Euler Class of the pencil of Fanos. Dim 5
#[F] = (18 c1^2 c2 + 9 c2^2) + h*( 6 c1^3 + 30 c2 c1 )
Fundamental_Class_of_Family=(18*c1**2*c2 + 9*c2**2) + h*(6*c1**3 + 30*c2*c1)

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


z=1
@cached_function
def I_ab(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)), w = hyperplane class on P1
  Res=0
  for d1 in range(0,d+1):
    d2=d-d1
    Res+=(-1)**d*(x1-x2 + (d1-d2)*z) * \
       prod(3*x1+w+k*z for k in range(1,3*d1+1)) * prod(2*x1+x2+w+k*z for k in range(1,2*d1+d2+1)) * prod(x1+2*x2+w+k*z for k in range(1,d1+2*d2+1)) * prod(3*x2+w+k*z for k in range(1,3*d2+1)) *\
       prod( (x1+k*z)**(-6) for k in range(1,d1+1) )* prod( (x2+k*z)**(-6) for k in range(1,d2+1) )

  Res2=Res/(x1-x2)
  return taylor(Res2, (x1,0), (x2,0), (w,0), 10)

@cached_function
def I_gr(d):
  #return I function of Fano inside Gr(2,6) viewed inside Gr(2,6) x P1. We need to change the variables x_1, x_2 to the Chern classes c_1, c_2.
  F=Rxw(I_ab(d))
  res=0
  for f,e in zip(F.coefficients(),F.exponents()):
    if e>=2:
      continue
    symf=SymmetricFunctions(QQ).from_polynomial(f)
    symf_e = SymmetricFunctions(QQ).elementary()(symf)
    symf_e = symf_e.restrict_parts(2)
    res+= sum(c*prod(Rc.gens()[i-1] for i in P) for P,c in symf_e)*h**e
  return res

#q=PowerSeriesRing(Rc, 'q').gen()
#qprecision=10
@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, 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,6)) #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),6)*c1**3/(c2**4*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)) 



KNOWN_LIST=[]
for d in range(1,10):
  print("--------------- Degree: ", d, " -------------------")
  print("Check:", GW1bps(d) == GWExp1(d), GW2bps(d) == GWExp2(d))
  print("GWNLRelation_refined(d)", GWNLRelation_refined(d))
  NewCases=[(m,a) for (m,a) in GWNLRelation_refined(d) if m>1 and (a > 0) and ((m,a) not in KNOWN_LIST)]
  if len(NewCases)==0:
    print("No new cases")
  elif len(NewCases)==1:
    KNOWN_LIST.append( NewCases[0] )
    print("Added:", NewCases[0])
  else:
    print("More than one new case:", NewCases)

  print("Currently Known cases", KNOWN_LIST)