(******************************************************************************
 *	LoopExp.m - expands one loop lattice integrals around a=0	      *
 ******************************************************************************)


BeginPackage["LoopExp`"];


(* define public functions *)
LoopExp::usage = "LoopExp[o, a, nc, bl, fl, fe, fp]
o: order in a,
a: lattice spacing,
nc: name for numerical constant,
bl: list of boson propagators,
fl: list of fermion propagators,
fe: list of extra fermion propagators in the denominator,
fp: name for Feynman parameters"

HardStart::usage = "HardStart[bs]
Specify boson index bs from where to start solving recurrence relations.
Must be called before SolveHard."

SolveHard::usage = "SolveHard[br, fr, nr]
sove recurrence relations over index range {br, fr, nr}"

SaveHard::usage = "SaveHard[] : save solutions for recurrence relations"
LoadHard::usage = "LoadHard[] : load solutions for recurrence relations"

(* recurrence relations and regulator are public *)
F;
\[Delta];


Begin["`Private`"];
Needs["Recurrence`"];


(* some small functions *)
listAnd[{i__}]:=And[i];
sq[x_]:=x.x;
PlusToList[x_]:=If[MatchQ[Head[x],Plus],List@@x,{x}];


(***
 * solve kontinuum Feynman integrals using Feynman parameters
 *)

(* main function, integrates a single summand *)
fint[m_,p_,e_,v_,fp_] := Module[
		{ret,i,esum,numer,denom,param,exp,np,x,dim,cl,v2,j,dir},

	np=Length[m];
	dim=Length[p[[1]]];
	esum=Sum[e[[i]],{i,np}];

	x=Array[dir[#]&,dim];
	param=Array[fp[#]&,np];
	param[[np]]=1-Sum[fp[i],{i,np-1}];
	denom=0;
	numer=1;
	For[i=1,i<=np,i++,
		denom+=param[[i]]*(sq[x+p[[i]]]+m[[i]]^2);
		numer*=param[[i]]^(e[[i]]-1);
	];
	cl=(Coefficient[denom,#,1]&/@x)/2;
	denom=Expand[denom-sq[x+cl]];
	If[denom==0,Return[0];];

	ret=PlusToList[Expand[Product[(x[[i]]-cl[[i]])^v[[i]],{i,dim}]]];
	For[j=1,j<=Length[ret],j++,
		v2=Exponent[ret[[j]],x];
		exp=esum-(Sum[v2[[i]],{i,dim}]+dim)/2;
		If[!listAnd[EvenQ/@v2],
			ret[[j]]=0;
			Continue[];
		];
		ret[[j]]=Together[ret[[j]]/Product[x[[i]]^v2[[i]],{i,dim}]]
			*FunctionExpand[Gamma[exp]
			/Product[Gamma[e[[i]]],{i, np}]
			*Product[Gamma[(1+v2[[i]])/2],{i,dim}]]
			*numer/denom^exp;
	];

	Return[Plus@@ret];
];

(* rules to split sums of integrals *)
fint/:fint[m_,p_,e1_,v_]^(e2_):=fint[m,p,e1*e2,v*e2];
fint/:fint[{m1__},{p1__},{e1__},v1_]*fint[{m2__},{p2__},{e2__},v2_]:=
		fint[{m1,m2},{p1,p2},{e1,e2},v1+v2];
intv/:intv[v_]^e_:=intv[v*e];
intv/:intv[v1_]*fint[m_,p_,e_,v2_]:=fint[m,p,e,v1+v2];
unit[i_,dim_]:=Module[{j},Table[Which[i==j,1,True,0],{j,dim}]];

(* integrate sums of Feynman intregrals using fint *)
fIntegrate[int_,prop_,x_,dim_,fp_]:=(int
		/.prop[m_,p_]:>fint[{m},{p},{-1},Array[0&,dim]]
		/.(x[[#]]->intv[unit[#,dim]]&/@Range[dim])
		/.fint[opt__]:>fint[opt,fp]);


(***
 * Functions to solve, save and restore recurrence relations
 *)

(* Wilson parameter *)
r=1;
(* dimension *)
dim=2;
(* regulator *)
reg=\[Delta];

HardStart[bs_]:=Module[{a1,a2,b,f,masters,recEq},
	recEq={	b[-1]+a1[1]+a2[1]-2,
		f[-1]+a1[1]+a2[1]-2+(a1[1]-1)^2+(a2[1]-1)^2
			-r^2*(a1[1]+a2[1]-2)^2,
		(f[-1]+a1[1]*(a1[1]-1)+a2[1]*(a2[1]-1))*b[2]-r^2,
		1+2f*a1[1]*(a1[1]-1)*f[1]*(2*a1[1]-1-2*r^2*(a1[1]+a2[1]-2))
			+2*a1*(a1[1]-1)+2*(b+reg)*a1[1]*(a1[1]-1)*b[1],
		1+2f*a2[1]*(a2[1]-1)*f[1]*(2*a2[1]-1-2*r^2*(a1[1]+a2[1]-2))
			+2*a2*(a2[1]-1)+2*(b+reg)*a2[1]*(a2[1]-1)*b[1]};

	masters={{-1,0,1,1},{0,0,1,1},{-1,1,1,1},{-1,1,1,2}};
	DefineRecurrence[F,recEq,{b,f,a1,a2},Priority->{{1},{3,4},{2}},
		Symmetric->{{3,4}},Masters->masters,
		Start->{bs,0,0,0},Directory->"."];
];

SolveHard[b_,f_,a_]:=SolveRecurrence[F,{b,a,f}];
SaveHard[]:=SaveRecurrence[F,"./F/F.m"];
LoadHard[]:=LoadRecurrence[F,"./F/F.m"];


(***
 * functions to save soft part on HD
 *)
sname="soft";
mySave[f_,e_]:=Put[e,Evaluate[f]];
myLoad[f_]:=Get[Evaluate[f]];
mySuffix=".m";
fileName[s_]:=ToFileName[{".",sname},StringJoin[s,mySuffix]];
saveFile[s_,ex_]:=mySave[fileName[s],ex];
loadFile[s_]:=myLoad[fileName[s]];
maxFilesPerDir=500;
fileName[s_,n_]:=ToFileName[{".",sname,StringJoin[s,
		ToString[Floor[n/maxFilesPerDir]*maxFilesPerDir]]},
		StringJoin[s,ToString[n],mySuffix]];
saveFile[s_,n_,ex_]:=Module[{dname,fname},
	fname=fileName[s,n];
	dname=DirectoryName[fname];
	If[!MatchQ[FileType[dname],Directory],
		CreateDirectory[dname];
	];
	mySave[fname,ex];
];
loadFile[s_,n_]:=myLoad[fileName[s,n]];


(***
 * own implementation of Series[], saves resources
 *)
softSeries[exp_,reg_,a_,amax_]:=Module[{l,i},
	l=PlusToList[Expand[exp]];
	For[i=1,i<=Length[l],i++,
		l[[i]]=Cancel[reg*l[[i]]/.reg->0]/reg
			+Cancel[D[reg*l[[i]],reg]/.reg->0];
	];
	l=CoefficientList[Plus@@l,a];
	l=Together[Take[l,Min[amax+1,Length[l]]]];
	For[i=Length[l],i>=1,i--,
		If[!MatchQ[l[[i]],0],Break[];];
		l=Delete[l,i];
	];
	Return[l];
];


(***
 * main function, expands soft and hard part
 *)
LoopExp[o_,a_,nc_,bl_,fl_,fe_,fp_] := Module[
		{coeff,nb,nf,x,p,ma,soft,hard,h,i,j,bprop,pos,coList,dexp,
		p1,p2,minpot,H,shift,numfp,div,masters,nil,DF,DB,sum,
		k,zeros,lsum,intc,hardtime,softtime,simptime},

	(* number of boson propagators *)
	nb=Sum[bl[[j,3]],{j,Length[bl]}];
	(* number of fermion propagators *)
	nf=Sum[fl[[j,3]],{j,Length[fl]}];
	(* number of 'extra' fermion propagators in denominator *)
	ne=Length[fe];
	(* leading order of hard part is a^shift *)
	shift=2*(nb+nf-ne)-dim;
	(* total leading order is a^minpot *)
	minpot=Min[0,shift];

	(* prepare some variables *)
	x=Array[Unique[]&,dim];
	coeff=Array[0&,o-minpot+1];
	saveFile["coeff",coeff];
	nil=Array[0&,dim];
	pi=1;
	setP[p_]:=a*p;

	If[o<minpot,Return[SeriesData[a,0,{},minpot,minpot,1]];];

	(* determine propagator for the regulator *)
	pos=Select[Join[bl,fl],#[[3]]>0&,1];
	If[Length[pos]==0,
		p=nil;
		ma=0;
	,
		p=setP[pos[[1,2]]];
		ma=a*pos[[1,1]]/2;
	];

	(***
	 * prepare soft and hard
	 *)
	(* introduce regulator and measure *)
	soft=Product[1/(1+h*x[[j]]^2),{j,dim}]*
		(h*Sum[(p[[j]]+x[[j]])^2/((1+h*x[[j]]^2)(1+h*p[[j]]^2))
			-(p[[j]]+x[[j]])^2,{j,dim}]
		+bprop[ma,p])^(-reg);
	hard=H[Join[{0,0},Array[1&,dim]]]*DB[ma,p,h]^(-reg);

	(* boson propagators *)
	For[i=1,i<=Length[bl],i++,
		p=setP[bl[[i,2]]];
		ma=a*bl[[i,1]]/2;
		If[bl[[i,3]]>0,
			soft*=(h*Sum[(p[[j]]+x[[j]])^2
					/((1+h*x[[j]]^2)*(1+h*p[[j]]^2))
					-(p[[j]]+x[[j]])^2,{j,dim}]
				+bprop[ma,p])^(-bl[[i,3]]);
		,
			soft*=(Sum[(p[[j]]+x[[j]])^2
					/((1+h*x[[j]]^2)*(1+h*p[[j]]^2))
					,{j,dim}]+ma^2)^(-bl[[i,3]])
		];
		hard*=DB[ma,p,h]^(-bl[[i,3]]);
	];

	(* fermion propagators *)
	For[i=1,i<=Length[fl],i++,
		p=setP[fl[[i,2]]];
		ma=a*fl[[i,1]]/2;
		If[fl[[i,3]]>0,
			soft*=(h*Sum[((x[[j]]+p[[j]]
				-h*x[[j]]^2*p[[j]]-x[[j]]*h*p[[j]]^2)
				/((1+h*x[[j]]^2)*(1+h*p[[j]]^2)))^2
				-(p[[j]]+x[[j]])^2,{j,dim}]
				+(ma+r*h*Sum[(p[[j]]+x[[j]])^2
				/((1+h*x[[j]]^2)(1+h*p[[j]]^2)),{j,dim}])^2
				+bprop[ma,p]-ma^2)^(-fl[[i,3]]);
		,
			soft*=(h*Sum[((x[[j]]+p[[j]]
				-h*x[[j]]^2*p[[j]]-x[[j]]*h*p[[j]]^2)
				/((1+h*x[[j]]^2)*(1+h*p[[j]]^2)))^2,
				{j,dim}]
				+(ma+r*h*Sum[(p[[j]]+x[[j]])^2
				/((1+h*x[[j]]^2)(1+h*p[[j]]^2)),{j,dim}])^2
				)^(-fl[[i,3]]);
		];
		hard*=DF[ma,p,h]^(-fl[[i,3]]);
	];

	(* 'extra' fermion propagators *)
	For[i=1,i<=ne,i++,
		p1=setP[fe[[i,1]]];
		p2=setP[fe[[i,2]]];
		soft*=Sum[(x[[j]]+p1[[j]]
				-h*x[[j]]^2*p1[[j]]-x[[j]]*h*p1[[j]]^2)
				/((1+h*x[[j]]^2)*(1+h*p1[[j]]^2))
			*(x[[j]]+p2[[j]]
				-h*x[[j]]^2*p2[[j]]-x[[j]]*h*p2[[j]]^2)
				/((1+h*x[[j]]^2)*(1+h*p2[[j]]^2)),
			{j,dim}];
		hard*=Sum[(x[[j]]+h*p1[[j]]
				-x[[j]]^2*h*p1[[j]]-x[[j]]*h^2*p1[[j]]^2)
				/((1+x[[j]]^2)*(1+h^2*p1[[j]]^2))
			*(x[[j]]+h*p2[[j]]
				-x[[j]]^2*h*p2[[j]]-x[[j]]*h^2*p2[[j]]^2)
				/((1+x[[j]]^2)*(1+h^2*p2[[j]]^2)),
			{j,dim}];
	];

	(***
	 * expand soft in h until no new terms of order o in a emerge
	 *)
	softtime=Timing[
	Print["Expanding soft part..."];
	intc=1;
	For[i=1,True,i++,
		saveFile["soft",soft];
		sum=soft/.h->0;
		soft=.;
		sum=PlusToList[Expand[sum]];
		lsum=Length[sum];
		For[j=1,j<=lsum,j++,
			saveFile["soft",j,sum[[j]]];
		];
		sum=.;
		zeros=0;
		For[j=1,j<=lsum,j++,
			coList=loadFile["soft",j];
			coList=fIntegrate[coList,bprop,x,dim,fp];
			(*saveFile["int",intc,coList];
			intc++;*)
			coList=softSeries[coList*a^(shift-minpot),
						reg,a,o-minpot];
			If[Length[coList]==0,
				zeros++;
				Continue[];
			];
			coeff=loadFile["coeff"];
			For[k=1,k<=Length[coList],k++,
				coeff[[k]]=Collect[coList[[k]]+coeff[[k]],
							reg,Together];
			];
			saveFile["coeff",coeff];
			Clear[coeff];
		];
		(*Print[intc-1];*)
		coList=.;
		If[zeros==lsum,Break[];];
		soft=loadFile["soft"];
		soft=D[soft,h]/i;
	];
	soft=.;
	coList=.;

	(* Laurent expansion in the regulator *)
	coeff=loadFile["coeff"];
	For[i=Max[0,shift]+1,i<=Length[coeff],i++,
		dexp=Exponent[coeff[[i]],reg,Min];
		If[dexp>=0,Continue[];];
		div=Together[Coefficient[coeff[[i]],reg,dexp]];
		numfp=Max[Cases[div,fp[_],\[Infinity]]
				/.fp->Identity];
		If[numfp<=0,Continue[];];
		coeff[[i]]-=div*reg^dexp;
		For[j=1,j<=numfp,j++,
			div=Integrate[div,{fp[j],0,1}];
		];
		coeff[[i]]=Collect[coeff[[i]]+div*reg^dexp,reg,Together];
	];
	][[1]];
	Print["    ",softtime];


	(***
	 * expand hard
	 *)
	hardtime=Timing[
	Print["Expanding hard part..."];
	If[o<shift,Goto[exit];];

	Derivative[0,nil,1][DB][m_,p_,h_]:=
		D[Sum[(h*p[[j]]+x[[j]])^2/((1+x[[j]]^2)(1+h^2*p[[j]]^2)),
			{j,dim}]+h^2*m^2,h];
	Derivative[0,nil,1][DF][m_,p_,h_]:=
		D[Sum[((x[[j]]+h*p[[j]]
			-x[[j]]^2*h*p[[j]]-x[[j]]*h^2*p[[j]]^2)
			/((1+x[[j]]^2)*(1+h^2*p[[j]]^2)))^2,{j,dim}]
		+(h*m+DB[0,p,h])^2,h];
	DB[_,_,0]=H[Join[{-1,0},nil]];
	DF[_,_,0]=H[Join[{0,-1},nil]];
	H/:H[i_]^e_:=H[i*e];
	H/:H[i1_]*H[i2_]:=H[i1+i2];

	For[i=0,i<=o-shift,i++,
		If[i!=0,hard=D[hard,h]/i;];
		coList=Expand[hard/.h->0];
		coList=coList/.(x[[#]]^e_:>(H[Join[{0,0},-unit[#,dim]]]-1)
				^Quotient[e,2]*x[[#]]^Mod[e,2]& 
					/@ Range[dim])
			/.(#->0& /@ x);
		coList=Expand[a^(shift-minpot)*coList];
		coList=Series[coList,{a,0,o-minpot}];
		coList=CoefficientList[coList,a]
				/.H[{b_,f__}]:>H[{b-reg,f}];
		For[j=1,j<=Length[coList],j++,
			coeff[[j]]+=coList[[j]];
		];
	];
	hard=.;
	coList=.;

	(* identify integrals that must be reduced to master integrals *)
	masters=Masters[F];
	hInd=Transpose[Join[Flatten[Indizes[H,coeff],1],masters]];
	Print["    Index Range: ",
		{Min[hInd[[1]]],Min[hInd[[2]]],Min[hInd[[3]]+hInd[[4]]]},
		" to ",
		{Max[hInd[[1]]],Max[hInd[[2]]],Max[hInd[[3]]+hInd[[4]]]}];
	coeff=coeff/.H->F; (* use recurrece relations *)

	(* expand master integrals in regulator *)
	dexp=Min[Exponent[coeff,reg,Min]];
	coeff=coeff/.Array[
		(F[masters[[#]]]->Sum[M[#,j]*reg^j,{j,0,-dexp}])&,
		Length[masters]];

	(* substitute known master integrals und relations between them *)
	M[1,0]=\[Pi]^2;
	M[2,0]=\[Pi]^2;
	M[1,1]=2*\[Pi]-\[Pi]^2+M[2,1];
	M[2,2]=\[Pi]^2+6*\[Pi]*(Log[2]-1)+M[1,2]+M[1,1];
	M[3,0]=4*\[Pi]^2/(3*Sqrt[3]);
	M[4,0]=(-18+(-27+32*Sqrt[3])*\[Pi])*\[Pi]/54;
	M[2,1]=nc+8/3*M[3,1]-2*M[4,1];

	(* try to eliminate regulator (the regulator MUST cancel everywhere) *)
	coeff=Normal[Series[coeff,{reg,0,0}]];
Label[exit];
	][[1]];
	Print["    ",hardtime];


	(***
	 * finally: simplify
	 *)
	simptime=Timing[
	Print["Simplify..."];
	(* TODO: p -> Tan[p*a/2]/a and expand *)
	coeff=SeriesData[a,0,
		Simplify[coeff/2^shift,Assumptions->{a>0}],
		minpot,o+1,1];
	][[1]];
	Print["    ",simptime];


	Return[coeff];
];


End[(*"`Private`"*)];
EndPackage[(*LatRec`*)];
