Expand a product of operators

Lets suppose I have a family of operators $f_i$, and unknown c-numbers $a,b,c$.
I want to expand such products: $(f_1+f_2+c)(f_3+f_2+b)$ into $b c+b f_2+b f_1+c f_2+c f_3+f_2^2+f_1 f_2+f_2 f_3+f_1 f_3$.
Equivalently, the action of the product on $x$:

 (f1+f2+c)@(f3+f2+b)[x]


should turn into

b*c*x + b*f2[x]+b*f1[x]+c*f2[x]+c*f3[x]+f2@f2[x]+f1@f2[x]+f2@f3[x]+f1@f3[x]


The problem I run into is making Mathematica assume that the symbols $a,b,c$ are numbers and not operators.

EDIT:

1) The operators $f_i$ should act linearly: f[a x]= a f[x].
2) Using NonCommutativeMultiply I have been able to get an expansion of the form:

b*c + b*f2+b*f1+c*f2+c*f3+f2**f2+f1**f2+f2**f3+f1**f3


After using a rule to replace ** with Composition I get:

 b*c + b*f2+b*f1+c*f2+c*f3+Composition[f2,f2]+Composition[f1,f2]+Composition[f2,f3]+Composition[f1,f3]


The problem is that the action of $(3a + 2b \, f1\circ f2)(x)$ on $x$

(3*a + 2*b*Composition[f1,f2])[x] //Through


gives

(3 a)[x] + (2 b Composition[f1, f2])[x]


I have tried multiplying the scalars by Identity:

(3*a Identity + 2*b*Composition[f1,f2])[x] //Through


But Mathematica still treats (3*a*Identity)[x] as a function.
What I want to get is this:

3*a*Identity[x] + 2*b*Composition[f1,f2][x]

• Composition is not commutative; I would have expected f1[b x] in the result instead of what you have. Jul 3, 2015 at 20:03
• You can't use _ for subscripts. Instead, use f[1][x], etc. Maybe you can use this?
– Jens
Jul 3, 2015 at 21:29
• @Guesswhoitis. I should have specified that the operators should be linear, thus f[a x]=a f[x]. Jul 4, 2015 at 14:17
• @Jens It was bad notation, fixed now. The solution proposed there doesn't work for me for reasons outlined in the edited post. Jul 4, 2015 at 14:48

First of all, the desired result Composition[f1, f2][x] auto-evaluates to something slightly different:

Composition[f1, f2][x]

(* ==> f1[f2[x]] *)


If that's OK with you, then I could simply re-use my answer to Having the derivative be an operator, without the first two lines in which I define the specific action of the operators as derivatives. So all you need is this:

multiplyOp[scalar_, op_] := Function[{f1}, scalar op[f1]];

addOps[ops__] := Function[{f1}, Total@Map[#[f1] &, {ops}]];

CircleDot[scalar_, op_] := multiplyOp[scalar, op];
CircleTimes[ops__] := Composition[ops];


The last line is optional, since there is also the symbol @* for Composition. But I prefer to have single operator symbols. The important thing is that the algebra becomes easier to program if we clearly distinguish between the operations that are involved. The addition of operators is not the same as Plus, so I define it as a special symbol; the same is true for multiplications. The shortcuts are:

• CirclePlus ⊕ typed as escc+esc
• CircleDot ⊙ typed as escc.esc
• CircleTimes ⊗ typed as escc*esc

With this, the example in the question can be entered as follows:

((3 a + 2 b)⊙f1⊗f2)[x]

(* ==> (3 a + 2 b) f1[f2[x]] *)