文章包括：组合数取模、Miller-Rabin 算法、Pollard-Rho 算法、威尔逊定理、类欧几里得算法、原根、BSGS 与 exBSGS。
相关内容：「算法笔记」基础数论、「算法笔记」CRT 与 exCRT（中国剩余定理 与 扩展中国剩余定理）。
一、组合数取模
如何求 (inom{n}{m}mod p)？组合数的求法。
- 杨辉三角。(inom{n}{m}=inom{n-1}{m-1}+inom{n-1}{m})（考虑最后一个元素是否选）。
- 直接暴力计算。(inom{n}{m}=frac{n!}{m!(n-m)!})
- 预处理阶乘和阶乘的逆元。先计算分子 (n!mod p)，再计算分母 (m! (n-m)! mod p) 的逆元，乘起来得到 (inom{n}{m}mod p)。
- Lucas 定理。适用于 (n,m) 较大，(p) 较小的情况（要求 (p) 是质数）。(inom{n}{m}equiv inom{nmod p}{mmod p} imes inom{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 1pmod p)。
一个自然的想法就是看一下 (p) 是否满足 (a^{p−1}equiv 1pmod p)。但是这样判断 (p) 是否为质数存在 反例。不满足这个条件的一定不是质数，满足这个条件的也不一定是质数。所以可以先把不满足这个条件的排除，然后在这个基础上打打补丁。
二次探测定理：若 (p) 为质数，(a^2equiv 1pmod p)，那么 (aequiv ±1pmod p)。
证明： (a^2equiv 1pmod p)，所以 (a^2-1equiv 0pmod p)，则 ((a+1)(a-1)equiv 0pmod p)。由于 (p) 是质数，所以 (p) 一定能被 (a+1) 和 (a-1) 中的其中一个整除，也就是 (a+1equiv 0pmod p) 或 (a-1equiv 0pmod p)，即 (aequiv ±1pmod p)。
(a^{p−1}equiv 1pmod p)，若 (p-1) 为偶数，就可以把 (p-1) 表示成 (2^k imes t) 的形式。则 (a^{2^k imes t}equiv 1pmod p)。根据二次探测定理，可以推断：(a^{2^{k-1} imes t}equiv ±1pmod p,a^{2^{k-2} imes t}equiv ±1pmod pcdots)
于是我们就可以用二次探测定理与费马小定理结合判定质数。
2. 具体流程
-
把 (p-1) 表示成 (2^k imes t) 的形式。随机选择一个 (a)，求出 (a^t mod p)。
-
不断地对 (a) 进行平方操作，同时用二次探测定理判断，若 (a^2equiv 1pmod p)，而 (a otequiv ±1pmod p)，则 (p) 不是质数，否则重复此步骤。当自乘 (k) 次时，此时 (a=p-1)，停止操作。
-
此时 (a=p-1)。用费马小定理判一下，若不满足 (a^{p−1}equiv 1pmod 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 ",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)mod n) 来生成随机数。其中 (c) 是一个随机的常数。先随机取一个数 (x_1)，然后对于所有 (1<ileq n)，(x_i=f(x_{i-1}))，即 (x_i=(x_{i-1}^2+c)mod n)。
2. Floyd 判环
由于生成的数列的某一项的值仅由前一项决定，且每一项可能的取值是有限的，因此 (x_i) 最终必定会形成环，而在形成环时就不能再继续操作了。举个栗子，当(n=50,c=2,x_1=1) 时，(f(x)) 生成的数据为 (1,3,11,23,31,11,23,31cdots)，数据在 (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(acmod b,b)=gcd(a,b))。其中 (c) 为正整数。
我们每过一段时间将这些差值进行 (gcd) 运算，令 (s=prod |x_0-x_j|mod 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 ",ans); //输出最大质因子 } return 0; }
四、威尔逊定理
威尔逊定理：((p-1)!equiv −1pmod p) 对任意质数 (p) 均成立，且当 (p) 不是质数时为 (0)，只在 (p=4) 时例外。
证明：
1. 若 (p) 为质数：因为 (p) 是质数，所以 (1sim p-1) 在模 (p) 意义下有逆元。
假如 (a^{-1}=b)。
-
当 (a=b) 时：(a=b) 说明 (a^2equiv 1pmod p)。在模数是质数的情况下，(a^2equiv 1pmod p) 意味着 (aequiv ±1pmod p)。
-
当 (a eq b) 时：在阶乘中让 (a) 和 (b) 配对，于是这两个数就消掉了。
两两配对，剩下的两个数为 (1) 和 (p-1)。因此，((p-1)!equiv -1pmod p)。
2. 若 (p) 不是质数：因为 (p) 不是质数，所以 (p) 能被分解成两个数的乘积，即 (p=ab)。
-
若 (p) 不是完全平方数，则存在 (p=ab,a eq b)。(a) 和 (b) 在 ((p-1)!) 中都出现过，因此，((p-1)!equiv 0pmod p)。
-
若 (p) 为完全平方数，(p=a^2)。当 (p>4) 时，(a<p,2a<p)，(a) 和 (2a) 在 ((p-1)!) 中都出现过，因此 ((p-1)!equiv 0pmod p)。举个栗子，当 (p=9) 时，((p-1)!) 中含有两个 (3)（一个 (3)，一个 (6)），所以 ((p-1)!equiv 0pmod p)。而在 (p=4) 时，((p-1)!) 中并没有含有两个 (2)，所以 (p=4) 时例外。
五、类欧几里得算法
后来知道可以从几何意义推导，并且几何意义的更容易理解，以前写的 从代数推导的 就没什么用了，所以不公开了 qwq。
从几何意义推导的看 黄队写的。
六、原根
1. 阶
定义：设 (gcd(a,m)=1)，则 (a) 在模 (m) 下的 阶 被定义为最小的使得 (a^xequiv 1pmod m) 的 (x)，记为 ( ext{ord}_m a)。
一些结论：
-
结论 1：( ext{ord}_m xleq varphi(m))，且 ( ext{ord}_m x) 一定存在。
证明：根据欧拉定理，当 (x=varphi(m)) 时 (a^xequiv 1pmod m) 成立，得证。
-
结论 2：(a^xequiv 1pmod mLeftrightarrow ext{ord}_mamid x)。
证明：设 ( ext{ord}_ma) 为 (y)。首先，(a^{ky}equiv 1pmod m)（(k) 为正整数）。其次，假设存在 (y mid x) 使得 (a^xequiv 1pmod m)，则有 (x=ky+r(0<r<y))，且有 (a^x=a^{ky+r}equiv 1pmod m) 成立。则 (a^requiv 1pmod m)。又 (r<y)，与已知矛盾。故结论成立。
-
结论 3：根据结论 2 有 ( ext{ord}_m amid varphi(m))。
-
结论 4： 设 ( ext{ord}_m a=x)，则 (a^0,a^1,cdots,a^{x-1}) 对模 (m) 两两不同余。
证明：反证法。假设存在 (0leq i<jleq x-1)，使得 (a^iequiv a^jpmod m)。由 (gcd(a,m)=1) 得，(frac{a^j}{a^i}=a^{j-i}equiv 1pmod m)。而 (0<j-i<x)，这与 ( ext{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) 的 (1simvarphi(m)) 次幂，复杂度为 (mathcal{O}(varphi(m)))。无法承受。
-
因为 ( ext{ord}_m amid varphi(m))，所以，要找 ( ext{ord}_ma)，只需枚举 (varphi(m)) 的约数。更进一步，若 (a^xequiv 1pmod m)，则 (a^{kx}equiv 1pmod m)（(k) 为正整数）。于是只需枚举 (varphi(m)) 的每个质因子 (p_i)，看 (frac{varphi(m)}{p_i}) 是否满足条件。如果存在 (frac{varphi(m)}{p_i}) 满足条件，说明 ( ext{ord}_ma<varphi(m))，(a) 必不是 (m) 的原根。否则，可以说明 (a) 必是 (m) 的原根。
定理 3：设 (g) 为 (m) 的一个原根，则集合 (S={g^smid 1leq sleq 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}} otequiv 1pmod m)。那么就有 ({(g^s)}^{frac{varphi(m)}{p}}equiv (g^{varphi(m)})^{frac{s}{p}}pmod m)。假设存在 (gcd(s,varphi(m)) eq 1) 使得 (g^s) 为原根，则存在一个 (p) 使得 (pmid s)。此时有 ((g^{varphi(m)})^{frac{s}{p}}equiv g^{kvarphi(m)}pmod m)（(k) 为正整数）。根据欧拉定理，(g^{varphi(m)}equiv 1pmod m)，得 ({(g^s)}^{frac{varphi(m)}{p}}equiv 1pmod m)，产生矛盾。故命题成立。
若 (s>varphi(m))，根据扩展欧拉定理，则 (g^sequiv g^{smod varphi(m)}pmod m)，转化为 (1leq sleq 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 ",tot); for(int i=1;i<=tot/d;i++) printf("%lld%c",ans[i*d],i==tot/d?' ':' '); if(tot/d==0) puts(""); } return 0; }
七、BSGS
详见 「算法笔记」BSGS。