一年一度的圣诞节快要来到了。每年的圣诞节小E都会收到许多礼物,当然他也会送出许多礼物。不同的人物在小E心目中的重要性不同,在小E心中分量越重的人,收到的礼物会越多。小E从商店中购买了n件礼物,打算送给m个人,其中送给第i个人礼物数量为wi。请你帮忙计算出送礼物的方案数(两个方案被认为是不同的,当且仅当存在某个人在这两种方案中收到的礼物不同)。由于方案数可能会很大,你只需要输出模P后的结果。
看到题之后,应该就能写出式子
式子看起来很麻烦,其实很好想
当我们给第一个人送礼物的时候,我们有\(n\)个礼物,要选\(w_1\)个,方案就有\(C_n^{w_{1}}\)种
然后给第二个人送礼物的时候,我们剩下\(n-w_1\)个礼物,要选\(w_2\)个,方案有\(C_{n-w_1}^{w_2}\)种
就这样一直送到第\(m\)个人,而根据乘法原理,就得到这个式子了
这道题似乎就做完了,但在算\(C_n^m(mod\ p_i^{c_i})\)时又遇到了瓶颈:\(p_i^{c_i}\)不是质数
既然不是质数,我们就不能用\(Lucas\)定理来求,这就需要用到\(ExLucas\)了
虽然是扩展的,但和\(Lucas\)完全沾不上边
我们观察到题目给了\(P=\prod _{i=1}^{t}p_{i}^{c_i}\)
而\(p_i\)是质数,那么所有的\(p_i^{c_i}\)都是互质的
那么我们只要对于每个\(p_i^{k=c_i}\)求出其\(C_n^m\ mod\ p_i^k\)的值,然后用中国剩余定理就可以算出最小正整数解
问题转化成了如何快速的求\(C_n^m\ mod\ p^k\)
把组合数展开,我们得到
那么只要快速求出模意义下的阶乘就好了
如果我们要求\(x!\ mod\ p^k\),假设\(x=17,p=2,k=3\),观察一下式子
化式子\(17!=\)
而在模\(p^k=2^3=8\)意义下\(1\times3\times5\times7=9\times11\times13\times15\)
所以式子就变成了
那么\(x!\)就变成了\(p^n\times t!\times(a_1\times a_2\times…a_m)^r\times a_{m+1}\times…a_{q}\)
可以看出\(n=t=\left \lfloor \frac{x}{p} \right \rfloor,r=\left \lfloor \frac{x}{p^k} \right \rfloor\),其中\(t!\)可以一直递归求解,然后后面的暴力算
这样这个题就做完啦QAQ
Code
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#define int long long
using namespace std;
int p[200000],n,m,P,cnt,z[200000],a[200000],ans;
void exgcd(int a,int b,int &x,int &y) //扩欧
{
if (!b)x=1,y=0;
else
{
exgcd(b,a%b,x,y);
int t=x;
x=y;
y=t-a/b*y;
}
}
int mypow(int a,int x,int p) //快速幂
{
int s=1;
while (x)
{
if (x&1)s=s*a%p;
a=a*a%p;
x>>=1;
}
return s;
}
int fac(int n,int a,int b) //求阶乘
{
if (!n)return 1;
int s=1;
for (int i=1;i<=b;i++) //处理a1*a2*…am
if (i%a!=0)
s=s*i%b;
s=mypow(s,n/b,b);
for (int i=1;i<=n%b;i++) //处理am+1*…aq
if (i%a!=0)
s=s*i%b;
return s*fac(n/a,a,b)%b; //递归求解
}
int inv(int a,int b) //求逆元
{
int x,y;
exgcd(a,b,x,y);
return (x%b+b)%b;
}
int C(int m,int n,int a,int b) //处理组合数
{
int nn=fac(n,a,b),mm=fac(m,a,b),nm=fac(n-m,a,b),po=0; //求阶乘
for (int i=n;i;i/=a) //处理n^p中的p
po+=i/a;
for (int i=m;i;i/=a)
po-=i/a;
for (int i=n-m;i;i/=a)
po-=i/a;
return nn*mypow(a,po,b)%b*inv(mm,b)%b*inv(nm,b)%b;
}
signed main()
{
cin>>n>>m>>P;
int pp=P;
for (int i=2;i*i<=P;i++) //分解
if (pp%i==0)
{
int x=i;
z[++cnt]=i;
while (pp%x==0)x*=i;
x/=i;
p[cnt]=x;
pp/=x;
}
if (pp!=1)p[++cnt]=pp,z[cnt]=pp;
for (int i=1;i<=cnt;i++)
a[i]=C(m,n,z[i],p[i]);
for (int i=1;i<=cnt;i++) //中国剩余定理
{
int pi=P/p[i],x,y;
exgcd(pi,p[i],x,y);
ans=((ans+pi*a[i]*x%P)+P)%P;
}
cout<<ans<<endl;
return 0;
}