两侧同时换到之前的修订记录 前一修订版 | |||
2020-2021:teams:legal_string:jxm2001:数论_5 [2020/10/16 14:25] jxm2001 |
2020-2021:teams:legal_string:jxm2001:数论_5 [2020/10/16 15:01] (当前版本) jxm2001 |
||
---|---|---|---|
行 70: | 行 70: | ||
==== 算法简介 ==== | ==== 算法简介 ==== | ||
- | $O(\sum p_i^k+\sum \log_{p_i^k}n)$ 计算 ${n\choose m}\bmod p$ 的算法,其中 $p=\prod p_i^k$。 | + | $O\left(\sum p_i^k+\sum \log_2 n\log_{p_i^k}n\right)$ 计算 ${n\choose m}\bmod p$ 的算法,其中 $p=\prod p_i^k$。 |
==== 算法实现 ==== | ==== 算法实现 ==== | ||
+ | |||
+ | 考虑计算 ${n\choose m}\equiv x_i\pmod {p_i^k}$,然后中国剩余定理合并答案。 | ||
+ | |||
+ | $$ | ||
+ | {n\choose m}\equiv \frac {n!}{m!(n-m)!}\pmod {p^k} | ||
+ | $$ | ||
+ | |||
+ | 分子不一定与模数互质,不能直接使用逆元计算。考虑提取出 $n!$ 所有的 $p$ 因子,记为 $g(n)$,同时记 $f(n)=\frac {n!}{p^{g(n)}}$。代入上式,有 | ||
+ | |||
+ | $$ | ||
+ | \cfrac{\frac {n!}{p^{g(n)}}}{\frac {m!}{p^{g(m)}}\frac {(n-m)!}{p^{g(n-m)}}}p^{g(n)-g(m)-g(n-m)}\equiv \cfrac{f(n)}{f(m)f(n-m)}p^{g(n)-g(m)-g(n-m)}\pmod {p^k} | ||
+ | $$ | ||
+ | |||
+ | 于是有 $(f(n),p^k)=1$,可以参与逆元计算。设 $n=ap^k+b$。 | ||
+ | |||
+ | $$ | ||
+ | n!=\prod_{i=1,p\mid i}^ni \prod_{i=1,p\not\mid i}^ni=(\lfloor\frac np\rfloor)!p^{\lfloor\frac np\rfloor}(\prod_{i=1,p\not\mid i}^{p^k}i)^a\prod_{i=1,p\not\mid i}^b i | ||
+ | $$ | ||
+ | |||
+ | 于是有 | ||
+ | |||
+ | $$ | ||
+ | f(n)=f(\lfloor\frac np\rfloor)(\prod_{i=1,p\not\mid i}^{p^k}i)^a\prod_{i=1,p\not\mid i}^b i | ||
+ | $$ | ||
+ | $$ | ||
+ | g(n)=g(\lfloor\frac np\rfloor)+\lfloor\frac np\rfloor | ||
+ | $$ | ||
+ | |||
+ | 经过 $O\left(\sum p_i^k\right)$ 预处理后即可 $O\left(\sum \log_2 n\log_{p^k}n\right)$ 递归计算答案。 | ||
+ | |||
+ | ==== 代码模板 ==== | ||
+ | |||
+ | [[https://www.luogu.com.cn/problem/P4720|洛谷p4720]] | ||
+ | |||
+ | <hidden 查看代码> | ||
+ | <code cpp> | ||
+ | int quick_pow(int a,int b,int p){ | ||
+ | int ans=1; | ||
+ | while(b){ | ||
+ | if(b&1)ans=1LL*ans*a%p; | ||
+ | a=1LL*a*a%p; | ||
+ | b>>=1; | ||
+ | } | ||
+ | return ans; | ||
+ | } | ||
+ | void exgcd(LL a,LL b,LL &tx,LL &ty){ | ||
+ | if(b==0){ | ||
+ | tx=1,ty=0; | ||
+ | return; | ||
+ | } | ||
+ | exgcd(b,a%b,ty,tx); | ||
+ | ty-=a/b*tx; | ||
+ | } | ||
+ | int inv(int x,int mod){ | ||
+ | LL t1,t2; | ||
+ | exgcd(x,mod,t1,t2); | ||
+ | return t1%mod; | ||
+ | } | ||
+ | namespace Lucas{ | ||
+ | const int MAXM=30; | ||
+ | vector<int> frac[MAXM]; | ||
+ | int p_cnt,prime[MAXM],y[MAXM],pod[MAXM]; | ||
+ | void Init(int p,int mod){ | ||
+ | frac[p_cnt].resize(mod+1); | ||
+ | frac[p_cnt][0]=1; | ||
+ | _rep(i,1,mod){ | ||
+ | if(i%p)frac[p_cnt][i]=1LL*frac[p_cnt][i-1]*i%mod; | ||
+ | else frac[p_cnt][i]=frac[p_cnt][i-1]; | ||
+ | } | ||
+ | prime[p_cnt]=p; | ||
+ | pod[p_cnt]=mod/p*(p-1); | ||
+ | y[p_cnt++]=mod; | ||
+ | } | ||
+ | void init(int mod){ | ||
+ | p_cnt=0; | ||
+ | int t=mod; | ||
+ | for(int p=2;p*p<=t;p++){ | ||
+ | if(t%p==0){ | ||
+ | int mod2=1; | ||
+ | while(t%p==0)t/=p,mod2*=p; | ||
+ | Init(p,mod2); | ||
+ | } | ||
+ | } | ||
+ | if(t!=1)Init(t,t); | ||
+ | } | ||
+ | int f(LL n,int i){ | ||
+ | int ans=1; | ||
+ | while(n){ | ||
+ | ans=1LL*ans*quick_pow(frac[i][y[i]],(n/y[i])%pod[i],y[i])%y[i]*frac[i][n%y[i]]%y[i]; | ||
+ | n/=prime[i]; | ||
+ | } | ||
+ | return ans; | ||
+ | } | ||
+ | int g(LL n,int i){ | ||
+ | int ans=0; | ||
+ | while(n)ans+=(n/=prime[i]); | ||
+ | return ans; | ||
+ | } | ||
+ | int cal(LL n,LL m,int i){ | ||
+ | int t=1LL*f(n,i)*inv(1LL*f(m,i)*f(n-m,i)%y[i],y[i])%y[i]; | ||
+ | t=1LL*t*quick_pow(prime[i],g(n,i)-g(m,i)-g(n-m,i),y[i])%y[i]; | ||
+ | return t; | ||
+ | } | ||
+ | int C(LL n,LL m){ | ||
+ | int mod=1,ans=0; | ||
+ | _for(i,0,p_cnt)mod*=y[i]; | ||
+ | _for(i,0,p_cnt){ | ||
+ | int x=cal(n,m,i); | ||
+ | ans=(ans+1LL*x*(mod/y[i])%mod*inv(mod/y[i],y[i]))%mod; | ||
+ | } | ||
+ | return (ans+mod)%mod; | ||
+ | } | ||
+ | } | ||
+ | int main() | ||
+ | { | ||
+ | LL n=read_LL(),m=read_LL(); | ||
+ | int mod=read_int(); | ||
+ | Lucas::init(mod); | ||
+ | enter(Lucas::C(n,m)); | ||
+ | return 0; | ||
+ | } | ||
+ | </code> | ||
+ | </hidden> |