Reformulating the Condition

Let the left-hand side of the given equality be \(L\coloneqq \min_{1\le i\le H}\max_{1\le j\le W}A_{i,j}\), and the right-hand side be \(R\coloneqq \max_{1\le j\le W}\min_{1\le i\le H}A_{i,j}\).

Choose an \(i\) attaining the \(\min\) on the left-hand side, and a \(j\) attaining the \(\max\) on the right-hand side. Then \(L\ge A_{i,j}\ge R\) holds. Thus, we always have \(L\ge R\), so the condition \(L=R\) can be restated as \(L\le R\).

For the \(i,j\) chosen above, the condition \(L=R\) implies

\[\max_{1\le q\le W}A_{i,q}=A_{i,j}=\min_{1\le p\le H}A_{p,j}\quad \cdots(\ast)\]

Conversely, the existence of such \(i,j\) is necessary and sufficient.

An Expression Using Inclusion–Exclusion

Consider counting the matrices for which no such \(i,j\) exist and subtracting this count from the total. We then expand the condition that \((\ast)\) fails for every \(i,j\) using inclusion–exclusion. Given several specified pairs \((i,j)\), we want to count the matrices \(A\) for which all these pairs satisfy \((\ast)\). We can verify that if \((i,j)\) and \((p,q)\) satisfy \((\ast)\), then \((i,q)\) and \((p,j)\) also satisfy \((\ast)\). Thus, what matters is not the set of specified pairs \((i,j)\) itself, but the pair consisting of the set of \(i\) values and the set of \(j\) values that occur in it, and ultimately the numbers of distinct \(i,j\) values.

Suppose that the numbers of distinct \(i,j\) values are \(h,w\ge 1\), respectively. The sum of the inclusion–exclusion coefficients can then be expressed as

\[ \sum_{a=0}^{h}\sum_{b=0}^{w}(-1)^a(-1)^b\binom{h}{a}\binom{w}{b}\sum_{k=0}^{(h-a)(w-b)}\binom{(h-a)(w-b)}{k}(-1)^k \]

Here, another application of inclusion–exclusion is used to enforce the condition that the numbers of distinct \(i,j\) values are \(h,w\), respectively. The coefficient from the original inclusion–exclusion is the final term \((-1)^k\). The \(3\)rd \(\Sigma\) is \(1\) only when \(a=h\) or \(b=w\), and is \(0\) otherwise. Note that \(0^0=1\). Consequently, this sum is \(-(-1)^{h+w}\). Now recall that we decided to subtract this count, including the case \(h=w=0\), from the total. Thus, as a contribution to the answer, the sum of the inclusion–exclusion coefficients is \((-1)^{h+w}\).

Fix an integer \(m\in [1,M]\) and consider the conditions on \(A\) for which \(L=R=m\). Suppose that \(h\) rows and \(w\) columns are specified as the sets of \(i\) and \(j\) satisfying \((\ast)\). At their intersections, we have \(A_{i,j}=m\). At positions belonging only to a specified row \(i\) or only to a specified column \(j\), the respective conditions are \(A_{i,q}\le m, A_{p,j}\ge m\). At positions belonging to neither a specified row nor a specified column, \(A_{p,q}\) may take any value. Therefore, the answer is

\[ \sum_{h=1}^{H}\sum_{w=1}^{W}\sum_{m=1}^{M}\binom{H}{h}\binom{W}{w}(-1)^{h+w}m^{h(W-w)}(M+1-m)^{(H-h)w}M^{(H-h)(W-w)} \]

We will compute this value efficiently.

Computation Using Lagrange Interpolation

By treating the occurrences of \(M\) inside the \(\sum\) as fixed constants and only the \(M\) at the upper limit of the \(\sum\) as a variable, we can regard the answer as a polynomial in \(M\). (Treating the occurrences of \(M\) inside the sum as variables as well would still give a polynomial, but would be inconvenient for computation.) Its degree is at most \(HW\). We therefore seek to determine this polynomial (implicitly), using Lagrange interpolation.

In the three nested sums \(\Sigma\) in the expression for the answer, move \(\sum_{m}\) to the outside and define the inner part as

\[ f(m)=\sum_{h=1}^{H}\sum_{w=1}^{W}\binom{H}{h}\binom{W}{w}(-1)^{h+w}m^{h(W-w)}(M+1-m)^{(H-h)w}M^{(H-h)(W-w)} \]

Using the binomial theorem, we can rewrite this as

\[ \begin{aligned} f(m)&=\sum_{h=1}^{H}\binom{H}{h}(-1)^h\sum_{w=1}^{W}\binom{W}{w}(-(M+1-m)^{H-h})^w(m^hM^{H-h})^{W-w}\\ &=\sum_{h=1}^{H}\binom{H}{h}(-1)^h\Big[(-(M+1-m)^{H-h}+m^hM^{H-h})^W-(m^hM^{H-h})^W\Big] \end{aligned} \]

Thus, for a given \(m\), we can compute \(f(m)\) in, for example, \(O(H\log(HW))\) time.

The answer is \(F(M)=\sum_{m=1}^{M}f(m)\), which we regard as a polynomial in \(M\). Once we have \(F(0),F(1),\dots, F(HW)\), we can compute \(F(M)\) using Lagrange interpolation. Computing these values ultimately reduces to computing \(f(1),f(2),\dots, f(HW)\). (Note that this argument would not work if \(f\) depended on \(M\) here.) Lagrange interpolation itself can be performed in time linear in the degree of the polynomial, after precomputing suitable factorial tables.

This gives an algorithm running in, for example, \(O(H^2W\log(HW))\) time.