Miller-Rabin素数判定算法

算法内容

在通常情况下,判断一个数是不是质数使用的是枚举因子的方式,这种方法的时间复杂度为

$O(\sqrt{n})$

当n的规模较小时比较常用.优点是正确率高(只要代码没写错就不可能判定出错).缺点是当n极大时,算法进行的无效判断实在太多了,导致时间效率难以接受.而这应该就是准确判断单个数是不是质数的时间下界了(暂时没找到严谨的证明).所以,为了追求效率,需要牺牲一定的准确性,这时就引入Miller-Rabin算法.

Miller-Rabin算法是基于费马小定理提出的.费马小定理的原版表述是:当p是一个质数时,对于任意的与p互质的数a,有

$a^{p-1} \equiv 1 (\mod p)$

但是,仅仅根据这一个定理来判断素数,正确率有些低.因为满足上面式子的合数实在是有点多.而且前面几个都不是很大,这些数有个统一的称呼是Carmichael数(OEIS的地址见这里).当a选到这些数时会就导致出现bug.所以还需要一些辅助手段帮助判断.

引入二次探测定理.当p是质数时,假如有

$a^{2} \equiv 1 (\mod p)$

则一定有

$a \equiv 1 (\mod p)$

或

$a \equiv -1 (\mod p)$

根据上面两个定理,就可以得到如下随机算法.Miller-Rabin的原理如下:

首先,对于整数N(显然,N是奇数),求出m和r,使得

$N=2^{r}m+1$

接下来进行随机测试,先rand出一个整数a,然后计算

$a^{m}$

之后,把上面的结果平方r次.每平方一次,就检查能否使用二次探测定理.平方完毕后,刚好计算出了

$a^{2^rm}=a^{p-1}$

再拿费马小定理判断一下.正确率基于以下结论:

$若p是一个合数,则p为合数需要的证据数量为\frac{p-1}{2}$

实际上,小于p的非证据很少,大概只有上面的一半.所以,判断q次的正确率可以达到

$1-\frac{1}{4^{q}}$

一般取q为10就可以了.实在是脸黑可以适当增加q.

代码模板

数据规模比较大,所以对乘法做了特殊处理.

inline LL Mul(LL a,LL b,LL MOD){
    LL Ans=a*b-(LL)((long double)a*b/MOD+0.5)*MOD;
    return (Ans<0)?Ans+MOD:Ans;
}

inline LL Pow(LL x,LL y,LL MOD){
    LL p=x,t=1,Res=1;
    for (;t<=y;(t&y)?Res=Mul(Res,p,MOD):0,p=Mul(p,p,MOD),t<<=1);
    return Res;
}

inline bool Miller_Rabin(LL x){
    int i=0,j=0,t=0;
    LL m=x-1,r=0,l=0;
    if (!(n&1)) return 0;
    while (!(m&1))
        t++,m>>=1;
    for (i=1;i<=10;i++){
        l=r=Pow(rand(),m,x);
        for (j=1;j<=t;j++){
            r=Mul(r,r,x);
            if (r==1&&l!=1&&l!=x-1) return 0;
            l=r;
        }
        if (r!=1)   return 0;
    }
    return 1;
}

关于整数N的选取,上面选择的是随机数,实际上可以钦定一些质数作为N,这样可以在一定程度上减少由于随机导致失误的概率,理论上下面的写法也会更快.

优化版:

inline LL Mul(LL a,LL b,LL MOD){
    LL Ans=a*b-(LL)((long double)a*b/MOD+0.5)*MOD;
    while (Ans<0)   Ans+=MOD;
    return (Ans<0)?Ans+MOD:Ans;
}

inline LL Pow(LL x,LL y,LL MOD){
    LL p=x,t=1,Res=1;
    for (;t<=y;(t&y)?Res=Mul(Res,p,MOD):0,p=Mul(p,p,MOD),t<<=1);
    return Res;
}

inline bool Miller_Rabin(LL x){
    int i=0,j=0,t=0;
    LL m=x-1,r=0,l=0;
    if (x==2)   return 1;
    if (!(x&1)) return 0;
    if (x<=800000)  return 1-v[x];
    while (!(m&1))  t++,m>>=1;
    //g[1] = 3;g[2] = 5;g[3] = 7;g[4] = 11;g[5] = 13;
    //g[6] = 17;g[7] = 19;g[8] = 23;
    for (i=1;i<=8;i++){
        l=r=Pow(g[i],m,x);
        for (j=1;j<=t;j++){
            r=Mul(r,r,x);
            if (r==1&&l!=1&&l!=x-1) return 0;
            l=r;
        }
        if (r!=1)   return 0;
    }
    return 1;
}

相关例题

[LOJ143]质数判定

完整代码:

#include <cstring>
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <vector>
#include <cmath>
#include <cstdlib>
#include <ctime>
#define N 100010
#define l(x) (x << 1)
#define r(x) ((x << 1) + 1)
#define LL long long
#define INF 0x3f3f3f3f
using namespace std;

int i;
LL n;
int v[N], p[N];

inline int Abs(int x) { return (x < 0) ? -x : x; }
inline void Swap(int &a, int &b) { a ^= b ^= a ^= b; }
inline int Min(int a, int b) { return (a < b) ? a : b; }
inline LL Max(LL a, LL b) { return (a > b) ? a : b; }

inline int read() {
    int p = 0;
    char c = getchar();
    while (c < 48 || c > 57) c = getchar();
    while (c >= 48 && c <= 57) p = (p << 1) + (p << 3) + c - 48, c = getchar();
    return p;
}

inline void Ready() {
    LL i = 0, x = 0;
    v[1] = 1;
    for (i = 2; i <= 100000; i++) {
        if (v[i])
            continue;
        p[++p[0]] = i;
        for (x = (LL)i * i; x <= 100000; x += i) v[x] = 1;
    }
}

inline LL Mul(LL a, LL b, LL MOD) {
    LL Ans = a * b - (LL)((long double)a * b / MOD + 0.5) * MOD;
    return (Ans < 0) ? Ans + MOD : Ans;
}

inline LL Pow(LL x, LL y, LL MOD) {
    LL p = x, t = 1, Res = 1;
    for (; t <= y; (t & y) ? Res = Mul(Res, p, MOD) : 0, p = Mul(p, p, MOD), t <<= 1)
        ;
    return Res;
}

inline bool Miller_Robin(LL x) {
    int i = 0, j = 0, t = 0;
    LL m = x - 1, r = 0, l = 0;
    if (!(n & 1))
        return 0;
    while (!(m & 1)) t++, m >>= 1;
    for (i = 1; i <= 10; i++) {
        l = r = Pow(rand(), m, x);
        for (j = 1; j <= t; j++) {
            r = Mul(r, r, x);
            if (r == 1 && l != 1 && l != x - 1)
                return 0;
            l = r;
        }
        if (r != 1)
            return 0;
    }
    return 1;
}

int main() {
    srand(time(NULL));
    Ready();
    while (~scanf("%lld", &n)) {
        if (n <= 100000) {
            printf((v[n]) ? "N\n" : "Y\n");
            continue;
        }
        if (Miller_Robin(n))
            printf("Y\n");
        else
            printf("N\n");
    }
    return 0;
}

优化版:

#include<cstring>
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<vector>
#include<cmath>
#include<cstdlib>
#include<ctime>
#define N 800010
#define l(x) (x<<1)
#define r(x) ((x<<1)+1)
#define LL long long
#define INF 0x3f3f3f3f
using namespace std;

LL T,x,Ans;
int v[N],p[N],s[N],g[10];

inline LL Abs(LL x){return (x<0)?-x:x;}
inline void Swap(int &a,int &b){a^=b^=a^=b;}
inline int Min(int a,int b){return (a<b)?a:b;}
inline int Max(int a,int b){return (a>b)?a:b;}
inline LL Gcd(LL a,LL b){return (!b)?a:Gcd(b,a%b);}

inline LL read(){
    LL p=0; char    c=getchar();
    while (c<48||c>57)  c=getchar();
    while (c>=48&&c<=57)    p=(p<<1)+(p<<3)+c-48,c=getchar();
    return p;
}

inline void Ready(){
    LL i=0,x=0;
    v[1]=1;
    for (i=2;i<=800000;i++){
        if (v[i])   continue;
        p[++p[0]]=i;
        for (x=(LL)i*i;x<=800000;x+=i)  v[x]=1;
    }
    g[1]=3; g[2]=5; g[3]=7; g[4]=11;    g[5]=13;
    g[6]=17;    g[7]=19;   g[8]=23;
}

inline LL Mul(LL a,LL b,LL MOD){
    LL Ans=a*b-(LL)((long double)a*b/MOD+0.5)*MOD;
    while (Ans<0)   Ans+=MOD;
    return (Ans<0)?Ans+MOD:Ans;
}

inline LL Pow(LL x,LL y,LL MOD){
    LL p=x,t=1,Res=1;
    for (;t<=y;(t&y)?Res=Mul(Res,p,MOD):0,p=Mul(p,p,MOD),t<<=1);
    return Res;
}

inline bool Miller_Rabin(LL x){
    int i=0,j=0,t=0;
    LL m=x-1,r=0,l=0;
    if (x==2)   return 1;
    if (!(x&1)) return 0;
    if (x<=800000)  return 1-v[x];
    while (!(m&1))  t++,m>>=1;
    for (i=1;i<=8;i++){
        l=r=Pow(g[i],m,x);
        for (j=1;j<=t;j++){
            r=Mul(r,r,x);
            if (r==1&&l!=1&&l!=x-1) return 0;
            l=r;
        }
        if (r!=1)   return 0;
    }
    return 1;
}

int main(){
    Ready();
    while (~scanf("%lld",&x)){
        if (x<=800000){
            printf((1-v[x])?"Y\n":"N\n");
            continue;  
        }
        (Miller_Rabin(x))?printf("Y\n"):printf("N\n");
    }
    return 0;
}

Pollard-Rho质因子分解算法

算法内容

在提升了判断质数的效率后,下一步自然是想办法提升质因数分解的效率.通常情况下,分解一个数的时间复杂度为

$O(\sqrt{n})$

当n的规模较大时,时间效率会变得极其低下.(这也是RSA加密安全性的保证之一).类似的,也可以利用随机算法通过牺牲一小部分正确性来提高效率.

很自然的一个想法是随机一个数x,判断是不是n的质数.但是这种朴素随机的命中的期望太小了.于是,有大神对这种方法进行了改进(网上说实际利用了生日悖论),得到了下面的随机方式:

$x_1=rand,c=rand$

$x_{i+1}=(x_{i}^2+c) \mod \ n$

然后检查

$d=gcd(Abs(x_{i+1}-x_{i},n))$

有人通过大量实验得出,d是n的约数的概率是

$n^{\frac{1}{4}}$

因为在计算的过程中要取模,所以这种随机方式在经过足够长的时间后就会出现环.关于判环,Floyd提出了一个很生动的方法:比如两个人在跑步,一个人的速度固定(每次只计算一次递归式),另一个人的速度是前一个人的速度的两倍(每次计算两回递归式),当两个人再次相遇时,第一个人一定已经跑了一圈(出现了环).画出来的图特别像希腊字母"rho".于是算法就叫Pollard-Rho算法.

当然,在实际写代码的时候,更为常见的是倍增判环的思想.每次让y记录x的取值,然后当x再经过当前迭代次数的一倍的后再与y作比较.反复记录,直到y与x相等.这样可以保证正确性,而且不会进行太多的冗余计算.

根据上面的分析,对一个数进行质因数分解的效率可以达到

$O(n^{\frac{1}{4}}logn)$

算导的说法是

$O(\sqrt{p}logn)$

p为n的因子.

威尔逊定理

存粹是翻书的时候恰好看到的,顺手记一下

威尔逊定理:

p是质数的充要条件是:

$(p-1)! \equiv -1 (\mod \ p)$

根据这个定理,可以构造一个以所有质数为零点的函数:

$f(x)=\sin(\pi \frac{(n-1)!+1}{n})$

声明: 本文采用 BY-NC-SA 协议进行授权，如无注明均为原创，转载请注明转自 Marvolo
本文地址: Miller-Rabin算法和Pollard-Rho算法学习笔记