https://www.51nod.com/onlineJudge/questionCode.html#!problemId=1303

所求=总矩阵对数-相离矩阵对数-相包含矩阵对数

总矩阵对数=${矩阵数 \choose 2}$

考虑如何求有多少个全1矩阵。

先枚举右下角的端点$(i,j)$，维护以$(i,j)$为右下角的全1矩阵个数

枚举行$i$，从左到右枚举列$j$

用单调栈维护$k=1..j$到$j$的$up$的$min$ ($up[i][j]$表示点$(i,j)$向上有连续多少个1)

以$(i,j)$为右下角的矩阵个数=$min_k$的和，可以维护。

考虑如何求相包含矩阵对数。

枚举较大的矩阵($n*m$)，则较小矩阵有${n+1 \choose 2} * {m+1 \choose 2}$个

所以用类似上面的方法，维护$1*f(m),n*f(m),n^2*f(m)$的和即可。($f(m)为{m+1 \choose 2} 的前缀和$)

考虑如何求相离矩阵对数。

即左右相离+上下相离-左上右下相离-右上左下相离

求出每个点作为四个方向的端点的矩阵个数，之后就轻松解决了。

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
#define pb push_back
#define rep(i,l,r) for(int i=l;i<=r;++i)
#define per(i,r,l) for(int i=r;i>=l;--i)

const int N=1e3+5,D=1e9+7;
ll mi(ll x,int y=D-2)
{
    ll ans=1;
    while(y)
    {
        if(y&1)ans=ans*x%D;
        x=x*x%D;y>>=1;
    }
    return ans;
}
int n,m;
char b[N][N];
int up[N][N],down[N][N];
ll S1[N],S2[N],f[N],f2[N];
//S1[n]=sigma i; S2[n]=sigma i^2; f[n]=sigma S[i]
ll a[N],h[N];
ll cnt[N][N],cnt1[N][N],cnt2[N][N],cnt3[N][N];
//cnt 右下角 ;cnt1 左下角 ;cnt2 右上角 ;cnt3 左上角 
ll c[N],c1[N];

int main()
{
    freopen("1.in","r",stdin);//freopen("2.out","w",stdout);
    cin>>n>>m;
    rep(i,1,n)scanf("%s",b[i]+1);

    int mn=max(n,m);
    rep(i,1,mn){S1[i]=S1[i-1]+i;f[i]=f[i-1]+S1[i];}
    rep(i,1,mn){S2[i]=S2[i-1]+i*i;(f2[i]=f2[i-1]+S2[i])%D;}

    rep(i,1,n)
    rep(j,1,m) 
    if(b[i][j]=='1')up[i][j]=up[i-1][j]+1;
    per(i,n,1)
    rep(j,1,m) 
    if(b[i][j]=='1')down[i][j]=down[i+1][j]+1;

    ll ans=0;
    //+C2(tot_cnt) 
    //-相包含
    //-横相离-纵相离+左上右下相离+右上左下相离

    //-相包含 
    rep(i,1,n)
    {
        int top=0,ncnt=0;
        //ncnt sigma h_k
        ll s0=0,s1=0,s2=0;
        //s0=sigma f[h_k]
        //s1=sigma (j-k+1)*f[h_k]
        //s2=sigma (j-k+1)^2*f[h_k]
        rep(j,1,m)
        {
            (s2+=2*s1+s0)%=D;
            (s1+=s0)%=D;
            while(h[top]>up[i][j])
            {
                ncnt-=h[top]*(a[top]-a[top-1]);
                (s0-=f[h[top]]*(a[top]-a[top-1]))%=D;
                (s1-=f[h[top]]*(S1[j-a[top-1]]-S1[j-a[top]]))%=D;
                (s2-=f[h[top]]*(S2[j-a[top-1]]-S2[j-a[top]]))%=D;
                --top;
            }
            h[++top]=up[i][j];
            a[top]=j;
            ncnt+=h[top]*(a[top]-a[top-1]);
            (s0+=f[h[top]]*(a[top]-a[top-1]))%=D;
            (s1+=f[h[top]]*(S1[j-a[top-1]]-S1[j-a[top]]))%=D;
            (s2+=f[h[top]]*(S2[j-a[top-1]]-S2[j-a[top]]))%=D;

            cnt[i][j]=ncnt;
            (ans+=s2+s1)%=D;
        }
    }
    (ans*=mi(2))%=D;
    ans=-ans;

        //get cnt2
    rep(i,1,n)
    {
        int top=0,ncnt=0;
        rep(j,1,m) 
        {
            while(h[top]>down[i][j])
            {
                ncnt-=h[top]*(a[top]-a[top-1]);
                --top;
            }
            h[++top]=down[i][j];
            a[top]=j;
            ncnt+=h[top]*(a[top]-a[top-1]);
            cnt2[i][j]=ncnt;
        }
    }

    a[0]=m+1;
    //get cnt1
    rep(i,1,n)
    {
        int top=0,ncnt=0;
        per(j,m,1) 
        {
            while(h[top]>up[i][j])
            {
                ncnt-=h[top]*(a[top-1]-a[top]);
                --top;
            }
            h[++top]=up[i][j];
            a[top]=j;
            ncnt+=h[top]*(a[top-1]-a[top]);
            cnt1[i][j]=ncnt;
        }
    }
    //get cnt3
    rep(i,1,n)
    {
        int top=0,ncnt=0;
        per(j,m,1)
        {
            while(h[top]>down[i][j])
            {
                ncnt-=h[top]*(a[top-1]-a[top]);
                --top;
            }
            h[++top]=down[i][j];
            a[top]=j;
            ncnt+=h[top]*(a[top-1]-a[top]);
            cnt3[i][j]=ncnt;
        }
    }

    //-横相离 
    rep(i,1,n)
    rep(j,1,m) c[j]+=cnt[i][j];
    rep(i,1,n)
    rep(j,1,m) c1[j]+=cnt1[i][j];
    ll sum=0;
    rep(i,1,m)sum+=c[i];
    sum%=D;
//-----  +C2(tot_cnt) 
    (ans+=(sum+1)*sum/2)%=D;
//-----    
    per(i,m,1)
    {
        (sum-=c[i])%=D;
        (ans-=sum*c1[i])%=D;
    }

    //-纵相离 
memset(c,0,sizeof(c));memset(c1,0,sizeof(c1));
    rep(i,1,n)
    rep(j,1,m) c[i]+=cnt[i][j];
    rep(i,1,n)
    rep(j,1,m) c1[i]+=cnt2[i][j];
    sum=0;
    rep(i,1,n)sum+=c[i];
    sum%=D;
    per(i,n,1)
    {
        (sum-=c[i])%=D;
        (ans-=sum*c1[i])%=D;
    }

    //+左上右下相离
    rep(i,1,n)
    rep(j,1,m)cnt[i][j]+=cnt[i][j-1];
    rep(i,1,n)
    rep(j,1,m)(cnt[i][j]+=cnt[i-1][j])%=D;
    rep(i,1,n)
    rep(j,1,m)(ans+=cnt3[i][j]*cnt[i-1][j-1])%=D;

    //+左下右上相离 
    rep(i,1,n)
    per(j,m,1)cnt1[i][j]+=cnt1[i][j+1];
    rep(i,1,n)
    rep(j,1,m)(cnt1[i][j]+=cnt1[i-1][j])%=D;
    rep(i,1,n)
    rep(j,1,m)(ans+=cnt2[i][j]*cnt1[i-1][j+1])%=D;

    cout<<(ans%D+D)%D;
}