这题 WC2018 时 wzj 搬过……
注意要预处理出 $1$ 到 $[(1000^2) imes 2]$ 的所有整数的阶乘。因为 work 函数中返回的 $x$ 和 $y$ 会达到 $1000^2$ 的级别，然后求组合数时还有一项为两个这么大的数相加。
一组卡上界的例子（稍微调整一下就可以使分母不是 $0$ 了，这里就设得极端一点方便理解）：Ax=Bx=0, Ay=By=500, x=500, y=0。
1 #include<bits/stdc++.h> 2 #define int long long 3 #define ll long long 4 #define N 505 5 #define L 1000005 6 #define mod 1000000007 7 using namespace std; 8 inline int read(){ 9 int x=0; bool f=1; char c=getchar(); 10 for(;!isdigit(c); c=getchar()) if(c=='-') f=0; 11 for(; isdigit(c); c=getchar()) x=(x<<3)+(x<<1)+(c^'0'); 12 if(f) return x; 13 return 0-x; 14 } 15 int n; ll jc[L+5],jcn[L+5],f[N]; 16 struct node{int x,y;}E,A,B,b[N]; 17 inline bool cmp(node a, node b){ 18 return a.x==b.x ? a.y<b.y : a.x<b.x; 19 } 20 ll Pow(ll x, int y){ 21 ll ret=1; 22 while(y){ 23 if(y&1) ret=ret*x%mod; 24 x=x*x%mod; 25 y>>=1; 26 } 27 return ret; 28 } 29 void work(int &x, int &y){ 30 int a1 = x*B.y - y*B.x, a2 = A.x*B.y - A.y*B.x; 31 int b1 = x*A.y - y*A.x, b2 = A.y*B.x - A.x*B.y; 32 if(!a2 || !b2 || a1/a2*a2!=a1 || b1/b2*b2!=b1){ 33 x=y=-1; return; 34 } 35 else 36 x=a1/a2, y=b1/b2; 37 } 38 ll C(int n, int m){ 39 return jc[n] * jcn[n-m] % mod * jcn[m] % mod; 40 } 41 ll cal(int x, int y){ 42 if(x<0 || y<0) return 0; 43 return C(x+y, y); 44 } 45 signed main(){ 46 //freopen("bzoj4767.in","r",stdin); 47 //freopen("1.out","w",stdout); 48 E.x=read(), E.y=read(), n=read(); 49 A.x=read(), A.y=read(), B.x=read(), B.y=read(); 50 work(E.x,E.y); if(E.x<0 || E.y<0){printf("0 "); return 0;}//goto faq;} 51 //printf("%d %d ",E.x,E.y); 52 for(int i=1; i<=n; ++i){ 53 b[i].x=read(), b[i].y=read(); 54 work(b[i].x, b[i].y); 55 if(b[i].x<0 | b[i].y<0 || b[i].x>E.x || b[i].y>E.y) --n, --i; 56 } 57 b[++n]=(node){E.x,E.y}; 58 sort(b+1,b+n+1,cmp); 59 jc[0]=jc[1]=1; 60 for(int i=2; i<=L; ++i) jc[i] = jc[i-1] * i % mod; 61 jcn[L] = Pow(jc[L], mod-2); 62 for(int i=L-1; i>=0; --i) jcn[i] = jcn[i+1] * (i+1) % mod; 63 for(int i=1; i<=n; ++i){ 64 f[i] = cal(b[i].x, b[i].y); 65 for(int j=1; j<i; ++j) 66 f[i] = ((f[i] - f[j] * cal(b[i].x-b[j].x, b[i].y-b[j].y) % mod) % mod + mod) % mod; 67 } 68 69 cout<<f[n]<<endl; 70 /*faq: 71 fclose(stdout); 72 cout<<f[n]<<endl; 73 */ 74 return 0; 75 }
别问为什么都是 zsy 的博客，zsy 天下第一（雾）（希望他进队后能继续带我这个蒟蒻啊）
看完他的博客学到了一个新科技
就是你会发现 crt 合并的时候，逆元部分并不用 exgcd，只需要预处理就可以了
这样原本看起来挺长的 crt 代码现在只有 2 行了
不过 zsy 写的是 $t^2$ 次 crt，我给改成只在最后做 $1$ 次 crt 了，结果只快了 4ms，还多用了些空间……由此可见 crt 跑得很快
1 #include<bits/stdc++.h> 2 #define ll long long 3 #define K 205 4 #define L 1000005 5 using namespace std; 6 inline ll read(){ 7 ll x=0; bool f=1; char c=getchar(); 8 for(;!isdigit(c); c=getchar()) if(c=='-') f=0; 9 for(; isdigit(c); c=getchar()) x=(x<<3)+(x<<1)+(c^'0'); 10 if(f) return x; 11 return 0-x; 12 } 13 14 const int p[]={1000003, 3, 5, 6793, 10007}; 15 ll n,m,mod; int k; 16 ll inv[5][L],jc[5][L],jcn[5][L],f[K]; 17 struct node{ll x,y;}b[K]; 18 inline bool cmp(node a, node b){ 19 return a.x==b.x ? a.y<b.y : a.x<b.x; 20 } 21 22 ll Pow(ll x, int y, int i){ 23 ll ret=1; 24 while(y){ 25 if(y&1) ret=ret*x%p[i]; 26 x=x*x%p[i]; 27 y>>=1; 28 } 29 return ret; 30 } 31 ll C(ll n, ll m, int i){ 32 if(n<m) return 0; 33 return jc[i][n] * jcn[i][n-m] % p[i] * jcn[i][m] % p[i]; 34 } 35 ll lucas(ll n, ll m, int i){ 36 if(m==0 || n==m) return 1; 37 return C(n%p[i], m%p[i], i) * lucas(n/p[i], m/p[i], i) % p[i]; 38 } 39 int cal(ll x, ll y){ 40 if(x<0 || y<0) return 0; 41 if(mod==p[0]) return lucas(x+y,y,0); 42 ll ret=0; 43 for(int i=1; i<=4; ++i){ 44 ll a=lucas(x+y,y,i), t=inv[i][mod/p[i]%p[i]]; 45 ret = (ret + a * (mod/p[i]) % mod * t % mod) % mod; 46 } 47 return ret; 48 } 49 int main(){ 50 for(int i=0; i<=4; ++i){ 51 jc[i][0]=jcn[i][0]=inv[i][0]=inv[i][1]=1; 52 for(int j=2; j<p[i]; ++j) inv[i][j] = (p[i]-p[i]/j) * inv[i][p[i]%j] % p[i]; 53 for(int j=1; j<p[i]; ++j) jc[i][j] = jc[i][j-1] * j % p[i], jcn[i][j] = jcn[i][j-1] * inv[i][j] % p[i]; 54 } 55 n=read(), m=read(), k=read(), mod=read(); 56 for(int i=1; i<=k; ++i) b[i].x=read(), b[i].y=read(); 57 b[++k]=(node){n,m}; sort(b+1,b+k+1,cmp); 58 for(int i=1; i<=k; ++i){ 59 f[i]=cal(b[i].x,b[i].y); 60 for(int j=1; j<i; ++j) f[i] = (f[i] - f[j] * cal(b[i].x-b[j].x, b[i].y-b[j].y) % mod + mod) % mod; 61 } 62 cout<<f[k]<<endl; 63 return 0; 64 } 65 /* 66 3 4 3 1000003 67 3 0 68 1 1 69 2 2 70 */ 71
1 #include<bits/stdc++.h> 2 #define ll long long 3 #define K 205 4 #define L 1000005 5 using namespace std; 6 inline ll read(){ 7 ll x=0; bool f=1; char c=getchar(); 8 for(;!isdigit(c); c=getchar()) if(c=='-') f=0; 9 for(; isdigit(c); c=getchar()) x=(x<<3)+(x<<1)+(c^'0'); 10 if(f) return x; 11 return 0-x; 12 } 13 14 const int p[]={1000003, 3, 5, 6793, 10007}; 15 ll n,m,mod; int k; 16 ll inv[5][L],jc[5][L],jcn[5][L],f[5][K]; 17 struct node{ll x,y;}b[K]; 18 inline bool cmp(node a, node b){ 19 return a.x==b.x ? a.y<b.y : a.x<b.x; 20 } 21 22 ll Pow(ll x, int y, int i){ 23 ll ret=1; 24 while(y){ 25 if(y&1) ret=ret*x%p[i]; 26 x=x*x%p[i]; 27 y>>=1; 28 } 29 return ret; 30 } 31 ll C(ll n, ll m, int i){ 32 if(n<m) return 0; 33 return jc[i][n] * jcn[i][n-m] % p[i] * jcn[i][m] % p[i]; 34 } 35 ll lucas(ll n, ll m, int i){ 36 if(m==0 || n==m) return 1; 37 return C(n%p[i], m%p[i], i) * lucas(n/p[i], m/p[i], i) % p[i]; 38 } 39 int cal(ll x, ll y, int i){ 40 if(x<0 || y<0) return 0; 41 return lucas(x+y,y,i); 42 } 43 int main(){ 44 for(int i=0; i<=4; ++i){ 45 jc[i][0]=jcn[i][0]=inv[i][0]=inv[i][1]=1; 46 for(int j=2; j<p[i]; ++j) inv[i][j] = (p[i]-p[i]/j) * inv[i][p[i]%j] % p[i]; 47 for(int j=1; j<p[i]; ++j) jc[i][j] = jc[i][j-1] * j % p[i], jcn[i][j] = jcn[i][j-1] * inv[i][j] % p[i]; 48 } 49 n=read(), m=read(), k=read(), mod=read(); 50 for(int i=1; i<=k; ++i) b[i].x=read(), b[i].y=read(); 51 b[++k]=(node){n,m}; sort(b+1,b+k+1,cmp); 52 if(mod==p[0]){ 53 for(int i=1; i<=k; ++i){ 54 f[0][i]=cal(b[i].x,b[i].y,0); 55 for(int j=1; j<i; ++j) f[0][i] = (f[0][i] - f[0][j] * cal(b[i].x-b[j].x, b[i].y-b[j].y, 0) % mod + mod) % mod; 56 } 57 } 58 else{ 59 for(int P=1; P<=4; ++P){ 60 for(int i=1; i<=k; ++i){ 61 f[P][i]=cal(b[i].x,b[i].y,P); 62 for(int j=1; j<i; ++j) f[P][i] = (f[P][i] - f[P][j] * cal(b[i].x-b[j].x, b[i].y-b[j].y, P) % p[P] + p[P]) % p[P]; 63 } 64 } 65 for(int i=1; i<=4; ++i) 66 f[0][k] = (f[0][k] + f[i][k] * (mod/p[i]) % mod * inv[i][mod/p[i]%p[i]] % mod) % mod; 67 } 68 cout<<f[0][k]<<endl; 69 return 0; 70 } 71 /* 72 3 4 3 1000003 73 3 0 74 1 1 75 2 2 76 */