矩阵颜色

Problem

一个 \(n\times m\) 的矩阵，第 \(i\) 行第 \(j\) 列元素有一个颜色 \(c_{i,j}\)，求所有子矩阵的颜色种类数的平均值。

\(n, m \le 100, c_{i, j} \le n \times m\)

Input

第一行两个正整数 \(n，m\)。
接下来 \(n\) 行，每行 \(m\) 个正整数，表示元素的颜色。

Output

输出子矩阵的平均颜色数，保留 \(9\) 位小数。

Sample

Input 1

2 3
1 2 1
2 1 2

Output 1

1.666666667

Solution

容易发现，颜色之间互不影响，可以分开考虑每种颜色对答案的贡献。

如下图所示，涂色的格子表示同一颜色，我们将其按照遍历的顺序编号。

为了避免重复计算，我们规定在统计包含编号为 \(i\) 的格子的矩阵数量时，矩阵不能够包含编号为 \(j(j < i)\) 的格子。比如当前我们在计算 \(8\) 号格的贡献，矩阵就不能包含编号在 \([1,7]\) 的格子。这些格子相当于障碍，我们将这些障碍标记出来。

这些障碍规定了矩阵的左右区间范围，且满足单调性，第 \(i\) 行的左右范围一定不大于第 \(i+1\) 行的左右范围。如下图，将左右范围补全。

由于矩阵要包含目标格子，坐标记为 \((x,y)\)，那么矩阵的上界一定在 \([1,x]\) 内。我们从第 \(x\) 行扫描到第 \(1\) 行，一边更新左右范围一边计算。

如下图，计算以第 \(5\) 行为上界时的矩阵范围（阴影部分）。

当扫到第 \(3\) 行时，范围缩小。

当上界为第 \(2\) 行时，范围再次缩小。

当上界为第 \(1\) 行时，发现上面有障碍阻拦，范围缩小至 \(0\)。

那么确定了左右范围，如何统计贡献呢？

先看看一维的情况。

如图，我们要选一个包含 \(x\)点 的区间，区间范围在 \((l, r)\) 之间（\(l,r\) 两点是障碍，不能取到）。左端点范围在 \((l, x]\)，有 \(x - l\) 种取值；右端点范围在 \([x, r)\)，有 \(r - x\) 种取值。符合条件的区间有 \((x-l) \times (r-x)\) 个。

回到矩阵上，一维情况求的是矩阵的长的取值。宽的上界由枚举确定，下界取值在 \([x,n]\) 内，有 \(n - x + 1\) 种取法。

因此，合法的矩阵总共有 \((n-x+1) \times (y-l) \times (r-y)\) 个。

最后，统计完一个格子的贡献后，要更新该格子对这种颜色的左右区间范围的限制。
代码：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <iomanip>

using namespace std;

const int kmax = 105;

int n, m, c[kmax][kmax];
int p[kmax][kmax * kmax];
int l[kmax], r[kmax];
long long res, resc;

void Calc(int x, int y, int c) {
    for (int i = 1; i <= n + 1; i++) {
        l[i] = 0, r[i] = m + 1; // 初始化左右区间
    }
    for (int i = 1; i < y; i++) {
        l[p[i][c]] = max(l[p[i][c]], i); // 计算左区间
    }
    for (int i = y + 1; i <= m; i++) {
        r[p[i][c]] = min(r[p[i][c]], i); // 计算右区间
    }
    for (int i = x; i > p[y][c]; i--) {
        l[i] = max(l[i], l[i + 1]), r[i] = min(r[i], r[i + 1]); // 边计算边更新
        //		cout << i << ' ' << l[i] << ' ' << r[i] << '\n';
        res += 1ll * (n - x + 1) * (y - l[i]) * (r[i] - y);
    }
}

int main() {
    cin >> n >> m;
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= m; j++) {
            cin >> c[i][j];
        }
    }
    resc = n * (n + 1) / 2 * m * (m + 1) / 2; // 总的矩阵数量
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= m; j++) {
            Calc(i, j, c[i][j]); // 计算贡献
            p[j][c[i][j]] = i; // 更新限制
        }
    }
    //	cout << res << '\n';
    cout << fixed << setprecision(9) << 1.0 * res / resc << '\n';
    return 0;
}