P13647 [NOISG 2016] Fabric 题解

· · 题解

下文中所说的 “高” 指所在行号更小，“低” 则相反。

首先考虑 K=1 的情况。考虑枚举矩形的右上端点，将行倒序枚举，不难发现此时的矩形会被每一列此时最高的洞给限制，且起到限制作用的洞高度一定单调递减。于是可以在扫描列的时候用单调栈更新起到限制作用的洞，计算答案不难手模一下得出。

很难推广到 K>1，主要原因还是在于一个去重的问题，即：难以找到一个能使计算答案不重不漏的式子。但是，K=1 的做法中，维护每列最高的洞以刻画限制应该是正确的。

注意到，每一列最高的洞的高度合起来组成了一个序列。若选定一个区间，提供限制的是这个区间中最高的洞。此时不难想到用单调栈维护出每个洞作为最高的洞的极长区间。把这个区间称为该点的管辖区间，那么对每个点只考虑当前枚举的行是最高行、在它的管辖区间内、且跨过它的矩形。这样做去重问题得到解决。

考虑计算答案。设 g(n,m) 表示在高度为 n、宽度为 m 的方格内，最高行是第一行、面积大于等于 K 的矩形数量。则对于上面的情况，设当前行列是 i,j，当前列最高的洞的高度是 b_j，当前列管辖区间是 L,R，该点的答案就是

g(b_j-i,R-L+1)-g(b_j-i,i-L)-g(b_j-i,R-i)

我们现在要在 O(NM) 的时间内递推算出所有的 g(n,m)。不好计算，于是先计算一个弱的问题 f(n,m) 表示在高度为 n、宽度为 m 的方格内面积大于等于 K 的矩形数量。显然 g(n,m)=f(n,m)-f(n-1,m)。这个可以推导：

\begin{aligned} f(n,m)&=\sum_{i=1}^n\sum_{j=1}^m[ij\ge K](n-i+1)(m-j+1)\\ &=\sum_{i=1}^n(n-i+1)\sum_{j=1}^m\bigg[j\ge\bigg\lceil\frac Ki\bigg\rceil\bigg](m-j+1) \end{aligned}

引入一个 h(n,m)=\sum\limits_{j=n}^m(m-j+1)。推导出递推公式：

\begin{aligned} h(n,m)&=\sum_{j=n}^m(m-j+1)\\ &=\sum_{j=n}^{m-1}(m-j+1)+1\\ &=\sum_{j=n}^{m-1}(m-j)+m-n+1\\ &=h(n,m-1)+m-n+1 \end{aligned}

那么有

\begin{aligned} f(n,m)&=\sum_{i=1}^n(n-i+1)h(\lceil K/i\rceil,m)\\ &=\sum_{i=1}^{n-1}(n-i+1)h(\lceil K/i\rceil,m)+h(\lceil K/n\rceil,m)\\ &=\sum_{i=1}^{n-1}(n-i)h(\lceil K/i\rceil,m)+\sum_{i=1}^{n-1}h(\lceil K/i\rceil,m)+h(\lceil K/n\rceil,m)\\ &=f(n-1,m)+\sum_{i=1}^nh(\lceil K/i\rceil,m) \end{aligned}

然后你会发现 f 不需要，可以直接由 h 推出 g。综上，递推预处理出 g,h，由上文所述统计答案即可。

::::success[码]

#include<bits/stdc++.h>
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<20,stdin),p1==p2)?EOF:*p1++)
#define il inline
using namespace std;

char buf[1<<20],*p1=buf,*p2=buf;
il int read(){
    int x=0;
    char ch=getchar();
    while(!isdigit(ch))ch=getchar();
    while(isdigit(ch))x=x*10+(ch^48),ch=getchar();
    return x;
}
using ll=long long;
constexpr int MAXN=2005;
int n,m,K,a[MAXN][MAXN],b[MAXN];
int stk[MAXN],L[MAXN],R[MAXN],top;
ll g[MAXN][MAXN],h[MAXN][MAXN],ans;
il void init(){
    for(int i=1;i<=m;i++)
        for(int j=i;j<=m;j++)
            h[i][j]=h[i][j-1]+j-i+1;
    for(int j=1;j<=m;j++){
        ll sm=0;
        for(int i=1;i<=n;i++){
            sm+=h[min((K+i-1)/i,m+1)][j];
            g[i][j]=sm;
        }
    }
}

int main(){
    n=read(),m=read(),K=read();
    init();
    for(int i=1;i<=n;i++)
        for(int j=1;j<=m;j++)
            a[i][j]=read();
    fill(b+1,b+m+1,n+1);
    for(int i=n;i;i--){
        for(int j=1;j<=m;j++)
            if(a[i][j])
                b[j]=i;
        stk[top=0]=0;
        for(int j=1;j<=m;j++){
            while(top&&b[stk[top]]>=b[j]) top--;
            L[j]=stk[top]+1;
            stk[++top]=j;
        }
        stk[top=0]=m+1;
        for(int j=m;j;j--){
            while(top&&b[stk[top]]>b[j]) top--;
            R[j]=stk[top]-1;
            stk[++top]=j;
        }
        for(int j=1;j<=m;j++)
            if(b[j]!=i)
                ans+=g[b[j]-i][R[j]-L[j]+1]-g[b[j]-i][j-L[j]]-g[b[j]-i][R[j]-j];
    }
    printf("%lld\n",ans);
    return 0;
}

::::

再介绍另一种比较类人的 O(N^{2.5}) 解法，来自机房同行，经过原作者同意在此记录。如果打算阅读下面内容，请扶稳坐好：

实际上，这个做法在常数优秀、块长合理的情况下，可以稳稳通过。

背后的故事：机房某骗分战神用正难则反+珂朵莉树在比赛时 A 了这道题，但是无法通过 luogu 数据。

感觉 S\ge k 的矩形面积没有 S<k 的矩形面积好算，考虑正难则反，计算矩形数量后转化为后面的问题。

计算矩形数量可以用单调栈简单计算，这里重点讲 S<k 的矩形计数。

考虑枚举右下角 (x,y)，然后计算一个数组 a_i，表示以当前点为右下角的高度为 i 的矩形中，面积 <k 的有几个。再设 b_i=\lfloor\frac{k-1}{i}\rfloor，则显然有 a_i\le b_i。同时，在当前列 y，设当前位置向上遇到第一个洞的横坐标为 x'，则 i\in[1,x-x'] 的 a_i 加一，i\in(x-x',n] 的 a_i 归零。

所以，我们现在需要一种 DS，可以实现前缀加一，后缀归零，以及在每一次操作后，数组整体对 b_i 取 \min。然后查询全局和。

发现这个东西是处理困难的，考虑大力分块。

由于枚举顺序是一行一行枚举的，所以 a_i 最大为 n，因此空间上允许我们对于每个块内维护 c_{i,j}，表示块 i 内 b_k' 值为 j 的数量，初始 b_i'=b_i。这样我们就可以通过 c 计算出每次给 \mathit{ans} 减去的值的增量。

由于有覆盖操作，可以只保留整块 +1，而不保留散块 +1，因此对每个整块维护 \text{tag}_i 表示这个块 +1 的次数。

覆盖时分散块归零和整块归零。散块归零暴力重构整个 c,b' 即可，方法是设修改第 K 块，对于一个位置 i，假如它没有被归零，则 b_i'\to b_i'-\text{tag}_K；否则 b_i'\to b_i。

整块归零时，考虑给每个块设 \text{vs}_i 表示有没有进行过散块归零操作。如果有，暴力重构；没有，则将 \text{tag} 归零即可。时间复杂度考虑均摊分析，共有 O(N) 次散块归零，因此整块归零时暴力重构次数也为 O(N)。

设块长为 B，则单次重构时间复杂度可以做到 O(B)，因此时间复杂度为 O(N^2(\frac{N}{B}+B))。由于暴力重构常数较大，所以当 B=20 左右时，代码常数最小。

::::success[码（原作者）]

#include<bits/stdc++.h>
using namespace std;
char buf[1<<20],*p1,*p2;
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<20,stdin),p1==p2)?0:*p1++)
template<typename T>
inline void read(T &x){
    bool f=1;x=0;char ch=getchar();
    while(ch<'0'||ch>'9') ch=getchar();
    while(ch>='0'&&ch<='9')
        x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
    x=(f?x:-x);
}
template<typename T>
inline void write(T x){
    if(x>9) write(x/10);
    putchar(x%10+'0');
}
const int N=2005,bl=20,M=105;
int n,m,k,pr[N],nx[N],ls[N],st[N],tp;
int b[N],num[N],tag[M],sm[M],c[M][N],gx[M],vs[M],snm[M],mx,sum;
long long ans;
inline void ret(int x,int cc=0){
    if(!vs[x]&&!cc){
        sum-=gx[x];
        tag[x]=gx[x]=0;
        sm[x]=c[x][0];
        return;
    }
    int rs=x*bl,ls=rs-bl+1;
    rs=min(rs,mx),sum-=gx[x];
    gx[x]=tag[x]=vs[x]=0;
    while(ls<=rs){
        c[x][num[ls]]--;
        c[x][num[ls]=b[ls]]++,ls++;
    }
    if(cc&&rs==mx) c[x][0]++;
    sm[x]=c[x][0];
}
inline void chg(int x){
    for(int i=1;i*bl-bl<x;i++)
        sum+=snm[i]-sm[i],gx[i]+=snm[i]-sm[i],sm[i]+=c[i][++tag[i]];
    int bid=x/bl+1,rs=min(bid*bl,mx);
    sum-=gx[bid],gx[bid]=0;
    for(int i=bid*bl-bl+1;i<=x;i++){
        c[bid][num[i]]--,num[i]=max(num[i]-tag[bid],0);
        c[bid][num[i]]++,gx[bid]+=b[i]-num[i];
    }
    for(int i=rs;i>x;i--) c[bid][num[i]]--,c[bid][num[i]=b[i]]++;
    tag[bid]=0,sum+=gx[bid];
    sm[bid]=c[bid][0],vs[bid]=1;
    for(int i=bid+1;i*bl-bl<mx;i++) ret(i);
}
int main(){
    read(n),read(m),read(k);
    for(int i=1;i<=n;i++,tp=0){
        for(int j=1,cc;j<=m;j++){
            read(cc);
            (cc?pr[j]=0:pr[j]++);
            while(tp&&pr[st[tp]]>pr[j]) nx[st[tp--]]=j;
            ls[j]=st[tp],nx[j]=m+1,st[++tp]=j;
        }
        mx=i;
        for(int j=1;j<=m;j++) ans+=1ll*pr[j]*(j-ls[j])*(nx[j]-j);
        for(int j=1;j<=i;j++) b[i]=min((k-1)/j,2001); 
        for(int j=1;j*bl-bl<i;j++) ret(j,1),snm[j]=min(j*bl,mx)-j*bl+bl;
        for(int j=1;j<=m;j++) chg(pr[j]),ans-=sum;
    }
    write(ans);
    return 0;
}

::::