(*** General d scalar conformal block from Zamolodchikov recursion relation (code by Joshua Sandor) ***)
hIn[l_,r_,η_,d_]:= ((l! GegenbauerC[l,((d-2)/2),η])/((((-2))^l)Pochhammer[d/2-1,l](((1-(r^2)))^((d-2)/2))(Sqrt[(((1+(r^2)))^2)-4(r^2)(η^2)])))
c[1,m_,l_,d_]:= -((m(((-2))^m))/(((m!))^2))((Pochhammer[((1-m)/2),m])^2)
c[2,m_,l_,d_]:= -((m l!Pochhammer[d/2+l,-m]Pochhammer[d/2+l-1,-m])/((((-2))^m)(((m!))^2)(l-m)!Pochhammer[d-2+l,-m]))((Pochhammer[((1-m)/2),m])^2)
c[3,m_,l_,d_]:= -((m(((-1))^(m/2))Pochhammer[(d-m)/2-1,m])/(2((((m/2)!))^2)Pochhammer[(d-m)/2+l-1,m]Pochhammer[(d-m)/2+l,m]))((Pochhammer[((2-(d+m)/2-l)/2),m/2])^2)((Pochhammer[(((d-m)/2+l)/2),m/2])^2)
Δ[1,m_,l_,d_]:= 1-l-m
Δ[2,m_,l_,d_]:= l+d-1-m
Δ[3,m_,l_,d_]:= (d-m)/2
ll[1,m_,l_,d_]:= l+m
ll[2,m_,l_,d_]:= l-m
ll[3,m_,l_,d_]:= l
MaxOrder = 30;
hInExpansion[d_,l_]:=hInExpansion[d,l]  = Table[SeriesCoefficient[ Series[hIn[l,r,η,d],{r,0,MaxOrder}],n],{n,0,MaxOrder}];
(*order needs to be even!*)
RecursiveBlock[0,Δb_, lb_,d_]:=RecursiveBlock[0,Δb, lb,d]= hInExpansion[d,lb][[1]]
RecursiveBlock[ord_,Δb_, lb_,d_]:=RecursiveBlock[ord,Δb, lb,d]= Module[{hinf,J1, J2, J3},
hinf = Sum[hInExpansion[d,lb][[ii]] (r^(ii-1)),{ii,1,ord+1}];
J1 = Sum[((c[1,2k,lb,d](((4r))^(2k)))/(Δb- Δ[1,2k,lb,d]))RecursiveBlock[ord-2k,Δ[1,2k,lb,d]+2k,ll[1,2k,lb,d],d],{k,1,ord/2}];
J3= Sum[((c[3,2k,lb,d](((4r))^(2k)))/(Δb-Δ[3,2k,lb,d]))(RecursiveBlock[ord-2k,Δ[3,2k,lb,d]+2k,ll[3,2k,lb,d],d]),{k,1,ord/2}];
If[lb>= 2,
J2 = Sum[((c[2,2k,lb,d](((4r))^(2k)))/(Δb- Δ[2,2k,lb,d]))RecursiveBlock[ord-2k,Δ[2,2k,lb,d]+2k,ll[2,2k,lb,d],d],{k,1,Min[ord/2,lb/2]}];
,J2 = 0];
Return[Simplify[hinf+J1+ J2+ J3]]
]
RecursiveBlockPoleExp[0,Δb_, lb_,d_]:=RecursiveBlockPoleExp[0,Δb, lb,d]= hInExpansion[d,lb][[1]]
RecursiveBlockPoleExp[ord_,Δb_, lb_,d_]:=RecursiveBlockPoleExp[ord,Δb, lb,d]= Module[{hinf,J1, J2, J3},
hinf = Sum[hInExpansion[d,lb][[ii]] (r^(ii-1)),{ii,1,ord+1}];
J1 = Sum[c[1,2k,lb,d](((4r))^(2k))Pole[Δb,Δ[1,2k,lb,d]]RecursiveBlock[ord-2k,Δ[1,2k,lb,d]+2k,ll[1,2k,lb,d],d],{k,1,ord/2}];
J3= Sum[c[3,2k,lb,d](((4r))^(2k))Pole[Δb,Δ[3,2k,lb,d]](RecursiveBlock[ord-2k,Δ[3,2k,lb,d]+2k,ll[3,2k,lb,d],d]),{k,1,ord/2}];
If[lb>= 2,
J2 = Sum[c[2,2k,lb,d](((4r))^(2k))Pole[Δb, Δ[2,2k,lb,d]]RecursiveBlock[ord-2k,Δ[2,2k,lb,d]+2k,ll[2,2k,lb,d],d],{k,1,Min[ord/2,lb/2]}];
,J2 = 0];
Return[Simplify[hinf+J1+ J2+ J3]]
]
RecursiveG[ord_,Δb_, lb_,d_]:= (((4r))^(Δb))RecursiveBlock[ord,Δb, lb,d]

(*** The general d conformal block with external weight factor stripped off, first for scalar intermediate primary, and then general internal spin recursively. This is slow but serves as a check for the Zamolodchikov recursion ***)
slowG[d_,Δ_,0,z_,zbar_,trunc_]:=Block[{u=z zbar,v=(1-z)(1-zbar)},u^(Δ/2) Sum[Pochhammer[Δ/2,m]^2 Pochhammer[Δ/2,(m+n)]^2/(m! n! Pochhammer[Δ+1-d/2,m] Pochhammer[Δ,2m+n]) u^m (1-v)^n,{m,0,trunc},{n,0,trunc}]]

beta[p_]:=p^2/(4(2p-1)(2p+1))

slowG[d_,Δ_,s_,z_,zbar_,trunc_]:=Block[{a=d/2-1},(1/2 (Δ+s-2a-2) (1/z+1/zbar-1) slowG[d,Δ+1,s-1,z,zbar,trunc]+(1-z)D[slowG[d,Δ+1,s-1,z,zbar,trunc],z]+(1-zbar) D[slowG[d,Δ+1,s-1,z,zbar,trunc],zbar]-(s-1)(Δ(Δ-2a+1)/((Δ-a)(Δ-a+1)) beta[1/2(Δ-s+2-2a)] slowG[d,Δ+2,s-2,z,zbar,trunc]+(Δ+s-1)/(s+a-1) slowG[d,Δ,s-2,z,zbar,trunc]))/((Δ-a)(s+2a-1)/(s+a-1)) ]/;s≥1

(*** Conversion between r, η and z, zbar ***)
subs={r->Sqrt[z zbar]/((1+Sqrt[1-z])(1+Sqrt[1-zbar])),
η->(z/(1+Sqrt[1-z])^2+zbar/(1+Sqrt[1-zbar])^2)/(2Sqrt[z zbar]/((1+Sqrt[1-z])(1+Sqrt[1-zbar])))};

(*** testing in 3D case ***)
block1=Assuming[x>0 && y>0 && ϵ>0,(RecursiveG[10,Δ,2,3]/(r^Δ)/.subs/.{z->ϵ x,zbar->ϵ y})+O[ϵ]^5];
block2=Assuming[x>0 && y>0 && ϵ>0,(slowG[3,Δ,2,z,zbar,6]/(r^Δ)/.subs//.{z->ϵ x,zbar->ϵ y})+O[ϵ]^5];
(*** compare Zamolodchikov recursion against spin recursion, in the spin 2 case (overall constant normalization factor is a matter of convention) ***)
block1/block2//Simplify

(*** 3D conformal bootstrap ***)
Crossing[Δ1_,Δ_,s_,x_,y_,trunc_]:=(((1-z)^Δ1(1-zbar)^Δ1 RecursiveG[20,Δ,s,3]/.subs/.{z->1/2+ϵ x,zbar->1/2+ϵ y})+O[ϵ]^trunc//Normal)/.ϵ->1

CrossingDeriv[Δ1_,Δ_,s_,deg_,trunc_]:=Block[{x,y,mat=CoefficientList[Crossing[Δ1,Δ,s,x,y,trunc],{x,y}]},Sum[A[m,n] mat[[m+1,n+1]],{m,0,deg},{n,m+1,deg-m,2}]//Flatten]

PrepareCoeff[Δ1_,gap_,weighttrunc_,weightstep_,spintrunc_,deg_,trunc_]:=Block[{h},Monitor[Join[{CrossingDeriv[Δ1,0,0,deg,trunc]},Table[CrossingDeriv[Δ1,h,spin,deg,trunc],{spin,0,spintrunc,2},{h,If[spin==0,gap,1+spin],weighttrunc,weightstep}]//Flatten],{h,spin}]]

vars[deg_]:=Table[A[m,n],{m,0,deg},{n,m+1, deg-m,2}]//Flatten;

derivatives=7;

maxweight=8;
stepweight=.2;
maxspin=6;

truncationlevel=9;

normalization=vars[derivatives][[1]];
normvars=Drop[vars[derivatives],1];

blockcoeffs[Δ1_,gap_]:=PrepareCoeff[Δ1,gap,maxweight,stepweight,maxspin,derivatives,truncationlevel]/.normalization->1

oracle[Δ1_,gap_]:=Block[{MM=blockcoeffs[Δ1,gap]},(normvars[[1]]/.LinearOptimization[MM[[1]],#≥0 & /@ MM,normvars])===Indeterminate]

oracle[.518,1.41]

oracle[.518,1.415]

oracle[.518,1.42]

oracle[.518,1.425]

oracle[.518,1.427]

oracle[.518,1.428]

oracle[.518,1.43]

PrepareOPEbound[Δ1_,gap_,weighttrunc_,weightstep_,spintrunc_,deg_,trunc_]:=Block[{h},Monitor[Table[CrossingDeriv[Δ1,h,spin,deg,trunc],{spin,0,spintrunc,2},{h,If[spin==0,gap,1+spin],weighttrunc,weightstep}]//Flatten,{h,spin}]]

derivatives=8;
maxweight=9;
stepweight=.1;
maxspin=8;
truncationlevel=9;

AllVars=vars[derivatives];

OPEboundcoeffs[Δ1_,gap_]:=PrepareOPEbound[Δ1,gap,maxweight,stepweight,maxspin,derivatives,truncationlevel]

NormalizationCondition[Δ1_]:=CrossingDeriv[Δ1,3,2,derivatives,truncationlevel];
MinimizationObjective[Δ1_]:=-CrossingDeriv[Δ1,0,0,derivatives,truncationlevel];

OPEbound[Δ1_,gap_]:=OPEbound[Δ1,gap]=Block[{MM=OPEboundcoeffs[Δ1,gap]},MinimizationObjective[Δ1]/.LinearOptimization[MinimizationObjective[Δ1],Join[#≥0 & /@ MM,{NormalizationCondition[Δ1]==1}],AllVars]]

UnnormalizedCTbound[Δ1_,gap_]:=Δ1^2/OPEbound[Δ1,gap]

UnnormalizedCTbound[.518,1.]
UnnormalizedCTbound[.5,1.]

UnnormalizedCTbound[.518,1.]/UnnormalizedCTbound[.5,1.]