[bzoj2142]礼物(扩展lucas定理+中国剩余定理)
题意：n件礼物，送给m个人，每人的礼物数确定，求方案数。
解题关键：由于模数不是质数，所以由唯一分解定理，
$\bmod = p_1^{{k_1}}p_2^{{k_2}}......p_s^{{k_s}}$
然后，分别求出每个组合数模每个$p_i^{{k_i}}$的值，这里可以用扩展lucas定理求解,（以下其实就是扩展lucas定理的简略证明）
关于$C_n^m\% {p^k}$，
$C_n^m = \frac{{n!}}{{m!(n - m)!}}$，
我们以$n=19,p=3,k=2$为例，
$\begin{array}{l}
19! = 1*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19\\
= (1*2*4*5*7*8*10*11*13*14*16*17*19)*36*(1*2*3*4*5*6)
\end{array}$
通过观察，我们可以将将上式分成三部分，
第一部分，3的幂，快速幂可以直接求解；
第二部分，$n!$项，可以递归求解；
第三部分，$(1*2*4*5*7*8*10*11*13*14*16*17*19)$，此项在模${3^2}$意义下是存在循环节${p^k}$的，可以暴力求出一个循环节，然后重复即可，最后一个循环节的长度一定小于${p^k}$，可以在不提升复杂度的基础上暴力。
那我们回归最初的问题，关于$\frac{{n!}}{{m!(n - m)!}}\bmod {p^k}$的求解，由于在模意义下牵扯到求逆元，而不互质是不存在逆元的，所以需将阶乘中与模数不互质的部分提取出来，而这一定是$p$的倍数。
$\frac{{n!}}{{m!(n - m)!}} = \frac{{\frac{{n!}}{{{p^{{k_1}}}}}*{p^{{k_1}}}}}{{\frac{{m!}}{{{p^{{k_2}}}}}*{p^{{k_2}}}*\frac{{(n - m)!}}{{{p^{{k_3}}}}}*{p^{{k_3}}}}} = \frac{{\frac{{n!}}{{{p^{{k_1}}}}}}}{{\frac{{m!}}{{{p^{{k_2}}}}}*\frac{{(n - m)!}}{{{p^{{k_3}}}}}}}*{p^{{k_1} - {k_2} - {k_3}}}\bmod {p^k}$，
而$\frac{{n!}}{{{p^{{k_1}}}}},\frac{{m!}}{{{p^{{k_2}}}}},\frac{{(n - m)!}}{{{p^{{k_3}}}}}$是与${p^k}$互质的，可以求逆元。
所以，我们只需求出每个阶乘的第二部分和第三部分，关于$p$的幂，直接将三个阶乘的结果求出即可。
这种方法可以扩展到任意阶乘模非质数的情况。
最后用中国剩余定理组合一下。
注意最后不同质因数之间是互质的，所以直接crt即可，不需扩展crt。
最终的解为$C_n^{n - w[1]}C_{n - w[1]}^{w[2]}C_{n - w[1] - w[2]}^{w[3]}......$
法一：组合数求解。
1 #include<algorithm> 2 #include<iostream> 3 #include<cstring> 4 #include<cstdio> 5 #include<cmath> 6 #include<cstdlib> 7 typedef long long ll; 8 using namespace std; 9 ll mod,n,m,w[10],ans,x,y,module[10002],piset[10002],r[10002],num; 10 11 ll mod_pow(ll x,ll n,ll p){ 12 ll res=1; 13 while(n){ 14 if(n&1) res=res*x%p; 15 x=x*x%p; 16 n>>=1; 17 } 18 return res; 19 } 20 21 ll extgcd(ll a,ll b,ll &x,ll &y){ 22 ll d=a; 23 if(b) d=extgcd(b,a%b,y,x),y-=a/b*x; 24 else x=1,y=0; 25 return d; 26 } 27 28 ll inv(ll t,ll mod){ extgcd(t,mod,x,y);return (x+mod)%mod;} 29 30 ll multi(ll n,ll pi,ll pk){//求非互质的部分 31 if (!n) return 1; 32 ll ans=1; 33 for (ll i=2;i<=pk;i++) if(i%pi) ans=ans*i%pk; 34 ans=mod_pow(ans,n/pk,pk); 35 for (ll i=2;i<=n%pk;i++) if(i%pi) ans=ans*i%pk; 36 return ans*multi(n/pi,pi,pk)%pk; 37 } 38 39 40 ll exlucas(ll n,ll m,ll pi,ll pk){//组合数 c(n,m)mod pk=pi^k 41 if(m>n) return 0; 42 ll a=multi(n,pi,pk),b=multi(m,pi,pk),c=multi(n-m,pi,pk); 43 ll k=0; 44 for(ll i=n;i;i/=pi) k+=i/pi; 45 for(ll i=m;i;i/=pi) k-=i/pi; 46 for(ll i=n-m;i;i/=pi) k-=i/pi; 47 return a*inv(b,pk)%pk*inv(c,pk)%pk*mod_pow(pi,k,pk)%pk;//组合数求解完毕 48 } 49 50 ll crt(int n,ll *r,ll *m){ 51 ll M=1,ret=0; 52 for(int i=0;i<n;i++) M*=m[i]; 53 for(int i=0;i<n;i++){ 54 ll w=M/m[i]; 55 ret+=w*inv(w,m[i])*r[i]; 56 ret%=M; 57 } 58 return (ret+M)%M; 59 } 60 61 ll fz(ll n,ll *m,ll *piset){//分解质因子 62 ll num=0; 63 for (ll i=2;i*i<=n;i++){ 64 if(n%i==0){ 65 ll pk=1; 66 while(n%i==0) pk*=i,n/=i; 67 m[num]=pk; 68 piset[num]=i; 69 num++; 70 } 71 } 72 if(n>1) m[num]=n,piset[num]=n,num++; 73 return num; 74 } 75 76 ll excomb(ll n,ll m){ 77 for(int i=0;i<num;i++){ 78 r[i]=exlucas(n,m,piset[i],module[i]); 79 } 80 return crt(num,r,module); 81 } 82 83 84 int main(){ 85 scanf("%lld",&mod); 86 scanf("%lld%lld",&n,&m); 87 ll sum=0; 88 for(int i=1;i<=m;i++) scanf("%lld",&w[i]),sum+=w[i]; 89 if(n<sum){ puts("Impossible");return 0;}//puts会自动换行 90 num=fz(mod,module,piset); 91 ans=1; 92 for(int i=1;i<=m;i++){ 93 n-=w[i-1]; 94 ll a1=excomb(n,w[i]); 95 ans=ans*a1%mod; 96 } 97 printf("%lld\n",ans); 98 return 0; 99 }
法二：多项式系数求解。
最终解为：
$\left( {\begin{array}{*{20}{c}}
n\\
{{w_1}{w_2}...{w_n}(n - sum)}
\end{array}} \right) = \frac{{n!}}{{{w_1}!{w_2}!...{w_m}!(n - sum)!}}$
1 #include<algorithm> 2 #include<iostream> 3 #include<cstring> 4 #include<cstdio> 5 #include<cmath> 6 #include<cstdlib> 7 typedef long long ll; 8 using namespace std; 9 ll mod,P,n,m,w[10],ans,x,y,module[10002],piset[10002],r[10002],num,jc[10]; 10 11 ll mod_pow(ll x,ll n,ll p){ 12 ll res=1; 13 while(n){ 14 if(n&1) res=res*x%p; 15 x=x*x%p; 16 n>>=1; 17 } 18 return res; 19 } 20 21 ll extgcd(ll a,ll b,ll &x,ll &y){ 22 ll d=a; 23 if(b) d=extgcd(b,a%b,y,x),y-=a/b*x; 24 else x=1,y=0; 25 return d; 26 } 27 28 ll inv(ll t,ll mod){ extgcd(t,mod,x,y);return (x+mod)%mod;} 29 30 ll multi(ll n,ll pi,ll pk){//求非互质的部分 31 if (!n) return 1; 32 ll ans=1; 33 for (ll i=2;i<=pk;i++) if(i%pi) ans=ans*i%pk; 34 ans=mod_pow(ans,n/pk,pk); 35 for (ll i=2;i<=n%pk;i++) if(i%pi) ans=ans*i%pk; 36 return ans*multi(n/pi,pi,pk)%pk; 37 } 38 39 40 ll exlucas(ll n,ll pi,ll pk){ 41 ll ans=multi(n,pi,pk); 42 for(int i=0;i<m;i++){ 43 jc[i]=multi(w[i+1],pi,pk); 44 ans=ans*inv(jc[i],pk)%pk; 45 } 46 ll k=0; 47 for(ll i=n;i;i/=pi) k+=i/pi; 48 for(int i=1;i<=m;i++) for(ll j=w[i];j;j/=pi) k-=j/pi; 49 return ans*mod_pow(pi,k,pk)%pk; 50 } 51 52 ll crt(int n,ll *r,ll *m){ 53 ll M=1,ret=0; 54 for(int i=0;i<n;i++) M*=m[i]; 55 for(int i=0;i<n;i++){ 56 ll w=M/m[i]; 57 ret+=w*inv(w,m[i])*r[i]; 58 ret%=M; 59 } 60 return (ret+M)%M; 61 } 62 63 ll fz(ll n,ll *m,ll *piset){//分解质因子 64 ll num=0; 65 for (ll i=2;i*i<=n;i++){ 66 if(n%i==0){ 67 ll pk=1; 68 while(n%i==0) pk*=i,n/=i; 69 m[num]=pk; 70 piset[num]=i; 71 num++; 72 } 73 } 74 if(n>1) m[num]=n,piset[num]=n,num++; 75 return num; 76 } 77 78 int main(){ 79 scanf("%lld",&mod); 80 scanf("%lld%lld",&n,&m); 81 ll sum=0; 82 for(int i=1;i<=m;i++) scanf("%lld",&w[i]),sum+=w[i]; 83 if(n<sum){ puts("Impossible");return 0;}//puts会自动换行 84 if (sum<n) w[++m]=n-sum; 85 num=fz(mod,module,piset); 86 for(int i=0;i<num;i++) r[i]=exlucas(n,piset[i],module[i]); 87 printf("%lld\n",crt(num,r,module)); 88 return 0; 89 }
