miller_rabin与pollard_rho

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 一些性质

1. \(a\) 选择不同，序列 \(f\) 也不同。

2. \(f\) 一开始不循环但最后一定会有循环节。

   证明：这是因为在 \(\mod n\) 的意义下，\(f\) 只有 \(n\) 中可能的取值。根据鸽巢原理，不难证明。

取模的函数一般会有循环节。

3. 若 \(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);
}

标签：return,gcd,1.2,v1,int,miller,pollard,v2,rho
来源： https://www.cnblogs.com/TianMeng-hyl/p/14811811.html