BZOJ4767 两双手（组合数学+容斥原理）

　　因为保证了两向量不共线，平面内任何一个向量都被这两个向量唯一表示。问题变为一张有障碍点的网格图由左上走到右下的方案数。

　　到达终点所需步数显然是平方级别的，没法直接递推。注意到障碍点数量很少，那么考虑容斥，即用总方案数减去经过障碍点的方案数。对每个障碍点计算其作为第一个经过的障碍点的方案数即可。

#include<iostream> 
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;
#define ll long long
#define N 510
#define P 1000000007
char getc(){char c=getchar();while ((c<'A'||c>'Z')&&(c<'a'||c>'z')&&(c<'0'||c>'9')) c=getchar();return c;}
int gcd(int n,int m){return m==0?n:gcd(m,n%m);}
int read()
{
    int x=0,f=1;char c=getchar();
    while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}
    while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return x*f;
}
int t,n,m,cnt,fac[N*N<<2],inv[N*N<<2],f[N],ans;
struct data
{
    int x,y;
    void get(){x=read(),y=read();}
    int operator *(const data&a) const
    {
        return x*a.y-y*a.x;
    }
    bool operator <(const data&a) const
    {
        return x+y<a.x+a.y;
    }
}e,a,b,ban[N],v[N];
int C(int n,int m){return 1ll*fac[n]*inv[m]%P*inv[n-m]%P;}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("bzoj4767.in","r",stdin);
    freopen("bzoj4767.out","w",stdout);
    const char LL[]="%I64d\n";
#else
    const char LL[]="%lld\n";
#endif
    e.get();t=read();a.get();b.get();
    for (int i=1;i<=t;i++) ban[i].get();
    if ((e*b)%(a*b)||(e*a)%(b*a)) {cout<<0;return 0;}
    n=(e*b)/(a*b),m=(e*a)/(b*a);
    for (int i=1;i<=t;i++)
    if ((ban[i]*b)%(a*b)||(ban[i]*a)%(b*a));
    else
    {
        int x=(ban[i]*b)/(a*b),y=(ban[i]*a)/(b*a);
        if (x>=0&&x<=n&&y>=0&&y<=m) v[++cnt]=(data){x,y};
    }
    sort(v+1,v+cnt+1);
    fac[0]=1;for (int i=1;i<=n+m;i++) fac[i]=1ll*fac[i-1]*i%P;
    inv[0]=inv[1]=1;for (int i=2;i<=n+m;i++) inv[i]=P-1ll*(P/i)*inv[P%i]%P;
    for (int i=2;i<=n+m;i++) inv[i]=1ll*inv[i]*inv[i-1]%P;
    ans=C(n+m,n);
    for (int i=1;i<=cnt;i++)
    {
        f[i]=1ll*C(v[i].x+v[i].y,v[i].x);
        for (int j=1;j<i;j++)
        if (v[i].x>=v[j].x&&v[i].y>=v[j].y)
        f[i]=(f[i]-1ll*f[j]*C(v[i].x-v[j].x+v[i].y-v[j].y,v[i].x-v[j].x)%P+P)%P;
        ans=(ans-1ll*f[i]*C(n-v[i].x+m-v[i].y,n-v[i].x)%P+P)%P;
    }
    cout<<ans;
    return 0;
}