数据范围：(n<=1000,m<=1000,|x_i|,|y_i|<=1000)，数据保证(a,b>=1,a+b<n)
Solution
（这种操作好像之前已经遇到过一次了但是并没有很好地掌握。。）
首先比较直接的一个想法当然是容斥，如果说我们有一种方法可以在一个比较合理的复杂度内计算出某个点集在(m)步之后走到同一个位置的方案数（记为(f(A))，其中(A)是一个点集），那么我们可以用简单的三元容斥得到答案：
所以现在的问题就是怎么算这个(f)
我们将所有的点都放在切比雪夫坐标系中（将原坐标系顺时针旋转45°并且单位长度缩小至原来的(frac{1}{sqrt 2})），直白一点就是原来在笛卡尔坐标系中的点((x,y))在这个新的坐标系中对应的点是((x+y,x-y))，这个时候我们走一步的操作从原来的在某一维上走(1)的长度变成了斜着走，反应到新坐标系中的坐标变化就是((xpm1,ypm1))，这就意味着走一步对(x)和(y)的影响是独立的，而不是像原来笛卡尔坐标系中的那样是((xpm 1,y))或者((x,ypm 1))，所以我们就可以将两维独立开来算了，最后的结果就是两维分别计算后的乘积
那么现在问题就由二维变成了一维：给你若干个点，求这些点移动(m)位最后走到同一个点的方案数
枚举最后走到的点(t)，那么我们要做的就是计算某个点(x)在(m)步移动之后到达(t)的方案数：首先(|t-x|)和(m)在(\% 2)意义下必须同余（或者说差必须是偶数），因为我们需要往(t)的方向走(|t-x|)步，剩下的(m-|t-x|)步需要两两抵消掉才能保证(m)步之后停在(t)这个位置，所以总的方案数就是在(m)步里面选(frac{m-|t-x|}{2})步往(t)的反方向走，也就是(inom m {frac{m-|t-x|}{2}})
注意因为新坐标系中的坐标是((x+y,x-y))，所以转化为一维问题后最后走到的位置应该(in[-3000,3000])，所以一次(f)计算的复杂度是(O(6000n))
代码大概长这个样子
#include<iostream>
#include<cstdio>
#include<cstring>
#include<vector>
#define vct vector<Dot>
using namespace std;
const int N=1010,MOD=1e9+7,V=3000;//just for debuging!!!
struct Dot{
int x,y;
void read(){
int x1,y1;
scanf("%d%d",&x1,&y1);
x=x1+y1; y=x1-y1;
}
}a[N];
int C[N][N],f[3];
int split[N];
int n,m,ans;
int plu(int x,int y){return (1LL*x+y)%MOD;}
int mul(int x,int y){return 1LL*x*y%MOD;}
int Abs(int x){return x<0?-x:x;}
vct operator + (const vct &x,const vct &y){
vct ret; ret.resize(x.size()+y.size());
copy(y.begin(),y.end(),copy(x.begin(),x.end(),ret.begin()));
return ret;
}
void prework(int n){
C[0][0]=1;
for (int i=1;i<n;++i){
C[i][0]=1; C[i][i]=1;
for (int j=1;j<i;++j)
C[i][j]=plu(C[i-1][j],C[i-1][j-1]);
}
}
int binom(int n,int m){return n<m||m<0?0:C[n][m];}
int calc(int x){
if ((m-x)&1) return 0;
return binom(m,(m-x)/2);
}
int dp(const vct &a){
int retx=0,rety=0,tmpx,tmpy,sz=a.size();
for (int goal=-V;goal<=V;++goal){
tmpx=tmpy=1;
for (int i=0;i<sz&&(tmpx||tmpy);++i){
tmpx=mul(tmpx,calc(Abs(goal-a[i].x)));
tmpy=mul(tmpy,calc(Abs(goal-a[i].y)));
}
retx=plu(retx,tmpx);
rety=plu(rety,tmpy);
}
return mul(retx,rety);
}
void solve(){
vct tmp[2];
int cnt,same;
for (int i=0;i<3;++i){
f[i]=dp(vct(a+split[i],a+split[i+1]));
cnt=0;
for (int j=0;j<3;++j)
if (i!=j)
tmp[cnt++]=vct(a+split[j],a+split[j+1]);
same=dp(tmp[0]+tmp[1]);
ans=plu(ans,MOD-mul(f[i],same));
}
ans=plu(ans,mul(f[0],mul(f[1],f[2])));
same=dp(vct(a+1,a+1+n));
ans=plu(ans,mul(2,same));
printf("%d
",ans);
}
int main(){
#ifndef ONLINE_JUDGE
freopen("a.in","r",stdin);
#endif
int x,y;
scanf("%d%d%d%d",&n,&x,&y,&m);
split[0]=1; split[1]=x+1; split[2]=split[1]+y; split[3]=n+1;
prework(N);
for (int i=1;i<=n;++i)
a[i].read();
solve();
}