(* ::Package:: *)

Print["MBsums v1.0.3 (14 Dec 2015) by Michal Ochman"];
Print["The author would like to thank Tord Riemann for many fruitful discussions."];
Print["Available at http://prac.us.edu.pl/~gluza/ambre/ -> Related and auxiliary Software"];
Print["Please cite as follows: M. Ochman, T. Riemann,"]; 
Print["``MBsums — a Mathematica package for the representation of Mellin–Barnes integrals by multiple sums'',"];
Print["Acta Physica Polon. 46 (2015) 2117, extended version: arXiv:1511.01323 "];
Print["User licence:"];
Print["CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html"];

BeginPackage["MBsums`"];

MBsumsInfo=True;

MBIntToSum;
MergeMBsums;
DoAllMBSums;

n;
MBsum;

k;
SimplifyMBsums;

Begin["`Private`"];

MergeMBsums[s1_]:=Block[{s2},

s2 = Plus @@ s1;

s2 = s2 //. {MBsum[x1__, y__, z__] + MBsum[x2__, y__, z__] :> 
     MBsum[x1 + x2, y, z], 
    a_ *MBsum[x2__, y__, z__] :> MBsum[a *x2, y, z]};
s2 = s2 /. {MBsum[x2__, y__, z__] :> MBsum[Factor[ExpandAll[x2]], y, z]};

s2 = PlusToList[s2];
Return[s2]
]

GammaArgs[vw_] := Block[{vvw,vvw2,x}, vvw = Expand[vw];
  vvw = If[ToString[Head[vvw]] == "Plus", List @@ vvw, {vvw}];
  vvw = Union[
	    Flatten[Map[Function[x1, Variables[Numerator[x1]]], vvw]]];
  vvw2 = Union[
	    Join[Cases[vvw, Gamma[x__] :> Expand[x]],
     Cases[vvw, PolyGamma[__, x__] :> Expand[x]]]];
  Return[vvw2]]

GammaArgsD[vw_] := Block[{vvw,vvw2,x}, vvw = Expand[vw];
  vvw = If[ToString[Head[vvw]] == "Plus", List @@ vvw, {vvw}];
  vvw = Union[
    Flatten[Map[Function[x1, Variables[Denominator[x1]]], vvw]]];
  vvw2 = Union[
    Join[Cases[vvw, Gamma[x__] :> Expand[x]], 
     Cases[vvw, PolyGamma[__, x__] :> Expand[x]]]];
  Return[vvw2]]


GammaArgsNeg[GammaArgsList_,VarZ_] := 
	Block[{imax,GammaArgsNegList,GammaArg},
	GammaArgsNegList={};
	imax=Length[GammaArgsList];
	
	Do[
		GammaArg=GammaArgsList[[i]];
		GammaArgsNegList = Append[GammaArgsNegList,If[Coefficient[GammaArg,VarZ]!=0,{GammaArg},{}]],
	{i,imax}];

	Return[Flatten[GammaArgsNegList]]]


NoZs[vw_, vVarZList_, imax_] :=  
	(* imax = Lenght[vVarZList] *)

	Block[{rN},

	rN=True;

	Do[
		If[Coefficient[vw,vVarZList[[i1]]]!=0,rN=False;Return[],0]
		,{i1,imax}];

Return[rN]
]


NoFrac[vw_] :=

	Block[{vvw},

	vvw=Expand[vw];
	vvw=If[ToString[Head[vvw]] == "Plus",List@@vvw,{vvw}];

	vvw=And@@Map[Function[x1,Abs[Denominator[x1]]==1],vvw];

Return[vvw]
]


LLAdd[l1_,l2_] := 
	Block[{imax,lout},
	imax=Length[l1];
	lout={};
	Do[lout=Append[lout,ToExpression[StringJoin[l1[[i]],l2[[i]]]]],{i,imax}];
	Return[lout]
]


nArgs1[NegArgsAll_, nTo0_] :=

	Block[{vvw},
(*	vvw = Select[NegArgsAll,Function[x,NoZs[x,vVarZList, imax]]];*)
	vvw = Select[NegArgsAll,Function[x,NumberQ[ExpandAll[x/.nTo0]]]];(* IntegerQ[n/.nTo0]*)
	vvw = Flatten[Map[PlusToList,vvw]];

	vvw = Select[vvw,Function[x,!NumberQ[x]]];
	vvw = Select[vvw,Function[x1,Abs[Denominator[x1]]=!=1]];

	vvw = Gather[vvw,Function[{x,y},Variables[x]===Variables[y]]];

Return[vvw]
]

PlusToList[vw_] :=

	Block[{vvw},

	vvw=Expand[vw];
	vvw=If[ToString[Head[vvw]] == "Plus",List@@vvw,{vvw}];

Return[vvw]
]

TimesToList[vw_] :=

	Block[{vvw},

	vvw=Expand[vw];
	vvw=If[ToString[Head[vvw]] == "Times",List@@vvw,{vvw}];

Return[vvw]
]

nArgs2[vw_] :=

	Block[{nv,dv},

	nv = Variables[vw[[1]]];
	dv = LCM @@Map[Function[x1,Abs[Denominator[x1]]],vw];

Return[{nv[[1]],dv}]
]




DecList[vw_] :=

	Block[{vvw,vPowList,vBase},

	vPowList=vw[[2]];
	vBase=vw[[1]]; 
	vvw=Map[Function[x1,{vBase,x1}],vPowList];

Return[vvw]
]


Coef1[vw_,vVarZ_] :=

	Block[{vc},

	vc=Coefficient[vw[[2]],vVarZ];

Return[{vw[[1]]^vc,vVarZ}]
]

DetCont2[IntXe2_,vVarZ_,vKin_] :=

	Block[{c},

	c = IntXe2;
	c=TimesToList[c];
	c=Cases[c,x_^y_:>{x,y}];
	c=Map[Function[x,{x[[1]],PlusToList[x[[2]]]}],c];
	c=Flatten[Map[DecList,c],1];
	c=Select[c,Function[x,Variables[x[[2]]]==={vVarZ}]];
	c=Map[Function[x,Coef1[x,vVarZ]],c];
	c=Map[Function[x,x[[1]]],c];
	c=Times@@c;
	c=c/.vKin;
	If[!NumberQ[N[c]],If[MBsumsInfo==True,
		Print["Found c = ",c," (not a number): please complete kinematic's list."],0];
		Return[{"R","L"}],0];

	If[Abs[c]>1,Return[{"L"}],0];
	If[Abs[c]<1,Return[{"R"}],0];
	If[Abs[c]==1,Return[{}],0]
]


DetCont[IntXe_,vVarZ_,vKin_] :=

	Block[{c},

	c = IntXe[[1]];
	c = PlusToList[c];
	Return[Flatten[Map[Function[x,DetCont2[x,vVarZ,vKin]],c]]]
]


CalcRes[IntXelem_,IntR_,VarZ_,CondVarZ_,nn_,nnRes_,ConditionsN_,nList_,NegArg_,NegArgsD_, nTo0_,epsTo0_] := 
	Block[{DoNotCont,DoCont,ConditionsRes,ConditionsRes0,ResSolved,rmax,w,endR,
			NegArgsDx,TuplesN,CasesN,CasesNv,wRes,ConditionsEv,x,y},

	(* Do not Simplify[any conditions] before Residue[], Expand[] is ok *)

	DoNotCont = NegArg[[2]]; 
	DoCont = NegArg[[1]];

	ConditionsRes=ConditionsN;
	ResSolved = ExpandAll[Solve[DoCont[[1]]==-nnRes,VarZ]];
	ConditionsRes0=ExpandAll[(CondVarZ/.ResSolved[[1]]/.Flatten[IntR]/.epsTo0)&&IntXelem[[2]]];
	
		rmax=Length[DoCont];

		endR=False;
		(* Create conditions for args. that do contribute *)
		Do[
			w=Expand[DoCont[[r]]/.ResSolved[[1]]];
			If[IntegerQ[ExpandAll[w/.nTo0]],
				ConditionsRes = ConditionsRes&& w<=0,endR=True;Break[]],
				(* &&Element[w, Integers] *)
			{r,rmax}];

		If[endR==True,Return[{{0,False}}],0];

		(* Create conditions for args. that do not contribute *)

		rmax=Length[DoNotCont];

		Do[
			w=Expand[DoNotCont[[r]]/.ResSolved[[1]]];
			If[IntegerQ[ExpandAll[w/.nTo0]],
				ConditionsRes = ConditionsRes&&w>=1,0],
			{r,rmax}];


		ConditionsRes =ExpandAll[ ConditionsRes];
		(* Check conditions *)
		If[ToString[FindInstance[ConditionsRes&&ConditionsRes0,nList, Integers]]==ToString[{}],
			Return[{{0,False}}],0];

		
		NegArgsDx=NegArgsD/.ResSolved[[1]];	
		NegArgsDx=Select[NegArgsDx,Function[w,IntegerQ[ExpandAll[w/.nTo0]]]];
		(* &&Element[w, Integers] *)

		rmax=Length[NegArgsDx];

		If[rmax==0,(*Print[DoCont," *** ",ResSolved[[1]]," *** ",ConditionsRes];*)		
			Return[{{Assuming[ConditionsRes&&Element[nList, Integers] ,
						Residue[IntXelem[[1]],{VarZ,ResSolved/.{{x_->y_}}:>y}]],
						ConditionsRes&&ConditionsRes0}}],
			(* ** *)
			TuplesN=Tuples[{">=1","<=0"},rmax];
			NegArgsDx=Map[Function[x,ToString[x,InputForm]],NegArgsDx];
			CasesN=Map[Function[t1,LLAdd[NegArgsDx,t1]],TuplesN];
			rmax=Length[CasesN];
			wRes={};
			Do[
				CasesNv=Expand[CasesN[[r]]];
				CasesNv=If[Length[CasesNv]>1,And@@CasesNv,CasesNv[[1]]];
				ConditionsEv=ExpandAll[ConditionsRes&&CasesNv];
				If[ToString[FindInstance[ConditionsEv&&ConditionsRes0,nList, 
					Integers]]==ToString[{}],{},
(*Print[DoCont," *** ",ResSolved[[1]]," *** (d) ",ConditionsEv];*)
				wRes = Append[wRes,{Assuming[ConditionsEv&&Element[nList, Integers] ,
						Residue[IntXelem[[1]],{VarZ,ResSolved/.{{x_->y_}}:>y}]],
						ConditionsEv&&ConditionsRes0}]],
			{r,rmax}];
			Return[wRes]
		]
	]

GGamma[x_,nTo0_]:=If[IntegerQ[x/.nTo0]&&NoFrac[x],(x-1)!,Gamma[x]]

GPolyGamma[k_,n_,nTo0_]:=If[IntegerQ[n/.nTo0]&&NoFrac[n], 
	If[k==0,
		HarmonicNumber[n-1,1]-EulerGamma,
		((-1)^k)*k!*(HarmonicNumber[n-1,k+1]-Zeta[k+1]) 
	],
	PolyGamma[k,n]
]

MBIntToSum[MBIntX_,vKin_,Intdesc_] :=
	Block[{Kin,IntX,IntR,IntR2List,noofR,VarLRList,VarZList,ConditionsN,
			ConditionsRes,kmax,IntXelem,ResSolved,ConditionsList,imax,nList,
			VarZ,nn,n,Result,jmax,x,y,IntR2ListEL,endR,CondVarZ,c1,c2,
			NegArgs,NegArgsD,NegArgsCa,NegArg,NegArg1,NegArg2,VarZCoef,
			FracCoef1,FracCoef2,a,a1,a2,resstr,epsTo0,nTo0,NegArgsAll,NegArgsDoC,
			rmax,ConditionsRes0},

	IntR = MBIntX[[2]]/.(x_->y_):>(x->Rationalize[y]) // ExpandAll;
	epsTo0=IntR[[1]];
	Kin=vKin/.(x_->y_):>(x->Rationalize[y]);

	IntR2List = IntR[[2]]/.(x_->y_):>{x,y};
	ConditionsList={};
	noofR=0;
	
		VarZList = Intdesc/.(x_->_):>x;
		VarLRList = Intdesc/.(x_->y_):>y;
		imax=Length[VarZList];

	If[vKin==={},
	(* Create conditions eg. z1 < 1/4 *)
		Do[
			IntR2ListEL=Select[IntR2List,Function[x,ToString[x[[1]]]==ToString[VarZList[[i]]]]][[1]];
			If[ToString[VarLRList[[i]]]=="L",
			ConditionsList = Append[ConditionsList,IntR2ListEL[[1]]<IntR2ListEL[[2]]],0];
			If[ToString[VarLRList[[i]]]=="R",
			ConditionsList = Append[ConditionsList,IntR2ListEL[[1]]>IntR2ListEL[[2]]];
			noofR=noofR+1,0] ,
		{i,imax}],0];

	IntX = {{MBIntX[[1]],True}};
	
	ConditionsN = True;
	endR=False;
	nList = {};
	nTo0 = {};
	Do[
		VarZ=VarZList[[i]];

		If[vKin==={},

			CondVarZ=ConditionsList[[i]];
			If[MBsumsInfo==True,
				Print[VarZ,"->",VarLRList[[i]]," ( Re ",ToString[CondVarZ,InputForm]," )"],0],

		(* Find contour based on vKin *)
			c1=Map[Function[x1,DetCont[x1,VarZ,Kin]],IntX];
			c1=Union[Flatten[c1]];
			If[Length[c1]>1,
				endR=True;
				If[MBsumsInfo==True,Print["Unable to find correct contour for ",VarZ,"."],0];
				Break[],0];
			If[c1==={},c2={ToString[VarLRList[[i]]]},c2=c1];
			IntR2ListEL=Select[IntR2List,Function[x,ToString[x[[1]]]==ToString[VarZList[[i]]]]][[1]];
			If[c2==={"L"},CondVarZ=IntR2ListEL[[1]]<IntR2ListEL[[2]],0];
			If[c2==={"R"},CondVarZ=IntR2ListEL[[1]]>IntR2ListEL[[2]]; noofR=noofR+1,0];
			If[MBsumsInfo==True,
				Print[VarZ,"->",c2[[1]]," ( Re ",ToString[CondVarZ,InputForm]," )"],0]];

		nn=ToExpression[StringJoin[{ToString[n],ToString[i]}]];
		ConditionsN = nn >=0;
		nList = Append[nList,nn];
		nTo0 =Append[nTo0,nn->0];
		Result={};
		kmax=Length[IntX];


		Do[
			IntXelem = IntX[[k]];
			(* Get arguments of Gamma and PolyGamma functions w.r.t. VarZ 
			in numerator and denominator *)
			NegArgs = GammaArgsNeg[GammaArgs[IntXelem[[1]]],VarZ]; 
			NegArgsD = GammaArgsNeg[GammaArgsD[IntXelem[[1]]],VarZ];
			NegArgsAll = Union[Join[NegArgs,NegArgsD]];

			NegArgsDoC={};
			rmax=Length[NegArgs];
			
			(* Select those arguments of Gamma and PolyGamma functions 
			that can give residue: NegArgsDoC *)
			Do[
				ConditionsRes=ConditionsN;
				ResSolved = Solve[NegArgs[[r]]==-nn,VarZ];
				ConditionsRes0=ExpandAll[(CondVarZ/.ResSolved[[1]]/.Flatten[IntR]/.epsTo0)
								&&IntXelem[[2]]];
				If[ToString[FindInstance[ConditionsRes&&ConditionsRes0,nList, 
					Integers]]==ToString[{}],0,NegArgsDoC=Append[NegArgsDoC,NegArgs[[r]]]],				  {r,rmax}];
	
			NegArgsCa=Reverse[Map[Function[x1, {x1, Complement[NegArgs, x1]}], Subsets[NegArgsDoC]]];
 			jmax=Length[NegArgsCa];
		(*	Print[k];*)
			Do[
				NegArg = NegArgsCa[[j]];
				If[NegArg[[1]]=={},Continue[],0];
				NegArg1 =NegArg[[1]];
				NegArg2 =NegArg[[2]];
				NegArg1 = SortBy[NegArg1, Function[x1, Abs[Coefficient[x1, VarZ]]]];
				NegArg={NegArg1,NegArg2};

				(* Correct for nn/M terms *)
				VarZCoef = Ceiling[Abs[Coefficient[NegArg1[[1]],VarZ]]];
				
				ResSolved = Solve[NegArg1[[1]]==-VarZCoef*nn,VarZ];
				FracCoef1 = Map[Function[x1, Denominator[Abs[Coefficient[x1, nn]]]],
						NegArgsAll/.ResSolved[[1]]];
				VarZCoef=LCM@@Append[FracCoef1,VarZCoef];


				(* *** *)
				ResSolved = Solve[NegArg1[[1]]==-VarZCoef*nn,VarZ];
	
				a=nArgs1[NegArgsAll/.ResSolved[[1]]/.epsTo0, nTo0];



				If[a==={},


					Do[
						Result = Append[Result,	CalcRes[IntXelem,IntR,VarZ,CondVarZ,
								nn,VarZCoef*nn+vzc-1,ConditionsN,nList,
								NegArg,NegArgsD,nTo0,epsTo0]],
						{vzc,VarZCoef}],

					(* *** *)

					a=Map[nArgs2,a]; (* {{n2, 2}, {n4, 6}} *)

					a1=Map[Function[x,x[[1]]->x[[1]]*x[[2]]+
							ToExpression["k"<>ToString[x[[1]]]]],a];
					a2=Map[Function[x,{ToExpression["k"<>ToString[x[[1]]]],0,x[[2]]-1}],a];

					a1=ToString[a1];
					a2=StringReplace[ToString[a2],{"{{"->",{","}}"->"}]"}];
					(* a2 = '',{kn2, 0, 1}, {kn4, 0, 5}]'' *)

					resstr="MBsums`Private`Result = Append[MBsums`Private`Result,
						MBsums`Private`CalcRes[
						MBsums`Private`IntXelem/."<>a1<>",
						MBsums`Private`IntR,
						MBsums`Private`VarZ,
						MBsums`Private`CondVarZ,
						MBsums`Private`nn,
						MBsums`Private`VarZCoef*MBsums`Private`nn+vzc-1,
						MBsums`Private`ConditionsN,
						MBsums`Private`nList,
						MBsums`Private`NegArg/."<>a1<>",
						MBsums`Private`NegArgsD/."<>a1<>",
						MBsums`Private`nTo0,MBsums`Private`epsTo0]]";

					ToExpression["Do["<>resstr<>",{vzc,1,MBsums`Private`VarZCoef}"<>a2]
				],
				{j,jmax}],

			{k,kmax}];

		IntX=Flatten[Result,1];
		IntX=Select[IntX,Function[x,ToString[x[[1]],InputForm]=!=ToString[0,InputForm]]];
		IntX=Map[Function[x,{Factor[ExpandAll[x[[1]]]],
				Assuming[Element[nList, Integers] ,Simplify[x[[2]]]]}],IntX],
		{i,imax}];
	
If[endR==True,Return[{}],0];

IntX=Map[Function[x,Append[x,nList]],IntX];
IntX=Map[Function[x,{Factor[x[[1]]]*(-1)^noofR,x[[2]],x[[3]]}],IntX];
IntX=Map[Function[x,MBsum@@{Factor[ExpandAll[ExpandAll[x[[1]]]/.{PolyGamma[k_,n_]:>GPolyGamma[k,n,nTo0],Gamma[n_]:>GGamma[n,nTo0]}]],
			x[[2]],x[[3]]}],IntX];

Return[IntX]
]


DoMBSum[x_, nmax_] :=
 	Block[{Result,cond,nList},
  
  	cond = x[[2]];
  	nList = x[[3]];
 If[nList==={}, If[cond==True, Return[x[[1]]], 0],
  	Result = Map[Function[x1, {x1, 0, nmax}], nList];
  	Result = Reverse[Append[Result, If[cond,ExpandAll[ x[[1]]], 0]]];
  
  Return[Sum @@ Result]]]


DoAllMBSums[x_,nmax_,SubstList_] := Plus@@ParallelMap[Function[y,DoMBSum[y/.SubstList,nmax]],x]// Expand



(* u.m *)

AndToList[vw_]:=Block[{vvw},vvw=LogicalExpand[vw];
vvw=If[ToString[Head[vvw]]=="And",List@@vvw,{vvw}];
Return[vvw]]

OrToList[vw_]:=Block[{vvw},vvw=LogicalExpand[vw];
vvw=If[ToString[Head[vvw]]=="Or",List@@vvw,{vvw}];
Return[vvw]]



f1[s1_] := Block[{conall},
  conall = Reduce[s1[[2]], s1[[3]], Integers];
  conall = OrToList[conall];
  Return[Map[Function[x, MBsum[s1[[1]], x, s1[[3]]]], conall]]]
 (* f1 : MBsum -> (Reduce) -> {MBsum[C],MBsum[C],...}*)



condRules[cond_,nList_] :=
 	Block[{d,sol,tosolve,tov},
  
	d=AndToList[cond];

 	tosolve = Cases[d, x__ == y__];

	(*{C[1] == n1, C[1] + C[2] == n2, 1 + C[2] + C[3] == n3}*)

	tov = tosolve /. x__ == y__ :> {x, y} // Flatten;
	sol=Solve[tosolve, Intersection[Variables[tov], nList]][[1]];
	(*If[Length[sol]>1,Print["Error"],0];*)

  Return[sol]]





f2[xv_] :=
 	Block[{Result,cond,nList,n2c},
  
  	cond = xv[[2]];
  	nList = xv[[3]];
  	n2c=condRules[cond,nList];
     	nList = Variables[nList/.n2c];
	cond = Refine[cond/.n2c,Element[nList,Integers]];

  Return[MBsum[xv[[1]]/.n2c,cond,nList]]] (* f2 : subst. n1->C[1]+1,... in MBsum in elem. of f1 *)



(* f3[{C[1], n2}] = {C[1] -> k1, n2 -> k2} *)

f3[nList_] := Block[{imax,res},
	imax=Length[nList];
	res={};
	Do[
		res = Append[res,nList[[i]]->ToExpression[StringJoin[{"k",ToString[i]}]]],
	{i,imax}];
Return[res]
]


f4[xv_] :=
 	Block[{nc2k},nc2k=f3[xv[[3]]];Return[xv/.nc2k] (* f4 : MBsum[C[1], n2] -> MBsum[k1,k2] *)
]

f42[x_] := f4[f2[x]]


SimplifyMBsum[vMBsum_] := Map[f42, f1[vMBsum]]


SimplifyMBsumAll[s1_]:=Flatten[Map[SimplifyMBsum, s1]]


SimplifyMBsums[s1_]:=Block[{s2},

s2 = Plus @@ SimplifyMBsumAll[s1];

s2 = s2 //. {MBsum[x1__, y__, z__] + MBsum[x2__, y__, z__] :> 
     MBsum[x1 + x2, y, z], 
    a_ *MBsum[x2__, y__, z__] :> MBsum[a *x2, y, z]};
s2 = s2 /. {MBsum[x2__, y__, z__] :> MBsum[Factor[ExpandAll[x2]], y, z]};

s2 = PlusToList[s2];
Return[s2]
]


End[];
EndPackage[];
