[bzoj2142]礼物(扩展lucas定理+中国剩余定理)

2023-02-26,,

题意: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]}......$

法一:组合数求解。

 #include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<cstdlib>
typedef long long ll;
using namespace std;
ll mod,n,m,w[],ans,x,y,module[],piset[],r[],num; ll mod_pow(ll x,ll n,ll p){
ll res=;
while(n){
if(n&) res=res*x%p;
x=x*x%p;
n>>=;
}
return res;
} ll extgcd(ll a,ll b,ll &x,ll &y){
ll d=a;
if(b) d=extgcd(b,a%b,y,x),y-=a/b*x;
else x=,y=;
return d;
} ll inv(ll t,ll mod){ extgcd(t,mod,x,y);return (x+mod)%mod;} ll multi(ll n,ll pi,ll pk){//求非互质的部分
if (!n) return ;
ll ans=;
for (ll i=;i<=pk;i++) if(i%pi) ans=ans*i%pk;
ans=mod_pow(ans,n/pk,pk);
for (ll i=;i<=n%pk;i++) if(i%pi) ans=ans*i%pk;
return ans*multi(n/pi,pi,pk)%pk;
} ll exlucas(ll n,ll m,ll pi,ll pk){//组合数 c(n,m)mod pk=pi^k
if(m>n) return ;
ll a=multi(n,pi,pk),b=multi(m,pi,pk),c=multi(n-m,pi,pk);
ll k=;
for(ll i=n;i;i/=pi) k+=i/pi;
for(ll i=m;i;i/=pi) k-=i/pi;
for(ll i=n-m;i;i/=pi) k-=i/pi;
return a*inv(b,pk)%pk*inv(c,pk)%pk*mod_pow(pi,k,pk)%pk;//组合数求解完毕
} ll crt(int n,ll *r,ll *m){
ll M=,ret=;
for(int i=;i<n;i++) M*=m[i];
for(int i=;i<n;i++){
ll w=M/m[i];
ret+=w*inv(w,m[i])*r[i];
ret%=M;
}
return (ret+M)%M;
} ll fz(ll n,ll *m,ll *piset){//分解质因子
ll num=;
for (ll i=;i*i<=n;i++){
if(n%i==){
ll pk=;
while(n%i==) pk*=i,n/=i;
m[num]=pk;
piset[num]=i;
num++;
}
}
if(n>) m[num]=n,piset[num]=n,num++;
return num;
} ll excomb(ll n,ll m){
for(int i=;i<num;i++){
r[i]=exlucas(n,m,piset[i],module[i]);
}
return crt(num,r,module);
} int main(){
scanf("%lld",&mod);
scanf("%lld%lld",&n,&m);
ll sum=;
for(int i=;i<=m;i++) scanf("%lld",&w[i]),sum+=w[i];
if(n<sum){ puts("Impossible");return ;}//puts会自动换行
num=fz(mod,module,piset);
ans=;
for(int i=;i<=m;i++){
n-=w[i-];
ll a1=excomb(n,w[i]);
ans=ans*a1%mod;
}
printf("%lld\n",ans);
return ;
}

法二:多项式系数求解。

最终解为:

$\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)!}}$

 #include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<cstdlib>
typedef long long ll;
using namespace std;
ll mod,P,n,m,w[],ans,x,y,module[],piset[],r[],num,jc[]; ll mod_pow(ll x,ll n,ll p){
ll res=;
while(n){
if(n&) res=res*x%p;
x=x*x%p;
n>>=;
}
return res;
} ll extgcd(ll a,ll b,ll &x,ll &y){
ll d=a;
if(b) d=extgcd(b,a%b,y,x),y-=a/b*x;
else x=,y=;
return d;
} ll inv(ll t,ll mod){ extgcd(t,mod,x,y);return (x+mod)%mod;} ll multi(ll n,ll pi,ll pk){//求非互质的部分
if (!n) return ;
ll ans=;
for (ll i=;i<=pk;i++) if(i%pi) ans=ans*i%pk;
ans=mod_pow(ans,n/pk,pk);
for (ll i=;i<=n%pk;i++) if(i%pi) ans=ans*i%pk;
return ans*multi(n/pi,pi,pk)%pk;
} ll exlucas(ll n,ll pi,ll pk){
ll ans=multi(n,pi,pk);
for(int i=;i<m;i++){
jc[i]=multi(w[i+],pi,pk);
ans=ans*inv(jc[i],pk)%pk;
}
ll k=;
for(ll i=n;i;i/=pi) k+=i/pi;
for(int i=;i<=m;i++) for(ll j=w[i];j;j/=pi) k-=j/pi;
return ans*mod_pow(pi,k,pk)%pk;
} ll crt(int n,ll *r,ll *m){
ll M=,ret=;
for(int i=;i<n;i++) M*=m[i];
for(int i=;i<n;i++){
ll w=M/m[i];
ret+=w*inv(w,m[i])*r[i];
ret%=M;
}
return (ret+M)%M;
} ll fz(ll n,ll *m,ll *piset){//分解质因子
ll num=;
for (ll i=;i*i<=n;i++){
if(n%i==){
ll pk=;
while(n%i==) pk*=i,n/=i;
m[num]=pk;
piset[num]=i;
num++;
}
}
if(n>) m[num]=n,piset[num]=n,num++;
return num;
} int main(){
scanf("%lld",&mod);
scanf("%lld%lld",&n,&m);
ll sum=;
for(int i=;i<=m;i++) scanf("%lld",&w[i]),sum+=w[i];
if(n<sum){ puts("Impossible");return ;}//puts会自动换行
if (sum<n) w[++m]=n-sum;
num=fz(mod,module,piset);
for(int i=;i<num;i++) r[i]=exlucas(n,piset[i],module[i]);
printf("%lld\n",crt(num,r,module));
return ;
}

[bzoj2142]礼物(扩展lucas定理+中国剩余定理)的相关教程结束。

《[bzoj2142]礼物(扩展lucas定理+中国剩余定理).doc》

下载本文的Word格式文档,以方便收藏与打印。