physics - Boson commutation relations

I have tried to adapt this answer to my problem of calculating some bosonic commutation relations, but there are still some issues.

The way I'm implementing the commutator is straightforward:

commutator[x_, y_] :=x**y-y**x

Example: if I want to compute $[b^\dagger b,ab^\dagger]$ I write

commutator[BosonC["b"]**BosonA["b"], BosonA["a"]**BosonC["b"]]

and the output is $ab^\dagger$ as it should be. However this fails when I compute $[a^\dagger a,ab^\dagger]$ (which should be $-ab^\dagger)$:

commutator[BosonC["a"]**BosonA["a"], BosonA["a"]**BosonC["b"]]
 Out: a†** a^2 ** b†- a ** b†** a†** a

How can I modify the code in this answer to have it work properly?

EDIT Building on the answers of @QuantumDot and @evanb, I came up with this solution. First I implement the commutator, with Distribute.

NCM[x___] = NonCommutativeMultiply[x];
SetAttributes[commutator, HoldAll]

NCM[] := 1
commutator[NCM[x___], NCM[y___]] := Distribute[NCM[x, y] - NCM[y, x]]
commutator[x_, y_] := Distribute[NCM[x, y] - NCM[y, x]]

Then I implement two tools, one for teaching Mathematica how to swap creation and annihilation operators and one is for operator ordering:

dag[x_] := ToExpression[ToString[x] ~~ "†"]

mode[x_] := With[{x†= dag[x]},
NCM[left___, x, x†, right___] := NCM[left, x†, x, right] + NCM[left, right]]

leftright[L_, R_] := With[{R† = dag[R], L† = dag[L]}, 
NCM[left___, pr : R | R†, pl : L | L†, right___] := NCM[left, pl, pr, right]]

Now I can use it like this: after evaluating the definitions I input (for instance)

mode[a]
mode[b]
leftright[a,b]

And finally I can evaluate commutators, for instance

commutator[NCM[a†,a] + NCM[b†,b], NCM[a,b†]]
(* 0 *)

Answer

The function NonCommutativeMultiply has too long of a name, so I make a short-hand version of it (NCP stands for non-commutative product):

NCP[x___] := NonCommutativeMultiply[x];

Now, here's the code

NCP[] := 1
NCP[left___, a, a†, right___] := NCP[left, a†, a, right] + NCP[left, right]
NCP[left___, b, b†, right___] := NCP[left, b†, b, right] + NCP[left, right]
NCP[left___, pl : a | a†, pr : b | b†, right___] := NCP[left, pr, pl, right]

Now your function:

SetAttributes[commutator, HoldAll]    
commutator[NCP[x___], NCP[y___]] := NCP[x, y] - NCP[y, x]

Let's give it a try:

commutator[NCP[b†, b], NCP[a, b†]]

b† ** a