Miller_Rabin算法&&Pollard-Rho算法

前述：近期在uva上刷到一个题，prime test，题意非常简单，就是让判断一个2^54以内的数是否为素数，年轻的我满脑子都是相关的暴力加优化，做之前我还去百度搜了一下判断一个数是否为素数的最佳方法，博主给的是素数的线性筛法，很自信的敲完代码，然后很自信的tle。做不了了，搜了一下题解，原来并不是一道简单题，涉及到本文提及的两个算法，Miller_Rabin算法&&Pollard-Rho算法，当然，一个是判断一个数是否为素数的，另外一个则是快速求一个数的质分解的，但是并非不搭边，前者是后者的基础，二者都是随机数算法，但前者失败的概率基本相当于中彩票，后者却十分高，但是相对的，效率也极快了，以前觉得随机数算法十分的神秘，现在发现并不是真的完全靠运气，而是用理论把失败的概率降到最低，很有趣的一类算法，当然这里原理的basic都不太精确，算是我对搜到的所有资料的总结吧，我也算是那种比较喜欢刨根问底的人了。

  Miller_Rabin算法：

  算法原理：

  这是一个基于费马小定理的算法，定理内容：若n是奇素数，a是任意正整数(1≤ a≤ n−1)，则 a^(n-1) ≡ 1 mod n。

 由费马小定理进行的一个推演（和许多博主一样，并未看懂，以后看懂了补写），得到相关的结论：如果n是一个奇素数，将n−1表示成2^s*r的形式，r是奇数，a与n是互素的任何随机整数，那么a^r ≡ 1 mod n或者对某个j (0 ≤ j≤ s−1, j∈Z) 等式a^(2jr) ≡ −1 mod n 成立。这里，为什么需要两个判断条件呢，因为少数非素数也满足费马定理的条件，这类数被称为Carmichael数。虽然这类数的个数是非常少的，但是仍然需要对这些数进行特殊的判定，也就是二次探测。

  也就是对于一个数n，首先，它如果是偶数，1，或者2，直接pass得到判断结果，剩下的就是>2的奇数，先将这个数进行分解，表示成为2^s*r的形式，其实就是得到相关的s和r，然后生成数个a,对每一个a进行以上两个条件的判定，若都满足就继续判断。

   看过一篇相关的文章说，这个对于a=2,满足条件的非素数最小是341，可见经过多次判断可以大大降低失败的概率，若判断次数为n，那么失败的概率约为（1/4）^n，是比较靠谱的一个算法了。

   代码模板：（qmul，qpow快速幂详解）

#include<iostream>

#include<cmath>

#include<cstring>

#include<algorithm>

#include<cstdio>

#include<stdlib.h>

#define times 20

using namespace std;

long long qmul(long long a,long long b,long long M){

    a%=M;

    b%=M;

    long long ans=0;

    while (b){

        if (b&1){

            ans=(ans+a)%M;

        }

        a=(a<<=1)%M;

        b>>=1;

    }

    return ans%M;

}

long long qpow(long long a,long long b,long long int M){

    long long ans=1;

    long long k=a;

    while(b){

        if(b&1)ans=qmul(ans,k,M)%M;

        k=qmul(k,k,M)%M;

        b>>=1;

    }

    return ans%M;

}

bool witness(long long a,long long n,long long x,long long sum){

    long long judge=qpow(a,x,n);

    if (judge==n-1||judge==1)return 1;

    while (sum--){

        judge=qmul(judge,judge,n);

        if (judge==n-1)return 1;

    }

    return 0;

}

bool miller(long long n){

    if (n<2)return 0;

    if (n==2)return 1;

    if ((n&1)==0)return 0;

    long long x=n-1;

    long long sum=0;

    while (x%2==0){

        x>>=1;

        sum++;

    }

    for (long long i=1;i<=times;i++){

        long long a=rand()%(n-1)+1;

        if (!witness(a,n,x,sum))return 0;

    }

    return 1;

}

int main(){

    long long n,m;

    while (scanf("%lld",&n)!=EOF){

        printf("%d\n",miller(n));

    }

}

   Pollard_rho算法

   算法思路：

这个算法并不是一个很说得上有basic的算法，但也不能说是一个很随意的算法，毕竟被称为一个算法。算法对于一个数n，先判断这个数是否是素数，若是素数，那么，记录为一个素因子，若不是，随机取x，y，当然并不是毫无规律，据说是经过大量实验，发现如下取法成功率最高：

先取一个随机数c,然后取另一个随机数x,令y=x，然后x=(x*x+c)%n,得到的x,和y是比较理想的判断数据。

然后取到的x和y求gcd(y-x,n)，若结果大于1，那么得到一个因子p，继续递归搜索n/p,

p的因子，直到搜到素数为止，若gcd结果是1，重复操作找随机数。Ps：中间有个据说是优化的代码，很多人都有这段，也有没有的，但是没有一个博主懂这段是做什么的。

这个算法因为在找寻随机数的过程中会出现成环的情况，有点类似于rho,因而得名。

  代码模板：（基于问题prime test,跟它的题意一样，这就是个裸的模板题，如下是ac代码）

#include<iostream>

#include<cmath>

#include<cstring>

#include<algorithm>

#include<cstdio>

#include<stdlib.h>

#include<time.h>

#define times 20

using namespace std;

long long total;

long long factor[110];

long long qmul(long long a,long long b,long long M){

    a%=M;

    b%=M;

    long long ans=0;

    while (b){

        if (b&1){

            ans=(ans+a)%M;

        }

        a=(a<<=1)%M;

        b>>=1;

    }

    return ans%M;

}

long long qpow(long long a,long long b,long long int M){

    long long ans=1;

    long long k=a;

    while(b){

        if(b&1)ans=qmul(ans,k,M)%M;

        k=qmul(k,k,M)%M;

        b>>=1;

    }

    return ans%M;

}

bool witness(long long a,long long n,long long x,long long sum){

    long long judge=qpow(a,x,n);

    if (judge==n-1||judge==1)return 1;

    while (sum--){

        judge=qmul(judge,judge,n);

        if (judge==n-1)return 1;

    }

    return 0;

}

bool miller(long long n){

    if (n<2)return 0;

    if (n==2)return 1;

    if ((n&1)==0)return 0;

    long long x=n-1;

    long long sum=0;

    while (x%2==0){

        x>>=1;

        sum++;

    }

    for (long long i=1;i<=times;i++){

        long long a=rand()%(n-1)+1;

        if (!witness(a,n,x,sum))return 0;

    }

    return 1;

}

long long gcd(long long a,long long b){

    return b==0?a:gcd(b,a%b);

}

long long pollard(long long n,long long c){

    long long x,y,d,i=1,k=2;

    x=rand()%n;

    y=x;

    while (1){

        i++;

        x=(qmul(x,x,n)+c)%n;

        d=gcd(y-x,n);

        if (d<0)d=-d;//注意，gcd可能为负数

        if (d>1&&d<n)return d;

        if (y==x)return n;

        if (i==k){

            y=x;

            k<<=1;

        }//就是这段，据说的优化，但是并不懂，自己也想了想，确实不知道在干嘛，有些代码里没这段

    }

}

void find(long long n){

    if (miller(n)){

        factor[++total]=n;

        return ;

    }

    long long p=n;

    while (p>=n){

        p=pollard(p,rand()%(n-1)+1);

    }

    find(n/p);

    find(p);

}

int main(){

    long long n,m,i,t;

    scanf("%lld",&t);

    while (t--){

        scanf("%lld",&n);

        if (miller(n)){

            printf("Prime\n");

        }

        else {

            memset(factor,0,sizeof(factor));

            total=0;

            find(n);

            sort(factor+1,factor+total+1);

            printf("%lld\n",factor[1]);

        }

    }

}