(*** Testing 4D scalar conformal blocks ***)(*** DD is the differential operator that, when acting on the 4D scalar block, gives one half of the SO(5,1) Casimir ***)
DD[f_]:=z^2(1-z) D[f,{z,2}]+((Δ1-Δ2-Δ3+Δ4)/2 -1) z^2 D[f,z]+(Δ1-Δ2)(Δ3-Δ4)/4 z f+2 z zbar/(z-zbar)(1-z) D[f,z]+zbar^2(1-zbar) D[f,{zbar,2}]+((Δ1-Δ2-Δ3+Δ4)/2 -1) zbar^2 D[f,zbar]+(Δ1-Δ2)(Δ3-Δ4)/4 zbar f+2 z zbar/(zbar-z)(1-zbar) D[f,zbar]

(*** DD0 is a new differential operator that splits into a holomorphic part and an anti-holomorphic part, whereas DD doesn't  ***)
DD0[f_]:=z^2(1-z) D[f,{z,2}]+((Δ1-Δ2-Δ3+Δ4)/2 +1) z^2 D[f,z]-2z D[f,z]+(Δ1-Δ2+2)(Δ3-Δ4-2)/4 z f+zbar^2(1-zbar) D[f,{zbar,2}]+((Δ1-Δ2-Δ3+Δ4)/2 +1) zbar^2 D[f,zbar]-2zbar D[f,zbar]+(Δ1-Δ2+2)(Δ3-Δ4-2)/4 zbar f

(*** checking that DD and DD0 are related by DD(1/(z-zbar)) = (1/(z-zbar))(DD0+2) ***)
DD[f[z,zbar]/(z-zbar)] -1/(z-zbar) (DD0[f[z,zbar]]+2f[z,zbar])//Simplify

(*** The 4D scalar conformal block ***)
Block4D[z_,zbar_]:=z^((Δ+s)/2+1) zbar^((Δ-s)/2)/(z-zbar) Hypergeometric2F1[(Δ+s)/2-(Δ1-Δ2)/2,(Δ+s)/2+(Δ3-Δ4)/2,Δ+s,z] Hypergeometric2F1[(Δ-s)/2-(Δ1-Δ2)/2-1,(Δ-s)/2+(Δ3-Δ4)/2-1,Δ-s-2,zbar]+zbar^((Δ+s)/2+1) z^((Δ-s)/2)/(zbar-z) Hypergeometric2F1[(Δ+s)/2-(Δ1-Δ2)/2,(Δ+s)/2+(Δ3-Δ4)/2,Δ+s,zbar] Hypergeometric2F1[(Δ-s)/2-(Δ1-Δ2)/2-1,(Δ-s)/2+(Δ3-Δ4)/2-1,Δ-s-2,z] 

(*** testing that the equation Casimir = Δ(Δ-4)+s(s+2) is obeyed by Taylor expansion in z, zbar ***)
(DD[Block4D[z, zbar]]/((Δ(Δ-4)+s(s+2))/2  Block4D[z, zbar])/.{Δ2->Δ1,Δ4->Δ3}/.{Δ->3,s->0}/.{z->ϵ z,zbar->ϵ zbar})+O[ϵ]^4//Simplify
(DD[Block4D[z, zbar]]/((Δ(Δ-4)+s(s+2))/2  Block4D[z, zbar])/.{Δ2->Δ1,Δ4->Δ3}/.{Δ->5,s->2}/.{z->ϵ z,zbar->ϵ zbar})+O[ϵ]^4//Simplify

(*** Testing the Casimir equation by plugging in random numerical values ***)
DD[Block4D[z,zbar]]-1/2 (Δ(Δ-4)+s(s+2)) Block4D[z,zbar]/.{Δ1->2,Δ2->2,Δ3->2,Δ4->2,Δ->3,s->0,z->.4+.5 I,zbar->.5-.3 I}

(*** testing numerical values of the conformal block itself ***)
Plot[(Block4D[x+.001I,x-.001I]//Re)/.{Δ1->2,Δ2->2,Δ3->2,Δ4->2,Δ->1.5,s->0},{x,0,1}]

(*** 4D conformal bootstrap ***)(*** After setting Δ1=Δ2=Δ3=Δ4, contribution to the crossing equation, (v^(Δ1))F[u,v]-(u^(Δ1))F[v,u], multiplied by (z-zbar) ***)
CrossingBlock[z_,zbar_]=(z-zbar)(1-z)^Δ1(1-zbar)^Δ1 Block4D[z,zbar]-(z-zbar)z^Δ1 zbar^Δ1 Block4D[1-z,1-zbar]/.{Δ2->Δ1,Δ3->Δ1,Δ4->Δ1}//Simplify

(*** Extracting the x^m xbar^n terms in the z=1/2+x expansion of CrossingBlock, for m+n=even (recall that we have multiplied by an extra (z-zbar) factor) ***)
CrossingDeriv[weight1_,weight_,spin_,deg_]:=Block[{mat=CoefficientList[((CrossingBlock[1/2+ϵ x,1/2+ϵ xbar]/.{Δ1->weight1,Δ->weight,s->spin})+O[ϵ]^(2 deg+1)//Normal//N)/.ϵ->1,{x,xbar}]},Sum[A[m,n] mat[[m+1,n+1]],{m,0,2 deg},{n,m+2,2 deg-m,2}]//Flatten]

(*** Collecting the Taylor-expanded crossing equations for spins up to spintrunc, weights from the unitarity bound or scalar gap up to weighttrunc, with sufficiently small weightstep ***)
PrepareCoeff[Δ1_,gap_,weighttrunc_,weightstep_,spintrunc_,deg_]:=Join[{CrossingDeriv[Δ1,0,0,deg]},Table[CrossingDeriv[Δ1,Δ,spin,deg],{spin,0,spintrunc,2},{Δ,If[spin==0,gap,2+spin],weighttrunc,weightstep}]//Flatten]

(*** Coefficients of the linear functional as variables for linear optimization ***)
vars[deg_]:=Table[A[m,n],{m,0,2 deg},{n,m+2,2 deg-m,2}]//Flatten;

(*** total derivative order on the crossing equation is effectively 2 times derivatives - 1 ***)
derivatives=4; 

(*** Numerical approximations are introduced here via cutoff on the spin, internal weight, and discretization of possible weight values; the need for cutoff and discretization of weights can be removed with "polynomial matrix programming", but we won't implement it here ***)
maxweight=12;
stepweight=.05;
maxspin=8;

(*** normalizing one of the linear functional coefficients so that the optimization problem can have a well-defined solution ***)
normalization=vars[derivatives][[1]];
normvars=Drop[vars[derivatives],1];

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

(*** oracle tells whether the gap is allowed (true) or not allowed (false)***)
oracle[Δ1_,gap_]:=Block[{MM=blockcoeffs[Δ1,gap]},(normvars[[1]]/.LinearOptimization[MM[[1]],#≥0 & /@ MM,normvars])===Indeterminate]

oracle[1.05,2.5]
oracle[1.05,2.51]
oracle[1.05,2.52]
oracle[1.05,2.53]
oracle[1.05,2.54]
oracle[1.05,2.55]