FeynCalc · vsht · Jan 2, 2021
diff --git a/FeynCalc/LoopIntegrals/FCFindLoopMomentumShifts.m b/FeynCalc/LoopIntegrals/FCFindLoopMomentumShifts.m
@@ -0,0 +1,123 @@
+(* ::Package:: *)
+
+(* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ *)
+
+(* :Title: FCFindLoopMomentumShifts											*)
+
+(*
+	This software is covered by the GNU General Public License 3.
+	Copyright (C) 1990-2020 Rolf Mertig
+	Copyright (C) 1997-2020 Frederik Orellana
+	Copyright (C) 2014-2020 Vladyslav Shtabovenko
+*)
+
+(* :Summary:  	Finds loop momentum shifts									*)
+
+(* ------------------------------------------------------------------------ *)
+
+FCFindLoopMomentumShifts::usage =
+"FCFindLoopMomentumShifts[source, target, {p1, p2, ...}] finds loop momentum shifts \
+that bring loop integrals or topologies in the list source to the form specified \
+in target. The integrals/topologies in intFrom and intTo are assumed to be equivalent \
+and their denominators must be properly ordered via FCToPakForm. Furthermore, target \
+must be provided as a list of FeynAmpDenominator objects, while intFrom is a list of \
+such lists.";
+
+FCFindLoopMomentumShifts::failmsg =
+"Error! FCFindLoopMomentumShifts has encountered a fatal problem and must abort the computation. \
+The problem reads: `1`"
+
+Begin["`Package`"]
+End[]
+
+Begin["`FCFindLoopMomentumShifts`Private`"]
+
+fcflsVerbose::usage = "";
+
+Options[FCFindLoopMomentumShifts] = {
+	FCI 						-> False,
+	FCVerbose 					-> False
+};
+
+FCFindLoopMomentumShifts[fromRaw_List, toRaw_List, lmoms_List, OptionsPattern[]] :=
+	Block[{	from,to,pakFormInts, res, time, x, fromPak, toPak, shifts},
+
+		If[	OptionValue[FCVerbose] === False,
+			fcflsVerbose = $VeryVerbose,
+			If[MatchQ[OptionValue[FCVerbose], _Integer],
+			fcflsVerbose = OptionValue[FCVerbose]];
+		];
+
+		If[OptionValue[FCI],
+			{from, to} = {fromRaw, toRaw},
+			{from, to} = FCI[{fromRaw, toRaw}]
+		];
+
+		FCPrint[1, "FCFindLoopMomentumShifts: Entering.", FCDoControl -> fcflsVerbose];
+		FCPrint[3, "FCFindLoopMomentumShifts: Entering with: ", from, FCDoControl -> fcflsVerbose];
+
+		If[!MatchQ[to,{__FeynAmpDenominator}|FCTopology[_,{__FeynAmpDenominator},__]],
+			Message[FCFindLoopMomentumShifts::failmsg,""];
+		];
+
+		shifts = findShifts[#,to,lmoms]&/@from;
+
+
+		FCPrint[1, "FCFindLoopMomentumShifts: Leaving.", FCDoControl -> fcflsVerbose];
+		FCPrint[3, "FCFindLoopMomentumShifts: Leaving with: ", shifts, FCDoControl -> fcflsVerbose];
+
+		shifts
+	];
+
+
+findShifts[from:{__FeynAmpDenominator},to:{__FeynAmpDenominator}, lmoms_List]:=
+	Block[{lhs, rhs, eq,mark, vars, sol},
+		lhs = MomentumCombine[from,FCI->True];
+		rhs = MomentumCombine[to,FCI->True];
+		{lhs, rhs} = {lhs, rhs} /. {
+			FeynAmpDenominator[PropagatorDenominator[Momentum[mom_, _], _]] :> mom,
+			FeynAmpDenominator[StandardPropagatorDenominator[Momentum[mom_, _], 0, _, {1, _}]] :> mom,
+			FeynAmpDenominator[CartesianPropagatorDenominator[CartesianMomentum[mom_, _], 0, _, {1, _}]] :> mom,
+			FeynAmpDenominator[StandardPropagatorDenominator[0, x_, _, {1, _}]]/; x=!=0 :> Unevaluated[Sequence[]],
+			FeynAmpDenominator[CartesianPropagatorDenominator[0, x_, _, {1, _}]]/; x=!=0 :> Unevaluated[Sequence[]]
+
+		};
+
+		FCPrint[3, "FCFindLoopMomentumShifts: Preliminary lhs: ", lhs, FCDoControl -> fcflsVerbose];
+		FCPrint[3, "FCFindLoopMomentumShifts: Preliminart rhs: ", rhs, FCDoControl -> fcflsVerbose];
+
+		If[	!FreeQ[{lhs,rhs},FeynAmpDenominator],
+			Message[FCFindLoopMomentumShifts::failmsg,"Failed to set up a proper system of equations."];
+			Abort[]
+		];
+
+		vars = mark/@(lmoms);
+		lhs = lhs /. Thread[Rule[lmoms,vars]];
+
+		FCPrint[3, "FCFindLoopMomentumShifts: Final lhs: ", lhs, FCDoControl -> fcflsVerbose];
+		FCPrint[3, "FCFindLoopMomentumShifts: Final rhs: ", rhs, FCDoControl -> fcflsVerbose];
+
+		eq = Thread[Equal[lhs,rhs]];
+
+		sol = Solve[eq,vars];
+
+		If[	sol==={},
+			Message[FCFindLoopMomentumShifts::failmsg,"Failed to find momentum shifts for one of the topologies."];
+			Abort[]
+		];
+
+		(First[sol] /. mark -> Identity)
+
+	]
+
+
+
+
+
+
+
+
+
+
+FCPrint[1,"FCFindLoopMomentumShifts.m loaded."];
+End[]
diff --git a/FeynCalc/LoopIntegrals/FCFindPakMappings.m b/FeynCalc/LoopIntegrals/FCFindPakMappings.m
@@ -0,0 +1,79 @@
+(* ::Package:: *)
+
+(* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ *)
+
+(* :Title: FCFindPakMappings														*)
+
+(*
+	This software is covered by the GNU General Public License 3.
+	Copyright (C) 1990-2020 Rolf Mertig
+	Copyright (C) 1997-2020 Frederik Orellana
+	Copyright (C) 2014-2020 Vladyslav Shtabovenko
+*)
+
+(* :Summary:  	Obtains a canonical (Pak) representation of the given
+				FeynCalc integral 											*)
+
+(* ------------------------------------------------------------------------ *)
+
+FCFindPakMappings::usage =
+"FCFindPakMappings[{int1, int2, ...}, {p1, p2, ...}] finds mappings between scalar
+multiloop-integrals int1, int2, ... that depend on the loop momenta p1, p2, ...
+using using the algorithm of Alexey Pak (arXiv:1111.0868).
+
+The current implementation is based on the FindEquivalents function from \
+FIRE 6 (arXiv:1901.07808)";
+
+FCFindPakMappings::failmsg =
+"Error! FCFindPakMappings has encountered a fatal problem and must abort the computation. \
+The problem reads: `1`"
+
+Begin["`Package`"]
+End[]
+
+Begin["`FCFindPakMappings`Private`"]
+
+fcfpmVerbose::usage = "";
+
+Options[FCFindPakMappings] = {
+	CharacteristicPolynomial	-> Function[{U,F}, U+F],
+	FCE 						-> False,
+	FCI 						-> False,
+	FCVerbose 					-> False,
+	FinalSubstitutions			-> {},
+	Function					-> Function[{U, F, pows, head, int}, {int, head[ExpandAll[U], ExpandAll[F], Transpose[pows]]}]
+};
+
+FCFindPakMappings[expr_, lmoms_List, OptionsPattern[]] :=
+	Block[{	pakFormInts, res, time, x},
+
+		If[	OptionValue[FCVerbose] === False,
+			fcfpmVerbose = $VeryVerbose,
+			If[MatchQ[OptionValue[FCVerbose], _Integer],
+			fcfpmVerbose = OptionValue[FCVerbose]];
+		];
+
+		FCPrint[1, "FCFindPakMappings: Entering.", FCDoControl -> fcfpmVerbose];
+		FCPrint[3, "FCFindPakMappings: Entering with: ", expr, FCDoControl -> fcfpmVerbose];
+
+		time=AbsoluteTime[];
+		FCPrint[1, "FCFindPakMappings: Calling FCToPakForm.", FCDoControl -> fcfpmVerbose];
+		pakFormInts = FCToPakForm[#, lmoms, FCI->OptionValue[FCI], FinalSubstitutions->OptionValue[FinalSubstitutions],
+			Check->False, Collecting->False, Names->x, CharacteristicPolynomial->OptionValue[CharacteristicPolynomial],
+			Function->OptionValue[Function]] & /@ expr;
+		FCPrint[1, "FCFindPakMappings: FCToPakForm done, timing: ", N[AbsoluteTime[] - time, 4], FCDoControl->fcfpmVerbose];
+
+		res = Reap[Sow [Sequence @@ #] & /@ pakFormInts, _, ##2 &][[2]];
+
+		If[	OptionValue[FCE],
+			res = FCE[res]
+		];
+
+		FCPrint[3, "FCFindPakMappings: Leaving.", FCDoControl -> fcfpmVerbose];
+		FCPrint[3, "FCFindPakMappings: Leaving with: ", res, FCDoControl -> fcfpmVerbose];
+
+		res
+	];
+
+FCPrint[1,"FCFindPakMappings.m loaded."];
+End[]
diff --git a/FeynCalc/LoopIntegrals/FCPakOrder.m b/FeynCalc/LoopIntegrals/FCPakOrder.m
@@ -0,0 +1,173 @@
+(* ::Package:: *)
+
+(* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ *)
+
+(* :Title: FCPakOrder														*)
+
+(*
+	This software is covered by the GNU General Public License 3.
+	Copyright (C) 1990-2020 Rolf Mertig
+	Copyright (C) 1997-2020 Frederik Orellana
+	Copyright (C) 2014-2020 Vladyslav Shtabovenko
+*)
+
+(* :Summary:  	Finds a canonical ordering of the Feynman parameters in
+				a polynomial using Pak's algorithm							*)
+
+(* ------------------------------------------------------------------------ *)
+
+FCPakOrder::usage =
+"FCPakOrder[poly, {x1, x2, ...}] determines a canonical ordering of the Feynman \
+parameters x1, x2, ... in the polynomial poly using the algorithm of Alexey Pak \
+(arXiv:1111.0868). Cf. also the PhD thesis of Jens Hoff (Hoff:2015kub, \
+10.5445/IR/1000047447) for a more detailed description of the algorithm.
+
+The current implementation is based on the PolyOrdering function from \
+FIRE 6 (arXiv:1901.07808)";
+
+FCPakOrder::failmsg =
+"Error! FCPakOrder has encountered a fatal problem and must abort the computation. \
+The problem reads: `1`"
+
+Begin["`Package`"]
+End[]
+
+Begin["`FCPakOrder`Private`"]
+
+fcpoVerbose::usage = "";
+fpIds::usage = "";
+
+Options[FCPakOrder] = {
+	FCVerbose		-> False,
+	MaxIterations	-> Infinity,
+	MonomialOrder	-> Lexicographic
+};
+
+
+FCPakOrder[poly_, fpars_List, OptionsPattern[]] :=
+	Block[{	coeffsList, polyGCD, res, matM, mPrefs, time, time0},
+
+		If[	OptionValue[FCVerbose] === False,
+			fcpoVerbose = $VeryVerbose,
+			If[	MatchQ[OptionValue[FCVerbose], _Integer],
+				fcpoVerbose = OptionValue[FCVerbose]
+			];
+		];
+
+		time0=AbsoluteTime[];
+		FCPrint[1, "FCPakOrder: Entering.", FCDoControl -> fcpoVerbose];
+		FCPrint[3, "FCPakOrder: Entering with: ", poly, FCDoControl -> fcpoVerbose];
+
+		time=AbsoluteTime[];
+		FCPrint[1, "FCPakOrder: Calculating coefficient lists.", FCDoControl -> fcpoVerbose];
+		coeffsList = GroebnerBasis`DistributedTermsList[poly, fpars, MonomialOrder -> OptionValue[MonomialOrder]];
+		FCPrint[1, "FCPakOrder: Done calculating coefficient lists, timing: ", N[AbsoluteTime[] - time, 4], FCDoControl->fcpoVerbose];
+
+		If[	coeffsList[[2]] =!= fpars,
+			Print[coeffsList[[2]]];
+			Print[fpars];
+			Message[FCPakOrder::failmsg, "Incorrect polynomial variables."];
+			Abort[]
+		];
+
+		FCPrint[3, "FCPakOrder: coeffsList: ", coeffsList, FCDoControl -> fcpoVerbose];
+
+		(*
+			The tranposed form of M is more convenient here, since we have
+			a direct access to the vector of the prefactors.
+		*)
+		time=AbsoluteTime[];
+		FCPrint[1, "FCPakOrder: Building up M^0.", FCDoControl -> fcpoVerbose];
+		matM = Transpose[Sort[Flatten /@ First[coeffsList]]];
+		mPrefs = Last [matM];
+		FCPrint[1, "FCPakOrder: Done building up M^0, timing: ", N[AbsoluteTime[] - time, 4], FCDoControl->fcpoVerbose];
+
+		time=AbsoluteTime[];
+		FCPrint[1, "FCPakOrder: Dividing by a common factor.", FCDoControl -> fcpoVerbose];
+		polyGCD = PolynomialGCD@@mPrefs;
+		If[	polyGCD =!= 1,
+			mPrefs = Cancel[mPrefs/polyGCD];
+			matM = Join[Most[matM], {mPrefs}]
+		];
+		FCPrint[1, "FCPakOrder: Done dividing by a common factor, timing: ", N[AbsoluteTime[] - time, 4], FCDoControl->fcpoVerbose];
+		FCPrint[3, "FCPakOrder: matM: ", matM, FCDoControl -> fcpoVerbose];
+
+		If[	!MatrixQ[matM],
+			Message[FCPakOrder::failmsg, "Failed to build up the matrix M^0."];
+			Abort[]
+		];
+
+		fpIds = Range[Length[fpars]];
+
+		time=AbsoluteTime[];
+		FCPrint[1, "FCPakOrder: Applying pakSort.", FCDoControl -> fcpoVerbose];
+		res = FixedPoint[pakSort[#, matM] &, {{}}, OptionValue[MaxIterations]];
+		FCPrint[1, "FCPakOrder: Done applying pakSort, timing: ", N[AbsoluteTime[] - time, 4], FCDoControl->fcpoVerbose];
+		If[	!FreeQ[res,pakSort],
+			Message[FCPakOrder::failmsg, "Failed to determine a canonical ordering."];
+			Abort[]
+		];
+
+
+		FCPrint[1, "FCPakOrder: Toal timing: ", N[AbsoluteTime[] - time0, 4], FCDoControl->fcpoVerbose];
+		FCPrint[1, "FCPakOrder: Leaving.", FCDoControl -> fcpoVerbose];
+		FCPrint[3, "FCPakOrder: Leaving with: ", res, FCDoControl -> fcpoVerbose];
+
+		res
+	];
+
+pakSort[_pakSort, ___] :=
+	(
+	Message[FCPakOrder::failmsg, "Invalid set of arguments submitted to pakSort."];
+	Abort[]
+	);
+
+pakSort[sigma_List, matM_?MatrixQ] :=
+Block[{	swappedRowsList, relVectors, maxVector, res},
+
+
+	FCPrint[4, "FCPakOrder: pakSort: Entering with sigma: ", sigma, FCDoControl -> fcpoVerbose];
+
+	(*	Generate all row swaps allowed in this iteration	*)
+	swappedRowsList = Function[xx, Sequence@@Map[Join[xx, {#}] &, Complement[fpIds, xx]]]/@sigma;
+
+	FCPrint[4, "FCPakOrder: pakSort: ", swappedRowsList, FCDoControl -> fcpoVerbose];
+	(*
+		Using row swaps from swappedRowsList we sort the first k rows together with the monomial
+		coefficients(!) canonically (or in any other definite way). Notice that the
+		original algorithm works on what would be M = matM^T, hence we need to transpose twice here:
+		first get to M, then sort and finally return to matM^T.
+	*)
+	relVectors = (matM[[Prepend[#, -1]]]//Transpose//Sort//Transpose//Last) & /@ swappedRowsList;
+	FCPrint[4, "FCPakOrder: pakSort: relVectors: ", relVectors, FCDoControl -> fcpoVerbose];
+	(*
+		relVectors is a list of vectors corresponding to the kth row obtained from matrices with
+		different row swaps. We need to select a "maximum" vector (using some canonical ordering)
+	*)
+
+	maxVector = relVectors[[Last[Ordering[relVectors]]]];
+	FCPrint[4, "FCPakOrder: pakSort: maxVector: ", maxVector, FCDoControl -> fcpoVerbose];
+
+	(* Keep all permutations that lead to matrices containing the maximum vector *)
+	res = MapThread[If[(#1 === maxVector), #2, Unevaluated[Sequence[]]] &, {relVectors,	swappedRowsList}];
+
+	FCPrint[4, "FCPakOrder: pakSort: Leaving with sigma: ", res, FCDoControl -> fcpoVerbose];
+
+	(* Final list of permuations aka sigma *)
+	res
+	] /; Length[First[sigma]] < Length[fpIds];
+
+
+pakSort[sigma_List, _?MatrixQ] :=
+	sigma /; Length[First[sigma]] === Length[fpIds];
+
+pakSort[sigma_List, _?MatrixQ] :=
+	(
+	Message[FCPakOrder::failmsg, "The list of permutations got too long."];
+	Abort[]
+	) /;
+	Length[First[sigma]] >= Length[fpIds];
+
+
+FCPrint[1,"FCPakOrder.m loaded."];
+End[]