(*graphs*)
Block[{$Messages={},$Output={},$Echo={},Print=(Null &)},
Quiet[Needs["gRaph`"]]  
]

With[{und=constructEdgeDesc[0,"UNDIRECTED"]},
With[{exts={
	constructCoupling[1,{constructEdgeDesc[1,"DIRECTED","INCOMING"]}],
	constructCoupling[1,{constructEdgeDesc[1,"DIRECTED","OUTGOING"]}],
	constructCoupling[2,{constructEdgeDesc[2,"DIRECTED","INCOMING"]}],
	constructCoupling[2,{constructEdgeDesc[2,"DIRECTED","OUTGOING"]}]},
	ints={
	constructCoupling[0,{und,und,und}],
	constructCoupling[1,{und,constructEdgeDesc[1,"DIRECTED","INCOMING"],constructEdgeDesc[1,"DIRECTED","OUTGOING"]}],
	constructCoupling[2,{und,constructEdgeDesc[2,"DIRECTED","INCOMING"],constructEdgeDesc[2,"DIRECTED","OUTGOING"]}]}},
	gs[1]=constructTower[exts,ints,1];
	gs[2]=constructTower[exts,ints,2];
]
]

(*formatting*)
dot/:MakeBoxes[dot[a_,b_],StandardForm]:=RowBox[{ToBoxes[a,StandardForm],"·",ToBoxes[b,StandardForm]}]
eta/:MakeBoxes[eta[lI[idx1_],lI[idx2_]],StandardForm]:=SuperscriptBox["η",RowBox[{SubscriptBox["μ",idx1],SubscriptBox["μ",idx2]}]]
propD/:MakeBoxes[propD[arg_],StandardForm]:=ToBoxes[arg^2-i0,StandardForm]
k/:MakeBoxes[k[idx_],StandardForm]:=SubscriptBox["k",idx]
k/:MakeBoxes[k[idx_][lI[idx2_]],StandardForm]:=SubsuperscriptBox["k",ToBoxes[idx,StandardForm],SubscriptBox["μ",idx2]]
delta/:MakeBoxes[delta[args_],StandardForm]:=RowBox[{"δ[",ToBoxes[args],"]"}]

toOutput[arg_,align_:Center]:=CellPrint@Cell[BoxData@ToBoxes@arg,"Output",Deployed->True,TextAlignment->align]

ginacToMMA[expr_,epOrder_:0]:=(expr
	//.Dispatch[{D->4-2ep,gamma->g,q[0]->q,b[0]->b,l[i_]:>k[i+1],m1->m[1],m2->m[2],idxd[a_,b_]->a[lI[b]],_idx->1,deltaP->delta}]
	/.Dispatch[{int[moms_,{},ds_,ps_,num_,es_]:>int[moms//Sort,es//Factor,ds//Sort,ps//Sort,Collect[Series[num,{ep,0,epOrder}],{ep,_flag},Factor]]}]
);

How not to compute gravitational observables at 6PM

The art of abstraction and optimization

Gregor Kälin

YOUNGST@RS - Lost in Translation: The language of Gravitational Waves
January 2026

Motivation

Impulse @4 PM (log scale!)

data = {
  "integrand" -> 52.9,
  "intermediate data" -> 424.7,
  "after IBPs" -> 15.8,
  "final result" -> 0.044
};
plot=BarChart[data[[All,2]],ChartLabels->Placed[data[[All,1]],Below],ScalingFunctions->"Log",ChartStyle->ColorData[97]/@Range[4],ImageSize->620,BarSpacing->0.5,Frame->True,FrameLabel->{None,"Size[MiB]"},LabelStyle->Directive[12],PlotRangePadding->Scaled[0.1]];
toOutput[plot]

({{, "Premature optimization is the root of all evil, or at least most evil in programming; 
the true issue is that programmers have worried about efficiency much too much, 
in the wrong places, and at the wrong times."  ---Donald Knuth}})

Overview

From model to integrands

Worldline formalism

Integrand construction: Feynman rules

Organization and representation

Feynman integral reduction

Integral families

Laporta & symbolic strategies

Integral parametrization: Having some fun with i0s

Assembling the final answer

Worldline formalism (PMEFT & WQFT)

toOutput[Row[{Style["1PM:", 30],Style[TraditionalForm[((Subscript[δS, eff])/(SubMinus[δx]))],30],Style[" = ",30],Show[Import["data/graphsRad1.svg","Graphics"],ImageSize->130]},Spacer[5]],Left]
toOutput[Row[{Style["2PM:", 30],Style[TraditionalForm[((Subscript[δS, eff])/(SubMinus[δx]))],30],Style[" = ",30],Show[Import["data/graphsRad2.svg","Graphics"],ImageSize->300]},Spacer[5]],Left]
toOutput[Row[{Style["3PM:", 30],Style[TraditionalForm[((Subscript[δS, eff])/(SubMinus[δx]))],30],Style[" = ",30],Show[Import["data/graphsRad3.svg","Graphics"],ImageSize->450]},Spacer[5]],Left]

toOutput[Show[Import["data/sourceRule1.png","Graphics"],ImageSize->400],Left]
toOutput[Show[Import["data/sourceRule2.png","Graphics"],ImageSize->1000],Left]

toOutput[Framed[Show[Import["data/Dp.png","Graphics"],ImageSize->900]]]

Purely mechanical!

(Note: Graph construction is very simple! No need to optimize!)

Integrand construction

Get["data/rules.dat"];

propDDim[{{lI[1],lI[2]},{lI[3],lI[4]}},1]

vGR[{{lI[1],lI[2]},{1,1}},{{lI[3],lI[4]},{1,2}},{{lI[5],lI[6]},{1,3}}]

vGR[{{lI[1], lI[2]}, {1, 1}}, {{lI[3], lI[4]}, {1, 2}}, {{lI[5], lI[6]}, {1, 3}},{{lI[7], lI[8]}, {1, 4}}]

vGR[{{lI[1], lI[2]}, {1, 1}}, {{lI[3], lI[4]}, {1, 2}}, {{lI[5], lI[6]}, {1, 3}}, {{lI[7], lI[8]}, {1, 4}},{{lI[9], lI[10]}, {1, 5}}]

Doesn' t  look  great, but: Systematic, parallelizable, predictable, straight forward!

Challenges, tricks, and optimization potential

Contractions lead to combinatorial blow up

Use dedicated programs/code: FORM/custom C++ implementation beats xAct/MMA by orders of magnitudes

order of contraction matters

Example: 3PM integrand. 
- Multiply everything together and contract: 1min+
- Glue and contract one-by-one, random order: 19s
- Glue and contract, starting from smallest factor, keep disconnected pieces disconnected, one-by-one: 6.5s

Feynman rules optimized for algebraic size

and not for the physics,

and not for integration.

Missing guiding principles?!

Integrand construction: The optimal integrand

Split into smallest possible (gauge invariant) building blocks

Minimal information, no duplication

Algebraically small

Optimal for integration

Optimize integrand for integration or optimize integration process for given integrand?

Finiteness

Reduce max power of loop momenta

Factorize loops?

Ideal for numerics?

Efficient to derive

IMHO: We are far away from optimality regarding any of these points!

Ideas & thoughts

toOutput[WordCloud[{{"CK duality & double copy",8},{"On-shell methods",10},{"Gauge invariance",10},{"Recursive construction",8},{"Self force",9},{"finiteness",8},{"loop momentum labelling",8},{"canonicalization",7},{"integrand database",8},{"low-level implementation",9}},ImageSize->600]]

Integrand: Organization and Representation

How we do it:

ginacToMMA[Get["data/Dp_WoLF.m"][[1]]]

How we should do it (abstract!):

toOutput[Row[{drawGraph[gs[1][[8]],"size"->200],Style[" x ",30],Style[" (numerator) ",30]},Spacer[5]]]

toOutput[Show[Import["data/y_tho_meme.jpg","Graphics",ImageSize->200]]]

Natural grouping (=?= gauge-invariant building blocks)

Graph isomorphism -> canonicalization -> symmetries & minimality

Ideal for IBPs

Challenges, tricks, and optimization potential

Need (causal) i0 prescription: non-trivial book keeping

graphs would “trivialize” that too!

What counts as same topology?

propagators vs. ISPs

subtopologies

Algebraic size becomes annoying at 5PM and challenging at 6PM

I repeat: Need better guiding principle!

Perhaps: Topology

Data representation

With millions of objects, data representation starts to matter.

Example:

Linear combination of 10 million fixed-size objects (64 bytes) with 64-bit integer coefficients (8 bytes)

Contiguous memory (e.g. C/C++ array/vector)

UnitConvert[Quantity[(10^7)(64+8),"Bytes"], "GigaBytes"]//N

Red-black tree (e.g. C++ std::map)

UnitConvert[Quantity[(10^7)(64+8+64(*overhead per element!*)),"Bytes"], "GigaBytes"]//N

Mathematica List / Association

{(10^54)//ByteCount,1//ByteCount}
rands=Table[RandomInteger[{0,(10^54)}],{(10^7)}];
UnitConvert[Quantity[rands->1//Thread//ByteCount,"Bytes"], "GigaBytes"]//N
UnitConvert[Quantity[rands->1//Thread//Association//ByteCount,"Bytes"], "GigaBytes"]//N
Clear[rands];

Integral families I

All-loop structure (for Δp):

toOutput[Show[Import["data/int.png","Graphics"],ImageSize->800]]

Want to map everything to as few families as possible. Not that trivial:

(*set of squared propagators at n loop, m'th family*)
sqPropMoms[1,1]={l[1],l[1]-q};
sqPropMoms[2,1]={l[1],l[2],l[1]+l[2]-q,l[1]-q,l[2]-q};
sqPropMoms[3,1]={l[1],l[2],l[3],l[1]-q,l[2]-q,l[3]-q,l[1]-l[2],l[2]-l[3],l[3]-l[1]};
sqPropMoms[3,2]={l[1],l[2],l[3],l[1]-q,l[2]-q,l[3]-q,l[1]+l[2]-q,l[2]+l[3]-q,l[1]+l[2]+l[3]-q};

sqPropsString[1,1]={"-l1^2", "-(l1-q)^2"};
sqPropsString[2,1]={"-l1^2", "-l2^2", "-(l1+l2-q)^2", "-(l1-q)^2", "-(l2-q)^2"};
sqPropsString[3,1]={"-l1^2", "-l2^2", "-l3^2", "-(l1-q)^2", "-(l2-q)^2", "-(l3-q)^2", "-(l1-l2)^2", "-(l2-l3)^2", "-(l3-l1)^2"};
sqPropsString[3,2]={"-l1^2", "-l2^2", "-l3^2", "-(l1-q)^2", "-(l2-q)^2", "-(l3-q)^2", "-(l1+l2-q)^2", "-(l2+l3-q)^2", "-(l1+l2+l3-q)^2"};

Options[toInt]={"PrintMap"->False};
toInt::usage = "toInt[expr_] will try to map all _int objects within the expr into the predefined integral families, i.e. the will become int[spec__][args__].
Note that it may use the u[1]<->u[2] symmetry.";
toInt::indexedErr = "Can only be used on scalar integrals, make sure to remove all tensorial quantities depending on the loop momenta.";
toInt::expErr = "Exponential doesn't have the correct form.";
toInt::deltaErr = "q deltas seem to be missing or have the wrong form.";
toInt::notFound = "No suitable topology found.";
SetAttributes[toInt,Listable];
toInt[0]=0;
toInt[expr_Plus,opts : OptionsPattern[]]:=toInt[#,opts]&/@expr
toInt[expr_Times,opts : OptionsPattern[]]:=toInt[#,opts]&/@expr
toInt[expr_?(FreeQ[#,_int]&),opts : OptionsPattern[]]:=expr
toInt[int[moms_,es_,ds_,ps_,arg_],OptionsPattern[]]:=Module[{tmp=prodToList@Factor[arg],deltas=ds,props=ps,coeff,i,kk,map},
	If[!FreeQ[arg,_?(MemberQ[DeleteCases[moms,q],#]&)[_lI]],Message[toInt::indexedErr];Abort[];];
	If[es=!=Exp[-I dot[q,b]],Message[toInt::expErr];Abort[];];
	If[Sort[Select[deltas,!FreeQ[#,q]&]]=!={delta[dot[q,u[1]]],delta[dot[q,u[2]]]},Message[toInt::deltaErr];Abort[];];
	(*pull as much as possible out of the integral*)
	coeff=Select[tmp,AllTrue[Append[moms,_log(*for nonlocal*)],Function[y,FreeQ[#,y]]]&];
	tmp=Complement[tmp,coeff];
	(*pull out delta[q,u]*)
	deltas=Select[deltas,FreeQ[#,q]&];
	(*pull out dot[q,q]*)
	coeff=Join[coeff,Cases[props,propR[q]^_.|propR[-q]^_.|propA[q]^_.|propA[-q]^_.|propF[q]^_.|propF[-q]^_.|propD[q]^_.|propD[-q]^_.]]/.{propR[q]|propR[-q]|propA[q]|propA[-q]|propF[q]|propF[-q]|propD[q]|propD[-q]->1/dot[q,q]};
	props=getPower/@DeleteCases[props,propR[q]^_.|propR[-q]^_.|propA[q]^_.|propA[-q]^_.|propF[q]^_.|propF[-q]^_.|propD[q]^_.|propD[-q]^_.];
	With[{
		ks=DeleteCases[moms,q],
		nn=Length[DeleteCases[moms,q]],
		(*don't try to change the numbers associated to the four propagators, they matter in the final steps when writing out the ints*)
		sqProps=Cases[props,{_propR|_propA|_propF|_propD,_}]/.{propR[a_]:>Sequence[a,1],propA[a_]:>Sequence[a,-1],propF[a_]:>Sequence[a,2],propD[a_]:>Sequence[a,-2]},
		linProps=Cases[props,{_propLin,_}]/.{propLin[a_,b_]:>Sequence[a,b/i0]}
		},
		map=findFamily[ks,First/@deltas,linProps[[All,1]],sqProps[[All,1]]];
		If[OptionValue["PrintMap"],Print[map[[3]]]];
		(*Jacobian/extra factor*)(map[[4]]/.{jac->Identity,extra[a_]:>(If[a=!=1,Print["Warning: unexpected factor of ",a,"appearing from delta function shifts."]];a)})
		(*-1's from linear props where we had to swap sign*)(-1)^Total[linProps[[First/@Position[map[[6]][[All,2]],-1],3]]]
		(*introducing a minus sign in all square props (WoLF doesn't have it)*)(-1)^Total[sqProps[[All,3]]]
		(Times@@coeff)*int[nn,map[[1]],map[[2]]][
			(*numerator*)
			(Times@@tmp/.Join[(*(dressIndex[#[[1]],lI[i_]]->dressIndex[#[[2]],lI[i]])&/@map[[3]],*)map[[3]]])/.If[
				map[[8]],
				{u[1]->u[2],u[2]->u[1](*do not swap uu[1] and uu[2] which only appear in indexed objects, they are technically outside of the int*)},{}
			],
			(*deltas*)
			1&/@deltas,
			(*linear props*)
			ReplacePart[Table[0,{nn}],#[[1,1]]->Total[#[[All,2]]]&/@GatherBy[
				MapIndexed[(#1[[1]]->If[#1[[2]]linProps[[First[#2],2]]===1,{linProps[[First[#2],3]],0,0},{0,linProps[[First[#2],3]],0}])&,map[[6]]],
				First
			]],
			(*square props*)
			ReplacePart[Table[0,{nProps[nn]-2nn}],#[[1,1]]->Total[#[[All,2]]]&/@GatherBy[
				MapIndexed[(#1[[1]]->Switch[#1[[2]]^sqProps[[First[#2],2]](*<- for feyn and dyson this sign doesn't matter, so it will be to power 2; clever, right?*)sqProps[[First[#2],2]],			
					1,ret[sqProps[[First[#2],3]]],
					-1,adv[sqProps[[First[#2],3]]],
					2,fey[sqProps[[First[#2],3]]],
					-2,dys[sqProps[[First[#2],3]]]])&,map[[7]]],
				First
			]]
		]
	]
]

findFamily::usage = "findFamily[loopMoms_,deltaArgs_,linProps_,sqProps_] returns {family, u's in deltas, loopMom relabelling, Jacobian, deltaArgs permutation with sign, lin prop permutation with sign, square prop permutation with sign, bool if u1<->u2 needs to be swapped in numerator}. 
props need to be given in WoLF language, i.e. propLin[_] or propR|propA|propF|propD[_] and with getPower applied to them, e.g. as a pair {propLin[_],power_}.";
findFamily::invalidUStructure = "Invalid u structure in deltas.";
findFamily::noLinIndepSqProps = "Couldn't find a basis of linearly independent square propagator arguments.";
findFamily[loopMoms_,deltaArgs_,linProps_,sqProps_]:=Module[
	{nn=Length[loopMoms],ls=Array[l,Length[loopMoms]],ansatz,us=Cases[#,_u,∞]&/@deltaArgs,usCompl,propsChosen={},deltaCnstrnts,extraFactor,targetDeltaArgs,targetLinArgs,swapu1u2=False},
	(*find the u sector*)
	If[Length[#]=!=1,Message[findFamily::invalidUStructure];Abort[]]&/@us;
	us=First/@us//Sort;
	If[
		Count[us,u[1]]<Count[us,u[2]],
		(*use u1<->u2 symmetry*)
		swapu1u2=True;
		us=us/.{u[1]->u[2],u[2]->u[1]}//Sort;
	];
	deltaCnstrnts=Join[{dot[q,u[1]]->0,dot[q,u[2]]->0},#->0&/@(dot@@@Transpose[{ls,us}])];
	us=First/@us;
	usCompl=us/.{1->2,2->1};
	targetDeltaArgs=dot@@@Transpose[{ls,u/@us}];
	targetLinArgs=dot@@@Transpose[{ls,u/@usCompl}];
	(*get the square (retarded) propagator momenta*)
	(*choose one (square) propagator momentum for each loop momentum to match afterwards*)
	Module[{loopMomsTmp=loopMoms,sqPropsTmp=SortBy[Union[canonicalizePoly/@sqProps](*getting rid of overall signs and duplicates due to different i0s*),LeafCount]},
		While[loopMomsTmp=!={},
			If[Length[sqPropsTmp]===0,Message[findFamily::noLinIndepSqProps];Abort[];];(*TODO: this might happen for other applications: maybe one loop momentum doesn't appear in any independent square prop?*)
			With[{lMom=SortBy[Function[{x},{Length[Select[sqPropsTmp,!FreeQ[#,x]&]],x}]/@loopMomsTmp,First][[1,2]]},
				With[{propPos=First@FirstPosition[sqPropsTmp,_?(!FreeQ[#,lMom]&)]},
					AppendTo[propsChosen,sqPropsTmp[[propPos]]];
					sqPropsTmp=Delete[sqPropsTmp,propPos];
				];
				(*step*)
				loopMomsTmp=DeleteCases[loopMomsTmp,lMom];
				(*check that they're linearly independent*)
				If[MatrixRank[Coefficient[propsChosen,#]&/@Complement[loopMoms,loopMomsTmp]]<Length[propsChosen],
					(*then: put only the lMom back, not the sqPropsTmp. TODO: maybe do better, the sqPropsTmp might still be usable for another lMom*)
					AppendTo[loopMomsTmp,lMom];
					propsChosen=Most[propsChosen];
				];
			];
		]
	];
	Catch[
		Do[
			(*make sure that every loop momentum appears at least once*)
			If[AnyTrue[ls,FreeQ[propChoice,#]&],Continue[];];
			doPermutation[
				propChoice,
				Function[{propChoiceX},
					Catch[
						With[{sol=Solve[#==0&/@(propChoiceX-propsChosen*signs),loopMoms]},
							If[Length[sol]=!=1||(*Jacobian*)Simplify[Abs[Det[Coefficient[sol[[1,All,2]],#]&/@ls]]]===0,Throw[1,"failed"]];
							With[{
								sqPs=sqProps/.sol[[1]]//Expand,
								linPs=linProps/.If[swapu1u2,{u[1]->u[2],u[2]->u[1]},{}]/.sol[[1]]/.deltaCnstrnts//Expand,
								ds=With[{deltaConstr=useDeltaConstraints[1,delta/@(deltaArgs/.If[swapu1u2,{u[1]->u[2],u[2]->u[1]},{}]/.sol[[1]]/.{dot[q,u[1]]->0,dot[q,u[2]]->0})]},
									If[deltaConstr===$Failed,{},With[{tmp=Times@@(deltaConstr[[2]])},extraFactor=tmp/.{_delta->1};First/@Cases[tmp,_delta,{0,∞}]]]
								]
								},
								(*try to find a map*)
								Throw[{
									family,
									us,
									sol[[1]],
									(*Jacobian&extraFactor*)jac[Det[Coefficient[sol[[1,All,2]],#]&/@ls]]^-4 extra[extraFactor],
									(*deltas*)
									Map[With[{tmp=FirstPosition[targetDeltaArgs,#]},If[tmp===Missing["NotFound"],(*not a correct mapping*)Throw[1,"failed"],First[tmp]]]&,ds],
									(*linear propagators*)
									Map[With[{tmp=pmPos[targetLinArgs,#]},If[tmp[[1]]===Missing["NotFound"],(*not a correct mapping*)Throw[1,"failed"],tmp]]&,linPs],
									(*square propagators*)
									Map[With[{tmp=pmPos[sqPropMoms[nn,family],#]},If[tmp[[1]]===Missing["NotFound"],(*not a correct mapping*)Throw[1,"failed"],tmp]]&,sqPs],
									swapu1u2
								},"final"];
							]
						]
					,"failed"];
				]
			]
		,{family,1,nFamilies[nn]},{propChoice,Subsets[sqPropMoms[nn,family],{nn}]},{signs,Tuples[{1,-1},{nn}]}];
		Print[{loopMoms,deltaArgs,linProps,sqProps}];Print["Could not map to any family."];Abort[];,
		"final"
	]
]

But: This is essentially just some (terrible implementation of) graph isomorphism.

Integral families II

With[{
n1="How to choose
integral families",
n2="If no match,
define new
family",
n3="Extend from
lower loops",
n4="symmetric",
n5="algebraically
small",
n6="topology",
n7="What about
ISPs?",
n8="minimal
building
blocks?",
nRand="random",
nFinal="reason for
inefficiency
of reduction?"},
Graph[{n1->n2,n1->n3,n1->n4,n1->n5,n1->n6,n6->n7,n7->n8,n2->nRand,n3->nRand,n4->nRand,n5->nRand,n7->nRand,nRand->nFinal},
VertexLabels->Placed["Name",Center],VertexStyle->{n8->Directive[LightRed,EdgeForm[Black]],n1->Directive[LightBlue,EdgeForm[Black]]},VertexSize->0.65,VertexLabelStyle->Directive[Black,14],GraphLayout->{"LayeredDigraphEmbedding"},ImageSize->900]
]//toOutput

Integral reduction: Laporta style

Efficient, automated implementations exist:

toOutput@Grid[{
  {Style["FIRE", Bold, 16], Style["Kira", Bold, 16]},
  {"automatic with good default choices", "manual seeding"},
  {Column[{"top-down:", "level-by-level", "sector-by-sector"}], Column[{"all at once:","forward solving","then back substitution"}]},
  {"read/write to database on disk", "all in memory"},
  {Column[{"finite field evaluation &","rational reconstruction"}],Column[{"finite field evaluation &","rational reconstruction"}]},
  {Column[{"rational arithmetics:","Fermat","Flint","Symbolica"}],Column[{"rational arithmetics:","Fermat","Flint"}]}
},
  Frame -> All,
  Background -> {None, {LightBlue, White, White, White}},
  Alignment -> Left,
  Dividers -> {None, {2 -> True}},
  Spacings -> {2, 1}
]

Challenges, tricks, and optimization potential

Both programs have parameters to fine tune.

We would like to have optimal seeding = smallest possible system

Problem-dependent or generic solution?

Optimal master integrals (MIs); depend on integral ordering

One-scale problem: Factorized denominators: rational reconstruction -> polynomial reconstruction

Factorial growth unavoidable

Integral reduction: LiteRed style II

(*Loading packages, setting directories*)
dataDir=NotebookDirectory[]<>"data/";
Needs["LiteRed2`"];
Needs["Fermatica`"];
$FermatCMD="/usr/bin/fer64";
SetDirectory[NotebookDirectory[]];

(*defining variables & kinematics*)
SetDim[d];
Declare[{l1,l2,q,u1,u2},Vector];
Declare[{x},Number]
SetConstraints[
sp[q,q]=-1;
sp[u1,q]=0;
sp[u2,q]=0;
sp[u1,u2]=(1+x^2)/(2 x);
sp[u1,u1]=1;
sp[u2,u2]=1;
];

(*define the integral family*)
NewDsBasis[
fam1s12,
{sp[l1,u1],sp[l2,u2],(*delta functions*)
sp[l1,u2],sp[l2,u1],(*linear propagators*)
-sp[l1,l1],-sp[l2,l2],-sp[l1+l2-q,l1+l2-q],-sp[l1-q,l1-q],-sp[l2-q,l2-q](*square propagators*)
},{l1,l2},
Directory->"fam1s12/",GenerateIBP->True];

(IBP[fam1s12]@@Array[n,9])[[1]]

(*Find zero sectors*)
AnalyzeSectors[fam1s12,CutDs->{1,1,0,0,0,0,0,0,0}]

ZeroSectors[fam1s12]

NonZeroSectors[fam1s12]

(*Find symbolic rules, sector-by-sector*)
NonZeroSectors[fam1s12][[1]]
SolvejSector[%,UseFermat->True,Depth->2]

jRules[fam1s12,1,1,0,0,0,0,1,1,1][[1;;5]]//TableForm

MIs[fam1s12]

(*perform a reduction*)
IBPReduce[j[fam1s12,1,3,0,0,0,0,3,2,3]]

Integral reduction: LiteRed style II

With[{
n1=Style["SpIRed",Bold],
n2="symbolic
rules",
n3="reusable",
n4="forward
solving:
trivial &
minimal",
n5="memory
efficient",
n6="finding rules
is hard",
n7="finite field
arithmetics",
n8="machine
learning",
n9="sparse
algorithms",
n13="...",
n11="trivially
parallelizable",
n12="identify
difficulties",
n10="syzygies"
},
Graph[{n1->n2,n2->n3,n2->n4,n1->n5,n1->n6,n6->n7,n6->n8,n6->n9,n6->n10,n1->n11,n1->n12,n6->n13},
VertexLabels->Placed["Name",Center],VertexStyle->{n1->Directive[LightBlue,EdgeForm[Black]],n2|n5|n6|n11|n12->Directive[LightGreen,EdgeForm[Black]],n7|n9|n10->Directive[LightRed,EdgeForm[Black]]},VertexSize->0.83,VertexLabelStyle->Directive[Black,14],GraphLayout->{"RadialEmbedding"},ImageSize->900]
]//toOutput

Integral parametrization: Having some fun with i0s

∫(d^d)l (1/((((l.n+i0))^(ν[1]))(((l.n-i0))^(ν[2]))((((((l.n))^2)))^(ν[3]))((((l^2)))^(ν[4]))((((((l-q))^2)))^(ν[5]))))
d=3-2ep, n={1,0,0},n.n==1, q.q==1, q.n==0

schwingerParam/:schwingerParam[coeff1_,exp1_]*schwingerParam[coeff2_,exp2_]:=schwingerParam[coeff1*coeff2,exp1+exp2]
dot[a___,b_Plus,c___]:=dot[a,#,c]&/@b
dot[a___,-b_,c___]:=-dot[a,b,c]
dot[a_,b_]/;Order[b,a]==1:=dot[b,a]
dot[a___,q*b_,c___]:=b dot[a,q,c]
dot[___,0,___]:=0

plusToList[expr_Plus]:=List@@expr
plusToList[expr_List]:=expr
plusToList[expr_]:={expr}

toFeynman[sum_Plus]:=toFeynman/@sum
toFeynman[schwingerParam[coeff_,exp_]]:=Module[{η},
	With[{m=Exponent[Simplify[coeff/.{x[i_]:>η x[i],i0->0},η>0],η]-1+Length[Union[Cases[{coeff,exp},x[_],∞]]]},
		coeff*Gamma[m+1]Exp[(I π)/2 (m+1)](exp/I+i0)^(-m-1)
	]
]
toFeynman[expr_Times]:=With[{exp=plusToList[FirstCase[expr,Exp[_]][[2]]//Expand]},toFeynman[schwingerParam[Exp[Total[Select[exp,FreeQ[#,x[_]]&]]](expr/.{Exp[_]:>1}),Total[Select[exp,!FreeQ[#,x[_]]&]]]]]

1. Parametrize

param=schwingerParam[Exp[-(I π ν[1])/2]/Gamma[ν[1]] x[1]^(ν[1]-1),I x[1](l[1][x]+i0)]*
schwingerParam[Exp[(I π ν[2])/2]/Gamma[ν[2]] x[2]^(ν[2]-1),-I x[2](l[1][x]-i0)]*
(*choose random i0*)schwingerParam[Exp[(I π ν[3])/2]/Gamma[ν[3]] x[3]^(ν[3]-1),-I x[3](l[1][x]^2-i0)]*
schwingerParam[Exp[(I π ν[4])/2]/Gamma[ν[4]] x[4]^(ν[4]-1),-I x[4](l[1][x]l[1][x]+dot[l[1],l[1]]-i0)]*
schwingerParam[Exp[(I π ν[5])/2]/Gamma[ν[5]] x[5]^(ν[5]-1),-I x[5](l[1][x]l[1][x]+dot[l[1]-q,l[1]-q]-i0)]/.{dot[q,q]->1}//Simplify;
param[[1]]*Exp[param[[2]]]

2. Perform the loop integrals

param2=With[{coeff=param[[1]],exp=param[[2]]/.{i0->0},nLoops=1},
	With[{
		aax=1/I Table[Coefficient[exp,l[ii][x]l[jj][x]],{ii,1,nLoops},{jj,1,nLoops}]//Simplify,
		bbx=1/I Table[1/2 Coefficient[exp,l[ii][x]],{ii,1,nLoops}]/.{l[_][x]:>0}//Simplify,
		ccx=exp/.{l[_][x]:>0}//Simplify
		},
		With[{
			aa=1/I Table[If[ii===jj,1,1/2]Coefficient[ccx,dot[l[ii],l[jj]]],{ii,1,nLoops},{jj,1,nLoops}]//Simplify,
			bb=1/I Table[1/2 Coefficient[ccx/.{dot[a___,b_l,c___]:>Times[a,b,c]},l[ii],1],{ii,1,nLoops}]/.{_l->0}//Simplify,
			cc=ccx/.{_l->0}//Simplify
			},
			schwingerParam[coeff*Exp[I d π/4] Product[(ev+i0)^(-1/2),{ev,Eigenvalues[aax]}] Product[(ev+i0)^(-(d-1)/2),{ev,Eigenvalues[aa]}],-I Sum[Inverse[aax][[ii,jj]]bbx[[ii]]bbx[[jj]],{ii,1,nLoops},{jj,1,nLoops}]-I Sum[Inverse[aa][[ii,jj]]dot[bb[[ii]],bb[[jj]]],{ii,1,nLoops},{jj,1,nLoops}]+cc]
		]
	]
]/.{dot[n,n]->1,dot[q,q]->1,dot[n,q]->0}//Simplify;

param2[[1]]*Exp[param2[[2]]]

3. Feynman trick

expr1=(toFeynman[param2]//Simplify)/.{
	(i0-x[4]-x[5])^(1/2-d/2)->(x[4]+x[5])^(1/2-d/2) Exp[I π (1/2-d/2)],
	(i0-x[3]-x[4]-x[5])^(-1/2)->(x[3]+x[4]+x[5])^(-1/2) Exp[I π (-1/2)],
	(i0+1/4 (-4 x[5]+(4 x[5]^2)/(x[4]+x[5])+(x[1]-x[2])^2/(x[3]+x[4]+x[5])))^exp_->(i0+(x[1]^2 (x[4]+x[5])-2 x[1] x[2] (x[4]+x[5])+x[2]^2 (x[4]+x[5])-4 x[4] x[5] (x[3]+x[4]+x[5])))^exp/(4^exp (x[4]+x[5])^exp (x[3]+x[4]+x[5])^exp)
}//Simplify

4. Perform x[1] and x[2] integral using a trick

tmp=HeavisideTheta[x[1]]HeavisideTheta[x[2]]Select[expr1,!FreeQ[#,x[1]]||!FreeQ[#,x[2]]&]/.{x[1]->x[1]+x[2]}//Simplify;
expr2Part1=Assuming[{x[1]>0},Integrate[tmp,{x[2],-∞,∞}]][[1]]
expr2Part1=Assuming[{Im[i0]>0,Re[i0]==0,x[3]>0,x[4]>0,x[5]>0},Integrate[expr2Part1,{x[1],0,∞}]][[1]]

expr2Part2=Assuming[{x[1]<0},Integrate[tmp,{x[2],-∞,∞}]][[1]]
expr2Part2=Assuming[{Im[i0]>0,Re[i0]==0,x[3]>0,x[4]>0,x[5]>0},Integrate[expr2Part2,{x[1],-∞,0}]][[1]]

expr2=(Select[expr1,FreeQ[#,x[1]]&&FreeQ[#,x[2]]&](expr2Part1+expr2Part2)//Simplify)/.{
	(1/(i0-4 x[4] x[5] (x[3]+x[4]+x[5])))^exp_->Exp[-I π exp]4^-exp x[4]^-exp x[5]^-exp (x[3]+x[4]+x[5])^-exp
}//Simplify

Takes too long to run, so here’s the result:

expr2=((2^(-1+ν[1]+ν[2])) (E^(-(1/2) I π (ν[1]-ν[2]))) Gamma[1-ν[1]-ν[2]] (Gamma[1-ν[1]] Gamma[ν[1]]+Gamma[1-ν[2]] Gamma[ν[2]]) Gamma[(1/2) (ν[1]+ν[2])] Gamma[(1/2) (-d+ν[1]+ν[2])+ν[3]+ν[4]+ν[5]] ((x[3])^(-1+ν[3])) ((x[4])^((1/2) (-2+d-ν[1]-ν[2]-2 ν[3]-2 ν[5]))) ((x[5])^((1/2) (-2+d-ν[1]-ν[2]-2 ν[3]-2 ν[4]))) (((x[4]+x[5]))^((1/2) (1-2 d+ν[1]+ν[2]+2 ν[3]+2 ν[4]+2 ν[5]))) (((x[3]+x[4]+x[5]))^((1/2) (-1+ν[1]+ν[2]))))/(Gamma[1-ν[1]] Gamma[ν[1]] Gamma[1-ν[2]] Gamma[ν[2]] Gamma[ν[3]] Gamma[ν[4]] Gamma[ν[5]]);

5. Integrate x[3]

expr3=Select[expr2,FreeQ[#,x[3]]&]Assuming[{x[4]>0,x[5]>0},Integrate[Select[expr2,!FreeQ[#,x[3]]&],{x[3],0,∞}][[1]]]//Simplify

6. Cheng-Wu Theorem: x[4] -> 1; integrate x[5]

expr3/.{x[4]->1}//Simplify;
finalRes=Select[%,FreeQ[#,x[5]]&]Integrate[Select[%,!FreeQ[#,x[5]]&],{x[5],0,∞}][[1]]//Simplify//FullSimplify

finalResBackup=(E^(-(1/2) I π (ν[1]-ν[2])) Csc[1/2 π (ν[1]+ν[2])] Gamma[1/2 (1-ν[1]-ν[2]-2 ν[3])] Gamma[1/2 (d-ν[1]-ν[2]-2 (ν[3]+ν[4]))] Gamma[1/2 (-d+ν[1]+ν[2])+ν[3]+ν[4]+ν[5]] Gamma[1/2 (d-ν[1]-ν[2]-2 (ν[3]+ν[5]))] (Sin[π ν[1]]+Sin[π ν[2]]))/(2 Sqrt[π] Gamma[ν[4]] Gamma[d-ν[1]-ν[2]-2 ν[3]-ν[4]-ν[5]] Gamma[ν[5]]);

7. Some interesting relations

(I^(+-)):=∫(d^d)l (1/((((l.n+i0))^(SubPlus[ν]))(((l.n-i0))^(SubMinus[ν]))(...))),
(I^(++:=))∫(d^d)l (1/((((l.n+i0))^((SubPlus[ν])+(SubMinus[ν])))(...)))
(I^(--:=))∫(d^d)l (1/((((l.n-i0))^((SubPlus[ν])+(SubMinus[ν])))(...)))
(I^(+-))==((Sin[π (SubPlus[ν])])/(Sin[π((SubPlus[ν])+(SubMinus[ν]))]))(e^(i π (SubMinus[ν])))(I^(++))+((Sin[π (SubMinus[ν])])/(Sin[π((SubPlus[ν])+(SubMinus[ν]))]))(e^(-i π (SubPlus[ν])))(I^(--))

finalRes-(((Sin[π ν[1]])/(Sin[π(ν[1]+ν[2])]))Exp[I π ν[2]](finalRes/.{ν[2]->0,ν[1]->ν[1]+ν[2]})+((Sin[π ν[2]])/(Sin[π(ν[1]+ν[2])]))Exp[-I π ν[1]](finalRes/.{ν[1]->0,ν[2]->ν[1]+ν[2]}))//FullSimplify

((((((l.n))^2)))^(ν[3]))->((((((l.n))^2)-(i0^2)))^(ν[3]))=(((l.n+i0))^(ν[3]))(((l.n-i0))^(ν[3]))

finalRes-(finalRes/.{ν[1]->ν[1]+ν[3],ν[2]->ν[2]+ν[3],ν[3]->0})//FullSimplify

Assembling: Proposition for an ideal pipeline

(*find symbolic reduction rules*)
{symbolicRules,MIs} = SpIRed[theFamily,neededSectors,ZeroSectors->True,Symmetries->True];

(*find the differential equation matrix*)
DE=applyRules[D[MIs,x],symbolicRules];

(*solve the differential equations*)
{toCanonicalBasis,canonicalBasisSolution,neededBCs}=ParallelTry[solveDE,{DE,#}&/@{Canonica,epsilon,fuchsia,INITIAL,Libra,askChristoph}];

(*compute the boundary conditions (via method of regions)*)
BCsSolution=computeBCs[neededBCs];

(*plug into integrand: reconstruct coefficients of special functions*)
FiniteFieldRationalReconstruct[
(Collect[#1,getSpecialFunctions[canonicalBasisSolution],Factor[Series[#,{ep,0,0}]]&]&)@*
(plugInBCsFF[#1,#2,BCsSolution]&)@*
(replaceFF[#1,#2,canonicalBasisSolution]&)@*
(transformFF[#1,#2,toCanonicalBasis]&)@*
(applyRulesFF[#1,#2,symbolicRules]&)
,integrand]

Collect[
Get["/home/gigi/Documents/pmeftdata/GR/3PM/Dp_res_pot.m"]+Get["/home/gigi/Documents/pmeftdata/GR/3PM/Dp_res_rad1.m"],
{_[lI[1]],_m,_Log},Factor]

Take home message

If you ever think about doing such computations, remember:

Data representation matters, design carefully

Abstract the problem:

Someone has probably already solved it!

Represent minimal information

Integral reduction requires fundamentally different approach