对于每个1，预处理出它头上最低的0距离它的相对高度 $h_{i,j}$ 。

然后一排一排的刷，对于每个1，统计出以当前这个1为右下角的矩形的方案数，那么最后的答案就是每个点的方案数之和。

对于每个统计的过程，向它左边找第一个 $h_{i,k}$ 使得 $h_{i,k}$ 严格小于 $h_{i,j}$ ，那么这个点的答案即为

$ans_{i,j}=ans_{i,k}+(j-k)*h_{i,j}$

然而这个找 $h_{i,k}$ 的过程就可以用单调栈维护

That's all

所以也就这些：

inline void solve(){
	reads(n),reads(m);
	for(int i=1;i<=n;i++){
		for(int j=1;j<=m;j++){
			reads(a[i][j]);
			h[i][j]+=a[i][j]*(1+h[i-1][j]);
		}
	}
	for(int i=1;i<=n;i++){
		tot=0;
		for(int j=1;j<=m;j++)sum[j]=0;
		for(int j=0;j<=m;j++){
			x=h[i][j];
			while(tot&&st[tot].y>x)tot--;
			st[++tot].y=x;
			st[tot].x=j;
			if(tot==1)st[tot-1].x=st[tot].x-1;
			if(x)sum[j]=sum[st[tot-1].x]+(st[tot].x-st[tot-1].x)*st[tot].y;
			ans+=sum[j];
	writes(ans);
}

最后附上一个完整code：

#include<bits/stdc++.h>
//#define int long long
//#define int __int128
#define endl "\n"
#define mkp(a,b) make_pair(a,b)
#define pii pair<int,int>
//#pragma GCC optimize(2)
using namespace std;
bool Test_MLE_start;
int __=1;
namespace FastIO{
//	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 reads(T &x){
		char c=getchar();
		int sum=0,f=1;
		while(!isdigit(c)){
			if(c=='-') f=-1;
			c=getchar();
		}
		while(isdigit(c)){
			sum=(sum<<3)+(sum<<1)+(c^48);
			c=getchar();
		}
		x=sum*f;
	}
	template<typename T>inline void reads_f(T &x){
		char c=getchar();
		int f=1,sum1=0,sum2=0,cnt=0;
		double sum3=0.00;
		while(!isdigit(c)){
			if(c=='-') f=-1;
			c=getchar();
		}
		while(isdigit(c)&&c!='.'){
			sum1=(sum1<<3)+(sum1<<1)+(c^48);
			c=getchar();
		}
		c=getchar();
		while(isdigit(c)){
			sum2=(sum2<<3)+(sum2<<1)+(c^48);
			cnt++;
			c=getchar();
		}
		sum3=(sum2*pow(10,-cnt)*1.0+sum1)*f;
		x=sum3;
	}
	template<typename T>inline void reads_bit(T &x){
		char c=getchar();
		int sum=0;
		while(!isdigit(c))c=getchar();
		while(isdigit(c)){
			sum=(sum<<1)+(c^48);
			c=getchar();
		}
		x=sum;
	}
	template<typename T>inline void writes(T x){
		if(x<0) putchar('-'),writes(-x);
		if(x>9) writes(x/10);
		putchar(x%10+48);
	}
}using namespace FastIO;
inline void files(){
	freopen("std.in","r",stdin);
	freopen("std.out","w",stdout);
}
const int N=3005;
int a[N][N];
int h[N][N];
int outs[N][N];
int sum[N];
long long ans;
int x;
struct block{
	int y,x;
}st[N];
int tot;
int n,m;
inline void clr(){
//	Don't forget!
}
inline void solve(){
	reads(n),reads(m);
	for(int i=1;i<=n;i++){
		for(int j=1;j<=m;j++){
			reads(a[i][j]);
			h[i][j]+=a[i][j]*(1+h[i-1][j]);
		}
	}
//	cout<<endl;
//	for(int i=1;i<=n;i++){
//		for(int j=1;j<=m;j++){
//			cout<<h[i][j]<<' ';
//		}
//		cout<<endl;
//	}
	for(int i=1;i<=n;i++){
		tot=0;
		for(int j=1;j<=m;j++)sum[j]=0;
		for(int j=0;j<=m;j++){
			x=h[i][j];
//			if(!x){
//				tot=0;
//				continue;
//			}
			while(tot&&st[tot].y>x)tot--;
			st[++tot].y=x;
			st[tot].x=j;
//			cout<<"("<<i<<","<<j<<"):";
			if(tot==1)st[tot-1].x=st[tot].x-1;
			if(x)sum[j]=sum[st[tot-1].x]+(st[tot].x-st[tot-1].x)*st[tot].y;
			ans+=sum[j];
//			outs[i][j]=sum[j];
//			cout<<endl;
		}
	}
//	for(int i=1;i<=n;i++){
//		for(int j=1;j<=m;j++){
//			cout<<outs[i][j]<<" ";
//		}
//		cout<<endl;
//	}
	writes(ans);
}
bool Test_MLE_end;
signed main(){
//	printf("%lf Mb\n",(&Test_MLE_end-&Test_MLE_start-1)/1024.0/1024.0);
//	files(); 
//	reads(__);
	cin.tie(0);cout.tie(0);
	ios::sync_with_stdio(false);
	while(__--) clr(),solve();
	return 0;
}
/*
6 11
1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 0 0 
1 1 0 1 1 1 1 0 1 1 1
1 0 1 0 0 1 1 1 1 1 1
0 1 1 1 1 1 0 1 1 1 1
1 1 1 1 1 0 1 1 1 1 1

363
*/