Miller-Rabin素数测试算法（POJ1811Prime Test）

时间：2023-03-08 17:16:04

题目链接：http://poj.org/problem?id=1811

题目解析：2<=n<2^54，如果n是素数直接输出，否则求N的最小质因数。

求大整数最小质因数的算法没看懂，不打算看了，直接贴代码，以后当模版用。

数据比较大，只能先用Miller_Rabin算法进行素数判断。
在用Pollard_rho分解因子。
#include <iostream>

#include <stdio.h>

#include <string.h>

#include <stdlib.h>

#include <math.h>

#include <time.h>

#include <algorithm>

typedef long long ll;

#define Time 15 //随机算法判定次数，Time越大，判错概率越小

using namespace std;

ll n,ans,factor[];//质因数分解结果（刚返回时是无序的）

ll tol;//质因数的个数，数组下标从0开始

//****************************************************************

// Miller_Rabin 算法进行素数测试

//速度快，而且可以判断 <2^63的数

//****************************************************************

long long mult_mod(ll a,ll b,ll c)//计算 (a*b)%c.   a,b都是ll的数，直接相乘可能溢出的

{

    a%=c;//                           利用二分思想减少相乘的时间

    b%=c;

    ll ret=;

    while(b)

    {

        if(b&)

        {

            ret+=a;

            ret%=c;

        }

        a<<=;

        if(a>=c)a%=c;

        b>>=;

    }

    return ret;

}

ll pow_mod(ll x,ll n,ll mod)//x^n%n

{

    if(n==)return x%mod;

    x%=mod;

    ll tmp=x;

    ll ret=;

    while(n)

    {

        if(n&) ret=mult_mod(ret,tmp,mod);

        tmp=mult_mod(tmp,tmp,mod);

        n>>=;

    }

    return ret;

}

//以a为基,n-1=x*2^t      a^(n-1)=1(mod n)  验证n是不是合数

//一定是合数返回true,不一定返回false

//二次探测

bool check(ll a,ll n,ll x,ll t)

{

    ll ret=pow_mod(a,x,n);

    ll last=ret;

    for(int i=; i<=t; i++)

    {

        ret=mult_mod(ret,ret,n);

        if(ret==&&last!=&&last!=n-) return true;//合数

        last=ret;

    }

    if(ret!=) return true;

    return false;

}

// Miller_Rabin()算法素数判定

//是素数返回true.(可能是伪素数，但概率极小)

//合数返回false;

bool Miller_Rabin(ll n)

{

    if(n<)return false;

    if(n==||n==||n==||n==)return true;

    if(n==||(n%==)||(n%==)||(n%==)||(n%==)) return false;//偶数

    ll x=n-;

    ll t=;

    while((x&)==)

    {

        x>>=;

        t++;

    }

    for(int i=; i<Time; i++)

    {

        ll a=rand()%(n-)+;//rand()需要stdlib.h头文件

        if(check(a,n,x,t))

            return false;//合数

    }

    return true;

}

//************************************************

//pollard_rho 算法进行质因数分解

//************************************************

ll gcd(ll a,ll b)

{

    if(a==)return ;

    if(a<) return gcd(-a,b);

    while(b)

    {

        long long t=a%b;

        a=b;

        b=t;

    }

    return a;

}

ll Pollard_rho(ll x,ll c)

{

    ll i=,k=;

    ll x0=rand()%x;

    ll y=x0;

    while()

    {

        i++;

        x0=(mult_mod(x0,x0,x)+c)%x;

        long long d=gcd(y-x0,x);

        if(d!=&&d!=x) return d;

        if(y==x0) return x;

        if(i==k)

        {

            y=x0;

            k+=k;

        }

    }

}

//对n进行素因子分解

void findfac(ll n)

{

    if(Miller_Rabin(n))//素数

    {

        factor[tol++]=n;

        return;

    }

    ll p=n;

    while(p>=n) p=Pollard_rho(p,rand()%(n-)+);

    findfac(p);//递归调用

    findfac(n/p);

}

int main()

{

    int T;

    //srand(time(NULL));加上RE不懂

    scanf("%d",&T);

    while(T--)

    {

        scanf("%lld",&n);//(n>=2)

        /*if(n==1)

        {

            printf("1\n");

            continue;

        }*/

        if(Miller_Rabin(n))

        {

            printf("Prime\n");

            continue;

        }

        tol=;

        findfac(n);//对n分解质因子

        ll ans=factor[];

        for(int i=; i<tol; i++)

            if(factor[i]<ans)

                ans=factor[i];

        /*for(int i=0;i<tol;i++)

        {

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

        }*/

        printf("%lld\n",ans);

    }

    return ;

}

算法解析：

由费马小定理可以知道，若p是素数且a是整数，则满足a^p==a(mod p)。若存在正整数a不满足a^p==a(mod p),那么p是合数。

定义：令a是一个正整数，若p是合数且满足a^p==a(mod p)，则p称为以a为基的伪素数。

Miller-Rabin素数测试算法原理： 假如p是素数，且(a,p)==1,(a为任意小于p的正整数），那么a^p-1==1(mod p)。如果a^p-1==1(mod p)，

则可认为n是素数，取多个底进行试验，次数越多，n为素数概率越大。（我的个人理解多次试验为p换基，使之成为伪素数的可能性大大减小）。

（Miller-Rabin测试：不断选取不超过n-1的基b(s次)，计算是否每次都有bn-1 ≡ 1(mod n)，若每次都成立则n是素数，否则为合数。）

转载：说Miller-Rabin测试以前先说两个比较高效的求a*b% n 和 ab %n 的函数，这里都是用到二进制思想，将b拆分成二进制，然后与a相加（相乘）

// a * b % n

//例如: b = 1011101那么a * b mod n = (a * 1000000 mod n + a * 10000 mod n + a * 1000 mod n + a * 100 mod n + a * 1 mod n) mod n 

ll mod_mul(ll a, ll b, ll n) {

    ll res = ;

    while(b) {

        if(b&)    res = (res + a) % n;

        a = (a + a) % n;//a=(a<<1)%n

        b >>= ;

    }

    return res;

}

这代码很棒，以后计算a*b时，如果里面有一个数很大，则可以选择上面的算法，(nlogn)的时间复杂度。

//a^b % n

//同理

ll mod_exp(ll a, ll b, ll n) {

    ll res = ;

    while(b) {

        if(b&)    res = mod_mul(res, a, n);

        a = mod_mul(a, a, n);

        b >>= ;

    }

    return res;

}

快速幂，没什么好说的。

核心代码：

开始程序时需加srand(time(NULL));

bool miller_rabin(ll n)

{

    for(int i=; i<=N; i++) //N为你打算测试的次数，N(10~20）

    {

        ll a=random(n-)+;//需头文件stdlib.h,random(X)产生0~X的随机数，+1产生1~n-1

        if(mod_exp(a,n-,mod)!=)

        {

            "合数";

        }

    }

}

注意，MIller-Rabin测试是概率型的，不是确定型的，不过由于多次运行后出错的概率非常小，所以实际应用还是可行的。（一次Miller-Rabin测试其成功的概率为3/4）

二次探测定理：（改进）

一个合数n,若对所有满足（b,n)=1的正整数b都有b^n-1==1(mod n)成立，（上面的反例，但出现这种数的几率不大），则称之为卡迈克尔数。

　二次探测　如果p是奇素数，则 x2 ≡ 1(mod p)的解为 x = 1 || x = p - 1(mod p);

可以利用二次探测定理在实现Miller-Rabin上添加一些细节，具体实现如下：

bool miller_rabin(ll n) {

    if(n ==  || n ==  || n ==  || n ==  || n == )    return true;

    if(n ==  || !(n%) || !(n%) || !(n%) || !(n%) || !(n%))    return false;

    ll x, pre, u;

    int i, j, k = ;

    u = n - ;    //要求x^u % n

    while(!(u&)) {    //如果u为偶数则u右移，用k记录移位数

        k++; u >>= ;

    }

    srand((ll)time());

    for(i = ; i < S; ++i) {    //进行S次测试

        x = rand()%(n-) + ;    //在[2, n)中取随机数

        if((x%n) == )    continue;

        x = mod_exp(x, u, n);    //先计算(x^u) % n，

        pre = x;

        for(j = ; j < k; ++j) {    //把移位减掉的量补上，并在这地方加上二次探测

            x = mod_mul(x, x, n);

            if(x ==  && pre !=  && pre != n-)    return false;    //二次探测定理，这里如果x = 1则pre 必须等于 1，或则 n-1否则可以判断不是素数

            pre = x;

        }

        if(x != )    return false;    //费马小定理

    }

    return true;

}

前边这个算法经过测试还是比较靠谱的，可以用作模板。

效率上，VC 10 RELEASE 模式下，采用三次循环 M - R，测试第 19999 个素数 224729 时，快除法快 而测试第 20000 个素数 224737 时，M - R 法快

因此，为保证最高效，测试大数 n 时，可以先对其使用前 19999 个素数进行快除法排除，而后再使用 M - R 测试。

AC_Von 原创，转载请注明出处：http://www.cnblogs.com/vongang/archive/2012/03/15/2398626.html

相关文章