首页 > 其他分享> 文章详细

miller_rabin与pollard_rho

2021-05-26 09:05:00 阅读：168 来源： 互联网

标签：return gcd 1.2 v1 int miller pollard v2 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 一些性质

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

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

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9