Loop integrals with odd number of loop momenta

[jruizvid]

Dear experts,

Mathematica is taking too long to handle relatively easy tasks. I'm not sure if it is related to my poorly efficient code (and I have a chance to improve it) or to Mathematica/Feyncalc itself.

After evaluating the trace and doing some expansion in my (tow-loop) amplitude I have ~150 terms. One simplification I want to do (before integrating) is removing the terms with an odd number of loop momenta in the numerator. I made a custom function called removeOddQ, but it takes ~40s for this amplitude.

To my surprise, the time it takes is not proportional to the number of terms. For instance, applying it to 100 and 50 terms,

Number of terms: 100 -> 72; time: 10.2s
Number of terms: 50 -> 33; time: 1.4s

Do you know why this happens?

I attach a simplified notebook showing this with 100 terms of my amplitude. (I had to put it in .txt to upload it.) removeOddK.txt

Another question: Is there a built-in function in FeynCalc that does this? I couldn't get FCMultiLoopTID to only remove the odd number of loop momenta. I would also appreciate if you see something clearly inefficient in my code. Here only the part with removeOddK

uncontractK[exp_] :=  exp // (Uncontract[#, k1, k2, Pair -> All] &) // (FCCanonicalizeDummyIndices[#, LorentzIndexNames -> {in1, in2, in3, in4}] &)

k1ork2[k_] := k === k1 || k === k2;

countks[exp_] :=  Count[exp, Pair[LorentzIndex[_, ___], Momentum[k_?k1ork2, ___]],   Infinity]

removeOddK[exp_] :=  FCLoopIsolate[uncontractK[exp], {k1, k2}, Head -> loopInt] /. loopInt[x_] :> If[OddQ[countks[x]], 0, x] // Contract // Simplify

[magnunm]

The only question here which really has to do with FeynCalc is whether there is built in functionality to do this, so I will address that first.
I don't think it has, and I don't think it needs it, as this is already covered quite conveniently by Mathematica itself.

If I understand you correctly you want to remove a term if and only if there is an odd number of the momenta k1, k2 in the numerator. If so, here is a possible solution: (test is your expression)

replaceAllRestricted[expr_, rules_List, heads_List] := 
  expr /. Join[a : Blank[#] :> a & /@ heads, rules];
replaceAllRestricted[
  test, 
  {Momentum[k1, D] -> a * Momentum[k1, D], Momentum[k2, D] -> a * Momentum[k2, D]},
  {FeynAmpDenominator}] /. a^n_Integer :> If[EvenQ[n], 1, 0]

For me this uses just 15 milliseconds for the 100 term expression.

Explanation:
The idea is to multiply all the momenta k1,k2 with a dummy constant a. Then all terms with an odd power of a has an odd number of momenta and vice versa. In the resulting expression we can replace a to an odd power with zero and a to an even power with one to get what we want. The only complication is that we want to exclude the denominator. This means we wish to replace all momenta except those contained in the head FeynAmpDenominator. That problem already has a good answer on the mathematica stackexchange, which is what I used above.

About your code:
That the time usage is not proportional to the number of terms is no surprise. That would only occur if the complexity of the algorithm was linear. Calling removeOddK a lot of things happen, you have calls then to: Uncontract, FCCanonicalizeDummyIndices, Count, FCLoopIsolate, Contract and Simplify. All of these are non-trivial functions with some cost in terms of time, you can test with AbsoluteTiming on each of the calls to see what takes the longest in this case.

[jruizvid]

Thank you @magnunm for your detailed answer. It works much faster with your code.

Just a small correction: the suggested solution uses the pattern a^n_Integer. This does not match to a^1. To solve it, since the only 'a's left are a^1, I just added the replacement rule a->0 at the end.

Also, for this function to work on D and 4 dimensions, one can use Momentum[k1, d___] -> a*Momentum[k1, d].

About the timing, the bottleneck was in FCCanonicalizeDummyIndices, probably my function could get by without it (with some fix), but I actually like better your solution. Thanks!

[jruizvid]

Some observation about the timing question: The time spent by FCanonicalizeDummyIndices grows as ~n^3 for >10 terms.

However, if you delete it completely (and use automatic indices), then Simplify takes much longer. So, overall, if you want to use Simplify for expressions with uncontracted indices, use FCanonicalizeDummyIndices first

[vsht]

I don't think that one needs FCCanonicalizeDummyIndices here. It's a rather expensive function that should be employed only when absolutely necessary. A quick solution that comes to my mind (one should see if I missed something here) goes along the lines of what was suggested by @magnunm. FCLoopIsolate alone is fairly fast and is not a deal breaker here.

(*Original code*)
out = test2 // removeOddK // AbsoluteTiming;
(*Newl code*)
AbsoluteTiming[
 aux1 = FCLoopIsolate[test2, {k1, k2}, Head -> loopInt];
 aux2 = aux1 /. {loopInt[a_FeynAmpDenominator b_] :> a loopInt[b]};
 loopIntList = Cases2[aux2, loopInt];
 loopIntListEval = 
  loopIntList /. {Momentum[k1, D] -> flag*Momentum[k1, D], 
      Momentum[k2, D] -> flag*Momentum[k2, D]} /. 
    flag^n_ /; EvenQ[n] :> 1 /. 
   loopInt[flag^_. b_] /; FreeQ[b, flag] :> 0;
 loopIntRule = 
  Thread[Rule[loopIntList, loopIntListEval /. loopInt -> Identity]];
 res = aux2 /. Dispatch[loopIntRule];]
(* {0.141442,Null} *)

out[[2]] - res // ExpandScalarProduct // ExpandAll // FDS
(*0*)