SingularType={"E12","E13","E14b","Z11","Z12","Z13b","W12","W13b","Q10","Q11","Q12b","S11","S12","U12","E14","Z13","W12b","W13","Q12","U12b","U12c","Z13T","Q12T","U12bT"};

DefiningEquation22[TT_List,TXYZ_List]:=Module[{tt=TT,aa,bb,cc,txyz=TXYZ,x,y,z},x=txyz[[1]];y=txyz[[2]];z=txyz[[3]];aa=tt[[1]];bb=tt[[2]];cc=tt[[3]];{bb* y^7+ aa* x^3+cc*z^2, aa* x*y^5+bb*x^3+cc*z^2,bb*y^8+aa*x^3+cc*z^2,
bb*y^5+aa*y*x^3+cc*z^2,bb*x*y^4+aa*y*x^3+cc*z^2,
bb*y^6+aa*y*x^3+cc*z^2,
bb*y^5+aa*x^4+cc*z^2,
aa y*x^4+bb*y^4+cc*z^2,
bb*y^4+cc*z^3+aa*y*x^2,
bb*y^3*z+cc*z^3+aa*y*x^2,bb*y^5+cc*z^3+aa*y*x^2,
cc*z^4+aa*x^2*y+bb*y^2*z, 
cc*z^3*x+aa*x^2*y+bb*y^2*z,
cc*z^4+bb*y^3+aa*x^3,cc*z^3+bb*y^4*x+aa*x^2,cc*z^3*y+bb*x*y^3+aa*x^2,aa*x^2*y+cc*z^5+bb*y^2,bb*y^2*x+cc*y*z^4+aa*x^2,aa*x^2*y+bb*x*y^3+cc*z^3,aa*x^2*y+bb*y^3+cc*z^4,aa*x^2*y+bb*x*y^2+cc*z^4,aa*x^3+bb*x*y^6+cc*z^2,aa*x^2+bb*x*y^5+cc*z^3,aa*x^2+bb*x*y^3+cc*z^4}];

JacWWBasis[TXYZ_List]:=Module[{txyz=TXYZ,x,y,z,res},x=txyz[[1]];y=txyz[[2]];z=txyz[[3]];res={{1,y,(y^2),x,(y^3),x y,(y^4),x (y^2),(y^5),x (y^3),x (y^4),x (y^5)},{1,y,(y^2),x,(y^3),x y,(y^4),x (y^2),(x^2),x (y^3),(x^2) y,(x^2) (y^2),(x^2) (y^3)},{1,y,(y^2),x,(y^3),x y,(y^4),x (y^2),(y^5),x (y^3),(y^6),x (y^4),x (y^5),x (y^6)},{1,y,x,(y^2),x y,(x^2),(y^3),x (y^2),(y^4),x (y^3),x (y^4)},{1,y,x,(y^2),x y,(x^2),(y^3),x (y^2),(x^2) y,x (y^3),(x^2) (y^2),(x^2) (y^3)},{1,y,x,(y^2),x y,(y^3),(x^2),x (y^2),(y^4),x (y^3),(y^5),x (y^4),x (y^5)},{1,y,x,(y^2),x y,(x^2),(y^3),x (y^2),(x^2) y,x (y^3),(x^2) (y^2),(x^2) (y^3)},{1,x,y,(x^2),x y,(y^2),(x^3),(x^2) y,x (y^2),(y^3),(x^2) (y^2),x (y^3),(x^2) (y^3)},{1,y,z,x,(y^2),y z,x z,(y^3),(y^2) z,(y^3) z},{1,y,z,x,(y^2),y z,(z^2),x z,(y^2) z,y (z^2),(y^2) (z^2)},{1,y,z,x,(y^2),y z,(y^3),x z,(y^2) z,(y^4),(y^3) z,(y^4) z},{1,z,x,y,(z^2),x z,y z,(z^3),x (z^2),y (z^2),y (z^3)},{1,z,x,y,(z^2),x z,y z,x y,x (z^2),y (z^2),x y z,x y (z^2)},{1,z,x,y,(z^2),x z,y z,x y,x (z^2),y (z^2),x y z,x y (z^2)},{1,y,x,(y^2),x y,(y^3),x (y^2),z,y z,x z,(y^2) z,x y z,(y^3) z,x (y^2) z},{1,y,z,(y^2),y z,x,(z^2),(y^2) z,x y,x z,x (y^2),x y z,x (y^2) z},{1,x,y,z,x z,y z,(z^2),x (z^2),y (z^2),(z^3),x (z^3),y (z^3)},{1,z,y,(z^2),y z,x,(z^3),y (z^2),x z,x y,x (z^2),x y z,x y (z^2)},{1,x,y,x y,(y^2),x (y^2),z,x z,y z,x y z,(y^2) z,x (y^2) z},{1,x,y,(y^2),z,x z,y z,(y^2) z,(z^2),x (z^2),y (z^2),(y^2) (z^2)},{1,x,y,x y,z,x z,y z,x y z,(z^2),x (z^2),y (z^2),x y (z^2)},{1,y,(y^2),(y^3),(y^4),(y^5),x,x y,x (y^2),x (y^3),x (y^4),(x^2),(x^2) y,(x^2) (y^2),(x^2) (y^3),(x^2) (y^4)},{1,y,(y^2),(y^3),(y^4),x,x y,x (y^2),x (y^3),z,y z,(y^2) z,(y^3) z,(y^4) z,x z,x y z,x (y^2) z,x (y^3) z},{1,y,(y^2),x,x y,z,y z,(y^2) z,x z,x y z,(z^2),y (z^2),(y^2) (z^2),x (z^2),x y (z^2)}};res];

DegtxyzOne[EQN_,TXYZ_List]:=Module[{eqn=EQN,txyz=TXYZ, tab,len,mat,vv,i},  tab=CoefficientRules[eqn,txyz];len=Length[tab];mat=Table[tab[[i,1]],{i,1,len}];vv=Table[1,{i,1,len}];Inverse[mat].vv];

(*Define a partial order among vectors. 
Input: uu, vv. 
Output: if uu≥ vv, then 1; otherwise, 0.*)
PartialOrder[UU_List, VV_List]:=Module[{uu=UU,vv=VV,len,res=0,i},len=Length[uu];If[Table[Sign[uu[[i]]-vv[[i]]+1],{i,1,len}]==Table[1,{i,1,len}],res=1];res];

(*Truncate monomials in variable varthat are of order larger than MM*)
PolyTruncate[FF_,VARTAB_,MM_]:=Module[{ff=FF,var=VARTAB,mm=MM,i,oldcoef,newcoef,len},oldcoef=CoefficientRules[ff,var];newcoef=CoefficientRules[ff,var];len=Length[oldcoef];
For[i=1,i≤len,newcoef[[i,2]]=Which[Total[oldcoef[[i,1]]]≤mm,oldcoef[[i,2]],True,0];i++];FromCoefficientRules[newcoef,var]];

(*Pick up monomials in variable varthat are of order larger than MM*)
PolyPickup[FF_,VARTAB_,MM_]:=Module[{ff=FF,var=VARTAB,mm=MM,i,oldcoef,newcoef,len},oldcoef=CoefficientRules[ff,var];newcoef=CoefficientRules[ff,var];len=Length[oldcoef];
For[i=1,i≤len,newcoef[[i,2]]=Which[Total[oldcoef[[i,1]]]==mm,oldcoef[[i,2]],True,0];i++];FromCoefficientRules[newcoef,var]];


(*Input: defining equation "EQN" and list "TXYZ" of variables.
  Output: the list of Jacobian ideal*)

JacIdeal[EQN_,TXYZ_List]:=Module[{eqn=EQN,i,txyz=TXYZ,len,res}, len=Length[txyz];
res= Table[D[eqn, txyz[[i]]],{i,len}] ;res];
 

(*The next two commands come from 
http://forums.wolfram.com/mathgroup/archive/2011/Mar/msg00362.html*)

moduleGroebnerBasis[polys_,vars_,cvars_,opts___]:=Module[{newpols,rels,len=Length[cvars],gb,j,k,rul},rels=Flatten[Table[cvars[[j]]*cvars[[k]],{j,len},{k,j,len}]];
newpols=Join[polys,rels];
gb=GroebnerBasis[newpols,Join[cvars,vars],opts];
rul=Map[(#:>{})&,rels];
gb=Flatten[gb/.rul];
Collect[gb,cvars]]

(*Input: Jacobian ideal "polys" and list "vars" of variables
Output: {Groebner basis, conversion matrix that express Groebner basis in terms of the original generators*)
conversionMatrix[polys_,vars_]:=Module[{aa,coords,pmat,len=Length[polys],newpolys,mgb,gb,convmat,fvar,rvars},coords=Array[aa,len+1];
fvar=First[coords];
rvars=Rest[coords];
pmat=Transpose[Join[{polys},IdentityMatrix[len]]];
newpolys=pmat.coords;
mgb=moduleGroebnerBasis[newpolys,vars,coords];
gb=mgb/.Join[{fvar->1},Thread[rvars->0]]/.0:>Sequence[];
convmat=Select[mgb,!FreeQ[#,fvar]&]/.fvar->0;
{gb,convmat/.Thread[rvars->Table[UnitVector[len,j],{j,len}]]}]

(*
(*The next command produces {list of degrees, list of a monomial basis}, when the number of  variables is equal to 3, 4, or 5. Such command is not necessary in the current calculations*)

JacWWOneBasis[EQN_,TXYZ_List]:=Module[{eqn=EQN,txyz=TXYZ,tab,len,monotab,i,tab22,tmax,s,j,k,max4=0,max5=0,res={},tab33={},tab44,degtab,t1,r,q,d1,len22,len44},tab= conversionMatrix[JacIdeal[eqn,txyz],txyz][[1]];len=Length[tab];degtab=DegtxyzOne[eqn,txyz];monotab=Table[Exponent[MonomialList[tab[[i]],txyz][[1]],txyz],{i,len}];
tmax=Map[Max,Transpose[monotab]];len22=Length[txyz];If[len22≥4,max4=tmax[[4]]];If[len22≥5,max5=tmax[[5]]];
For[i=0,i≤tmax[[1]],For[j=0,j≤tmax[[2]],For[k=0,k≤tmax[[3]],
For[r=0,r≤max4,For[q=0,q≤max5, 
If[len22==3,t1={i,j,k},If[len22==4,t1={i,j,k,r},t1={i,j,k,r,q}]];d1=Total[Table[PartialOrder[t1,monotab[[s]]],{s,len}]];
If[d1==0,AppendTo[tab33,{degtab.t1,t1}]];q++];r++];k++];j++];i++];tab44=Transpose[Sort[tab33]];len44=Length[tab44[[2]]];{tab44[[1]],Table[FromCoefficientRules[{tab44[[2]][[i]]->1},txyz],{i,len44}]}];

*)


(*Input: defining equation EQN, variables TXYZ, a given monomial basis JAC, a Groebner basis GROBAS.
Output: {the matrix of residue pairing ((Subscript[ϕ, i]), (Subscript[ϕ, j])), dual  basis {(ϕ^1),(ϕ^2),..., (ϕ^μ)}*)
Metric[EQN_,TXYZ_List,JAC_List,GROBAS_List]:=Module[{eqn=EQN,txyz=TXYZ,jac=JAC, i,j,monocw, Groebner=GROBAS,degtab,degxyz,len,mat,cw,jacGB},degxyz=DegtxyzOne[eqn,txyz];len=Length[jac];  
degtab=Table[Exponent[jac[[i]],txyz].degxyz,{i,len}];cw=degtab[[len]];monocw=(PolynomialReduce[jac[[len]],Groebner,txyz])[[2]];
mat=Table[0,{i,len},{j,len}];For[i=1,i≤len,For[j=1,j≤len,If[degtab[[i]]+degtab[[j]]==cw,mat[[i,j]]=Coefficient[(PolynomialReduce[jac[[i]]*jac[[j]],Groebner,txyz])[[2]],monocw]];j++];i++];{mat,Inverse[mat].jac}];

 

(*Input: a polynomial FF, variables TXYZ, a given monomial basis JAC={(Subscript[ϕ, 1]),(Subscript[ϕ, 2]),..., (Subscript[ϕ, μ])} of Jac(W), a Groebner basis GROBAS,   the dual basis {(ϕ^1),(ϕ^2),..., (ϕ^μ)}, and the dimension μ of Jac(W).
Output:  *)
FindRemainer[FF_,TXYZ_List,JAC_List,GROBAS_List,DUALJAC_List,LEN_]:=Module[{f=FF, txyz=TXYZ,jac=JAC,Groebner=GROBAS,dualjac=DUALJAC,len=LEN,monocw,tab,i},
monocw=(PolynomialReduce[jac[[len]],Groebner,txyz])[[2]];
tab=Table[Coefficient[(PolynomialReduce[ dualjac[[i]]*f,Groebner,txyz])[[2]],monocw],{i,len}];tab.jac];

(* Input: a polynomial FF, together with other data. 
Output: reduction of FF in the cohomology *)

CohomologyReduction[FF_,TXYZ_List,JAC_List,GROBAS_List,DUALJAC_List,LEN_,CONV_List,LENXYZ_]:=Module[{f=FF, txyz=TXYZ,jac=JAC,Groebner=GROBAS,dualjac=DUALJAC, lenxyz=LENXYZ,conv=CONV,monocw,tab,i,len=LEN,mm=0,tab11={},polyrr,poly,red11},poly=f;While[mm==0,polyrr=FindRemainer[poly, txyz,jac,Groebner,dualjac,len];AppendTo[tab11,polyrr];red11=(PolynomialReduce[poly-polyrr,Groebner,txyz][[1]]).conv;poly=Expand[-t*Total[Table[D[red11[[i]],txyz[[i]]],{i,lenxyz}]]];If[ToString[poly]==ToString[0],mm=1]];Expand[Total[tab11]]];
  

 
(* Output the  primitive form up to order MM. 
 *)
PrimitiveForm[EQN_,TXYZ_List,JAC_List,MM_,VAR_]:=Module[{eqn=EQN,txyz=TXYZ,mm=MM,var=VAR,jac=JAC,len,lenxyz,tab,Groebner,conv,para,sca,k,psi,expotab,polyrow,red11,dualjac,psirow,len22,pos22,i,eone,vv,res11={},primitive},len=Length[jac];lenxyz=Length[txyz];tab=conversionMatrix[JacIdeal[eqn,txyz],txyz];Groebner=tab[[1]];conv=tab[[2]]; para=Table[(Subscript[var, j]),{j,1,len}]; dualjac=Metric[eqn,txyz,jac,Groebner][[2]];sca=Sum[(para.jac)^k/(k! t^k),{k,1,mm}]; psi=Table[0,{i,1,len},{j,1,len}];expotab=Table[Exponent[jac[[i]],txyz],{i,1,len}];For[i=1,i≤len,polyrow=jac[[i]]*sca; red11=Coefficient[t^mm*CohomologyReduction[polyrow,txyz,jac,Groebner,dualjac,len,conv,lenxyz],t^mm];psirow=CoefficientRules[red11,txyz]; len22=Length[psirow];
For[k=1,k≤len22,pos22=Position[expotab,psirow[[k]][[1]]][[1,1]]; psi[[i,pos22]]=psirow[[k]][[2]];k++]; i++]; eone=IdentityMatrix[len][[1]];vv=eone;AppendTo[res11,1];For[i=1,i≤mm,vv=-(vv.psi);AppendTo[res11,PolyTruncate[vv.jac,para,mm]];i++];primitive=Total[res11]]; 

(*Output the 4-pt function*)
FourPointFunction[EQN_,TXYZ_List,JAC_List,MM_,VAR_,VARVV_]:=
 Module[{eqn=EQN,txyz=TXYZ,jac=JAC,mm=MM,var=VAR,varvv=VARVV,len,para,paravv,expotab,sca,k,lenxyz,GrConv,Groebner,conv,MatDual,basismat,res11,dualjac,red,flatvar,coe22,flatvartab,uutab,vvtab,pos22,i,j,coe33,tab33,res},
 len=Length[jac];para=Table[(Subscript[var, j]),{j,1,len}];paravv=para/.var->varvv;
expotab=Table[Exponent[jac[[i]],txyz],{i,1,len}];sca=Sum[(para.jac)^k/(k! t^k),{k,mm}]; 
lenxyz=Length[txyz];GrConv=conversionMatrix[JacIdeal[eqn,txyz],txyz];Groebner=GrConv[[1]];conv=GrConv[[2]];  MatDual=Metric[eqn,txyz,jac,Groebner];dualjac=MatDual[[2]];basismat=MatDual[[1]];
red=CohomologyReduction[sca,txyz,jac,Groebner,dualjac,len,conv,lenxyz]; 
 flatvar=PolyTruncate[Coefficient[t^mm*red,t^(mm-1)],para,2];
coe22=PolyTruncate[Coefficient[t^mm*red,t^(mm-2)],para,3];
flatvartab=CoefficientRules[flatvar,txyz];
uutab=Table[0,{i,1,len}];For[i=1,i≤len,pos22=Position[expotab,flatvartab[[i,1]]][[1,1]];uutab[[pos22]]=flatvartab[[i,2]];i++];vvtab=(2para-uutab)/.var->varvv;  
res11=CoefficientRules[PolyTruncate[Expand[coe22/.Table[para[[i]]->vvtab[[i]],{i,1,len}]],paravv,3],txyz]; 
tab33=Table[0,{i,1,len}]; For[i=1,i≤len,coe33=res11[[i]];pos22=Position[expotab,coe33[[1]]][[1,1]]; 
tab33[[pos22]]=PolyPickup[coe33[[2]],paravv,3];i++];
res=Expand[-1*(paravv.basismat).tab33/4];res];

 
(*Print out final ouput. SCAL is used to take a multiplication of the 4-point function*)

Output[TT_List,POS_,TXYZ_List,MM_,VAR_,VARVV_]:= Module[{mm=MM, tt=TT,k,pos=POS,var=VAR, varvv=VARVV,mat,i,j, len,para, txyz=TXYZ,jac,eqn,deg},jac=JacWWBasis[txyz][[pos]];eqn=DefiningEquation22[tt,txyz][[pos]];len=Length[jac];deg=DegtxyzOne[eqn,txyz];para=Table[(Subscript[var, j]),{j,1,len}];i=pos;Print["Type: ",SingularType[[i]]];Print["Defining Equation is: ",eqn];
Print["Degree of variables", txyz, " : ", deg];
Print["Monomial Basis= ", jac]; Print["Degree of the monomial basis: " , Table[Exponent[jac[[i]],txyz].deg,{i,len}]]; Print["Parameters= ",para]; Print["deltaWW= ", para.jac];Print["Up to order ", mm, ", the primitive form is: "];Print[" "];Print[ PrimitiveForm[eqn,txyz,jac,mm,var]];Print[" "];  Print["The four point function multiplying with (-1)  is given by: "];Print[" "];Print[ Expand[FourPointFunction[eqn,txyz,jac,mm,var,varvv]]]];

(*tt={a, b, c} denotes the normalization. t is used to generate the flat coordinates here *)
tt={1,1,1};txyz={x, y,z};
For[pos=1,pos≤24,
Output[tt,pos,txyz, 3,u,t];pos++];

,

(*Output the 5-pt function for U12c*)
FivePointFunction[EQN_,TXYZ_List,JAC_List,MM_,VAR_,VARVV_,PRIM_]:=
 Module[{eqn=EQN,txyz=TXYZ,jac=JAC,mm=MM,var=VAR,varvv=VARVV,prim=PRIM,len,para,paravv,expotab,sca,k,lenxyz,GrConv,Groebner,pr22,conv,MatDual,basismat,res11,dualjac,red,flatvar,coe22,flatvartab,uutab,vvtab,pos22,i,j,coe33,tab33,res,uut22,uut33,uv22,uv33,uv25},
 len=Length[jac];para=Table[(Subscript[var, j]),{j,1,len}];paravv=para/.var->varvv;
expotab=Table[Exponent[jac[[i]],txyz],{i,1,len}];pr22=PolyPickup[prim,para,2]+1;sca=Sum[(para.jac)^k*pr22/(k! t^k),{k,mm}]; 
lenxyz=Length[txyz];GrConv=conversionMatrix[JacIdeal[eqn,txyz],txyz];Groebner=GrConv[[1]];conv=GrConv[[2]];  MatDual=Metric[eqn,txyz,jac,Groebner];dualjac=MatDual[[2]];basismat=MatDual[[1]];
red=CohomologyReduction[sca,txyz,jac,Groebner,dualjac,len,conv,lenxyz]; 
 flatvar=PolyTruncate[Coefficient[t^mm*red,t^(mm-1)],para,3];
coe22=PolyTruncate[Coefficient[t^mm*red,t^(mm-2)],para,4];
flatvartab=CoefficientRules[flatvar,txyz];
uutab=Table[0,{i,1,len}];For[i=1,i≤len,pos22=Position[expotab,flatvartab[[i,1]]][[1,1]];uutab[[pos22]]=flatvartab[[i,2]];i++];uut22=Table[PolyPickup[uutab[[i]],para,2],{i,1,len}];
uut33=Table[PolyPickup[uutab[[i]],para,3],{i,1,len}];vvtab=(para-uut22)/.var->varvv;  
uv22=uut22;uv25=Expand[uv22/.Table[para[[i]]->vvtab[[i]],{i,1,len}]];uv33=Table[PolyPickup[uv25[[i]],paravv,3],{i,1,len}];
vvtab=((para-uut22-uut33)/.var->varvv)-uv33;
res11=CoefficientRules[PolyTruncate[Expand[coe22/.Table[para[[i]]->vvtab[[i]],{i,1,len}]],paravv,4],txyz]; 
tab33=Table[0,{i,1,len}]; For[i=1,i≤len,coe33=res11[[i]];pos22=Position[expotab,coe33[[1]]][[1,1]]; 
tab33[[pos22]]=PolyPickup[coe33[[2]],paravv,4];i++];
res=Expand[(paravv.basismat).tab33/5];res];


(*Print out final ouput. SCAL is used to take a multiplication of the 4-point function*)

Output55[TT_List,POS_,TXYZ_List,MM_,VAR_,VARVV_]:= Module[{mm=MM, tt=TT,k,pos=POS,var=VAR, varvv=VARVV,mat,i,j, len,para, txyz=TXYZ,jac,eqn,deg,prim},jac=JacWWBasis[txyz][[pos]];eqn=DefiningEquation22[tt,txyz][[pos]];len=Length[jac];deg=DegtxyzOne[eqn,txyz];para=Table[(Subscript[var, j]),{j,1,len}];i=pos;Print["Type: ",SingularType[[i]]];Print["Defining Equation is: ",eqn];
Print["Degree of variables", txyz, " : ", deg];
Print["Monomial Basis= ", jac]; Print["Degree of the monomial basis: " , Table[Exponent[jac[[i]],txyz].deg,{i,len}]]; Print["Parameters= ",para]; Print["deltaWW= ", para.jac];Print["Up to order ", 3, ", the primitive form is: "];Print[" "];
prim=PrimitiveForm[eqn,txyz,jac,3,var];Print[ prim];Print[" "];  Print["The four point function multiplying with (-1)  is given by: "];Print[" "];Print[ Expand[FourPointFunction[eqn,txyz,jac,3,var,varvv]]];Print["The five point function is given by: "];Print[" "];Print[ Expand[FivePointFunction[eqn,txyz,jac,mm,var,varvv,prim]]]];

(*tt={a, b, c} denotes the normalization. t is used to generate the flat coordinates here *)
tt={1,1,1};txyz={x, y,z};
For[pos=21,pos≤21,
Output55[tt,pos,txyz, 4,u,t];pos++];