矩阵颜色

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;
}