公式解説「包除原理による表式」までは同様です。

\(\displaystyle f(h,w)=\binom{H}{h}\binom{W}{w}(-1)^{h+w}M^{(H-h)(W-w)}\) として、

\(\displaystyle\sum_{h=1}^{H}\sum_{w=1}^{W}\sum_{m=1}^{M}f(h,w)m^{h(W-w)}(M+1-m)^{(H-h)w}\) を計算するのが目的です。\(M\) のシグマを外して、これを \(m\) の多項式としたものを陽に求め、最後に \(m\) に \(1,2,...,M\) を代入することを考えます。

\(h,w\) は全探索してしまいます。これを固定したときに影響する項は \((H-h)w\) 項で、二項定理を用いて容易に加算できます。この部分の計算量は \(O(H^2W^2)\) ですが、定数倍に \(\frac{1}{4}\) がかかり高速です。適宜 \(m\) と \((M+1-m)\) を入れ替えて \(\min(h(W-w),(H-h)w)\) 項に加算することでさらに高速になります。

こうして得られた \(m\) の \(O(HW)\) 次式の \(m\) に \(1,2,\dots,M\) を代入して総和を求めたいです。これは、\(\displaystyle \sum_{n=1}^{M}n^k\) を \(k=1,2,...,O(HW)\) について列挙する問題になります。

これはよく知られた問題（リンク先はmaspyさんのブログ）で、\(O(HW\log HW)\) の計算量で求めることができます。