差别

这里会显示出您选择的修订版和当前版本之间的差别。

--- 2020-2021:teams:legal_string:jxm2001:数论_5 [2020/10/16 10:20]
jxm2001
+++ 2020-2021:teams:legal_string:jxm2001:数论_5 [2020/10/16 15:01] (当前版本)
jxm2001
@@ 行 5: / 行 5: @@
 ==== 算法简介 ====
-$O(p+\log_pn)$ 计算 ${n\choose m}\bmod p$ 的算法。
+$O(p+\log_pn)$ 计算 ${n\choose m}\bmod p$ 的算法，$p$ 为素数。
 ==== 算法实现 ====
@@ 行 67: / 行 67: @@
 ===== 拓展卢卡斯 =====
+==== 算法简介 ====
+$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>

CVBB ACM Team

用户工具

站点工具

差别

页面工具