What is the cleverest way to calculate $[\hat{a}^{M},\hat{a}^{\dagger N}]$ when $[\hat{a},\hat{a}^{\dagger}]=1$? What are some elegant ways to calculate
$$[\hat{a}^{M},\hat{a}^{\dagger N}]\qquad\text{with} \qquad[\hat{a},\hat{a}^{\dagger}]=1,$$
other than brute force calculation?

(EDIT) I got the same result as Qmechanic. I think Prathyush's suggestion should be equivalent to mine of the correspondence up to a canonical transformation. Here is my calculation:
$\begin{array}{c}
\mbox{representation of }\left(\hat{a},\hat{a}^{\dagger}\right)\mbox{ on polynomial space }span\left\{ \frac{x^{n}}{\sqrt{n!}}\right\} _{n\ge0}\\
\hat{a}\left[f\left(x\right)\right]=\frac{d}{dx}f\left(x\right)\;;\;\hat{a}^{\dagger}\left[f\left(x\right)\right]=xf\left(x\right)\;;\;\left[\hat{a},\hat{a}^{\dagger}\right]\left[f\left(x\right)\right]=id\left[f\left(x\right)\right]\\
\left|0\right\rangle \sim 1\;;\;\left|n\right\rangle \sim x^{n}/\sqrt{n!}
\end{array}$
$\begin{array}{c}
\mbox{calculate the normal ordering }\left[\hat{a}^{M},\hat{a}^{\dagger}{}^{N}\right]\mbox{:}\\
\sim\left[\frac{d^{M}}{dx^{M}},x^{N}\right]=\frac{d^{M}}{dx^{M}}\left(x^{N}\star\right)-x^{N}\frac{d^{M}}{dx^{M}}\left(\star\right)\\
\sim\left\{ \overset{min\left\{ M,N\right\} }{\underset{k=0}{\sum}}\frac{N!}{\left(N-k\right)!}C_{M}^{k}\left(\hat{a}^{\dagger}\right)^{N-k}\left(\hat{a}\right)^{M-k}\right\} -\left(\hat{a}^{\dagger}\right)^{N}\left(\hat{a}\right)^{M}\\
\end{array}$

One comment on 02-12-2012:
The representation I was using is actually related to Bergmann representation with the inner product for Hilbert space (polynomials) being:
$$\left\langle f\left(x\right),g\left(x\right)\right\rangle :=\int dxe^{-x^{2}}\overline{f\left(x\right)}g\left(x\right)\,,x\in\mathbb{R}\,,\, f,g\in\mathbb{C}\left[x\right]$$
 A: You can compute this via purely combinatorial 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 terms obtained in the process.
The nontrivial 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 of functions to take an expression, systematically apply the commutation rules, and graph the resulting "tree 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

A: @YunlongLian points out an interesting way to do this in comments, stating that the commutator $[\partial_x, x]=1$ naturally satisfies this if we ignore that $x$ here is not self-adjoint, which turns out to be an irrelevant detail.
From there, we look at $[\partial_x^n, x^m]$ and the first option is that its derivative gets taken once: the derivative is $mx^{m-1}$ and there are $n$ derivative operators that can do it. Then we want two derivatives, then three: so if in commuting the derivative is taken $k$ times, then the sum over all $k$ must be,
$$
[\partial_x^n, x^m]=\sum_{k=1}^{\min(m,n)}
{n\choose k}\frac{m!}{(m-k)!} x^{m-k} \partial_x^{n-k},
$$where the construction $m!/(m-k)! = m (m-1)\dots(m-k+1)$ properly gets the leading term of the derivative.
When I found this it seemed arrestingly visually straightforward and I felt it deserves to be posted as a full-fledged answer.
A: 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.  


*

*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$.

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

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

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

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

*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}.
$$

*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}.
$$
