1.1 miller_rabin
一个检查素数的算法。
是一个概率算法，并不能保证绝对。
1.1.1一些定理
- 如果 \(n\) 为素数，取 \(n-1=d\times 2^r\)，\(\forall a<n,a\in Z^+\) ，有 \(a^d\equiv 1 \mod n\) 或者 \(\exist 0\leq i<r,s.t.a^{d\times 2^i}\equiv -1\mod n\) （费马小定理推论）
1.1.2 思路
我们找任意个数小于 \(n\) ，对 \(n\) 进行检验，看上面两个性质是否都满足，如果都不满足，则 \(n\) 不是素数；否则 \(n\) 有很大概率是一个素数。
我们要选若干个 \(a<n\) ，一般我们的算法是若干个小质数，例如 \(2,3,5,11,17,19,23,37\)，这些数在 \(longlong\) 范围内已经可以判断出全部的素数了。
1.1.3 代码：
const int INF=0x3f3f3f3f;
inline int read(){
int x=0,f=1;
char ch=getchar();
while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
return x*f;
}
inline void ksm(int a,int b,int mod){
int res=1;
while(b){
if(b&1) (res*=a)%=mod;
a=a*a%mod;
b>>=1;
}
return res;
}
inline bool miller_rabin(int n,int a){
int d=n-1,r=0;
while(!(d&1)) d>>=1,r++;
int x=ksm(a,d,n);
if(x==1) return 1;
for(int i=0;i<=r-1;i++){
if(x==n-1) return 1;
x=1ll*x*x;
}
return 0;
}
const int prime_list[]={0,2,3,5,11,27,19,23,37};
inline bool is_prime(int n){
if(n<2) return 0;
for(int i=1;i<=8;i++){
if(n==prime_list[i]) return 1;
if(!miller_rabin(n,prime_list[i])) return 0;
}
return 1;
}
1.2 质因数分解Pollard_Rho
对 \(n\) 作质因数分解。
1.2.1 思路1
把 \(n\) 分成是否质数两类，目前要解决的问题是要找到一个 \(a,s.t.a|n\) 。
为了降低复杂度，考虑随机化，我们不断随机 \(a\)，直到 \(a|n\) 。
目前要解决的问题是随机的次数太多了，导致复杂度很高。
考虑更改条件，我们只需要找到 \(a,s.t.\gcd(a,b)\not=1\) ，这样就能找到一个 \(n\) 的质因子。
我们需要期望随机 \(\dfrac{n}{n-\phi(n)}\) 次。
虽然复杂度降了一点，但是仍然不能接受。
1.2.2思路2
我们考虑随便选 \(k\) 个数：\(a_1,a_2,...a_k\) 然后两两作差：\(b_1,b_2,...b_{k^2}\)，然后考虑这样的一个问题：我们给定一个数 \(c\) ，所有的 \(a_i,b_i\) 在 \(\mod c\) 意义下 ，那么 \(\nexists i,s.t.b_i=\) ，的概率0是多少？
我们大约估算一下，其实 \(b_i\) 也可以看做是从 \(0\) 到 \(c-1\) 中的随机数，所以概率是 \((\dfrac{c-1}{c})^{k^2}\)，考虑到 \((\dfrac{x-1}{x})^x \approx \dfrac{1}{e}\)，所以当 \(k\) 取 \(10\sqrt c\) 时才可以吧概率降到一个比较小的范围。
我们回到以前的问题，通过把 \(c\) 换成 \(n\) 的一个因子，我们发现复杂度仍然是不可接受的，尤其是对于一些大质数的乘积。
接下来我们考虑怎样随机选这 \(k\) 个数。
1.2.3 伪随机函数
\(f_i=(f_{i-1}^2+a)\mod n\)
1.2.3.1 一些性质
-
\(a\) 选择不同，序列 \(f\) 也不同。
-
\(f\) 一开始不循环但最后一定会有循环节。
证明：这是因为在 \(\mod n\) 的意义下，\(f\) 只有 \(n\) 中可能的取值。根据鸽巢原理，不难证明。
取模的函数一般会有循环节。
-
若 \(g|f_i-f_j\) 则 \(g|f_{i+1}-f_{j+1}\)
证明：把 \(f_{i+1},f_{j+1}\) 表示出来用平方差即可证明。
我们考虑用这个方法生成 \(k\) 个数，根据性质 \(3\) 我们不用检验 \(k^2\) 个差，两个下标之差为定值，则有共同因子。
1.2.4 双指针法
第一个指针指在第一个位置，第二个指针指在第二个位置，第一个指针跳一步，第二个指针跳两步，若后面这个指针赶上了前面这个指针，则我们认为已经遍历完了循环节。
不难发现，我们实际上是在遍历下标之差。当我们发现两个下标所代表数相同时，通过画图我们可以知道，这等价于两个下标相同了，因为他们已经到达了同一位置，可以认为下标相同，这个时候我们没有必要在进行下去。
1.2.5 代码1
我们同时用 miller_rabin 判素数，因为如果是一个素数，用 pollard_rho 就会在 \(29\) 行死循环。
inline int f(int a,int n,int c){
return (1ll a*a+c)%n;
}
inline int gcd(int a,int b){
return b==0?a:gcd(b,a%b);
}
inline void pollard_rho(int n){
int c=rand()%(n-1)+1;
int v1=c;
int v2=f(v1,n,c);
while(v1!=v2){
int g=gcd(n,abs(v1-v2));
if(g!=1) return g;
v1=f(v1,n,c);
v2=f(v2,n,c);
v2=f(v2,n,c);
}
return n;
}
inline void factorize(int n){
if(miller_rabin(n)){
z[++cnt]=n;
return;
}
int g=n;
while(g==n) g=pollard_rho(n);
factorize(n/g);
factorize(g);
}
1.2.6 玄学优化
我们发现上面这个代码的时间复杂度取决于循环节的长度，在实际实践中发现实际上并不长。
我们做了许多次没有用的 \(gcd\) ，考虑优化这个过程。
我们算许多次 \(gcd\) 相当于我们算一个 \(\gcd((f_j-f_i)\times (f_{j+1}-f_i)\times ...,n)\)，我们同时采取一个倍增的方法去写，每次处理一个倍增区间中的数，把差乘起来在取 \(gcd\)，这是在许多人写的时候总结出来的可以这么写，没有必要追究为什么要这么写。
实际在写的时候我们通常每 \(128\) 作一次 \(gcd\) ，因为万一这个倍增区间很长，而答案可能早就出现了，这当然很亏，所以我们有上述做法。
我们每次从 \(v1\) 开始，取一下倍增区间内每一个 \(v2\)，把差连乘并取 \(gcd\)，
如果没有找到，\(v2\) 现在已经是 \(v1\) 进行 \(t\) 步之后了，我们更换起点，继续上述过程。
实际应用中，下面这份代码确实比上面那个快。
需要注意的一点是，我们不需要在最后 return n ，原因是如果 \(v1\) 和 \(v2\) 是一样的话，\(mul=0\) ，而 \(\gcd(0,n)=n\)，这是直接会 return n 。
1.2.7代码2
inline int gcd(int a,int b){
return b==0?a:gcd(b,a%b);
}
inline void pollard_rho(int n){
int c=rand()%(n-1)+1;
int v1=0;//c
for(int s=1,t=2;s<<=1,t<<=1){//枚举每次处理数的个数
int v2=v1,mul=1;
for(int a=s,step=1;a<t;a++,step++){
v2=f(v2,n,c);
mul=1ll*mul*abs(v1-v2)%n;
if(step==127){
step=0;
int v=gcd(mul,n);
if(v>1) return v;
}
}
int v=gcd(mul,n);
if(v>1) return v;
v1=v2;
}
}
inline void factorize(int n){
if(miller_rabin(n)){
z[++cnt]=n;
return;
}
int g=n;
while(g==n) g=pollard_rho(n);
factorize(n/g);
factorize(g);
}