Factor out prefixes when expanding dummy indices

[eschnett]

Factor out prefixes that do not contain the dummy index when expanding dummy indices. This reduces code size slightly.

For example, when expanding

a b[ui] c[li]

then we want to generate

a (b1 c1 + b2 c2 + b3 c3)

instead of distributing "a" over all terms. This reduces code size when "a" is a large expression.

[eschnett]

Here is the patch:

$ git diff
diff --git a/Tools/CodeGen/TensorTools.m b/Tools/CodeGen/TensorTools.m
index 50d108c..65630ab 100644
--- a/Tools/CodeGen/TensorTools.m
+++ b/Tools/CodeGen/TensorTools.m
@@ -565,19 +565,52 @@ makeSum[x_] :=
   ThrowError["Expression " <> ToString[x] <> 
     " is not recognized, and tensor indices will not be expanded"];

-sumComponentsOfDummyIndex[x_, i_] := 
+sumComponentsOfDummyIndex1[x_, i_] := 
   Apply[Plus,listComponents[x, i, toggleIndex[i]]];

+sumComponentsOfDummyIndex[x_, i_] := 
+  Module[{contains, indexName, hasIndex, termsWithoutIndex, termsWithIndex},
+         (* We can't use MatchQ because we want to look into all
+            subexpressions as well *)
+         contains[expr_, form_] := Count[expr, form, Infinity] != 0;
+         indexName[TensorIndex[n_,_]] := n;
+         hasIndex[j_, y_] := contains[y, TensorIndex[indexName[j],_]];
+         termsWithoutIndex = 1;
+         termsWithIndex = x;
+         (* Do not try to split a term if it has only a single factor *)
+         If[Head[x]==Times || Head[x]==TensorProduct,
+            termsWithoutIndex = Select[x, !hasIndex[i,#]&];
+            termsWithIndex = Select[x, hasIndex[i,#]&]];
+         (termsWithoutIndex *
+          Apply[Plus, listComponents[termsWithIndex, i, toggleIndex[i]]])];
+
 (* Given an expression and a list of lower indices which are dummies
    in that expression, expand the dummy indices in the expression into
    components. *)
-makeSumOverDummies[x_] := 
+makeSumOverDummies1[x_] := 
   Module[{is},
     is = dummiesIn[x];
     If[is == {}, 
       x, 
       makeSumOverDummies[sumComponentsOfDummyIndex[x, First[is]]]]];

+makeSumOverDummies[x_] :=
+  Module[{count, indexName, countIndex, is},
+         count[expr_, form_] := Count[expr, form, Infinity];
+         indexName[TensorIndex[n_,_]] := n;
+         countIndex[j_, y_] := count[y, TensorIndex[indexName[j],_]];
+         (* find all dummies *)
+         is = dummiesIn[x];
+         (* count how often they occur *)
+         is = Map[{#, countIndex[#, x]}&, is];
+         (* sort by this count *)
+         is = SortBy[is, Last];
+         Print["msod x ",x];
+         Print["msod is ",is];
+         (* drop counts *)
+         is = Map[First, is];
+         Fold[sumComponentsOfDummyIndex, x, is]];
+
 (* -------------------------------------------------------------------------- 
    Rules for converting expressions into different forms
    -------------------------------------------------------------------------- *