What is the cleverest way to calculate $[\hat{a}^{M},\hat{a}^{\dagger N}]$ when $[\hat{a},\hat{a}^{\dagger}]=1$?

The standard way is to use generating functions (in this case a la coherent states). Usually one would like the resulting formula to be normal-ordered.

1. Recall the following version $$\tag{1} e^Ae^B~=~e^{[A,B]}e^Be^A$$ of the Baker-Campbell-Hausdorff formula. The formula (1) holds if the commutator $[A,B]$ commutes with both the operators $A$ and $B$.

2. Put $A=\alpha a $ and $B=\beta a^{\dagger}$, where $\alpha,\beta\in\mathbb{C}$.

3. Let $[a, a^{\dagger}]=\hbar {\bf 1}$, so that the commutator $[A,B]=\alpha\beta\hbar {\bf 1}$ is a $c$-number.

4. Now Taylor-expand the exponential factors in eq. (1).

5. For fixed orders $n,m\in \mathbb{N}_0$, consider terms in eq. (1) proportional to $\alpha^n\beta^m$.

6. Deduce that the the antinormal-ordered operator $a^n(a^{\dagger})^m$ can be normal-ordered as $$\tag{2} a^n(a^{\dagger})^m~=~\sum_{k=0}^{\min(n,m)} \frac{n!m!\hbar^k}{(n-k)!(m-k)! k!}(a^{\dagger})^{m-k}a^{n-k}. $$

7. Finally, deduce that the normal-ordered commutator is $$\tag{3} [a^n,(a^{\dagger})^m]~=~\sum_{k=1}^{\min(n,m)} \frac{n!m!\hbar^k}{(n-k)!(m-k)! k!}(a^{\dagger})^{m-k}a^{n-k}. $$

You can compute this via purely combinatoric arguments.

The trick is to observe that $aa^\dagger=a^\dagger a +1$ means that you can think of the process of going from $a^{m}a^{\dagger n}$ to $a^{\dagger n}a^m$ as one composed of a number of steps, in each one of which a choice is made as to whether have $a$ "act" on $a^\dagger$, thus having them both disappear, or to have them switch places. You thus build a "tree of possibilities" whose leaves are the normally ordered obtained in the process. The hard part is counting the number of resulting factors of each type.

The possible factors in the final expression are $a^{\dagger (n-k)}a^{(m-k)}$ for $k=0,...,\min(n,m)$. This is easy to see: every time you apply the commutation relation you either maintain all operators or annihilate one pair of them. The final expression must therefore have the form $$a^m a^{\dagger n}= \sum_{k=0}^{\min(n,m)}c_k a^{\dagger (n-k)}a^{(m-k)}$$ for some integer coefficients $c_k$. These coefficients are derived by counting the number of ways in which a given term can be obtained. To do this, we simply need to find the number of ways in which we can pair subsets of size $k$ of $a^\dagger$ and $a$.

For example, $c_0=1$ follows easily from the fact that there is only one way to not pair anything with anything else. At the same time, $c_1=nm$, because each of the $m$ annihilation operators can interact once with each of the $n$ creation ones. More generally, we have $$c_k = k! \binom{n}{k}\binom{m}{k}.$$ To see it, notice that the number of ways to single out a subset of $k$ operators $a$ and $k$ operators $a^\dagger$ (ignoring ordering) is $\binom{n}{k}\binom{m}{k}$. Then, the number of ways in which these $k$ operators can all annihilate each others is $k!$: the first $a$ can interact with any of the $k$ terms on the right, then the second $a$ can interact with any of the remaining $k-1$ terms, etc.

To obtain the final form of the commutator we just remove the $k=0$ term from the RHS: $$[a^m,a^{\dagger n}]=\sum_{k=1}^{\min(n,m)}k!\binom{n}{k}\binom{m}{k} a^{\dagger(n-k)}a^{(m-k)}.$$

---

To better visualise what's going with these calculations, I wrote a couple functions to take an expression, systematically apply the commutation rules, and graph the resulting "three of possibilities". Here is for example the result for $a^2 a^{\dagger 2}$ (I use $b\equiv a^\dagger$ in the graph for better readability):

Here is the Mathematica code generating this figure (it requires having MaTeX installed):

Get["MaTeX`"]
distribute[args_] := args //. 
    {HoldPattern[nc[l___, Plus[m__], r___]] :> 
      Total[(nc[l, #1, r] & ) /@ {m}], 
     nc[l___, (c_)*nc[m__], r___] :> c*nc[l, m, r], 
     nc[l___, nc[m__], r___] :> nc[l, m, r], nc[-(a_), b_] :> 
      -nc[a, b], nc[a_, -(b_)] :> -nc[a, b], nc[nc[l__], r_] :> 
      nc[l, r], nc[l_, nc[r__]] :> nc[l, r], nc[a_] :> a, nc[] -> 1}; 
singleStepExpand[expr_] := (Map[distribute])[
    expr /. {nc[l___, a, b, r___] :> {nc[l, b, a, r], nc[l, r]}}]; 
stepExpand[expr_] := (If[Length[#1] > 1, stepExpand /@ #1, #1] & )[
    singleStepExpand[expr]]; 
stepExpandFullStory[expr_] := 
   (If[Head[#1] === List, Append[{expr}, stepExpandFullStory /@ #1], 
      #1] & )[singleStepExpand[expr]]; 
firstIfList[expr_] := If[Head[expr] === List, First[expr], expr]; 
nestedListToListOfEdges[expr_] := Cases[expr, 
    {l:nc[__], {first_, second_}} :> Sequence @@ 
      {DirectedEdge[l, firstIfList[first]], DirectedEdge[l, 
        firstIfList[second]]}, All]; 
groupPowers[args_] := args //. {nc[l___, a_, a_, r___] :> 
      nc[l, a^2, r], nc[l___, (a_)^(n_), a_, r___] :> 
      nc[l, a^(n + 1), r], nc[l___, a_, (a_)^(n_), r___] :> 
      nc[l, a^(n + 1), r]}; 
makeCsToBrackets[expr_] := expr //. 
    {c[a_, b_] :> StringJoin["[", ToString[makeCsToBrackets[a]], ",", 
       ToString[makeCsToBrackets[b]], "]"], 
     (a_)^(n_) :> StringJoin[ToString[makeCsToBrackets[a]], "^", 
       ToString[n]]}; 
beautify[expr_] := makeCsToBrackets[expr] //. 
    {nc[args__] :> MaTeX[StringJoin @@ ToString /@ {args}, 
       Magnification -> 1.5], s_String :> 
      MaTeX[s, Magnification -> 1.5]}; 
edgesToGraphWithNiceLabels[edges_] := 
   Graph[edges, VertexLabels -> 
     (#1 -> beautify[groupPowers[#1]] & ) /@ DeleteDuplicates[
       Flatten[edges, Infinity, DirectedEdge]], 
    GraphLayout -> "LayeredDigraphEmbedding"];

nc[a, a, b, b] // stepExpandFullStory // nestedListToListOfEdges // edgesToGraphWithNiceLabels