「算法笔记」基础数论 2
文章包括：组合数取模、Miller-Rabin 算法、Pollard-Rho 算法、威尔逊定理、类欧几里得算法、原根、BSGS 与 exBSGS、莫比乌斯反演、杜教筛。
前置内容：「算法笔记」基础数论、「算法笔记」CRT 与 exCRT。
upd：写这篇文章时处于什么也不懂的状态，后来重写了 MR 和 PR，想看可找我 QwQ。
一、组合数取模
如何求 \(\binom{n}{m}\bmod p\)？组合数的求法。
- 杨辉三角。\(\binom{n}{m}=\binom{n-1}{m-1}+\binom{n-1}{m}\)（考虑最后一个元素是否选）。
- 直接暴力计算。\(\binom{n}{m}=\frac{n!}{m!(n-m)!}\)
- 预处理阶乘和阶乘的逆元。先计算分子 \(n!\bmod p\)，再计算分母 \(m! (n-m)! \bmod p\) 的逆元，乘起来得到 \(\binom{n}{m}\bmod p\)。
- Lucas 定理。适用于 \(n,m\) 较大，\(p\) 较小的情况（要求 \(p\) 是质数）。\(\binom{n}{m}\equiv \binom{n\bmod p}{m\bmod p}\times \binom{n/p}{m/p} \pmod p\)。也就是把 \(n,m\) 表示成 \(p\) 进制，对 \(p\) 进制下的每一位分别计算组合数，最后再乘起来。
- \(p\) 不是质数：扩展 Lucas 定理。代码。
二、Miller-Rabin 算法
Miller-Rabin 算法可以 \(\mathcal{O}(\log_2 n)\) 判断一个数是否是质数。
1. 基本思想
费马小定理：若 \(p\) 为质数，且 \(\gcd(a,p)=1\)，则 \(a^{p−1}\equiv 1\pmod p\)。
一个自然的想法就是看一下 \(p\) 是否满足 \(a^{p−1}\equiv 1\pmod p\)。但是这样判断 \(p\) 是否为质数存在 反例。不满足这个条件的一定不是质数，满足这个条件的也不一定是质数。所以可以先把不满足这个条件的排除，然后在这个基础上打打补丁。
二次探测定理：若 \(p\) 为质数，\(a^2\equiv 1\pmod p\)，那么 \(a\equiv ±1\pmod p\)。
证明： \(a^2\equiv 1\pmod p\)，所以 \(a^2-1\equiv 0\pmod p\)，则 \((a+1)(a-1)\equiv 0\pmod p\)。由于 \(p\) 是质数，所以 \(p\) 一定能被 \(a+1\) 和 \(a-1\) 中的其中一个整除，也就是 \(a+1\equiv 0\pmod p\) 或 \(a-1\equiv 0\pmod p\)，即 \(a\equiv ±1\pmod p\)。
\(a^{p−1}\equiv 1\pmod p\)，若 \(p-1\) 为偶数，就可以把 \(p-1\) 表示成 \(2^k\times t\) 的形式。则 \(a^{2^k\times t}\equiv 1\pmod p\)。根据二次探测定理，可以推断：\(a^{2^{k-1}\ \times t}\equiv ±1\pmod p,a^{2^{k-2}\ \times t}\equiv ±1\pmod p\cdots\)
于是我们就可以用二次探测定理与费马小定理结合判定质数。
2. 具体流程
-
把 \(p-1\) 表示成 \(2^k\times t\) 的形式。随机选择一个 \(a\)，求出 \(a^t \bmod p\)。
-
不断地对 \(a\) 进行平方操作，同时用二次探测定理判断，若 \(a^2\equiv 1\pmod p\)，而 \(a\not\equiv ±1\pmod p\)，则 \(p\) 不是质数，否则重复此步骤。当自乘 \(k\) 次时，此时 \(a=p-1\)，停止操作。
-
此时 \(a=p-1\)。用费马小定理判一下，若不满足 \(a^{p−1}\equiv 1\pmod p\)，则 \(p\) 不是质数。
-
通过了所有的测试，则返回 True。
正确性：每一次通过测试的数不是质数的概率为 \(\frac{1}{4}\)，则测试 \(t\) 次，错误的概率为 \(\frac{1}{4^t}\)。所以 \(a\) 取 \(2,3,5,7,11,13,17,19,23,29\) 进行判断错误率趋近于 \(0\)。对于不同范围数的 Miller-Rabin 检验可以取不同的 \(a\)。
//HDU 2138 #include<bits/stdc++.h> #define int long long using namespace std; int n,x,a[15]={2,3,5,7,11,13,17,19,23,29},ans; int mul(int x,int n,int mod){ int ans=mod!=1; for(x%=mod;n;n>>=1,x=x*x%mod) if(n&1) ans=ans*x%mod; return ans; } bool solve(int p){ if(p==1) return 0; int t=p-1,k=0; while(t%2==0) k++,t/=2; //将 p-1 分解成 (2^k)*t 的形式 for(int i=0;i<10;i++){ if(p==a[i]) return 1; int x=mul(a[i],t,p),nxt=x; //计算出 a^t%p if(x==1) continue; //小优化，如果此时结果为 1，那么无论如何自乘也为 1 for(int j=1;j<=k;j++){ nxt=x*x%p; //不断自乘 if(nxt==1&&x!=1&&x!=p-1) return 0; //若 a^2%p=1，而 a!=1 且 a!=p-1，则 p 不是质数 x=nxt; //继续 } if(x!=1) return 0; //此时 a=p-1。用费马小定理判断 } return 1; } signed main(){ while(~scanf("%lld",&n)){ ans=0; for(int i=1;i<=n;i++){ scanf("%lld",&x); if(solve(x)) ans++; } printf("%lld\n",ans); } return 0; }
三、Pollard-Rho 算法
Pollard-Rho 算法可以 \(\mathcal{O}(\sqrt[4]{n} \log_2 n)\) 找到一个数 \(n\) 的一个因子。
1. 优化随机算法
首先，一种想法是通过随机的方法，猜测一个数是否为 \(n\) 的因数。
考虑优化这个随机算法。
由于 \(\gcd(x,n)\mid n\)，所以只要选适当的 \(x\) 满足 \(1<\gcd(x,n)<n\) ，就可以求得一个约数 \(\gcd(x,n)\)。满足这样的 \(x\) 不少，\(x\) 有若干个质因子，每个质因子及其倍数都是可行的。
根据生日悖论，取多组数据作差能够优化取数的精确性。
不妨取一组数 \(x_1,x_2,\cdots,x_n\)，若有 \(1<\gcd(|x_i-x_j|,n)<n\)，那么我们就找到了 \(n\) 的一个因子。
生成 \(x_i\) 的方法：通过 \(f(x)=(x^2+c)\bmod n\) 来生成随机数。其中 \(c\) 是一个随机的常数。先随机取一个数 \(x_1\)，然后对于所有 \(1<i\leq n\)，\(x_i=f(x_{i-1})\)，即 \(x_i=(x_{i-1}^2+c)\bmod n\)。
2. Floyd 判环
由于生成的数列的某一项的值仅由前一项决定，且每一项可能的取值是有限的，因此 \(x_i\) 最终必定会形成环，而在形成环时就不能再继续操作了。举个栗子，当 \(n=50,c=2,x_1=1\) 时，\(f(x)\) 生成的数据为 \(1,3,11,23,31,11,23,31\cdots\)，数据在 \(3\) 以后都在 \(11,23,31\) 之间循环。
假设两个人在赛跑，A 的速度快，B 的速度慢，经过一定时间后，A 一定会和 B 相遇，且相遇时 A 跑过的总距离减去 B 跑过的总距离一定是圈长的 \(n\) 倍。
设 \(a=f(1),b=f(f(1))\)，每一次更新 \(a=f(a),b=f(f(b))\)，只要检查在更新过程中 \(a\) 和 \(b\) 是否相等，如果相等了，那么就出现了环。
我们每次对 \(d=\gcd(|x_i-x_j|,n)\)，判断 \(d\) 是否满足 \(1<d<n\)，若满足则可直接返回 \(d\)。形成环时直接返回 \(n\) 本身，并且在后续操作里调整随机常数 \(c\)，重新分解。
//基于 Floyd 判环的 Pollard-Rho 算法 int f(int x,int c,int n){ return (x*x+c)%n; } int find(int n){ int c=rand()%(n-1)+1,a=f(0,c,n),b=f(f(0,c,n),c,n); while(a!=b){ int d=__gcd(abs(a-b),n); if(d>1) return d; a=f(a,c,n),b=f(f(b,c,n),c,n); } return n; //返回 n，在后续操作里调整随机常数 c }
3. 倍增优化
考虑到对于每个 \(x_i\) 都要计算一遍 \(\gcd\)，显然很慢。可以通过乘法累积来减少求 \(\gcd\) 的次数。
如果 \(1<\gcd(a,b)\)，则有 \(1<\gcd(ac,b)\)，并且有 \(1<\gcd(ac\bmod b,b)=\gcd(a,b)\)。其中 \(c\) 为正整数。
我们每过一段时间将这些差值进行 \(\gcd\) 运算，令 \(s=\prod |x_0-x_j|\bmod n\)，如果某一时刻得到 \(s=0\) 那么表示分解失败，退出并返回 \(n\) 本身。每隔 \(2^k\) 个数，计算是否满足 \(1<\gcd(s,n)<n\)。
Pollard-Rho 算法模板：（也可以不用 __int128 用快速乘。代码。）
//Luogu P4718 #include<bits/stdc++.h> #define int long long #define LLL __int128 using namespace std; int t,n,x,a[15]={2,3,5,7,11,13,17,19,23,29},ans; int mul(int x,int n,int mod){ int ans=mod!=1; for(x%=mod;n;n>>=1,x=(LLL)x*x%mod) if(n&1) ans=(LLL)ans*x%mod; return ans; } bool solve(int p){ //Miller-Rabin 算法判定质数 if(p==1) return 0; int t=p-1,k=0; while(t%2==0) k++,t/=2; for(int i=0;i<10;i++){ if(p==a[i]) return 1; int x=mul(a[i],t,p),nxt=x; if(x==1) continue; for(int j=1;j<=k;j++){ nxt=(LLL)x*x%p; if(nxt==1&&x!=1&&x!=p-1) return 0; x=nxt; } if(x!=1) return 0; } return 1; } int f(int x,int c,int n){ return ((LLL)x*x+c)%n; } int find(int n){ //求 n 的一个因数 if(n%2==0) return 2; if(solve(n)) return n; int x=0,y=0,c=rand()%(n-1)+1,s=1,d; for(int k=1;;k<<=1,x=y,s=1){ for(int i=1;i<=k;i++){ y=f(y,c,n),s=(LLL)s*abs(y-x)%n; if(i%127==0){d=__gcd(s,n);if(d>1) return d;} } d=__gcd(s,n);if(d>1) return d; } } void work(int x){ if(x<=ans||x<2) return ; if(solve(x)){ans=max(ans,x);return ;} //用 Miller-Rabin 算法判断是否出现质因子，并更新答案 int p=x; while(p>=x) p=find(x); //用 Pollard-Rho 算法找一个因子 p while(x%p==0) x/=p; //将 x 除去因子 p work(x),work(p); //分别递归分解 x 和 p } signed main(){ scanf("%lld",&t); while(t--){ srand(time(0)); scanf("%lld",&n),ans=0; if(solve(n)) puts("Prime"); else work(n),printf("%lld\n",ans); //输出最大质因子 } return 0; }
四、威尔逊定理
威尔逊定理：\((p-1)!\equiv −1\pmod p\) 对任意质数 \(p\) 均成立，且当 \(p\) 不是质数时为 \(0\)，只在 \(p=4\) 时例外。
证明：
1. 若 \(p\) 为质数：因为 \(p\) 是质数，所以 \(1\sim p-1\) 在模 \(p\) 意义下有逆元。
假如 \(a^{-1}=b\)。
-
当 \(a=b\) 时：\(a=b\) 说明 \(a^2\equiv 1\pmod p\)。在模数是质数的情况下，\(a^2\equiv 1\pmod p\) 意味着 \(a\equiv ±1\pmod p\)。
-
当 \(a\neq b\) 时：在阶乘中让 \(a\) 和 \(b\) 配对，于是这两个数就消掉了。
两两配对，剩下的两个数为 \(1\) 和 \(p-1\)。因此，\((p-1)!\equiv -1\pmod p\)。
2. 若 \(p\) 不是质数：因为 \(p\) 不是质数，所以 \(p\) 能被分解成两个数的乘积，即 \(p=ab\)。
-
若 \(p\) 不是完全平方数，则存在 \(p=ab,a\neq b\)。\(a\) 和 \(b\) 在 \((p-1)!\) 中都出现过，因此，\((p-1)!\equiv 0\pmod p\)。
-
若 \(p\) 为完全平方数，\(p=a^2\)。当 \(p>4\) 时，\(a<p,2a<p\)，\(a\) 和 \(2a\) 在 \((p-1)!\) 中都出现过，因此 \((p-1)!\equiv 0\pmod p\)。举个栗子，当 \(p=9\) 时，\((p-1)!\) 中含有两个 \(3\)（一个 \(3\)，一个 \(6\)），所以 \((p-1)!\equiv 0\pmod p\)。而在 \(p=4\) 时，\((p-1)!\) 中并没有含有两个 \(2\)，所以 \(p=4\) 时例外。
五、类欧几里得算法
后来知道可以从几何意义推导，并且几何意义的更容易理解，以前写的 从代数推导的 就没什么用了，所以不公开了 qwq。
从几何意义推导的看 黄队写的。
六、原根
1. 阶
定义：设 \(\gcd(a,m)=1\)，则 \(a\) 在模 \(m\) 下的 阶 被定义为最小的使得 \(a^x\equiv 1\pmod m\) 的 \(x\)，记为 \(\text{ord}_m a\)。
一些结论：
-
结论 1：\(\text{ord}_m x\leq \varphi(m)\)，且 \(\text{ord}_m x\) 一定存在。
证明：根据欧拉定理，当 \(x=\varphi(m)\) 时 \(a^x\equiv 1\pmod m\) 成立，得证。
-
结论 2：\(a^x\equiv 1\pmod m\Leftrightarrow\text{ord}_ma\mid x\)。
证明：设 \(\text{ord}_ma\) 为 \(y\)。首先，\(a^{ky}\equiv 1\pmod m\)（\(k\) 为正整数）。其次，假设存在 \(y\nmid x\) 使得 \(a^x\equiv 1\pmod m\)，则有 \(x=ky+r(0<r<y)\)，且有 \(a^x=a^{ky+r}\equiv 1\pmod m\) 成立。则 \(a^r\equiv 1\pmod m\)。又 \(r<y\)，与已知矛盾。故结论成立。
-
结论 3：根据结论 2 有 \(\text{ord}_m a\mid \varphi(m)\)。
-
结论 4： 设 \(\text{ord}_m a=x\)，则 \(a^0,a^1,\cdots,a^{x-1}\) 对模 \(m\) 两两不同余。
证明：反证法。假设存在 \(0\leq i<j\leq x-1\)，使得 \(a^i\equiv a^j\pmod m\)。由 \(\gcd(a,m)=1\) 得，\(\frac{a^j}{a^i}=a^{j-i}\equiv 1\pmod m\)。而 \(0<j-i<x\)，这与 \(\text{ord}_m a=x\) 的定义矛盾。所以定理成立。
2. 原根
定义：\(\gcd(g,m)=1\)，若 \(g\) 模 \(m\) 的阶为 \(\varphi(m)\)，则称 \(g\) 为 \(m\) 的 原根。
定理 1：根据阶的结论 4 可得，若 \(g\) 是模 \(m\) 的原根，则 \(g^0,g^1,\cdots,g^{\varphi(m)-1}\) 构成模 \(m\) 的简化剩余系。
定理 2：只有 \(1,2,4,p^a,2p^a\) 存在原根，其中 \(p\) 是奇素数，\(a\) 是任意正整数。证明略。
检验原根：判断一个数 \(a\) 是否为 \(m\) 的原根。
-
朴素算法：枚举 \(a\) 的 \(1\sim\varphi(m)\) 次幂，复杂度为 \(\mathcal{O}(\varphi(m))\)。无法承受。
-
因为 \(\text{ord}_m a\mid \varphi(m)\)，所以，要找 \(\text{ord}_ma\)，只需枚举 \(\varphi(m)\) 的约数。更进一步，若 \(a^x\equiv 1\pmod m\)，则 \(a^{kx}\equiv 1\pmod m\)（\(k\) 为正整数）。于是只需枚举 \(\varphi(m)\) 的每个质因子 \(p_i\)，看 \(\frac{\varphi(m)}{p_i}\) 是否满足条件。如果存在 \(\frac{\varphi(m)}{p_i}\) 满足条件，说明 \(\text{ord}_ma<\varphi(m)\)，\(a\) 必不是 \(m\) 的原根。否则，可以说明 \(a\) 必是 \(m\) 的原根。
定理 3：设 \(g\) 为 \(m\) 的一个原根，则集合 \(S=\{g^s\mid 1\leq s\leq \varphi(m),\gcd(s,\varphi(m))=1\}\) 给出 \(m\) 的全部原根。因此，模 \(m\) 意义下的原根有 \(\varphi(\varphi(m))\) 个。根据阶的结论 4，这 \(\varphi(\varphi(m))\) 个原根模 \(m\) 两两不同余。
略证：由原根的检验得，若 \(g\) 为 \(m\) 的原根，则对于 \(\varphi(m)\) 的质因子 \(p\)，有 \(g^{\frac{\varphi(m)}{p}}\not\equiv 1\pmod m\)。那么就有 \({(g^s)}^{\frac{\varphi(m)}{p}}\equiv (g^{\varphi(m)})^{\frac{s}{p}}\pmod m\)。假设存在 \(\gcd(s,\varphi(m))\neq 1\) 使得 \(g^s\) 为原根，则存在一个 \(p\) 使得 \(p\mid s\)。此时有 \((g^{\varphi(m)})^{\frac{s}{p}}\equiv g^{k\varphi(m)}\pmod m\)（\(k\) 为正整数）。根据欧拉定理，\(g^{\varphi(m)}\equiv 1\pmod m\)，得 \({(g^s)}^{\frac{\varphi(m)}{p}}\equiv 1\pmod m\)，产生矛盾。故命题成立。
若 \(s>\varphi(m)\)，根据扩展欧拉定理，则 \(g^s\equiv g^{s\bmod \varphi(m)}\pmod m\)，转化为 \(1\leq s\leq \varphi(m)\) 的问题。
求原根：
- 求一个原根：原根密度较大，大约是 \(n^{0.25}\)。那么可从 \(2\) 开始枚举 \(g\)，并进行检验，可找到最小的原根。也可以直接随机一个数并进行检验。复杂度均约为 \(\mathcal{O}(n^{0.25})\)。
- 求所有原根：首先找到 \(m\) 最小的原根 \(g\)，则根据定理 3，\(m\) 的每一个原根都形如 \(g^k\) 的形式，且满足 \(\gcd(k,\varphi(m))=1\)。
//Luogu P6091 #include<bits/stdc++.h> #define int long long using namespace std; const int N=2e6+5; int t,n=1e6,d,x,cnt,p[N],phi[N],g,tot,tmp,ans[N]; bool vis[N],f[N]; vector<int>v; int mul(int x,int n,int mod){ int ans=mod!=1; for(x%=mod;n;n>>=1,x=x*x%mod) if(n&1) ans=ans*x%mod; return ans; } bool solve(int a){ //判断 a 是否为 n 的原根 if(mul(a,x,n)!=1) return 0; for(int i=0;i<(int)v.size();i++){ int p=v[i]; if(mul(a,x/p,n)==1) return 0; //看 φ(m)/p 是否满足条件 } return 1; } signed main(){ scanf("%lld",&t),vis[0]=vis[1]=1,phi[1]=1; for(int i=2;i<=n;i++){ if(!vis[i]) p[++cnt]=i,phi[i]=i-1; for(int j=1;j<=cnt&&i*p[j]<=n;j++){ vis[i*p[j]]=1; if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j];break;} phi[i*p[j]]=phi[i]*phi[p[j]]; } } f[1]=f[2]=f[4]=1; //f[i]=1 表示 i 存在原根 for(int i=2;i<=cnt;i++) for(int j=p[i];j<=n;j*=p[i]) f[j]=f[2*j]=1; while(t--){ scanf("%lld%lld",&n,&d); if(!f[n]){puts("0"),puts("");continue;} vector<int>().swap(v),x=phi[n],g=tot=0,tmp=1; for(int i=2;i<=sqrt(x);i++){ //求出 φ(m) 的质因子 if(x%i) continue; v.push_back(i); while(x%i==0) x/=i; } if(x>1) v.push_back(x); x=phi[n]; for(int i=1;i<=n;i++) if(solve(i)){g=i;break;} //找到 n 最小的原根 for(int i=1;i<=x;i++){ tmp=tmp*g%n; if(__gcd(i,x)==1) ans[++tot]=tmp; } sort(ans+1,ans+1+tot),printf("%lld\n",tot); for(int i=1;i<=tot/d;i++) printf("%lld%c",ans[i*d],i==tot/d?'\n':' '); if(tot/d==0) puts(""); } return 0; }
七、BSGS
八、莫比乌斯反演
详见 「算法笔记」莫比乌斯反演。
九、杜教筛
详见 「算法笔记」杜教筛。
