BZOJ2142: 礼物

Posted 2020-08-31

传送门

组合数模终极版！

扩展lucas搞似乎不会T掉，我用的是另一种方法，顺便吐槽一下某人的辣鸡题解，他的一个公式打错了，我也没去验证那个公式的正确性就直接套用导致小数据轻松过大数据狂wa。

怎么把组合数列出来不说了，大家都能看出来。

问题是怎么把一个组合数对一个合数模，我们知道这样一个公式：$C_N^M=\frac{N!}{M! \times (N-M)! }$。那么我们现在就转化成了怎么将一个数的阶乘对一个合数取余。

这里指的当然是把含有模数的因子给分离出来，因为只要把因子给分离出来我们组合数就可以很好的分子分母上下化简。

我们首先考虑阶乘的$\prod \nolimits^{i \in [1,p_i^j-1]} i$这一段，我们可以考虑把每一段$p_i$的倍数提出来，因为$p_i$是一个质数，所以其他数的模肯定不会出现0的情况。

然后我们把提出来的数放到一块，从这里再提出来一个$p_i$，然后递归处理即可。

对于左边的，因为这一段的答案和以后同余，作为循环节快速幂即可。

//BZOJ 2142
//by Cydiater
//2017.2.20
#include <iostream>
#include <queue>
#include <map>
#include <ctime>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <string>
#include <algorithm>
#include <cstdio>
#include <iomanip>
#include <bitset>
#include <set>
#include <vector>
#include <complex>
using namespace std;
#define ll long long
#define up(i,j,n)	for(ll i=j;i<=n;i++)
#define down(i,j,n)	for(ll i=j;i>=n;i--)
#define cmax(a,b)	a=max(a,b)
#define cmin(a,b)	a=min(a,b)
#define FILE		"nt2011_gift"
inline ll read(){
	char ch=getchar();ll x=0,f=1;
	while(ch>‘9‘||ch<‘0‘){if(ch==‘-‘)f=-1;ch=getchar();}
	while(ch>=‘0‘&&ch<=‘9‘){x=x*10+ch-‘0‘;ch=getchar();}
	return x*f;
}
ll mod,M,N,a[10],cnt,w[10],sum=0,ans=1;
struct Factor{
	ll p,pc;
}fac[10];
struct Split{
	ll base,ind,P;
};
void exgcd(ll a,ll b,ll &x,ll &y){
	if(!b){x=1;y=0;return;}
	exgcd(b,a%b,y,x);
	y-=a/b*x;
}
ll inv(ll num,ll P){
	ll x,y;
	exgcd(num,P,x,y);
	x=((x%P+P)%P+P)%P;
	return x;
}
Split operator * (Split x,Split y){
	ll base=x.base*y.base%x.P,ind=x.ind+y.ind;
	return (Split){base,ind,x.P};
}
Split operator / (Split x,Split y){
	ll base=x.base*inv(y.base,y.P)%y.P,ind=x.ind-y.ind;
	return (Split){base,ind,x.P};
}
namespace solution{
	ll quick_pow(ll base,ll ind,ll P){
		ll tmp=1;
		while(ind){
			if(ind&1)(tmp*=base)%=P;
			(base*=base)%=P;ind>>=1;
		}
		return tmp;
	}
   	Split split(ll num,ll k){
        	Split res;
      		ll p=fac[k].p,pc=fac[k].pc;
        	res.base=1;res.ind=0;res.P=pc;
        	ll tmp=1;
        	up(i,2,pc-1)if(i%p!=0)(tmp*=i)%=pc;
        	(tmp=quick_pow(tmp,(num/pc),pc))%=pc;
        	if(num>p)res=split(num/p,k);
        	up(i,1,num%pc)if(i%p!=0)(tmp*=(i+pc))%=pc;
        	(res.base*=tmp)%=pc;
        	res.ind+=num/p;
        	return res;
    	}
	ll CRT(){
		ll a1=a[1],m1=fac[1].pc;
		up(i,2,cnt){
			ll a2=a[i],m2=fac[i].pc,x,y;
			exgcd(m1,m2,x,y);
			x=(x*(a2-a1)%m2+m2)%m2;
			a1+=x*m1;
			m1*=m2;
		}
		return a1;
	}
	ll C(ll A,ll B){
		if(A<B)return 0;
		up(i,1,cnt){
			Split tmp=split(A,i)/(split(B,i)*split(A-B,i));
			a[i]=tmp.base*quick_pow(fac[i].p,tmp.ind,fac[i].pc)%fac[i].pc;
		}
		return CRT();
	}
	void Prepare(){
		mod=read();M=read();N=read();
		up(i,1,N)sum+=(w[i]=read());		
		ll tmp=mod,lim=sqrt(1.0*tmp);
		up(i,2,lim){
			if(!(tmp%i))fac[++cnt]=(Factor){i,1};
			while(tmp&&!(tmp%i)){
				fac[cnt].pc*=i;
				tmp/=i;
			}
		}
		if(tmp>1)fac[++cnt]=(Factor){tmp,tmp};
	}
	void Solve(){
		if(sum>M){
			puts("Impossible");
			return;
		}
		up(i,1,N){
			(ans*=C(M,w[i]))%=mod;
			M-=w[i];
		}
		cout<<ans<<endl;
	}
}
int main(){
	freopen(FILE".in","r",stdin);
	freopen(FILE".out","w",stdout);
	//freopen("input.in","r",stdin);
	using namespace solution;
	Prepare();
	Solve();
	return 0;
}