BZOJ2142 礼物 EXlucas

全场最佳...毒瘤题

一眼得出公式：C(n,sum)*C(sum,a[1])*C(sum-a[1],a[2])*....

然而数据范围n<=1e9

显然可以用lucas定理算，但是模数tmd不一定是个质数

于是需要用到ex_lucas......

把p拆成若干个质因子幂乘积的形式，然后对于每一个质因子幂次求解，最后中国剩余定理乘起来即可

具体过程：

1、拆

得到底数q与q^t，t表示q因子的数量

这种情况下

我们可以对C求解

原因是：这种情况下，我们只需要把所有q因子提出来，就可以求逆元

如何求解：

2、求 C(i,j)%(q^t)

我们可以吧含有q的项拆出来，方法如下：

例如q=3,t=2;

求20!在模域q^t下的值

20!=1*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19*20

=1*2*4*5 * 7*8*10*11 * 13*14*16*17  * 19*20 *3^6*(1+2+3+4+5+6)

=(1*2*4*5)^3 *3^t *19*20 * 6!

于是我们可以对6!继续此操作 最终求出20!

这个地方解释一下为什么1*2*4*5=7*8*10*11

7*8*10*11=(1+6)*(2+6)*(4+6)*(5+6) 显然，展开后只有1*2*4*5这一项不含因子6

于是我们就可以在log*某较大常数  的时间内求出k!的值了

然后暴力算一下i!中q的数量-j!和(i-j)!中q的数量，乘进答案里

余下的部分则可以用上述算法求出阶乘后用exgcd求逆元，暴力算出组合数

3、中国剩余定理

求关于方程

x%p1=a1

x%p2=a2

x%p3=a3

...

所有p都互质

的解

我们设π p=P

P/p[i]=V[i]

V[i]在p[i]模域下的逆元为REV[i]

为什么这样？

我们可以发现一个很有用的性质

REV[I]*V[i]这个东西

模p[i]=1;模p[j] （j!=i) =0

因为V[i]是除了p[i]以外所有p的乘积

而REV是V[i]在p[i]意义下的逆元

所以我们的一个解ANS=sigma(a[i]*REV[i]*V[i])%P

其他解。。如果要用的话，ANS'=ANS+K*P;

以上就是ex_Lucas的所有要用到的东西了

CODE:

#include<iostream>
#include<cstdio>
#include<queue>
#include<cmath>
#include<algorithm>
#include<map>
#include<cstring>
#define For(i,j,k) for(ll i=j;i<=k;++i)
#define Dow(i,j,k) for(ll i=k;i>=j;--i)
#define pb push_back
#define inf 1e18
#define maxn 500001
#define ll long long
using namespace std;
inline ll read()
{
	ll t=0,f=1;char c=getchar();
	while(c<'0'||c>'9')	{if(c=='-')	f=-1;c=getchar();}
	while(c<='9'&&c>='0')	t=t*10LL+c-48LL,c=getchar();
	return t*f;
}
ll mo,n,m,sum,x[200001];
inline ll ksm(ll x,ll y,ll m){ll sum=1;for(;y;y/=2){if(y&1)	sum=sum*x%m;x=x*x%m;}	return sum;}
inline void exgcd(ll a,ll b,ll &x,ll &y)
{
	if(b==0)	{x=1;y=0;}
		else	exgcd(b,a%b,y,x),y-=(a/b)*x;
}
inline ll rev(ll x,ll pr)
{
	ll tx=0,ty=0;
	exgcd(x,pr,tx,ty);
	return (tx%pr+pr)%pr;
}
inline ll mul(ll x,ll p,ll pr)
{
	if(x==0)	return 1;
	ll tmp=1;
	For(i,1,pr)	if(i%p)	tmp*=i,tmp%=pr;
	tmp=ksm(tmp,x/pr,pr);
	ll tep=x%pr;
	For(i,2,tep)	if(i%p)	tmp*=i,tmp%=pr;
	tmp=tmp*mul(x/p,p,pr);
	return tmp%pr;	
}
inline ll crt(ll x,ll y,ll p,ll pr)
{
	ll t1=mul(x,p,pr),t2=mul(y,p,pr),t3=mul(x-y,p,pr);
	ll tot=0;
	ll tmp=x;while(tmp)	tot+=tmp/p,tmp/=p;
	tmp=y;while(tmp)	tot-=tmp/p,tmp/=p;
	tmp=x-y;while(tmp)	tot-=tmp/p,tmp/=p;
	ll ans=t1*rev(t2,pr)%pr*rev(t3,pr)%pr*ksm(p,tot,pr)%pr;
	return ans*(mo/pr)%mo*rev(mo/pr,pr)%mo;
}
inline ll lucas(int x,int y)
{
	ll tmp=mo,ans=0;
	For(i,2,tmp)
	{
		ll mul=1;
		while(tmp%i==0)	mul*=i,tmp/=i;
		if(mul!=1)	ans=(ans+crt(x,y,i,mul));	
	}
	return ans;
}
int main()
{
	mo=read();
	n=read();m=read();
	For(i,1,m)	x[i]=read(),sum+=x[i];
	if(sum>n)	{puts("Impossible");return 0;}
	ll ans=lucas(n,sum)%mo;
	For(i,1,m)
	{
		ans=ans*lucas(sum,x[i]),sum-=x[i],ans%=mo;
	}
	printf("%lld\n",ans);
}
  
  

   
   优化后的EX_LUCAS
  
  
  
  

   
   
#include<iostream>
#include<cstdio>
#include<queue>
#include<cmath>
#include<algorithm>
#include<map>
#include<cstring>
#define For(i,j,k) for(ll i=j;i<=k;++i)
#define Dow(i,j,k) for(ll i=k;i>=j;--i)
#define pb push_back
#define inf 1e18
#define maxn 500001
#define ll long long
using namespace std;
inline ll read()
{
	ll t=0,f=1;char c=getchar();
	while(c<'0'||c>'9')	{if(c=='-')	f=-1;c=getchar();}
	while(c<='9'&&c>='0')	t=t*10LL+c-48LL,c=getchar();
	return t*f;
}
ll mo,n,m,sum,x[200001],rest,blo;
inline ll ksm(ll x,ll y,ll m){ll sum=1;for(;y;y/=2){if(y&1)	sum=sum*x%m;x=x*x%m;}	return sum;}
inline void exgcd(ll a,ll b,ll &x,ll &y)
{
	if(b==0)	{x=1;y=0;}
		else	exgcd(b,a%b,y,x),y-=(a/b)*x;
}
inline ll rev(ll x,ll pr)
{
	ll tx=0,ty=0;
	exgcd(x,pr,tx,ty);
	return (tx%pr+pr)%pr;
}
ll pre[2000001];
inline ll mul(ll x,ll p,ll pr)
{
	if(x==0)	return 1;
	ll tmp=pre[pr];
	tmp=ksm(tmp,x/pr,pr);
	tmp=tmp*pre[x%pr];
	tmp=tmp*mul(x/p,p,pr);
	return tmp%pr;	
}
inline ll crt(ll x,ll y,ll p,ll pr)
{
	pre[0]=1;
	For(i,1,pr) pre[i]=pre[i-1]*(i%p?i:1)%pr;	
	ll t1=mul(x,p,pr),t2=mul(y,p,pr),t3=mul(x-y,p,pr);
	ll tot=0;
	ll tmp=x;while(tmp)	tot+=tmp/p,tmp/=p;
	tmp=y;while(tmp)	tot-=tmp/p,tmp/=p;
	tmp=x-y;while(tmp)	tot-=tmp/p,tmp/=p;
	ll ans=t1*rev(t2,pr)%pr*rev(t3,pr)%pr*ksm(p,tot,pr)%pr;
	return ans*(mo/pr)%mo*rev(mo/pr,pr)%mo;
}
inline ll lucas(ll x,ll y)
{
	ll tmp=mo,ans=0;
	For(i,2,tmp)
	{
		ll mul=1;
		while(tmp%i==0)	mul*=i,tmp/=i;
		if(mul!=1)	ans=(ans+crt(x,y,i,mul))%mo;	
	}
	return ans;
}