Miller Rabin算法的改进

BelowLuminous edited 6 年,9 月前

Millar-Rabin质数检验方法:
根据费马小定理,如果p是素数,a<p,那么有a^(p-1) mod p = 1。
直观想法我们直接取若干个a,如果都有一个不满足,那么p就是合数。
遗憾的是,存在Carmichael数:你无论取多少个a,有一个不满足,算我输。
比如:561 = 11X51就是一个Carmichael数。
那么就很江了啊。。我们需要改进算法。
首先有:如果p是素数,x是小于p的正整数,且x^2 mod p = 1,那么要么x=1,要么x=p-1
(这个废话,x=p-1模意义下等于x=-1)
然后我们可以展示下341满足2^340 mod 341 = 1,却不是素数(341=3111)的原因:
2^340 mod 341 = 1
2^170 mod 341 = 1
2^85 mod 341 = 32
(32这个数很py啊怎么不等于340也不等于1啊。。这明显有交易嘛)
那么就能说明这个数不是素数。
如果是素数,一定是从p-1变到1,或是把所有2的次幂去除完,本来就等于1(这样平方完就一直是1了)
所以要么把所有2的次幂去除完,本来就等于1,要么存在某一个次幂=p-1(这样就正常多了)
这就是Millar-Rabin素数验证的二次探测。
应该来说Millar-Rabin算法也是挺好写的
其中mul(a,b,c)表示a
b%c(因为a*b会爆longlong,所以用快速加)
namespace Millar_Rabin {
const int Prime[14] = {0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41};
const int PN = 13;

inline bool witness(int pr, ll res, int times, ll n) {
    ll p = pwr2((ll)pr, res, n);
    if(p == 1) return 0;   
    for (int i=0; i<times; ++i) {
        if(p == n-1) return false;
        if(p == 1) return false;
        p = mul(p, p, n);
    }
    return true;
}

inline bool main(ll n) {
    for (int i=1; i<=PN; ++i) {
        if(n == Prime[i]) return 1;
        if(n % Prime[i] == 0) return 0;
    }
    ll p = n-1;
    int times = 0;
    while(!(p&1)) {
        ++times;
        p >>= 1;
    }
    for (int i=1; i<=PN; ++i) 
        if(witness(Prime[i], p, times, n)) return false;
    return true;
}

}

Comments