程序员

Miller–Rabin

2020-08-04  本文已影响0人  dachengquan

Miller–Rabin算法用于检测大素数,时间复杂度是O(k\log^3 n)
使用费马小定理和二次探测检测一个数是否是质数。a^{p-1} \equiv 1 \pmod{p}当p是质数时,这个公式成立。但是这个公式成立p不一定是质数。因为有些合数,对于特定的a同样得到这个结果。其中有一类数无论a是什么都满足这个公式,例如561。如果不存在这样的数选择多次不同的a进行验证,使用费马小定理就可以检测出一个数是否是质数。
二次探测定理:如果p是素数 ,x^2\equiv1\pmod p的解是x\equiv 1 \pmod porx\equiv -1 \pmod p
证明过程就是
p\mid x^2-1
p\mid (x-1)(x+1),p|(x-1)p|(x+1)都能使这个式子成立。而x在1到p范围内的解只有1和p-1。

Miller-Rabin算法流程

给定一个N判断是否是整数,将N-1分解为2^t *d的形式。
我们要判断对于N来说是否能使费马小定理成立。随机选择一个a,1<a<p
a^{N-1} = a^{d*2^t} = (a^d)^{2^t}\equiv 1\pmod p费马小定理告诉我们如果N是质数最后结果模p同余1。计算a^{N-1}的过程中,可以分为计算a^d,然后对a^d进行t次平方,第二步的过程我们可以使用二次探测定理验证。

举例子:一个质数561 计算N-1 = 560

分解560 = 2^4 *35
2^{35} \equiv 263 \pmod {561}
(2^{35})^2 \equiv 166 \pmod {561}
(2^{35})^4 \equiv 67 \pmod {561}
(2^{35})^8 \equiv 1 \pmod {561}
根据二次探测定理可知p不是质数。
67^2 \equiv 1 \pmod {561}
67^2 -1\equiv 0 \pmod {561}
(67-1)(67+1)\equiv 0 \pmod {561}
66*68\equiv 0 \pmod {561}
561|66*68说明66*68拥有561的全部质因数,66<561并且68<561,所以561的质因数必定来自66的一部分和68的一部分。

a如何选取?

通常应该取随机值,但是有人测试过在下面的范围取固定的几个a就可以了。
如果n <2,047,则测试a = 2 就足够了;
如果n <1,373,653,则测试a = 2和3 就足够了;
如果n <9,080,191,则测试a = 31和73 就足够了;
如果n <25,326,001,则足以测试a = 2、3和5;
如果n <3,215,031,751,测试a = 2、3、5 和7 就足够了;
如果n <4,759,123,141,测试a = 2、7和61 就足够了;
如果n <3,474,749,660,383,测试a = 2、3、5、7、11 和13 就足够了;
如果n <341,550,071,728,321,则足以测试a = 2、3、5、7、11、13 和17。
如果n <3,825,123,056,546,413,051,则足以测试a = 2、3、5、7、11、13、17、19 和23
随机选择a,为了确保正确率通常选择的k比较高,我们知道有以上结论,每次测试直接选择前几个质数就可以了。当我们的质数选择10时,已经可以测试long long范围内的数了。
Miller-Rabin时间复杂度是O(k\log^3 n)k是选取a的次数。

#include<bits/stdc++.h>
#define ll long long 
#define i128 __int128
using namespace std;
ll n;
int primes[19]={2,3,5,7,11,13,17,19,23,29,31,37};
ll poww(ll a,ll b,ll c){
    ll ans=1;i128 base=a;
    while(b){
        if(b&1)ans=ans*base%c;
        base=base*base%c;
        b>>=1;
    }
    return ans;
}
bool miller_rabin(ll n){
    for(int i=0;i<=8;i++){
        if(primes[i]==n)return 1;
        else if(primes[i]>n)return 0;
        ll t=poww(primes[i],n-1,n),x=n-1;
        if(t!=1)return 0;
        while(t==1&&x%2==0){
            x/=2;
            t=poww(primes[i],x,n);
            if(t!=1&&t!=n-1)return 0;
        }
    }
    return 1;
}
int main(){
    ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
    while(cin>>n){
         if(miller_rabin(n))puts("Y");else puts("N");
    }
    return 0;
}

使用了__int128省去了自己写快速乘法的时间。
判定1e18的质数

上一篇 下一篇

猜你喜欢

热点阅读