网格图上有三个矩形
A:(x_1,y_1)-(x_2,y_2),B:(x_3,y_3)-(x_4,y_4),C:(x_5,y_5)-(x_6,y_6)
求 a \to b \to c(a \in A,b \in B,c \in C) 的方案数。
其中 a \to b 表示 a 向下向右走到 b 的路径。两个方案不同当且仅当 a,b,c 不同，或 a \to b,b \to c 不同。
考虑枚举中间节点 b，需要计算的是一个点到一个矩形的方案数。将矩形差分，只需要计算点 (0,0) 到矩形 (0,0)-(x,y) 的方案数。
\begin{aligned}
\sum_{i=0}^{x} \sum_{j=0}^{y} \binom{i+j}{i} &= \sum_{i=0}^{x} \binom{y+i+1}{i+1} \\
&= \sum_{i=-1}^{x} \binom{y+i+1}{y} - \binom{y}{y} \\
&= \binom{x+y+2}{y+1} - 1
\end{aligned}
差分的四项求和之后 -1 就被消去了，而 \binom{x+y+2}{y+1} 恰为 (0,0) \to (x+1,y+1) 的方案数。记 a \to b 表示 a 到 b 路径的方案数，那么 b 到 A 的方案数为
b \to (x_1-1,y_1-1) + b \to (x_2,y_2) - b \to (x_1-1,y_2) - b \to (x_2,y_1-1)
通过差分相当于将 $A,B$ 分别转化为 $4$ 个点，而枚举 $4 \times 4$ 个点对是可以接受的，考虑枚举起点 $a$ 和终点 $c$，此时对答案贡献的系数为两个点对应系数的乘积。
一条 $a \to c$ 的路径产生的方案数为路径与 $B$ 的交点个数（每个交点都可以充当 $b$）。可以枚举路径进入和离开 $B$ 的点计算方案数，但这样仍是 $O(n^2)$ 的。若矩形包含 $a$ 那就不用枚举进入的点了。继续差分，将 $B$ 差分为四个以 $a$ 为左上角的矩形。枚举路径离开矩形的点计算方案数即可。
时间复杂度 $O(n)$，带一个比较大的常数。
```cpp
#include<bits/stdc++.h>
using namespace std;
#define ar(x) array<int,x>
#define pb push_back
const int N=2e6+100,MOD=1e9+7;
inline int qm(int x,int y=MOD-2){ int res=1; for(;y;y>>=1,x=1ll*x*x%MOD) if(y&1) res=1ll*res*x%MOD; return res; }
inline void mod(int &x){ if(x>=MOD) x-=MOD; if(x<0) x+=MOD; }
#define gc getchar()
#define rd read()
inline int read(){
int x=0,f=0; char c=gc;
for(;c<'0'||c>'9';c=gc) f|=(c=='-');
for(;c>='0'&&c<='9';c=gc) x=(x<<1)+(x<<3)+(c^48);
return f?-x:x;
}
int jc[N],fjc[N];
void init(int n){
jc[0]=1; for(int i=1;i<=n;++i) jc[i]=1ll*jc[i-1]*i%MOD;
fjc[n]=qm(jc[n]); for(int i=n-1;~i;--i) fjc[i]=1ll*fjc[i+1]*(i+1)%MOD;
}
inline int CM(int n,int m){ return n<0||m<0||m>n?0:1ll*jc[n]*fjc[m]%MOD*fjc[n-m]%MOD; }
inline int calc(int x_1,int y_1,int x_2,int y_2){ return CM(abs(x_2-x_1)+abs(y_2-y_1),abs(x_2-x_1)); }
int x[7],y[7],ans=0;
vector<ar(3)> A,B,C;
void sol(ar(3) a,ar(3) b,ar(3) c){
int res=0;
for(int i=a[1];i<=b[1];++i) mod(res+=1ll*calc(b[0],i,a[0],a[1])*calc(b[0]+1,i,c[0],c[1])%MOD*(b[0]-a[0]+i-a[1]+1)%MOD);
for(int i=a[0];i<=b[0];++i) mod(res+=1ll*calc(i,b[1],a[0],a[1])*calc(i,b[1]+1,c[0],c[1])%MOD*(i-a[0]+b[1]-a[1]+1)%MOD);
mod(ans+=a[2]*b[2]*c[2]*res);
}
int main(){
init(N-1);
for(int i=1;i<=6;++i) x[i]=rd;
for(int i=1;i<=6;++i) y[i]=rd;
A.pb({x[2],y[2],1}),A.pb({x[1]-1,y[1]-1,1}),A.pb({x[2],y[1]-1,-1}),A.pb({x[1]-1,y[2],-1});
B.pb({x[4],y[4],1}),B.pb({x[3]-1,y[3]-1,1}),B.pb({x[4],y[3]-1,-1}),B.pb({x[3]-1,y[4],-1});
C.pb({x[6]+1,y[6]+1,1}),C.pb({x[5],y[5],1}),C.pb({x[6]+1,y[5],-1}),C.pb({x[5],y[6]+1,-1});
for(auto a:A) for(auto b:B) for(auto c:C) sol(a,b,c);
printf("%d\n", ans);
return 0;
}
```