【JZOJ 5270】神奇的矩阵

Description

Solution

把b排序，从小到大插入，绝对值就没用了，
所以就可以容易算出它在答案中乘上的系数，
设f[i][j]表示以i,j作为左上角，这个k*k的矩形里有多少个1（有多少个数已被加入），
那么点x,y作为+出现的次数就是以x-k+1,y-k+1作为左上角的k*k矩形内f的和，

这个用线段树套线段树很容易想到，但常数太大了，
我用的是树状数组（没错，用树状数组区间加区间查询）
想了解这个黑科技，点这里

复杂度：$O(n^2\log(n)^2)$

Code

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fod(i,a,b) for(int i=a;i>=b;i--)
#define min(q,w) ((q)>(w)?(w):(q))
#define max(q,w) ((q)<(w)?(w):(q))
#define NX(q) ((q)&(-(q)))
using namespace std;
const int N=505,mo=1e4+7;
typedef long long LL;
int read(int &n)
{
    char ch=' ';int q=0,w=1;
    for(;(ch!='-')&&((ch<'0')||(ch>'9'));ch=getchar());
    if(ch=='-')w=-1,ch=getchar();
    for(;ch>='0' && ch<='9';ch=getchar())q=q*10+ch-48;n=q*w;return n;
}
int m,n,m1;
LL ans;
struct qqww
{
    int x,y,v;
}a[N*N];
int  f1[N][N],f[N][N];
int f2[N][N],f3[N][N];
int fd,fd1,fd2,fd3;
bool PX(qqww q,qqww w){return q.v<w.v;}
void change(int x,int y,int e)
{
    if(x>n||y>m)return;
    for(int i=x;i<=n;i+=NX(i))
        for(int j=y;j<=m;j+=NX(j))
        {
            f[i][j]+=e;
            f1[i][j]=(f1[i][j]+e*(n-x+1)*(m-y+1))%mo;
            f2[i][j]=(f2[i][j]+e*(n-x+1))%mo;
            f3[i][j]=(f3[i][j]+e*(m-y+1))%mo;
        }
}
void find(int x,int y)
{
    fd1=fd2=fd3=fd=0;
    for(int i=x;i;i-=NX(i))
        for(int j=y;j;j-=NX(j))
        {
            fd+=f[i][j];
            fd1=(fd1+f1[i][j])%mo;
            fd2=(fd2+f2[i][j])%mo;
            fd3=(fd3+f3[i][j])%mo;
        }
}
void add(int x,int y,int x1,int y1)
{
    change(x,y,1);
    change(x,y1+1,-1);
    change(x1+1,y,-1);
    change(x1+1,y1+1,1);
}
int Gans(int x,int y,int x1,int y1)
{
    LL s,s1,s2,s3;
    find(x-1,y1);s=fd,s1=fd1,s2=fd2,s3=fd3;
    find(x-1,y-1);s-=fd,s1-=fd1,s2-=fd2,s3-=fd3;
    LL t=0;
    t=(s3-s*(m-y1))%mo*(x1-x+1)%mo;
    find(x1,y-1);s=fd,s1=fd1,s2=fd2,s3=fd3;
    find(x-1,y-1);s-=fd,s1-=fd1,s2-=fd2,s3-=fd3;
    t=(t+(s2-s*(n-x1))%mo*((y1-y+1))%mo)%mo;
    find(x-1,y-1);s=fd,s1=fd1,s2=fd2,s3=fd3;
    t=(t+s*(y1-y+1)%mo*(x1-x+1))%mo;

    find(x1,y1);s=fd,s1=fd1,s2=fd2,s3=fd3;
    find(x-1,y1);s-=fd,s1-=fd1,s2-=fd2,s3-=fd3;
    find(x1,y-1);s-=fd,s1-=fd1,s2-=fd2,s3-=fd3;
    find(x-1,y-1);s+=fd,s1+=fd1,s2+=fd2,s3+=fd3;
    s1=s1-s3*(n-x1)-s2*(m-y1);s1=(s1%mo+mo)%mo;
    s=s%mo;
    s1=(s1+(n-x1)*(m-y1)%mo*s)%mo;
    return (s1+t)%mo;
}
int main()
{
    int q,w,e;
    read(n),read(m),read(m1);
    fo(i,1,n)fo(j,1,m)read(a[i*m+j-m].v),a[i*m+j-m].x=i,a[i*m+j-m].y=j;
    sort(a+1,a+1+n*m,PX);
    fo(i,1,n*m)
    {
        q=a[i].x,w=a[i].y;
        int x=max(1,q-m1+1),y=max(1,w-m1+1);
        int x1=min(q,n-m1+1),y1=min(w,m-m1+1);
        e=Gans(x,y,q,w);
        a[i].v=a[i].v%mo;
        ans=(ans+(e*2-(x1-x+1)*(y1-y+1)%mo*(m1*m1-1))%mo*a[i].v)%mo;
        add(x,y,x1,y1);
    }
    printf("%lld\n",(ans+mo)*2%mo);
    return 0;
}