
(*  File is available at 
    http://www.us.edu.pl/~gluza/ambre

   In case of any use, please, refer to
   J. Gluza, K. Kajda, T.Riemann, Computer Physics Communications 177(2007)879
   and
   J.Gluza, K. Kajda, T. Riemann, V. Yundin, arXiv:1010.1667
*)

(*------------------------------------------------------------*)
(* AMBREnLOOP - Automatic Mellin Barnes REpresentation N loop *)
(* by K. Kajda                                                *)
(*------------------------------------------------------------*)         

Print["AMBRE by K.Kajda   ver: 2.1.1 \nlast modified Aug 2017"]

(* List of last minute changes, 05.10.2010
1. BarnesLemma[repres_,num_] -> BarnesLemma[repres_,num_,powersOfPropagators_]
 powersOfPropagators, see BarnesLemma::usage below 
2. misprint:
OptimizeResult -> True changed to OptimizedResult -> False
*)

BeginPackage["AMBRE`"]


MBrepr::usage = "MBrepr[{numerator},{propagators},{internal momenta}] is a function for input integrals"
invariants::usage = "A list of invariants e.g. {p1^2->s}"
momentum::usage = "Momentum for a given subloop"
BarnesLemma::usage = "BarnesLemma[repr,1,powersOfPropagators] tries to apply 1st Barnes lemma on a given representation. 
\n BarnesLemma[repr,2,powersOfPropagators] is used to apply 2nd Barnes lemma.\n\n These two functions automatically apply lemmas on every M-B integral variable, only if a given variable appears in Gamma functions but not in functions related to kinematic variables. powersOfPropagators: e.g. {n1,n2,...}"
PR::usage = "PR[momenta,mass,power of propagator] is a propagator. Mass must be a symbol."
INT::usage = "Function that displays intermediate output: INT[numerator, MB representation, propagators1, propagators2].\n -numerator: it has the following form e.x. {k2[mu1],k2[mu2]}-->k2[mu1]k2[mu2] and is exactly the same as the one implemented in the old version of AMBRE i.e. ver: 1.2.\n -MB representation: Mellin-Barnes representation obtained during calculation of current subloop/iteration part. It is multiplied by Mellin-Barnes representation which was calculated in the previous step/iteration.\n -propagators1: propagators extracted out of F polynomial in the current iteration and merged by the powers with diagram propagators.\n -propagators2: keeps propagators in same form as propagators2 does. The role of propagators2 is as follows. Sometimes extracted propagators doesn't fit in the next iteration because of a different internal momentum, instead it/they will work fine during one of the following iterations."
eps::usage = "Epsilon variable"
X::usage = "Feynman parameter"
FX::usage = "FX[X[1]+X[2]]^2 = (X[1]+X[2])^2"
z::usage = "M-B integral variable, written in a form z1,z2,..."
g::usage = "Metric tensor"
d::usage = "d - dimension"
Contract::usage=""


Powers::error="Use non-numeric powers of propagators"

(***********************************************)
(*               Options                       *)
(***********************************************)
Options[MBrepr]={Text->True, OptimizedResult->False, BarnesLemma1->True, BarnesLemma2->False};
Options[BarnesLemma]={Text->False,Shifts->True};


Begin["`Private`"]



(***********************************************)
(*               Functions                     *)
(***********************************************)
SFred[x_]:=StyleForm[x, FontColor -> RGBColor[1, 0, 0]];
SFblue[x_]:=StyleForm[x, FontColor -> RGBColor[0, 0, 1]];
SFgray[x_]:=StyleForm[x, FontColor -> RGBColor[0.5, 0.5, 0.5]];
ToExprStr[x_,y_]:=ToExpression[ToString[x]<>ToString[y]];
UnsUnApp[x_,y_]:=Reap[Sow[1,Append[x,y]],_,#1&][[2]];
UnsortedUnion[x_]:=Reap[Sow[1, x], _, #1 &][[2]];
ToTake[x_]:=If[Length[x]>iteration,Take[x,iteration],x];
ToList[x_]:=If[Head[x]===Plus,List @@ x,{x}];

IFprint[str__,opt_]:=If[Text/.opt,Print[str]];

IFPrint[x_]:=If[False,Print[x]];
IFPrint[x_,y_]:=If[False,Print[x,y]];

OptChange[fun_,new___Rule]:=
  Module[{i},
     Map[If[Length[i=Position[{new},First[#]->_,1,1]]===0, 
	    #,{new}[[Sequence @@ First[i]]]]&, 
	 Options[fun]] ]

PropRules={PR[x_,y_,a_]PR[x_, y_, b_]->PR[x,y,a+b],
           PR[s1__,a_,n1_]PR[s2__,a_,n2_]:>PR[s1,a,n1+n2]/;Expand[(s1)^2-(s2)^2]==0};

(*set up global variables*)
d = 4 - 2 eps;
(***********************************************)
(*           M-B representation                *)
(***********************************************)
MBrepr[numerator_,propagators_,momenta_,options___Rule]:=
   Module[{OptimizedResult, option,loop=Length[momenta],
	   powersOfPropagators,powersCheck,numeratorTableStructure,propagatorsTableStructure,
	   momentum,subloopNumerator,subloopPropagator,subLoop,subloopRepresentation,
	   newNumerators={Num[1]},newPropagators={PRs[1]},newRepresentation={Repr[1]},
	   dotInvariants,result},

	(* Merges parts with common gamma functions *)
	OptimizedResult[result_]:=
	   Module[{noGammaResult, gammaResult, indexOld=1, indexNew, newResult={}},

		noGammaResult = (result /. Gamma[__] -> 1);
		gammaResult = Split[result/noGammaResult];
		
		Do[
 			indexNew = indexOld + Length[gammaResult[[i]]];
 			newResult = Append[newResult, 
					Take[noGammaResult, {indexOld, indexNew - 1}]*gammaResult[[i]]];
 			indexOld = indexNew, 
		{i, 1, Length[gammaResult]}];
		newResult = Map[Factor[Plus @@ #]&, newResult];

	  Return[newResult]];


	option=OptChange[MBrepr,options];

	(* Checks if_ powers of propagators are Symbols *)
	powersOfPropagators=Apply[List, First[propagators]] /.PR[_,_,n_]->n;
	powersCheck = 
 		Map[Head[#]&, powersOfPropagators];
	If[MemberQ[powersCheck, Integer | Real], {
		Message[Powers::error];
		result="error";
		Goto["Abort"];
	}];


	(* Prepares numerator structure & propagators to relevant iteration *)
	numeratorTableStructure=DoNumeratorTable[numerator, loop, momenta];
	propagatorsTableStructure=DoPropagatorsTable[List @@ First[propagators], loop, momenta];

	(* External momenta *)
	external=numeratorTableStructure[[2]];
	IFprint[SFred[">>External momenta = "],If[external=!=1,external,"N/A"],option];

	(* LoopByLoop/Re-insertion calculation *)
	IFprint[SFred[">>Starting LoopByLoop calculation"],option];
	Do[
		momentum=momenta[[iteration]];
		IFprint[SFblue["--iteration nr: "],iteration,SFblue[" with momentum: "],momentum,option];
		If[iteration==1,{znumb={{1,0}}; PropList=propagators}];       (* TBD later!!! *)
		
		(* Reads current iteration numerator/propagator *)
		subloopNumerator=numeratorTableStructure[[1,iteration]];
		subloopPropagator=propagatorsTableStructure[[iteration]];

		(* Adds numerator/representation/propagator from previous subloop result *)
		subloopNumerator=newNumerators /.Num[x_]:>Num[{x*Apply[Times, subloopNumerator] /.Times->Sequence}];
		subloopRepresentation=newRepresentation;
		subloopPropagator=newPropagators /.PRs[x_]:>PRs[x*subloopPropagator //.PropRules];
		subloopPropagator=subloopPropagator /.PRs[x_]:>ShiftPropagators[x, momentum];

		(* Calculates current subloop part *)
		subLoop=Transpose[{subloopNumerator,subloopRepresentation,subloopPropagator}];
		IFprint[SFgray["  Run "],"?INT",SFgray[" to see description of below output \n"],subLoop /.{Num[x_],Repr[y_],PRs[z1_,z2_]}:>INT[x,y,z1,If[z2=!=1,z2,"N/A"]],option];
		subLoop=subLoop	/.{Num[x_],Repr[y_],PRs[z1_,z2_]}->Int[x,y,z1,z2];
		subLoop=subLoop /.Int[subloopNumerator_,subloopRepresentation_,subloopPropagator_,laterPropagator_]:>
				 SubLoop[iteration,subloopNumerator,subloopRepresentation,subloopPropagator,laterPropagator,momenta];

		IFprint[SFgray["  F polynomial during this iteration \n"],fupc,option];

		(* New numerators/representation/propagators parts *)
		newNumerators=Map[#[[1]]&, subLoop]//Flatten;
		newRepresentation=Map[#[[2]]&, subLoop]//Flatten;
		newPropagators=Map[#[[3]]&, subLoop]//Flatten,
	{iteration,1,loop}];



	(* Contracting and finalizing output *)
	IFprint[SFred[">>Contracting and finalizing output"],option];
	result=newRepresentation /.Repr[x_]->x;

	IFprint[SFblue["--contracting..."],option];
	If[external=!=1,{
		dotInvariants=invariants /.(p_^2->x__)->(p.p->x) /.(p1_*p2_->x__)->(p1.p2->x);
		IFPrint["dotInvariants ",dotInvariants];
		
		result=Contract[newRepresentation /.Repr[x_]->x, external] /.dotInvariants;
	}];
	result=DeleteCases[result,0];

	(* Optimalization of M-B representation *)
	IFprint[SFblue["--finalizing output..."],option];
	If[OptimizedResult /.option, {
		If[Length[result]>1, result=Apply[List, Plus @@ result]]; (* Sums same parts of M-B representation *)
		result = OptimizedResult[result] }];

	(* Automatic application of Barnes Lemmas *)
	If[BarnesLemma1 /.option, {
		IFprint[SFred[">>Checking Barnes 1-st lemma..."],option];
		result=BarnesLemma[#, 1, powersOfPropagators]& /@ result}];
	If[BarnesLemma2 /.option, {
		IFprint[SFred[">>Checking Barnes 2-nd lemma..."],option];
		result=BarnesLemma[#, 2, powersOfPropagators]& /@ result}];


	Label["Abort"];

 Return[ Map[Factor[#]&, result] ]]



(***********************************************)
(*            Do numerator table               *)
(***********************************************)
DoNumeratorTable[numerator_, loop_, momenta_]:=
  Module[{ModifiedPartition, numeratorTable, external, numerator2, indices, internalMomentaMatrix, seqRul, gmnRul},


ModifiedPartition[x_]:=If[EvenQ[Length[x]], Partition[x,2], Append[Partition[x,2], Last[x]] ];

(* Is it scalar or tensor numerator? *)
If[numerator==={1}, {								
	numeratorTable=Table[{1}, {i,1,loop}];
	indices=0;
	external=1;
}, {
	numerator2=numerator /.Times->Dot /.x_^_->x.x;

	(* Is it numerator with external momenta or not? *)
	If[ Union[Map[Head[#]&,numerator2]]=!={Dot} , {				
		indices=Map[# /.k_[mu_]->mu &, numerator2];    (* Without Map it fails for_one dim list *)
 		external=1}, {


		(* Will change internal momenta scalars (if_any) into sequence inside a list *)
		internalMomentaMatrix=Flatten[Table[momenta[[i]].momenta[[j]], {i,loop},{j,loop}]];

       		seqRul=internalMomentaMatrix /.ka_ .kb_->(ka.kb->Sequence[ka,kb]);	
 		gmnRul=internalMomentaMatrix /.ka_ .kb_->({ka[mu_],kb[nu_]}->g[mu*nu]);
 		numerator2=numerator2 /.seqRul;
		(**)

 		indices=Table[ToExpression["mu"<>ToString[i]], {i,1,Length[numerator2]}];

		(* Generates external projector *)
	 	external=Flatten[numerator2 /.Dot[x_,y_]:>DeleteCases[List[x,y],Alternatives @@ momenta]];
 		external=Map[# /.{k_,mu_}->k[mu]&, 
				Apply[List,Transpose[{external, indices}]] ];
		kExternal=Cases[external, Apply[Alternatives, Map[# /.x_->x[_]&, momenta]]];
		pExternal=Times @@ Complement[external, kExternal];
		kExternal=Times @@ Flatten[Map[#/.gmnRul&, ModifiedPartition[kExternal]]];
		external=pExternal*kExternal;

 		numerator2=Cases[numerator2 /. Dot->Sequence, Alternatives @@ momenta];
 		numerator2=Apply[NonCommutativeMultiply,Transpose[{numerator2, indices}]] 
					/.{k_,mu_}->k[mu] /.NonCommutativeMultiply->List
		}
	];
	numeratorTable=Map[Cases[numerator2,#[_]]&, momenta] /.{}->{1};
}];


 Return[{numeratorTable, external}]]



(***********************************************)
(*           Do propagators table              *)
(***********************************************)
DoPropagatorsTable[propagators_, loop_, momenta_]:=
  Module[{kiCasePropagators, deleteWithMomenta={}, propagatorsTable={}},

	Do[
		kiCasePropagators=Cases[propagators, PR[a_.*Part[momenta,i]+x_.,__]];
		kiCasePropagators=DeleteCases[kiCasePropagators, Alternatives @@ deleteWithMomenta];
		deleteWithMomenta=Append[deleteWithMomenta, kiCasePropagators] //Flatten;
		propagatorsTable=Append[propagatorsTable, kiCasePropagators],
	{i,1,loop}];

	propagatorsTable = Map[(Times @@ #)&, propagatorsTable];

 Return[propagatorsTable]]


(***********************************************)
(*           Propagators used now/later        *)
(***********************************************)
ShiftPropagators[propagators_, momentum_]:=
  Module[{kiCasePropagators, nokiCasePropagators},

	kiCasePropagators=Times @@ Cases[propagators, PR[a_.*momentum+x_.,__]];
	nokiCasePropagators=propagators/kiCasePropagators;

 Return[PRs[kiCasePropagators, nokiCasePropagators]]]


(***********************************************)
(*        Main for given subloop               *)
(***********************************************)
SubLoop[iteration_,subloopNumerator_,subloopRepresentation_,subloopPropagator_,laterPropagator_,momenta_]:= 
  Module[{currentMomentum, ipow, Nup, tensorial, newNumerators, newPropagators, 
	  afterXintegration, newRepresentation},	 


	currentMomentum=momenta[[iteration]];
	fupc = Expand[FUpolyfun[subloopPropagator,currentMomentum,momenta,invariants,option]//Simplify];
	ipow = List @@ subloopPropagator /. PR[__,nu_]->nu;
	Nup = Expand[Total[ipow] - d/2(* - 2 + eps *)]; (*dim*)

	tensorial = SubloopNumerator[iteration,subloopNumerator,currentMomentum,momenta,subloopPropagator,Nup];

	afterXintegration=Xintegrate[tensorial[[3]], subloopPropagator, ipow];
	newRepresentation=(-1)^(Expand[Nup + d/2(* + 2 - eps *)])/Apply[Times,Gamma[#]& /@ ipow]*afterXintegration*Part[tensorial,2]; (*dim*)

	newNumerators=Num[#]& /@ Flatten[Part[tensorial,1]];
	newRepresentation=Repr[#]& /@ Flatten[newRepresentation*subloopRepresentation];
	newPropagators=PRs[#]& /@ Flatten[Part[tensorial,4]*laterPropagator];
	
 Return[{newNumerators,newRepresentation,newPropagators}]]



(***********************************************)
(*            Mellin-Barnes formula            *)
(***********************************************)
MBarnes[iteration_,fupc_,integral_,Nup_]:= 
  Module[{imin,fuplist,poly=fupc,nfu,
	  MB1,pow=Nup,MB2,barnes,
	  FXcheck,aux,temp={}},
	   
	
    znumb = Take[znumb,iteration];
    imin = znumb[[iteration,2]];
    
    Label[FXstart];
    fuplist = ToList[Expand[poly]];
 
    nfu = Length[fuplist];
    MB1 = Sum[ToExprStr[z,i+imin],{i,1,nfu-1}]+pow;
    MB2 = Product[fuplist[[i]]^ToExprStr[z,i+imin]*
		  Gamma[-ToExprStr[z,i+imin]],{i,1,nfu-1}];
    barnes = MB2*fuplist[[nfu]]^(-MB1)*Gamma[MB1]/Gamma[pow];
 
    znumb = UnsUnApp[znumb,{iteration+1,imin+nfu-1}];
    FXcheck = Cases[PowerExpand[barnes*aux],FX[_]^_];
    
    If[Length[FXcheck] != 0,
       {
	 {poly,pow} = FXcheck /.{FX[x_]^a_}->{x,-a};
	 temp = Append[temp,barnes/. FX[__]->1];
	 imin = znumb[[iteration+1,2]];
	 znumb = Delete[znumb,-1];
         Goto[FXstart]
       }];

    If[Length[temp]!=0,barnes=First[temp]*barnes];


 Return[barnes]]



(***********************************************)
(*            Integration over  x              *)
(***********************************************)
Xintegrate[fnu_,integral_,ipow_] := 
  Module[{tomult,expr,xonly,numG,denG},
     
    tomult = Product[X[i]^ipow[[i]], {i,1,Length[ipow]}];
    expr = PowerExpand[tomult*fnu];

    xonly = expr/(expr /.X[_]->1);	     

    numG = xonly /. Power[_,x__]:>Gamma[Expand[x]];
    denG = numG //. Gamma[a_]*Gamma[b_]:>Gamma[Expand[a+b]];
           
 Return[numG/denG]]



(***********************************************)
(*         Vector & Tensor parts               *)
(***********************************************)
SubloopNumerator[iteration_,numerator_,currentMomentum_,momenta_,subloopPropagator_,Nup_] :=
    Module[{rank=Length[numerator],FF,PQ,AP,P,Qin,APobject,indices,newNumerators,
	    gams,Allgams,XandPR,Xonly,newPropagators,NoIntMom,NoXnIntMom,NoXPRgams,result},
    

    FF[i_]:=MBarnes[iteration,fupc,subloopPropagator,Nup-i];
    PQ[i_]:=Qin /.a_.*p_*X[n_]:>a*p[i]X[n];

    AP[ind_,rank_]:= 
      Module[{poss,rul1,rul2,rul3,rul4,A,result},
    
      	  poss = Select[Permutations[ind],Signature[#]==1&];
       	  rul1 = A[x_]*P[y_]:>Map[A[Take[#,x]]*P[Take[#,-y]]&,poss];
       	  rul2 = A[x_]:>(Times @@ Map[g[Times @@ #]& ,Partition[x,2]]);
       	  rul3 = P[x_]:>(Times @@ Map[P,x]);

       	  result =
	      Map[UnsortedUnion[#]&,
	          Table[A[r]P[rank-r],{r,0,rank,2}] /.rul1 /.rul2 /.rul3];

    Return[result]];


    If[numerator==={1}, 
    {newNumerators={1};
     gams={{Gamma[Nup]*FF[0]}};
     Allgams=gams /.X[_]->1 /.PR[_,_]->1;
     XandPR=PowerExpand[gams]/PowerExpand[Allgams] /.PR[k_,m_]^n_->PR[k,m,-n];
     Xonly=XandPR /.PR[_,_,_]->1;
     newPropagators=XandPR /.X[_]->1;	

    Goto[NumkEND]}];

    indices=Map[(# /.k_[mu_]->mu)&,Flatten[numerator]];
    Qin=ToPoly[subloopPropagator, currentMomentum][[3]];
    APobject=Expand[AP[indices,rank] /.P->PQ]/.Plus->List;
    NoIntMom=APobject /.Map[#[_]->1&,momenta];
    newNumerators=APobject/NoIntMom;
    NoXnIntMom=NoIntMom /.X[_]->1;

    gams=Table[Gamma[Nup-r/2]/(-2)^(r/2)*FF[r/2],{r,0,rank,2}];	
    NoXPRgams=gams /.PR[_,_]->1 /.X[_]->1;
    XandPR=PowerExpand[gams]/PowerExpand[NoXPRgams] /.PR[k_,m_]^n_->PR[k,m,-n];	
    Allgams=NoXnIntMom*NoXPRgams;
    Xonly=(NoIntMom/NoXnIntMom)*(XandPR /.PR[_,_,_]->1);
    newPropagators=Xonly*XandPR /.X[_]->1;

     Label[NumkEND];


   Return[{newNumerators,Allgams,Xonly,newPropagators,rank}]]  



(***********************************************)
(*              Contract                       *)
(***********************************************)
Contract[repr_,external_]:=
      Module[{rulA,rulB,rulC,rulD,rulE,result},
    	
    rulA = g[mu1_*mu3_]*g[mu2_*mu3_]:>g[mu1*mu2];
    rulB = g[__]^2->d;
    rulC = g[mu1_*mu2_]*p_[mu2_] -> p[mu1];
    rulD = p1_[x_]*p2_[x_] :> p1.p2/;p1=!=g/;p1=!=X/;p1=!=Gamma;
    rulE = p1_[mu1_]^2 :> p1.p1/;p1=!=X /;p1=!=Gamma;

    result = (repr*external) 
	 //.rulA //.rulB //.rulC //.rulD //.rulE;

  Return[result]]


(***********************************************)
(*           Qin, jpol, tt                     *)
(***********************************************)
ToPoly[integral_, momentum_] := 
  Module[{jpol, ttpol, Qin, denom, integrand},

    (*Print["#ToPoly: {integral,momentum} = ", {integral,momentum}];*)
    (*denom = List @@ integral;*) 
    denom = List[List @@ (integral)] // Flatten; (*updated 16.04.2014*)
    integrand =
      Sum[X[i]*denom[[i]] /.PR[k_,m_,_]->k^2-m^2,{i,Length[denom]}];
    ttpol = Coefficient[integrand, momentum*momentum];
    jpol = integrand - momentum*ttpol*momentum;
    Qin = Expand[-1/2*Coefficient[jpol, momentum]];
    jpol = Expand[jpol + 2*momentum*Qin];

 Return[{ttpol,jpol,Qin}]]


(***********************************************)
(*            Finds F function                 *)
(***********************************************)
FUpolyfun[integral_,momentum_,momenta_,invariants_,option_]:= 
  Module[{jpol, ttpol, Qin,(*denom, integrand,*)fpoly,fpolyToList},

    {ttpol,jpol,Qin}=ToPoly[integral,momentum]; (*Print["#FUpolyfun: {integral,momentum} = ", {integral,momentum}];*)
    Upoly = ttpol;
    (*fpoly = Expand[Simplify[ttpol*(Expand[Qin*1/ttpol*Qin] - jpol)]];*)
    fpoly = Expand[Simplify[Expand[Qin*Qin - ttpol*jpol]]]; (* Print["#FUpolyfun: {ttpol,jpol,Qin}=",{ttpol,jpol,Qin}]; *)

    If[momentum == Last[momenta],
	fpoly = Expand[fpoly /. invariants],
	{},
	fpoly = Expand[PRfind[fpoly,integral,momenta,invariants] /.invariants]
	];

    masslist = Union[List @@ integral /. PR[_, m_, _] -> m];
    masslist = DeleteCases[masslist, 0];
    Do[
  	caslist = Cases[fpoly, masslist[[i]]^2  X[n_]^2];
  	If[Length[caslist] >= 2,
	  {casexpr = caslist /. a_^2 -> a;
      	   expr = Simplify[Plus @@ casexpr];
      	   FXexpr = expr /. x_.*(y__) -> x^2*FX[y]^2;
      	   FXrule = {First[caslist]->FXexpr-(expr)^2+First[caslist]};
      	   fpoly = Expand[fpoly /. FXrule]; 
      	  },{}],
    {i, 1, Length[masslist]}];
    
    (* Print["Fpoly -> ", fpoly];  *)

    (***********improvements***************)
    (*negative masses*)
    cleanmassPR = - m_^2 X[i1_] X[i2_] - PR[p_, m_] X[i1_] X[i2_] -> - PR[p, 0] X[i1] X[i2];
    fpoly = fpoly //. cleanmassPR; 
    (*U polynomial*)
    If[(* Length[fpoly/.{X[_]->1, -1->1}] > 1 *) Head[fpoly]==Plus,{
    (*IF*)

    fxlst = Cases[(List@@fpoly), __ FX[__]^2]; (* Print["fxlst = ", fxlst]; *)
    invlist = DeleteDuplicates[Complement[(List@@fpoly), fxlst] /. {X[_]->1}]; (* Print["invlist = ", invlist]; *)
    xlist = Coefficient[Plus@@Complement[(List@@fpoly), fxlst], invlist]; (* Print["xlist = ", xlist]; *)
    
    templist = xlist;
    Do[ If[ Head[templist[[i]]]==Plus,
            templist[[i]] = templist[[i]] /. {Plus->List},{}, 
            templist[[i]] = List[ templist[[i]] ] ],
       {i,1,Length[templist]}]; (* Print["templist1 = ", templist]; *)
    
    Do[ If[ Length[ templist[[i]] ] >= Length[ List @@ Upoly ], 
            If[ Length[ templist[[i]] ] == Length[ List @@ Upoly ],
                {isec = Intersection @@ (templist[[i]] /. {Times->List, X[idx_]^2->List[X[idx], X[idx]]}),
                 If[isec != {}, templist[[i]] = isec]},
                {}
               ] 
           ],
      {i,1,Length[ templist ]}];  (* Print["templist2 = ", templist]; *)

    Do[templist[[i]] = templist[[i]] /. List -> Plus, {i,1,Length[templist]}];
    
    fpoly = Expand[Plus @@ (invlist*templist) + Plus @@ fxlst];
    (*END IF*)
    }];
    (*U polynomial with FX[]*)
    If[(* Length[fpoly/.{X[_]->1, -1->1}] > 1 *) Head[fpoly]==Plus,
    (*IF*){
        fxlst = Cases[(List@@fpoly), __ FX[__]^2], (* Print["fxlst = ", fxlst], *)
        If[Length[fxlst] >= 1, (*!!! >= 1 and Do loop for fxlst[[i]] !!!*)
        (*IF*)Do[{
            fxvar = fxlst[[i]] /. FX[_] -> 1, (* Print["fxvar = ", fxvar], *)
            fxarg = fxlst[[i]] /. __ FX[arg_]^2 -> arg, (* Print["fxarg = ", fxarg], *)
            tofactor = Coefficient[fpoly,fxvar] /. FX[_] -> 0, (* Print["tofactor = ", tofactor], *)
            If[Head[tofactor] == Plus,
            (*IF*){
                factorexpr = Factor[tofactor],
                If[Head[factorexpr] == Times,
                (*IF*){
                    tlst = List@@(factorexpr), FXcheck = False,
                    If[Count[tlst, fxarg] == 1 && Plus@@tlst == Upoly, FXcheck = True]
                }(*ENDIF*)]
            }(*ENDIF*)],
            If[FXcheck, fpoly = (fpoly /. fxvar -> 0) + Expand[fxvar*fxarg]] 
        },{i,1,Length[fxlst]}](*ENDIF*)] 
    }(*ENDIF*)];
(*     If[Length[fpoly/.{X[_]->1, -1->1}] > 1,{
    (*IF*)
    fxlst = Cases[(List@@fpoly), __ FX[__]^2]; Print["fxlst = ", fxlst];
    (* invlist = DeleteDuplicates[Complement[(List@@fpoly), fxlst] /. {X[_]->1}]; Print["invlist = ", invlist]; *)
    invlist = DeleteDuplicates[(List@@fpoly) /. {X[_]->1, FX[_]->1}]; Print["invlist = ", invlist];
    
    (* xlist = Coefficient[Plus@@Complement[(List@@fpoly), fxlst], invlist]; Print["xlist = ", xlist];
    fxlist = Cases[#,FX[__]^2]&/@(List@@##&/@Coefficient[fpoly, invlist]); Print["fxlist = ", fxlist]; *)
    xlist = Coefficient[Plus@@Complement[(List@@fpoly), fxlst], invlist]; Print["xlist = ", xlist];
    fxlist = Coefficient[fpoly, invlist]; Print["fxlist = ", fxlist];

    templist = xlist;
    Do[ If[ Head[templist[[i]]]==Plus,
            templist[[i]] = templist[[i]] /. {Plus->List},{}, 
            templist[[i]] = List[ templist[[i]] ] ],
       {i,1,Length[templist]}]; Print["templist1 = ", templist];
                                      
    Do[ If[ fxlist[[i]] != {}, 
            If[ Length[ templist[[i]] ] == Length[ (fxlist[[i]] /. FX[arg_]^2 -> arg)[[1]] ],
                {tvr = templist[[i]] /. {Times->List, X[idx_]^2->List[X[idx], X[idx]]}, (* Print["tvr = ", tvr], *)
                 isec = Intersection @@ tvr, (* Print["isec = ", isec], *)
                 If[isec != {}, {tocheck = (Plus@@(Complement[#,isec]&/@tvr))[[1]], (* Print["tocheck =", tocheck], *)
                    If[tocheck == (fxlist[[i]] /. FX[arg_]^2 -> arg)[[1]] && tocheck + isec[[1]] == Upoly, 
                       {templist[[i]] = tocheck, fxlist[[i]] = {}}]
                  }]},
                {}
               ] 
           ],
      {i,1,Length[ templist ]}];  

    Do[templist[[i]] = templist[[i]] /. List -> Plus, {i,1,Length[templist]}]; (* Print["templist2 = ", templist]; *)
    Do[fxlist[[i]] = fxlist[[i]] /. {{} -> 0, {var_} -> var}, {i,1,Length[templist]}]; (* Print["fxlist = ", fxlist]; *)

    fpoly = Expand[Plus @@ (invlist*templist) + Plus @@ (invlist*fxlist)];
    (*END IF*)
    }]; *)
    (*linear FX*)
    (*????????*)
    (**************************************)

 Return[fpoly]]


(*************************************************)
(*          Creates PRs in F polynomial          *)
(*************************************************)
PRfind[fpoly_,integral_,momenta_,invariants_]:=
  Module[{intmom,masrul,nomass,preprop,p1,p2,ints,intrul,F,
	  PRnom,nom,PRmass,rul1,rul2,nm,PRadd,aux},
   
     
    intmom = Alternatives @@ (a_. momenta);
    masrul = DeleteCases[Union[List @@ integral /.PR[_,a_,_]->a],0];
    nomass =
      Simplify[#]& /@ Collect[fpoly/.(#->0& /@ masrul),X[a_]X[b_]];
    preprop = Plus @@ Cases[nomass+aux,(intmom)^2*_|((intmom)+__)^2*_];
	   
    p1 = preprop /.(x_)^2->PR[x];
    p2 = Expand[fpoly-preprop];

    ints= PropList[[1]];
    intrul = ints /.PR[a_,b_,_]->(PR[a]->PR[a,b]+b^2);
    intrul = DeleteCases[List @@ (intrul*(aux->0)),aux->0];

    F = Expand[p1+p2 /.intrul /.invariants];       
    
    PRnom = Cases[F,x_. PR[a_]];
    If[Length[PRnom] == 0,{},
    {nom = PRnom /. w_. PR[a_]X[n_]X[m_]->X[n]X[m];
     PRMass = Cases[F,Alternatives @@ (a_. m_^2*nom)];
     rul1 = PR[x_]-m_^2->PR[x,m];
     rul2 =-PR[x_]+m_^2->-PR[x,m];
     nm = Plus @@ Join[PRnom, PRMass];
     PRadd = Collect[nm,nom] /.{rul1,rul2};
     F=F-nm+PRadd}];

   
 Return[Expand[F /. PR[x_]->PR[x,0]]]]


(*************************************************)
(*          Small Barnes 1 Lemma                 *)
(*************************************************)
BarnesLemma[repres_, num_, powersOfPropagators_, options___Rule] :=
    Module[ {option, 
      ZETS, RSC, Zcheck, 
      BL1, BL2, BL3, BL4, 
      B2a,
      Bmes1, Bmes2, Bmes3, 
      Zpowers, ShiftRul, powlist, gamlist, 
      varlist, dimrepr,reprList, BLrulToApply, 
      before, after,
        counter = 0, BList = {}, zel, denR, ifden, numR, ifnum, result},
        
        option = OptChange[BarnesLemma,options]; 

        (* powersOfPropagators = Apply[List, First[integral]] /. PR[_,_,n_]->n; *) (* Print[powersOfPropagators]; *)
        
        ZETS[zlst_] :=
            Module[ {FML, rslt},
                FML[x_] :=
                    Flatten[Map[(List @@ #)&, x]];
                rslt = DeleteCases[FML[FML[zlst]], _Integer] // Union;
                Return[rslt]
            ];
        
        RSC[x_,y_] :=
            ToExprStr[z, #]& /@ 
            Reverse[Sort[
            ToExpression[StringDrop[ToString[#], 1]]& /@ 
            Complement[x, y, {}]]];
        
        Zcheck[x_,z_] :=
            Cases[x,Gamma[a_.*z+__]|Gamma[a_.*z]];
       

     (* Using 1st or 2nd Barnes Lemma *)
        Switch[num,
           1,
           
        BL1[z_] :=
            Gamma[z + (l1_.)]*Gamma[z + (l2_.)]*Gamma[ - z + (l3_.)]*Gamma[ - z + (l4_.)]
            -> Gamma[l1 + l3]*Gamma[l1 + l4]*Gamma[l2 + l3]*
            (Gamma[l2 + l4]/Gamma[l1 + l2 + l3 + l4]);
        BL2[z_] :=
            Gamma[z + (l1_.)]^2*Gamma[ - z + (l3_.)]*Gamma[ - z + (l4_.)] 
            -> Gamma[l1 + l3]^2*Gamma[l1 + l4]^2/Gamma[2*l1 + l3 + l4];
        BL3[z_] :=
            Gamma[z + (l1_.)]*Gamma[z + (l2_.)]*Gamma[ - z + (l3_.)]^2 
            -> Gamma[l1 + l3]^2*Gamma[l2 + l3]^2/Gamma[l1 + l2 + 2*l3];
        BL4[z_] :=
            Gamma[z + (l1_.)]^2*Gamma[ - z + (l3_.)]^2
            -> Gamma[l1 + l3]^4/Gamma[2*l1 + 2*l3];
        
        Bmes1 = ">> Barnes 1st Lemma will be checked for: ";
        Bmes2 = ">> Representation after 1st Barnes Lemma:  <<";
        Bmes3 = "   1st Barnes Lemma was applied for: ";
        
        (* BLrul = Dispatch[{BL1[zel], BL2[zel], BL3[zel], BL4[zel]}], *)
        BLrul[z_] := {BL1[z],BL2[z],BL3[z],BL4[z]},
        (*----------------------------------------*)          
           2,
        
        B2a[z_] :=
            Gamma[(l1_.) + z]*Gamma[(l2_.) + z]*Gamma[(l3_.) + z]*
            Gamma[(l4_.) - z]*Gamma[(l5_.) - z]/Gamma[(l6_.) + z]
            :> Gamma[l1 + l4]*Gamma[l2 + l4]*Gamma[l3 + l4]*Gamma[l1 + l5]*
            Gamma[l2 + l5]*Gamma[l3 + l5]/(Gamma[l1 + l2 + l4 + l5]*
            Gamma[l1 + l3 + l4 + l5]*Gamma[l2 + l3 + l4 + l5]) 
            /;Expand[l1 + l2 + l3 + l4 + l5 - l6] === 0;
            
        Bmes1 = ">> Barnes 2nd Lemma will be checked for: ";
        Bmes2 = ">> Representation after 2nd Barnes Lemma:  <<";
        Bmes3 = "   2nd Barnes Lemma was applied for: ";
        
        (* BLrul = Dispatch[{B2a[zel]}] *)
        BLrul[z_] := {B2a[z]}
        ];

     (* Preparing representation *)
        reprList = ToList[Expand[repres /. Gamma[x_] :> Gamma[Expand[x]]]];
        (* Print["#B: reprList -> ", reprList]; *)    

     (* Creating list of zs NOT to be used *)
        Zpowers = LSTofZ[reprList,1,powersOfPropagators];
        (* Print["#B: Zpowers -> ", Zpowers]; *)      

     (*---------------------------*)
     (*      Shifting part...     *) (* Print["!!!!",(Shifts /. option)];Print["!!!!", option]; Print["!!!!",(Shifts)]; *)
        If[ (Shifts /. option),
            ShiftRul = SHIFT[Zpowers];
            If[ Length[ShiftRul] != 0,
                IFprint[">>", SFred[" Shifting: "], ShiftRul, option];
                reprList = (# /. ShiftRul /. Gamma[x_] :>
                      Gamma[Expand[x]])& /@ reprList;
                Zpowers = LSTofZ[reprList, 1, powersOfPropagators],
                {}
            ],
            {}
        ];

     (*---------------------------*)
     
        powlist = ZETS[Zpowers];


     (* Creating list of zs which appear in Gammas *)
        gamlist = ZETS[LSTofZ[reprList,2,powersOfPropagators]];
     
            
     (* Creating list of zs to check with B-L*)
        varlist = RSC[gamlist,powlist];
        dimrepr = Length[gamlist];
        (* Print["#B: Zpowers -> ", Zpowers];
        Print["#B: powlist -> ", powlist];
        Print["#B: gamlist -> ", gamlist];
        Print["#B: varlist -> ", varlist]; *)

     (* Beggining of checking *)
        IFprint[Bmes1, varlist, " <<\n   Starting with dim=", 
              dimrepr, " representation...", option]; 
        result = reprList;
        Do[
        zel = varlist[[i]];
        If[ Text /. option,
            WriteString[$Output, "\n", i, ". ", "Checking ", zel]
        ];
        denR = Denominator[First[result]];
        ifden = Length[Zcheck[denR,zel]];
        numR = Numerator[First[result]];
        ifnum = Length[Zcheck[numR,zel]];
        If[ num == 1 && ifden != 0,
            Goto["break"]
        ];
        If[ num == 1 && ifnum > 4,
            Goto["break"]
        ];
        If[ num == 2 && ifden == 0,
            Goto["break"]
        ];
        If[ num == 2 && ifnum > 5,
            Goto["break"]
        ];

        (* before = result;
        result = (# /. BLrul)& /@ result; *)
        before=result; BLrulToApply = Dispatch[BLrul[zel]];
        result=(# /. BLrulToApply)& /@ result; (* Print["BLrulToApply = ", BLrulToApply]; *)
        after = result;

        If[ before =!= after,
            {
              counter = ++counter;
              BList = Prepend[BList, zel];
              If[ Text /. option,
                  WriteString[$Output, "...Barnes Lemma was applied."]
              ]
            },
            {}
        ];
        Label["break"],
        {i,1,Length[varlist]}];
        
        IFprint[Bmes2,option];
        
        If[ Length[BList] == 0,
            IFprint[SFred["   Could not apply Barnes-Lemma"], option],
            {
            IFprint[SFred[Bmes3], BList,
              "\n   Obtained representation has: dim=",
              dimrepr - counter, option]
            }
        ];
        Return[Plus @@ result (* ShiftRul must be also here *)]
    ]

LSTofZ[xli_,contr_,lstpowers_] :=
    Module[ {pow,result},
      
        powers = Map[Abs[#]&, lstpowers] /. Abs[x_] -> x;
        powers = # -> 0& /@ DeleteCases[(List @@ powers), _Integer];
        
        Switch[contr,
           1,
          pow = First[xli] /. Gamma[_] -> 1 /. d -> 0 /. eps -> 0 /. powers,
               2,
          rul = {1/Gamma[x__]^a_. -> x, Gamma[x_]^a_. -> x};
          pow = Cases[First[xli], Gamma[__]^n_.] /. rul /. d -> 0 /. eps -> 0 /.powers (*dim*)
        ];
        result = DeleteCases[
            If[ Head[pow] === Times || Head[pow] === List,
                List @@ pow,
                {pow}
            ] /. _^x_ :> 
                 Expand[x], _Integer | _Real | _Rational, 2];
        If[ result == 1,
            result = {},
            {},
            {}
        ];
        Return[result]
    ]

SHIFT[zlist_] :=
    Module[ {SD, prep, crosTab, comel, nocomel, ord, result},
        
        SD[x_] :=
            ToExpression[StringDrop[ToString[x], 1]];
        
        prep = DeleteCases[zlist, _Symbol];
        crosTab = DeleteCases[INT[prep, -1] // Flatten, 0];
        
        If[ Length[crosTab] == 1,
            GOTO["direct1"],
            {}
        ];
        (* Print["#S: crosTab - > ", crosTab]; *)
        
        comel = DeleteCases[INT[crosTab, 1] // Flatten, 0];
        If[ Length[comel] == 0,
            GOTO["direct1"],
            {}
        ];
        (* Print["#S: comel - > ", comel]; *)
        
        nocomel = DeleteCases[
                     DeleteCases[#, Alternatives @@ comel] & /@ crosTab, _Symbol];
        (* Print["#S: nocomel - > ", nocomel]; *)
        ord = Sort[#, SD[#1] < SD[#2] &] & /@ nocomel;
        (* Print["#S: ord - > ", ord]; *)
        
        Label["direct1"];
        result = (First[#] -> First[#] - (Take[#, 1 - Length[#]])) & /@ ord;
        result = result /. {Rule[x_, y_] :> Rule[-x, -y] /; Head[x] == Times};
        (* Print["#S result -> ", result]; *)
        
        Return[result]
    ];

(* INT for SHIFT *)
    INT[x_, a_] :=
        Module[ {ToL, ToS, temp, res = {}},
        
            ToL[s_] :=
                {List @@ (s)} // Flatten;
            ToS[s_] :=
                Plus @@ s;
        
            Do[ If[ i == j,
                    res = Append[res, {0}],
                    {temp = Intersection[ToL[Part[x, i]], (a)*ToL[Part[x, j]] ],
                     If[ Length[temp] == 1,
                         res = Append[res, {0}],
                         res = Append[res, ToS[temp]]
                     ] }
                ],  
            {i, 1, Length[x]}, {j, 1, i}];
            (* Print["#INT res ->", res]; *)
        
            Return[res]
        ]


End[]


EndPackage[]
