P2216 [HAOI2007] 理想的正方形 题解

user100566 · 2024-12-03 14:04:01 · 题解

前言

这道题是一道很好的深入理解 ST 算法的例题，通过优化，ST 算法成功地取得了这道题的总 222ms，最大 32ms,8.04MB 的最优解，与次优解拉开断层差距。

本题解讲解如何从基础的 ST 算法开始优化，在这之前你需要掌握基本的 ST 算法。

题目分析

题目描述很直接，考虑用 ST 算法实现，初始化 ST 表后，暴力枚举正方形区域（共 O(ab) 个），然后用 ST 查出最值并更新答案即可。主要需要考虑空间大小，初始化时间和单次查询时间。

基础的 ST 做法

二维 ST 表的定义

用 A 表示原数组，首先我们来看一维 ST 算法定义的数组 B[i][p](1\le i\le n, 0\le p\le k)，最小的 k 是多少？基础二维 ST 表空间较大，有必要分析需要分配的长度。
ST 的原理是用两个可重叠的区间合并为要查询的区间，只要两个最长区间长度的和不小于最大查询长度即可，不妨设最大查询长度为 m，那么 k 满足 2\times 2^k\ge m\Rightarrow k=\lceil\log_2m\rceil -1，注意 p 从 0 开始，数组长度需设为 k+1。

接下来考虑二维 ST 表，参考一维下 B[i][p] 表示从 A[i] 开始的长度为 2^p 的线段的最值，我们可以定义 B[i][j][p][q] 表示以 A[i][j] 为左上角，大小为 2^p\times 2^q 的矩形的最值。

因为查询长度 n\le100，块长需达到 2^6=64，则 p, q 维的长度定义为 7。

求解 ST 表

求解这个 ST 表比较简单，设现在是在求最大值 ST 表，参考一维 ST 表的倍增转移方程，得到下面的初始化方法：

注意，并不需要同时通过两个方程转移，实现时可以先从小到大枚举 p，然后从小到大枚举 q，若 q=0 则通过第一个方程转移，否则通过第二个方程转移。

处理查询

接下来考虑查询操作，设要查询的矩形为 A[i_l\sim i_r][j_l\sim j_r]，先在 i 上划分出两个长度为二的幂次的块，例如 [10, 20] 划分为 [10, 17] 和 [13, 20]，然后在 j 上继续划分即可，总计需要查询 4 个块。下面的图展示了如何查询 A[10\sim 20][3\sim 8]：

可以看出查询被划分为 4 个 8\times4 的矩形 [10\sim17][3\sim6], [10\sim17][5\sim8], [13\sim20][3\sim6], [13\sim20][5\sim8]，在 4 个子最值中取出最值即为查询结果。

复杂度分析

ST 表的大小即空间复杂度为 O(ab(\log n)^2)，初始化复杂度 O(ab(\log n)^2)，查询复杂度 O(ab)（忽略求子块长的 O(\log n) 复杂度，这可以通过初始化规避），总时间复杂度 O(ab(\log n)^2)。

然而我们计算程序需要的空间大小，最大最小 ST 表总计需要内存 1000\times1000\times7\times7\times4\text{B}\times2=3.92\times10^8\text{B}\approx373.84\text{MB}\gg128\text{MB}，MLE，必须优化。

正方形优化

ST 表压维

注意到查询的矩形始终是正方形的，所以查询时的子块也是正方形的，这意味着我们没有必要定义两个 p, q 分别表示两个方向上的长度，只需要一个 p 表示正方形的大小就可以了。
定义新的压维 ST 表 B[i][j][p] 表示以 A[i][j] 为左上角，边长为 2^p 的正方形的最值。

压维 ST 表的求解

一个大正方形可以由 4 个小正方形组合，倍增转移方程如下：

处理查询

同基础 ST 表的查询，由于查询的是正方形，一定能找到 4 个子正方形，注意如果不保证查询为正方形，这一步无法进行，例如上面举例的 A[10\sim 20][3\sim 8] 就无法划分为 4 个子正方形。

复杂度分析

容易知道时空复杂度均为 O(ab\log n)，现在的最大最小 ST 表总计需要内存 1000\times1000\times7\times4\text{B}\times2=5.6\times10^7\text{B}\approx53.4\text{MB}\ll128\text{MB}，已经可以通过本题。

自我滚动优化

ST 算法的本质是倍增 DP，既然是 DP，那么就可以使用滚动数组优化空间复杂度，倍增长度这一维是可以滚掉的，不过为什么普通的 ST 算法不使用滚动优化呢？
原因显然：对于不同的查询长度，需要不同长度的块来合并为查询区间，例如查询 [1\sim 9] 必须用长度为 8 的子块合并，查询 [1\sim 63] 必须用长度为 32 的子块合并，这两种不同长度的子块在这种情况下无法相互替代，因此必须存储倍增过程中的中间计算结果用于查询。
然而本题询问的正方形边长 n 为定值！例如 n=100 时，只需要用到 64\times64 的子正方形的最值，即查询时只用到了 B[\ ][\ ][6]，这是我们使用滚动数组优化的依据。

滚动 ST 表

我们把 p 也压维，定义滚动 ST 表 B[i][j] 表示外层循环 p 意义下，以 A[i][j] 为左上角，边长为 2^p 的正方形的最值。

滚动 ST 表的求解

同压维 ST 表的求解，需要注意的是外层 p 的枚举次数需要随着 n 变动，如果锁定为一个大值，那么无法处理 n 较小的查询。
同时滚动时注意 i, j 的取值顺序。

处理查询

查询正方形的大小是固定的，对于给定的正方形左上角坐标，要查询的 4 个子正方形的左上角坐标是可以预处理出 \text{offset} 的，详见代码。

复杂度分析

容易知道这一步优化后空间复杂度变为 O(ab)，最好可以只定义两个 a\times b 的 ST 表，总内存仅为存储输入数据所需内存的两倍！

递推优化

设想如果要求查询的矩形不是正方形，那么无法使用正方形优化，但是依然可以使用自我滚动优化：先滚动 p，再滚动 q，同样可以求出滚动 ST 表，观察两种方法的转移方程：

1. B[i][j]=\max\{B[i][j], B[i][j+2^{p-1}], B[i+2^{p-1}][j], B[i+2^{p-1}][j+2^{p-1}]\}
2. B[i][j]=\max\{B[i][j], B[i+2^{p-1}][j]\}, B[i][j]=\max\{B[i][j], B[i][j+2^{q-1}]\}

将 1. 中的 \max 拆成二元 \max，会拆出 3 个 \max，相比之下 2. 只有 2 个，后者优于前者。
实际上并不一定真的先滚动 p，再滚动 q，可以轮流滚动，由于查询是正方形，滚动会同时结束，可以在一个循环内同时完成 p, q 的滚动。

当然，查询过程中也可以使用这样的优化。

图形解释

观察下面一张格点有向图：

注意，这张图省略了主对角线外的斜有向边，例如 (2,1)\rightarrow(4,2)\rightarrow(8,4)\rightarrow(16,8) 就被省略，事实上可以通过这些隐藏的边转移。

其中，每个节点都可视为 ST 表的状态（仅含 i, j 维），点上的 (x, y) 表示对应状态下的矩形的长为 x 宽为 y（x=2^p, y=2^q）。

查询时，需要根据查询长度在 i, j 上找到合适长度，例如 A[10\sim 20][3\sim 8] 需要查找图中的节点 (8,4)，A[1\sim 31][1\sim 31] 需要查找图中的节点 (16, 16) 。

回顾基础 ST 表 B[i][j][p][q]，我们发现这个表包括了整张图的全部节点，根据上面给出的思路（优先通过第二个方程转移）可以得到它的转移图：

这张图共 O(\log_2n\times\log_2n) 个节点，O(\log_2n\times\log_2n) 条边，据此可知时空复杂度为 O(ab\times (\log_2n)^2)。

对于正方形优化中的压维 ST 表 B[i][j][p]，我们发现它只包括了对角线上的节点，通过向右下的有向边转移，其转移图为：

对角线上的节点数共 \log_2n 个，据此可知时空复杂度为 O(ab\times \log_2n)。

滚动 ST 表利用了询问只用到了右下角的节点这一特性，只存储一个节点的信息，然后通过自我滚动转移到新的节点上。

边上的数值表示了转移所需的时间，前面已经提到了，先向下再向右由于直接向右下，递推优化的转移图如下：

上面演示的是先滚动 p 再滚动 q 的转移图，我们发现路径的长度减少了。

由于询问是正方形的，到达最终节点向下和向右的边数相等，所以可以同步滚动，转移图如下：

这仅仅是为了将转移写在一个循环内，两者没有性能上的明显差距。

请读者借助以以上图片思考以下问题：

1. 如果需要找出的是 n\times m 的矩形，如何实现？
2. 如果正方形的边长可以是给定的 k 个 \le n 的值（假设题目为此缩小了数据范围或者延长了时间限制），能否用 O(ab) 的空间复杂度实现？
3. 如果找出的矩形的大小可以是给定的 k 个 (n, m)，并且满足 n_i\le n_{i+1}, m_i\le m_{i+1}，能否用 O(ab) 的空间复杂度实现？
4. 跳出本题，考虑二维数组 RMQ，每次查询一个 n\times m 的矩阵，其中 n 是一个在所有询问前给出的定值，m 与查询有关，那么最好的时空复杂度分别为多少？如果强制在线呢？

代码

关于取值细节，代码中的注释有详细解释，其它优化如快读，预处理等请参见代码。

云剪贴板存档

#include <cstdio>
#include <algorithm>
using namespace std;
// 快读 
#ifdef ONLINE_JUDGE
#define getchar() getchar_unlocked()
#endif
inline int read(){
    int c = getchar();
    while(c<48||c>57) c = getchar();
    int x = 0;
    while(48<=c&&c<=57){
        x = (x<<3)+(x<<1)+(c^48);
        c = getchar();
    }
    return x;
}

int a, b, n;
// dp1 维护最大值 dp2 维护最小值 
int dp1[1000][1000], dp2[1000][1000];

int main(){
    a = read(); b = read(); n = read();
    // 即使去掉这个特判也是正确的 
    if(n==1){
        putchar('0');
        return 0;
    }
    for(int i=0; i<a; ++i){
    for(int j=0; j<b; ++j){
        // 不需要存原矩阵，直接作为 dp 初值 
        dp2[i][j] = dp1[i][j] = read();
    }
    }

    int maxT = 0;
    while((2<<maxT)<n) ++maxT;
    int offset, maxi, maxj;
    for(int t=0; t<maxT; ++t){
        offset = 1<<t;
        // 在 i 方向上倍增 
        // i 方向上区间长是 2*offset，区间为 [i, i+2*offset-1]
        // i+2*offset-1<a => i<a-2*offset+1 => i<=a-2*offset
        maxi = a-(offset<<1);
        // j 方向上区间长是 offset，区间为 [j, j+offset-1]
        // j+offset-1<b => j<b-offset+1 => j<=b-offset
        maxj = b-offset;
        for(int i=0; i<=maxi; ++i){
        for(int j=0; j<=maxj; ++j){
            // std::max, std::min 快得不止一点 
            //if(dp1[i+offset][j]>dp1[i][j]) dp1[i][j] = dp1[i+offset][j];
            //if(dp2[i+offset][j]<dp2[i][j]) dp2[i][j] = dp2[i+offset][j];
            dp1[i][j] = max(dp1[i][j], dp1[i+offset][j]);
            dp2[i][j] = min(dp2[i][j], dp2[i+offset][j]);
        }
        }
        // 在 j 方向上倍增 
        // maxi = a-(offset<<1);
        maxj = b-(offset<<1);
        for(int i=0; i<=maxi; ++i){
        for(int j=0; j<=maxj; ++j){
            dp1[i][j] = max(dp1[i][j], dp1[i][j+offset]);
            dp2[i][j] = min(dp2[i][j], dp2[i][j+offset]);
        }
        }
    }

    offset = n-(1<<maxT);
    // 在 i 方向上倍增 
    // 区间长是 n，区间为 [i, i+n-1]
    // i+n-1<a => i<a-n+1 => i<=a-n
    maxi = a-n;
    // 区间长是 2^maxT，区间为 [j, j+2^maxT-1]
    // j+2^maxT-1<b => j<b-2^maxT+1 => j<=b-2^maxT
    maxj = b-(1<<maxT);
    for(int i=0; i<=maxi; ++i){
    for(int j=0; j<=maxj; ++j){
        dp1[i][j] = max(dp1[i][j], dp1[i+offset][j]);
        dp2[i][j] = min(dp2[i][j], dp2[i+offset][j]);
    }
    }
    // 在 j 方向上倍增 
    // maxi = a-n;
    maxj = b-n;
    int ans = 1000000000;
    for(int i=0; i<=maxi; ++i){
    for(int j=0; j<=maxj; ++j){
        // 直接更新 ans，不用更新 dp 
        ans = min(ans, max(dp1[i][j], dp1[i][j+offset])-min(dp2[i][j], dp2[i][j+offset]));
    }
    }
    printf("%d", ans);
    return 0;
}