# Multiplying two expressions with a specific function, efficiently

Paraphrasing, the problem I am trying to solve is the following. I am given two "matrix valued" sums

expr1 = Sum[c[i] M[i], {i, range1}]
expr2 = Sum[d[j] M[j], {j, range2}]


where the cs and ds are some coefficients and the Ms are some unspecified (symbolic) matrices. I would like to multiply these two sums, with the multiplication of two such matrices understood by me, i.e. specified by some function, let's call it "f". I am trying to do this efficiently. I do not know beforehand which matrices will appear, only that they have head M.

Defining an explicitly linear matrix multiplication works, but seems to be quite inefficient, at least in the way that I managed to come up with (my coefficients are further generic expressions).

Using Collect and Pattern matching instead, I managed to come up with the following

desiredfunction[a_, b_] :=
Module[{c = Collect[a, M[_]], d = Collect[b, M[_]]},
c[] = List; d[] = List;
Outer[#1[[1 ;; -2]]*#2[[1 ;; -2]]*f[#1[[-1]], #2[[-1]]] &, c, d]//Flatten//Tr]


i.e. I explicitly turn the sums into lists of terms, and generate the outer product of the two lists, multiplying the coefficients and combine the matrices using "f", and finally adding everything together.

This does what I want it to, but I cannot shake the feeling that this is such an elementary kind of operation (combining two expressions written in the form coefficients*objects with a rule to combine the objects and coefficients) that there must be a better way to do this, maybe via in-built functions. Am I overlooking something obvious here?

For copy-paste into mathematica:

expr1 = Sum[c[i] M[i], {i, 1, 3}];
expr2 = Sum[d[i] M[i], {i, 2, 6}];

desiredfunction[a_, b_] :=
Module[{c = Collect[a, M[_]], d = Collect[b, M[_]]},
c[] = List;
d[] = List;
Outer[#1[[1 ;; -2]]*#2[[1 ;; -2]]*f[#1[[-1]], #2[[-1]]] &, c, d] //
Flatten // Tr]

desiredfunction[expr1, expr2]


ClearAll[desiredF]

  True