[SDOI2013]方程
Description
给定方程
X1+X2+. +Xn=M
我们对第l..N1个变量进行一些限制：
Xl < = A
X2 < = A2
Xn1 < = An1
我们对第n1 + 1..n1+n2个变量进行一些限制：
Xn1+l > = An1+1
Xn1+2 > = An1+2
Xnl+n2 > = Anl+n2
求：在满足这些限制的前提下，该方程正整数解的个数。
答案可能很大，请输出对p取模后的答案，也即答案除以p的余数。
Input
输入含有多组数据，第一行两个正整数T，p。T表示这个测试点内的数据组数，p的含义见题目描述。
对于每组数据，第一行四个非负整数n，n1，n2，m。
第二行nl+n2个正整数，表示A1..n1+n2。请注意，如果n1+n2等于0，那么这一行会成为一个空行。
Output
共T行，每行一个正整数表示取模后的答案。
Sample Input
3 1 1 6
3 3
3 0 0 5
3 1 1 3
3 3
Sample Output
6
0
【样例说明】
对于第一组数据，三组解为(1，3，2)，(1,4,1),(2,3,1)
对于第二组数据，六组解为(1，1，3)，(1,2,2),(1,3,1),(2,1,2),(2,2,1),(3,1,1)
HINT
n < = 10^9 , n1 < = 8 , n2 < = 8 , m < = 10^9 ,p<=437367875
对于l00%的测试数据： T < = 5，1 < = A1..n1_n2 < = m，n1+n2 < = n
出题人很贱，故意使模数为合数
首先，是简单的排列组合
对于n2个条件，直接m-=ai-1就行了
对于满足n1条件，可以用容斥，枚举一个不满足减去，两个加上，以此类推
不满足条件直接把m-=ai，就是说直接分ai给i使i不符合
用隔版法，把m个物品分成n个非空子集：C(n-1,m-1)
由于n，m很大，所以只能lucas，但因为模数是合数，所以不可以
于是就有了拓展lucas
http://www.cnblogs.com/Przz/p/5409566.html
先分解模数Mod的质因数，Mod=p1^t1*p2^t2*.....
对于每个p=p1^t1取模，这样就可以用中国剩余定理(要求模数互质)
这样是为了方便快速阶乘取模
C(x,y)=y!*(x!)-1*(y-x!)-1
说明一下快速阶乘取模，因为p只含一个因数,拿p=3^2,即pi=3，ti=2，n=19举例
n！=(1*2×4×5×7×8×10×11×13×14×16×17×19)*(3*6*9*12*15*18)
=(1*2×4×5×7×8×10×11×13×14×16×17)*3^6*(1*2*3*4*5*6)*19
因为1×2×3×4×5×6×7×8≡10×11×12×13×14×15×16×17 (mod 9)
n!=(fac[9-1]^(19/9))*fac[19%9]*((19/3)!%p)
后面部分递归处理，至于红色部分，在阶乘前就算出x!,y!,(y-x)!中pi的含量a,b,c
用f(x)=x/pi+f[x/pi]递归求出
将答案乘pi^(a-b-c)
还有求出了阶乘，还要算逆元，由于p=pi^ti是合数
所以费马和线性模逆元公式都无法使用，但庆幸的是可以用拓展欧几里德
ax≡1(mod p) 的条件是a与p互质
意味着a只要含有pi就没有逆元
很简单，把所有pi的倍数事先挑出来(上面红色部分就是为了现在)
在阶乘前缀积中不予计算，这样算出来的阶乘就可以套用拓展欧几里德求逆元了
至于中国剩余定理这里不做说明
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 using namespace std; 6 typedef long long ll; 7 ll fac[100011],tot,B[10001],d[10001],pri[10001],A[10001],c[10001],n,n1,n2,m; 8 ll a[10001],Mod,ans; 9 ll qpow(ll x,ll y,ll p) 10 { 11 ll res=1; 12 while (y) 13 { 14 if (y&1) res=res*x%p; 15 x=x*x%p; 16 y/=2; 17 } 18 return res; 19 } 20 ll exgcd(ll a,ll b,ll &x,ll &y) 21 { 22 if (!b) 23 { 24 x=1;y=0; 25 return a; 26 } 27 ll d=exgcd(b,a%b,x,y); 28 int t=x;x=y;y=t-(a/b)*y; 29 return d; 30 } 31 ll reverse(ll a,ll b)//求逆元 32 { 33 ll x,y; 34 exgcd(a,b,x,y); 35 return (x%b+b)%b; 36 } 37 ll calfac(ll x,ll p,ll t)//快速阶乘取模 38 { 39 if (x<t) return fac[x]; 40 ll s=qpow(fac[p-1],x/p,p); 41 s=s*fac[x%p]%p; 42 s=s*calfac(x/t,p,t)%p; 43 return s; 44 } 45 ll calc(ll x,ll y,ll p,ll t)//组合数取模(拓展lucas) 46 {int i; 47 ll ap=0,bp=0,cp=0; 48 fac[0]=1; 49 for (i=1;i<=p-1;i++) 50 if (i%t==0) fac[i]=fac[i-1]; 51 else fac[i]=fac[i-1]*i%p; 52 for (i=y;i;i/=t) ap+=i/t; 53 for (i=x;i;i/=t) bp+=i/t; 54 for (i=y-x;i;i/=t) cp+=i/t; 55 ap=ap-bp-cp; 56 ll s=qpow(t,ap,p); 57 ap=calfac(y,p,t);bp=calfac(x,p,t);cp=calfac(y-x,p,t); 58 s=(((s*ap%p)*reverse(bp,p)%p)*reverse(cp,p)%p); 59 return s; 60 } 61 ll C(ll x,ll y,ll p)//中国剩余定理 62 {int i; 63 if (x>y) return 0; 64 for (i=1;i<=tot;i++) 65 B[i]=calc(x,y,d[i],pri[i]); 66 for (i=1;i<=tot;i++) 67 A[i]=reverse(c[i],d[i]); 68 ll s=0; 69 for (i=1;i<=tot;i++) 70 { 71 s+=((B[i]*A[i]%p)*c[i]%p); 72 s%=p; 73 } 74 return s; 75 } 76 void dfs(ll id,ll m,ll cnt)//容斥原理 77 { 78 if (id>n1) 79 { 80 ans=(ans+C(n-1,m-1,Mod)*cnt+Mod)%Mod;//隔版法求方案 81 return; 82 } 83 dfs(id+1,m,cnt); 84 if (m>=a[id]) dfs(id+1,m-a[id],-cnt); 85 } 86 int main() 87 {ll T,x,i; 88 cin>>T>>Mod; 89 x=Mod; 90 for (i=2;i*i<=Mod;i++) 91 if (x%i==0) 92 {int cnt=0; 93 d[++tot]=1; 94 pri[tot]=i; 95 while (x%i==0) 96 { 97 ++cnt; 98 x/=i; 99 d[tot]*=i; 100 } 101 102 } 103 if (x!=1) 104 { 105 tot++; 106 d[tot]=x; 107 pri[tot]=x; 108 } 109 for (i=1;i<=tot;i++) c[i]=Mod/d[i]; 110 while (T--) 111 { 112 cin>>n>>n1>>n2>>m; 113 for (i=1;i<=n1;i++) 114 scanf("%lld",&a[i]); 115 for (i=1;i<=n2;i++) 116 { 117 scanf("%lld",&x); 118 if (x) m-=x-1; 119 } 120 if (m<n1) 121 cout<<0<<endl; 122 else 123 { 124 ans=0; 125 dfs(1,m,1); 126 cout<<ans<<endl; 127 } 128 } 129 }
