Consider the following defined commutation relations:

$$[\hat a,\hat a^{\dagger}]=1$$ 

$$[\hat b,\hat b^{\dagger}]=1$$

$$[\hat a,\hat b]=0$$

(where the usual algebra of commutators holds)

Let us now consider the operator:

$$\hat H=\hat a^{\dagger}\hat a+\hat b^{\dagger}\hat b+\hat a^{\dagger}\hat a(b^{\dagger}+b)+(a^{\dagger}+a)(b^{\dagger}+b)$$

I now want to evaluate the following using Mathematica:
$$[\hat A,\hat H]+[\hat A,[\hat A,\hat H]]+[\hat A,[\hat A,[\hat A,\hat H]]].....$$ 

where,
$\hat A=(\hat a^{\dagger}-\hat a+\hat b^{\dagger}-\hat b)$

As can be seen, this is a really difficult task to perform by hand and I want to use Mathematica to do it. Since I am new at Mathematica, I don't even know how to begin.

Maybe I would have to somehow define the above relations as rules but I don't then know how to make Mathematica follow the commutator algebra rules.

I am not looking for a complete solution (though it's okay if it comes as an answer since I can always use it to match with my code:), just a guide on how to start this problem in Mathematica.

**EDIT:**
From the link in the comment by evanb I realized that it's rather hard to do this in Mathematica. Has anybody used the SNEG package who can guide me on how to install or use it in Mathematica?

I am hoping that that should make my work easier. Since explicitly writing the code might take a very long time. I even have to change $\hat A$ in the above expression and try for many such $\hat A$'s.