ClearAll[pmax,lmax,lmax2,h1,GVTab,PTTabnew,DTTabnew,PTTabGv,h1fun,h1int,h1ser];
SetDirectory[NotebookDirectory[]];
(* choose your favorite CICY *)
<<X42.m
{CY,chi,kappa,c2}

(* Castelnuovo bound on the D0-charge for PT invariants *)
MMin[Q_]:=Ceiling[-((Q^2)/(2kappa))-(Q/2)];
(* Castelnovo bound on the genus of GV invariants *)
gMax[Q_]:=Floor[((Q^2)/(2 kappa))+(Q/2)+1];
(* action of the spectral flow *)
mun[mu_,n_]:=mu+n kappa;
mn[mu_,M_,n_]:=M-n mu-(kappa/2)((n^2)+n);
(* spectral flow transformation bringing D2 brane charge to the interval [0,κ/2] *)
SFmu[κ_,μ_]:=Module[{x= Round[μ/κ]},If[κ x > μ,κ x-μ ,μ-κ x]];
SFM[κ_,μ_,M_]:=Module[{x= Round[μ/κ]},If[κ x > μ,M+μ x+(κ/2)x(1-x)-(κ x-μ),M+μ x+(κ/2)x(1-x)]];

ChiO1[mu_,M_,n_]:=(kappa/2)n(n-1)+(n-1) mu-M+(kappa/6)+(c2/12);

(* rules for substitution of invariants *)
repGV={Gv[g_,m_]:>If[g>gMax[m]||m==0,0,If[g+1<=Length[GVTab],If[m<=Length[GVTab[[g+1]]],GVTab[[g+1,m]],Gv[g,m]],Gv[g,m]]]};
repPT={PT[mu_,n_]:>If[n<MMin[mu],0,If[mu+1<=Length[PTTabnew],If[n-MMin[mu]+2<=Length[PTTabnew[[mu+1]]],PTTabnew[[mu+1,n-MMin[mu]+2]],PT[mu,n]],PT[mu,n]]]};repDT={DT[mu_,n_]:>If[n<MMin[mu],0,If[mu+1<=Length[DTTabnew],If[n-MMin[mu]+2<=Length[DTTabnew[[mu+1]]],DTTabnew[[mu+1,n-MMin[mu]+2]],DT[mu,n]],DT[mu,n]]]};

repPTGv={PT[mu_,n_]:>If[n<MMin[mu],0,If[mu+1<=Length[PTTabGv],If[n-MMin[mu]+2<=Length[PTTabGv[[mu+1]]],PTTabGv[[mu+1,n-MMin[mu]+2]],PT[mu,n]],PT[mu,n]]]};

JTab=CoefficientList[h1int,q];
repJ={J[mu_,n_]:>If[n<MMin[mu],0,If[mu+1<=Length[JTab]&& n+1<=Length[JTab[[mu+1]]]+MMin[mu],JTab[[mu+1,n-MMin[mu]+1]],J[mu,n]]]};

(* valid only for m close to Castelnuovo bound *)
repPTApprox={PT[b_,m_]:>Sum[Binomial[2g-2,g-1-m]Gv[g,b],{g,1,Floor[b^2/(2kappa)+b/2+1]}]/.repGV};
(* same but do not evaluate GV invariants *)
repPTApprox0={PT[b_,m_]:>Sum[Binomial[2g-2,g-1-m]Gv[g,b],{g,1,Floor[b^2/(2kappa)+b/2+1]}]};

(* FULL VERSION OF THEOREM *)
(* Determination of the minimal value of the spectral flow parameter n *)
fcon[x_]:=If[x<1, x+(1/2),If[x<(15/8),Sqrt[2x+(1/4)],If[x<(9/4),(2/3)x+(3/4),If[x<3,(1/3)x+(3/2),(1/2)x+1
]]]];
Con[beta_,m_]:=((2 beta)/3)fcon[(beta/kappa)]+m;
nMin2 [mu_,M_]:=Module[{n},n=1;
While[Con[mun[mu,n],mn[mu,M,n]]>=0,n++];
n
];
NTab2=Array[nMin2 [#1-1,#2-1+MMin[#1-1]]&,{IntegerPart[kappa/2]+1,10}];
(* adds the first column encoding the shifts of the D0-charge *)
 NTab2=Transpose[Join[{Table[-MMin[mu],{mu,0,IntegerPart[kappa/2]}]},Transpose[NTab2]]];
MatrixForm[NTab2]

(* limits of summation *)
ULbeta[beta_,m_]:=beta+((3 kappa m)/(2 beta))+(kappa/2);
ULm[beta_,m_,betap_]:=((((beta-betap))^2)/(2 kappa))+(1/2)(beta+betap)+m;
LLm[betap_]:=-((((betap))^2)/(2 kappa))-(1/2)betap;
(* answer from Soheyla's Thm 2 *)
S2[mu_,M_,n_]:=Module[{Mun,Mn},Mun=mun[mu,n];Mn=mn[mu,M,n];((((-1))^(ChiO1[mu,M,n]+1))/(ChiO1[mu,M,n]))(PT[Mun,Mn]-Sum[Sum[(((-1))^(ChiO1[mu,M,n]-mup+mp+1))(ChiO1[mu,M,n]-mup+mp)PT[mup,mp] J[SFmu[kappa,Mun-mup],SFM[kappa,Mun-mup,Mn-mp+mup]],{mp,Ceiling[LLm[mup]],Floor[ULm[Mun,Mn,mup]]}],{mup,1,IntegerPart[ULbeta[Mun,Mn]]}])
];

Ch[mu_,M_,n_]:=Module[{Mun,Mn},Mun=mun[mu,n];Mn=mn[mu,M,n];If[(-((3 Mn)/(2 Mun)))^2>2((Mun/kappa)),{(Mun/kappa),-((3 Mn)/(2 Mun))},{}]];
Rankm1WC[{x_,alpha_}]:=Module[{b1,b2,rmax,Li,r01max,xmin},
(* follows reasoning in proof of Lemma 2.1 *)
If[alpha^2<2x,Print["The BMT line does not intersect the parabola"]];
If[alpha^2==2x,Print["The BMT line is tangent to parabola"]];
b1=alpha-Sqrt[alpha^2-2x];
b2=alpha+Sqrt[alpha^2-2x];
Print["(b1,b2)=",N[{b1,b2}]];
If[alpha>fcon[x],Print["Satisfies Thm 1.1"]];
rmax=Floor[b1/(b2-b1)];
r01max=Floor[b2/(b2-b1)];
Li=Drop[Flatten[Table[Table[{r,c},{c,Ceiling[b2 r],Floor[b1(1+r)]}],{r,0,rmax}],1],1];
Print["Allowed (r,c):", Li];
Table[
Print["(r,c)=",Li[[i]],": E0:",Li[[i]],", E1:",{-1,0}-Li[[i]]];
Do[If[r0+r1>=1&& r0+r1<=r01max&&r1-1-Li[[i,1]]>=0,
xmin=-b2^2/2+1/2(b2^2-b1^2)(r0+r1);
Print["{{r0',r0},{r1',r1}=",{{Li[[i,1]]+r0,r0},{r1-1-Li[[i,1]],r1}},", x<=",N[xmin]]];
,{r0,0,r01max},{r1,0,r01max}];
,{i,Length[Li]}];
];
Rankm1WC0[{x_,alpha_}]:=Module[{b1,b2,rmax,Li,r01max,xmin},
(* follows reasoning in proof of Lemma 2.1 *)
b1=alpha-Sqrt[alpha^2-2x];
b2=alpha+Sqrt[alpha^2-2x];
rmax=Floor[b1/(b2-b1)];
r01max=Floor[b2/(b2-b1)];
Li=Drop[Flatten[Table[Table[{r,c},{c,Ceiling[b2 r],Floor[b1(1+r)]}],{r,0,rmax}],1],1]
];

(* coefficients for the minimal value of n *)
DTs=Table[Table[S2[mu,M,nMin2[mu,M]],{M,MMin[mu],7+MMin[mu]}],{mu,0,kappa/2}]/.repPT;
Charges=Table[Table[Ch[mu,M,nMin2[mu,M]],{M,MMin[mu],7+MMin[mu]}],{mu,0,kappa/2}]/.repPT;
Table[Table[S2[mu,M,nMin2[mu,M]-1],{M,MMin[mu],7+MMin[mu]}],{mu,0,kappa/2}]//MatrixForm
DTs//MatrixForm
DTs/.repJ//MatrixForm
Charges//MatrixForm

(* This is for X42 *)
DTs/.repJ/.{J[0,0]->-6,J[0,1]->880,J[1,0]->0,J[1,1]->-5120,J[2,-1]->0,J[2,0]->0,J[3,-2]->0,J[3,-1]->0,J[3,0]->7680,J[4,-3]->0}//MatrixForm

(* try if n-1 already produces the correct result *)
DTs1=Table[Table[S2[mu,M,nMin2[mu,M]-1],{M,MMin[mu],7+MMin[mu]}],{mu,0,kappa/2}]/.repPT;
Charges1=Table[Table[Ch[mu,M,nMin2[mu,M]-1],{M,MMin[mu],7+MMin[mu]}],{mu,0,kappa/2}]/.repPT;
Table[Table[S2[mu,M,nMin2[mu,M]-1],{M,MMin[mu],7+MMin[mu]}],{mu,0,kappa/2}]//MatrixForm
DTs1//MatrixForm
DTs1/.repJ//MatrixForm
Charges1//MatrixForm

(* This is for X42 *)
DTs1/.repJ/.{J[0,0]->-6,J[0,1]->880,J[1,0]->0,J[1,1]->-5120,J[2,-1]->0,J[2,0]->0,J[3,-2]->0,J[3,-1]->0,J[3,0]->7680,J[4,-3]->0}//MatrixForm

(* This is for X2222 *)
DTs1/.repJ/.{J[0,0]->-8,J[1,0]->0,J[2,-1]->0,J[2,0]->0,J[3,-1]->0,J[3,0]->0,J[4,-2]->0,J[4,-1]->0,J[5,-3]->0,J[5,-2]->0,J[5,-1]->0,J[6,-4]->0,J[6,-3]->0,J[6,-2]->0,J[7,-5]->0,J[7,-4]->0,J[7,-3]->0,J[8,-7]->0,J[8,-6]->0}//MatrixForm

(* allowed values of (r,c) for each (x,\alpha) *)
Charges1/.{x_,a_}:>Rankm1WC0[{x,a}]//MatrixForm

Rankm1WC[{17/16,3/2}]