Lflame

【bzoj 2142】 大组合数取模 --- 分段递归处理

    题意:求 m 个 C(a,b) 乘起来模 P 的答案。

    (m <= 5 , a、b <=10^9 , P 可能为非质数,P唯一分解后 1≤pi^ci≤10^5)


    这种做法在很久之前就写过一遍了,然而当时调了很久很久....

    于是再写了一遍,顺带也总结一下吧。


    考虑在模 pi^ci 下求 C(a,b) , 再用中国剩余定理合并起来。

    要怎么求呢 ? 由于有除法,且对一个数求逆元需要与模数互质,我们考虑将 C(a,b) 中的 pi 提出来,即写成 pi^k * t 的形式。

    这样 t 就可以求逆元了 。所以我们求出分子的 k1 和 t1 ,以及分母的 k2 和 t2,那么 C(a,b) 的 k 就是 k1-k2 , t 就是 t1/t2。得到 k 和 t 过后就可以知道 C(a,b) 是多少了。

    

    好的,现在问题转化成了怎么求分子和分母的 k、t,即需要求 a! 中的 k 和 t。令 MOD=pi^ci。我们将这 a 个数中是 pi 的倍数的数找出来,而其它的数直接累乘进 t。找到的这些数即为 pi、pi*2、pi*3… 然后将这些数都提出一个 pi 的因数,提出后这些数变成了 1、2、3...[a/p],也即 [a/p]!,递归处理即可。递归层数为 log 级别。

    但由于 a 会达到 10^9 级别,由题意知 MOD<=10^5,那么我们将 a! 按每 MOD 个数分成一段,每一整段的贡献是相同的,这样每一层最多处理O(MOD) 个数。

    整个算法就是这样了。


#include<iostream>

#include<cstdio>

#include<cstring>

#include<algorithm>

using namespace std;


int m ;

int que[1005] , qnum , ci[1005] ;

long long P , mm[1005] , aa[1005] , n , w[6] ;


long long pow_mod(long long a,long long b,long long c)

{

long long ret=1 ;

while(b){

if(b&1)ret=ret*a%c ;

a=a*a%c , b>>=1 ;

}

return ret ;

}


long long ex_gcd(long long a,long long b,long long &x,long long &y)

{

if(b==0){

x=1 , y=0 ;

return a;

}

long long ret=ex_gcd(b,a%b,x,y) ;

long long tx=x ;

x=y , y=tx-a/b*y ;

return ret;

}


long long get_rev(long long a,long long b)//求 a 在模 b 下的逆元 (保证存在)

{

long long x,y ;

ex_gcd(a,b,x,y) ;

return (x%b+b)%b ;

}


long long xnum , xx ;

long long get_jc(long long a,long long p,long long k,long long MOD,long long &num)

{

if(a<=1)return 1;

long long ret=1 , tmp=xx , tnum=xnum ;

long long t=a/MOD ;

tnum*=t , num=num+tnum , ret=ret*pow_mod(tmp,t,MOD)%MOD ;

t=a-t*MOD ;

for(int i=1;i<=t;++i){

if(i%p==0)tnum ++ , num ++ ;

elseret=ret*i%MOD ;

}

ret=ret*get_jc(tnum,p,k,MOD,num)%MOD ;

return ret;

}


long long get_C(long long a,long long b,long long p,long long k,long long MOD)

{

xnum=0 , xx=1 ;

for(int i=1;i<=MOD;++i) {

if(i%p==0)xnum++ ;

elsexx=xx*i%MOD ;

}

long long ret=1 , num=0 , tnum=0 , tmp=1 ;

ret=ret*get_jc(a,p,k,MOD,num)%MOD ;

tmp=tmp*get_jc(b,p,k,MOD,tnum)%MOD*get_jc(a-b,p,k,MOD,tnum)%MOD ;

num-=tnum , ret=ret*get_rev(tmp,MOD) ;

if(num>=k)return 0 ;

for(int i=1;i<=num;++i)ret=ret*p%MOD ;

return ret;

}


long long C(long long a,long long b)

{

long long ret=0 ;

for(int i=1;i<=qnum;++i){

long long t=get_C(a,b,que[i],ci[i],mm[i]) , tmp=P/mm[i] ;

ret=(ret+t*tmp%P*get_rev(tmp,mm[i])%P)%P ;

}

return ret;

}


int main()

{

scanf("%lld%lld%d",&P,&n,&m) ;

long long sum=0 ;

for(int i=1;i<=m;++i)scanf("%lld",&w[i]) , sum+=w[i] ;

if(sum>n){

printf("Impossible\n") ;

return 0 ;

}

long long tmp=P ;

for(int i=2;(long long)i*i<=P;++i){

while(tmp%i==0){

tmp/=i ;

if(que[qnum]!=i)que[++qnum]=i ;

ci[qnum]++ ;

}

}

if(tmp>1)que[++qnum]=(int)tmp , ci[qnum]=1 ;

for(int i=1;i<=qnum;++i){

mm[i]=1 ;

for(int j=1;j<=ci[i];++j)mm[i]*=que[i] ;

}

long long ans=1 , res=n ;

for(int i=1;i<=m;++i){

ans=ans*C(res,w[i])%P ;

res-=w[i] ;

}

if(!ans){

printf("Impossible\n") ; return 0 ;

}

printf("%lld\n",ans) ;

return 0 ;

}

评论