#include<iostream>#include<vector>classBinomModPrime{intp;std::vector<int>fa,ifa;// Calculate binom(n, k) mod p for n, k < p.intcalc(intn,intk){if(n<k)return0;longlongres=fa[n];res=(res*ifa[k])%p;res=(res*ifa[n-k])%p;returnres;}public:BinomModPrime(intp):p(p),fa(p),ifa(p){// Factorials mod p till p.fa[0]=1;for(inti=1;i<p;++i){fa[i]=(longlong)fa[i-1]*i%p;}// Inverse of factorials mod p till p.ifa[p-1]=p-1;// Wilson's theorem.for(inti=p-1;i;--i){ifa[i-1]=(longlong)ifa[i]*i%p;}}// Calculate binom(n, k) mod p.intbinomial(longlongn,longlongk){longlongres=1;while(n||k){res=(res*calc(n%p,k%p))%p;n/=p;k/=p;}returnres;}};intmain(){intt,p;std::cin>>t>>p;BinomModPrimebm(p);for(;t;--t){longlongn,k;std::cin>>n>>k;std::cout<<bm.binomial(n,k)<<'\n';}return0;}
该实现的时间复杂度为 ,其中, 为询问次数。
exLucas 算法
Lucas 定理中对于模数 要求必须为素数,那么对于 不是素数的情况,就需要用到 exLucas 算法。虽然名字如此,该算法实际操作时并没有用到 Lucas 定理。它的关键步骤是 计算素数幂模下的阶乘。上文的第二个证明指出了它与 Lucas 定理的联系。
#include<iostream>#include<vector>// Extended Euclid.voidex_gcd(inta,intb,int&x,int&y){if(!b){x=1;y=0;}else{ex_gcd(b,a%b,y,x);y-=a/b*x;}}// Inverse of a mod m.intinverse(inta,intm){intx,y;ex_gcd(a,m,x,y);return(x%m+m)%m;}// Coefficient in CRT.intcrt_coeff(intm_i,intm){longlongmm=m/m_i;mm*=inverse(mm,m_i);returnmm%m;}// Binominal Coefficient Calculator Modulo Prime Power.classBinomModPrimePower{intp,a,pa;std::vector<int>f;// Obtain multiplicity of p in n!.longlongnu(longlongn){longlongcount=0;do{n/=p;count+=n;}while(n);returncount;}// Calculate (n!)_p mod pa.longlongfact_mod(longlongn){boolneg=p!=2||pa<=4;longlongres=1;while(n>1){if((n/pa)&neg)res=pa-res;res=res*f[n%pa]%pa;n/=p;}returnres;}public:BinomModPrimePower(intp,inta,intpa):p(p),a(a),pa(pa),f(pa){// Pretreatment.f[0]=1;for(inti=1;i<pa;++i){f[i]=i%p?(longlong)f[i-1]*i%pa:f[i-1];}}// Calculate Binom(n, k) mod pa.intbinomial(longlongn,longlongk){longlongv=nu(n)-nu(n-k)-nu(k);if(v>=a)return0;autores=fact_mod(n-k)*fact_mod(k)%pa;res=fact_mod(n)*inverse(res,pa)%pa;for(;v;--v)res*=p;returnres%pa;}};// Binominal Coefficient Calculator.classBinomMod{intm;std::vector<BinomModPrimePower>bp;std::vector<longlong>crt_m;public:BinomMod(intn):m(n){// Factorize.for(intp=2;p*p<=n;++p){if(n%p==0){inta=0,pa=1;for(;n%p==0;n/=p,++a,pa*=p);bp.emplace_back(p,a,pa);crt_m.emplace_back(crt_coeff(pa,m));}}if(n>1){bp.emplace_back(n,1,n);crt_m.emplace_back(crt_coeff(n,m));}}// Calculate Binom(n, k) mod m.intbinomial(longlongn,longlongk){longlongres=0;for(size_ti=0;i!=bp.size();++i){res=(bp[i].binomial(n,k)*crt_m[i]+res)%m;}returnres;}};intmain(){intt,m;std::cin>>t>>m;BinomModbm(m);for(;t;--t){longlongn,k;std::cin>>n>>k;std::cout<<bm.binomial(n,k)<<'\n';}return0;}