LG的数学计划之素数测试(Miller Rabin算法)

我们知道对于素数的一些处理方法,

1.暴枚求质因数,

效率比较。。。。不排除大力出奇迹的情况
代码简单这里不给出

2.素数筛选法,

这个算法比起第一种效率好多了, 是许多题都会考察到的,
算法的思想是枚举每一种可能的质因数, 来求出大范围的素数,
即第一遍我们筛去所有2的倍数, 在筛去所有3的倍数, 5, 7。。。。。。。
最后把 sqrt(N)也筛一遍就好了, 那么可知剩下的就都是素数了,
给出不严格证明:
对于任意合数x, 必定能分解成若干质数的积, 而在
1~N数列中, 最大的情况必定是sqrt(N), 那么划去所有质数能做积得到的数, 剩下的必然是无法分解出这些质因数;(仅仅是理解性证明)

void Prime_pre() {
    for(int i = 2; i*i <= N; ++i)
        if(!p[i])
            for(int j = i<<1; j <= N; j += i) p[j] = true;
}

给出的代码(!) 为质数

3.下面给出本篇博客的主角素数测试算法,

因为虽然第二种算法对于处理区间的算法时很实用但是有时我们不需要这么大的区间, 并且有时候题目给出的数可能相当大, 空间是不够的, 那么就需要我们能够快速判定一个数是否为素数, —Miller Rabin随机算法
首先给出这个算法的理论依据
我们有费马小定理

a^(p-1) ≡ 1 % p(p为质数, a, p互质)
显然对于这个定理的逆定理是显然不成立的, 当我们可以说如果对于1 < a < p
如果有a^(p-1) % p != 1 那么p一定不是质数
所以我们在(1~p-1)中多随机几个a, 都来判定一下这些a满不满足这些鬼
多试几个就能得出正确率比较高的结论;

其实有这个结论就够了, 但我们要更高的效率
二次探测定理

0< a < p (p是素数) a^2 % p = 1 则a1 = 1, a2 = p-1;
证明略过(不是太会, 应该是用类似环?的方式来求证, 网上没有好的证明, 有读者会请不吝赐教)

有了这个算法, 就能优化了
一开始对于每一个a我们都要求出a^(p-1)
但其实我们可以进行分解

我们可以先把p−1分解成2t×u(u∈{x|x=2∗k+1,k∈N})的形式
这样原式就是(a^u)^2t, 就是a^u的t次平方注意是平方, 那么我们就可以用二次探测定理了, 我们把x0 = a^u, 设x[i]为x0的i次平方, 那么一旦出现x[i] % p == 1, 那么我们就判定x[i-1]是否是1或者p-1, 不是的话p一定不是质数, 是不是觉得没什么区别, 其实区别很大, 你用第一种办法时,测试p0是否为素数 ,假如我们选20个a得到p0是素数, 那么这20个a可能出现a^(p0-1) % p0==1 但只是碰巧, p0是个合数
如虽然2的340次方除以341余1,但341=11*31
这时采用第一种方法我们就必须多试50个, 甚至100个a才可以排除这种情况
冲突几率增大, 那么就浪费了时间,
那么我们用二次探测定理判断条件就更加苛刻, 可能就能排除这些最后结果满足, 中间却不满足的情况, 只需20次仍能得出p0是素数的正确结论
论毕, 给出代码

//a*b 采用分治方法来做, 加速防爆
ll pls(ll a, ll b, ll p) {
    ll res = 0;
    for(; b; b >>= 1, a = (a+a)%p)
        if(b & 1) res = (res+a)%p;
    return res%p;
}

ll pow(ll a, ll b, ll p) {
    ll res = 1;
    for(; b; b >>= 1, a = pls(a, a, p))
        if(b & 1) res = pls(res, a, p)%p;
    return res % p;
}
//p是要判断的数
bool M(ll p) {
    ll x[60] = {0}, s = 20;
    ll rest = p-1, t = 0;
    //分解p-1
    while(!(rest%2)) {
        t ++;
        rest >>= 1;
    }
//s就是我们要试几个a来最后拿结论
    while(s--) {
        ll a = rand()%(p-1)+1;
        x[0] = pow(a, rest, p);
        rep(i, 1, t) {
            x[i] = pow(x[i-1], 2, p)%p;
            if(x[i] == 1) 
                if((x[i-1] != 1) && (x[i-1] != p-1)) return false;
        }
        if(x[t] ^ 1) return false;
    }
    return true;
}

写的好长啊, 看完的点个赞吧。。。。

原文地址:https://www.cnblogs.com/pbvrvnq/p/8530161.html