(* 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.0 \nlast modified 18 Jun 2010"] (* 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->False}; 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}; (***********************************************) (* 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 = Total[ipow]-2+eps; tensorial = SubloopNumerator[iteration,subloopNumerator,currentMomentum,momenta,subloopPropagator,Nup]; afterXintegration=Xintegrate[tensorial[[3]], subloopPropagator, ipow]; newRepresentation=(-1)^(Nup+2-eps)/Apply[Times,Gamma[#]& /@ ipow]*afterXintegration*Part[tensorial,2]; 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}, denom = List @@ integral; 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}, {ttpol,jpol,Qin}=ToPoly[integral,momentum]; Upoly = ttpol; fpoly = Expand[Simplify[ttpol*(Expand[Qin*1/ttpol*Qin] - jpol)]]; 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]}]; 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, BLru, pow, Zpowers, ShiftRul, powlist, gamlist, varlist, dimrepr,reprList, before, after, counter=0, BList={},zel, denR, ifden, numR, ifnum, result}, option=OptChange[BarnesLemma,options]; ZETS[zlst_]:= Module[{FML,result}, FML[x_]:=Flatten[Map[(List @@ #)&, x]]; result=DeleteCases[FML[FML[zlst]],_Integer]//Union; Return[result]]; 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]}], (*----------------------------------------*) 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]}] ]; (* Preparing representation *) reprList=ToList[Expand[repres/.Gamma[x_]:>Gamma[Expand[x]]]]; (* Creating list of zs NOT to be used *) Zpowers=LSTofZ[reprList,1,powersOfPropagators]; (*---------------------------*) (* Shifting part... *) 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]; (* 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; 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]] 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 /.eps->0/.powers, 2, rul={1/Gamma[x__]^a_.->x,Gamma[x_]^a_.->x}; pow=Cases[First[xli],Gamma[__]^n_.]/.rul /.eps->0 /.powers ]; 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]]; INT[x_, a_] := Table[If[i == j, {0}, Intersection[Part[x, i], Expand[(a)*Part[x, j]]]], {i, 1, Length[x]}, {j, 1, i}]; prep = DeleteCases[zlist, _Symbol]; crosTab = DeleteCases[INT[prep, -1] // Flatten, 0 | -1, 3]; If[Length[Tab] == 1, GOTO["direct1"], {}]; comel = DeleteCases[INT[crosTab, 1] // Flatten, 0]; If[Length[comel] == 0, GOTO["direct1"], {}]; nocomel = DeleteCases[ DeleteCases[#, Alternatives @@ comel] & /@ crosTab, _Symbol]; ord = Sort[#, SD[#1] < SD[#2] &] & /@ nocomel; Label["direct1"]; result = (First[#] -> First[#] - (Take[#, 1 - Length[#]])) & /@ ord; Return[result]] End[] EndPackage[]