zoukankan      html  css  js  c++  java
  • Miller-Rabin

    引子 

    一个数是素数(也叫质数),当且仅当它的约数只有两个——1和它本身。规定这两个约数不能相同,因此1不是素数。对素数的研究属于数论范畴,你可以 看到许多数学家没事就想出一些符合某种性质的素数并称它为某某某素数。整个数论几乎就围绕着整除和素数之类的词转过去转过来。对于写代码的人来说,素数比 想像中的更重要,Google一下BigPrime或者big_prime你总会发现大堆大堆用到了素数常量的程序代码。平时没事时可以记一些素数下来以备急用。

    素数的几个性质   

    1. 素数的个数无限多(不存在最大的素数)
      证明:反证法,假设存在最大的素数P,那么我们可以构造一个新的数$2 * 3 * 5 * 7 * … * P + 1$(所有的素数乘起来加$1$)。显然这个数不能被任一素数整除(所有素数除它都余1),这说明我们找到了一个更大的素数。

    2. 存在任意长的一段连续数,其中的所有数都是合数(相邻素数之间的间隔任意大)
      证明:当$0<a<=n$时,$n!+a$能被$a$整除。长度为$n-1$的数列$n!+2, n!+3, n!+4, …, n!+n$中,所有的数都是合数。这个结论对所有大于$1$的整数$n$都成立,而$n$可以取到任意大。

    3. 所有大于$2$的素数都可以唯一地表示成两个平方数之差。
       证明:大于$2$的素数都是奇数。假设这个数是$2n+1$。由于$(n+1)^2=n^2+2n+1$,$(n+1)^2$和$n^2$就是我们要找的两个平方数。下面证明这个方案是唯一的。如果素数$p$能表示成 $a^2-b^2$,则$p=a^2-b^2=(a+b)(a-b)$。由于p是素数,那么只可能$a+b=p$且$a-b=1$,这给出了a和b的唯一解。

    4. 当$n$为大于2的整数时,$2^n+1$和$2^n-1$两个数中,如果其中一个数是素数,那么另一个数一定是合数。
       证明:$2^n$不能被$3$整除。如果它被$3$除余$1$,那么$2^n-1$就能被$3$整除;如果被3除余$2$,那么$2^n+1$就能被$3$整除。总之,$2^n+1$和$2^n-1$中至少有一个是合数。

    这是大数学家$Fermat$证明的,叫做$Fermat$小定理($Fermat's Little Theorem$)(费马小定理)。$Euler$对这个定理进行了推广,叫做$Euler$定理。$Euler$一生的定理太多了,为了和其它的“$Euler定理$”区别开来,有些地方叫做$Fermat$小定理的$Euler$推广。

    $Euler$定理中需要用一个函数$f(m)$,它表示小于$m$的正整数中有多少个数和$m$互素(两个数只有公约数$1$称为互素)。为了方便,我们通常用记号$φ(m)$来表示这个函数(称作$Euler$函数)。$Euler$指出,如果$a$和$m$互素,那么$a^φ(m) ≡ 1 (mod m)$。可以看到,当$m$为素数时,$φ(m)$就等于$m-1$(所有小于$m$的正整数都与$m$互素),因此它是$Fermat$小定理的推广。

    定理的证明和$Fermat$小定理几乎相同,只是要考虑的式子变成了所有与m互素的数的乘积:$m_1 * m_2 … m_{φ(m)} ≡ (a * m_1)(a * m_2) … (a * m_{φ(m)}) (mod m)$。

    ​ 素数研究的历史

    谈到$Fermat$小定理,数学历史上有很多误解。很长一段时间里,人们都认为$Fermat$小定理的逆命题是正确的,并且有人亲自验证了 $a=2, p<300$的所有情况。国外甚至流传着一种说法,认为中国在孔子时代就证明了这样的定理:如果$n$整除$2^{n-1}-1$,则$n$就是素数。后来某个英国学者进行考证后才发现那是因为他们翻译中国古文时出了错。1819年有人发现了$Fermat$小定理逆命题的第一个反例:虽然$2$的$340$次方除以$341$余$1$,但$341=11*31$。后来,人们又发现了$561$,$ 645$, $1105$等数都表明$a=2$时$Fermat$小定理的逆命题不成立。虽然这样的数不多,但不能忽视它们的存在。于是,人们把所有能整除$2^{n-1}-1$的合数$n$叫做伪素数($pseudoprime$),意思就是告诉人们这个素数是假的。


    ​ 不满足$2^{n-1} mod n = 1$的$n$一定不是素数;如果满足的话则多半是素数。这样,一个比试除法效率更高的素性判断方法出现了:制作一张伪素数表,记录某个范围内的所有伪素数,那么所有满足$2^{n-1} mod n = 1$且不在伪素数表中的$n$就是素数。之所以这种方法更快,是因为我们可以使用二分法快速计算$2^{n-1} mod n$ 的值,这在计算机的帮助下变得非常容易;在计算机中也可以用二分查找有序数列、$Hash$表开散列、构建$Trie$树等方法使得查找伪素数表效率更高。


    ​ 有人自然会关心这样一个问题:伪素数的个数到底有多少?换句话说,如果我只计算$2^{n-1} mod n$的值,事先不准备伪素数表,那么素性判断出错的概率有多少?研究这个问题是很有价值的,毕竟我们是$OIer$,不可能背一个长度上千的常量数组带上考场。 统计表明,在前$10$亿个自然数中共有$50847534$个素数,而满足$2^{n-1} mod n = 1$的合数$n$有$5597$个。这样算下来,算法出错的可能性约为$0.00011$。这个概率太高了,如果想免去建立伪素数表的工作,我们需要改进素性判断的算 法。

    ​ 最简单的想法就是,我们刚才只考虑了$a=2$的情况。对于式子$a^{n-1} mod n$,取不同的$a$可能导致不同的结果。一个合数可能在$a=2$时通过了测试,但$a=3$时的计算结果却排除了素数的可能。于是,人们扩展了伪素数的定义,称满足 $a^{n-1} mod n = 1$的合数$n$叫做以$a$为底的伪素数($pseudoprime to base a$)。前$10$亿个自然数中同时以$2$和$3$为底的伪素数只有$1272$个,这个数目不到刚才的$frac{1}{4}$。这告诉我们如果同时验证$a=2$和$a=3$两种情况,算法出错 的概率降到了0.000025。容易想到,选择用来测试的$a$越多,算法越准确。通常我们的做法是,随机选择若干个小于待测数的正整数作为底数$a$进行若干次 测试,只要有一次没有通过测试就立即把这个数扔回合数的世界。这就是$Fermat$素性测试。
    ​ 人们自然会想,如果考虑了所有小于$n$的底数$a$,出错的概率是否就可以降到$0$呢?没想到的是,居然就有这样的合数,它可以通过所有$a$的测试。$Carmichael$第一个发现这样极端的伪素数,他把它们称作$Carmichael$数。你一定会以为这样的数一定很大。错。第一个$Carmichael$数小得惊人,仅仅是一个三位数,$561$。前$10$亿个自然数中$Carmichael$数也有$600$个之多。$Carmichael$数的存在说明,我们还需要继续加强素性判断的算法。

    ​ Miller_Rabin

    $Miller$和$Rabin$两个人的工作让$Fermat$素性测试迈出了革命性的一步,建立了传说中的$Miller-Rabin$素性测试算法。 新的测试基于下面的定理:**如果$p$是素数,$x$是小于$p$的正整数,且$x^2 mod p = 1$,那么要么$x=1$,要么$x=p-1$。这是显然的,因为$x^2 mod p = 1$相当于$p$能整除$x^2-1$,也即$p$能整除$(x+1)(x-1)$。由于$p$是素数,那么只可能是$x-1$能被$p$整除(此时$x=1$)或$x+1$能被$p$整除(此时 $x=p-1$)。


    ​ 我们下面来演示一下上面的定理如何应用在$Fermat$素性测试上。前面说过$341$可以通过以$2$为底的$Fermat$测试,因为$2^{340} mod 341=1$。如果$341$真是素数的话,那么$2^{170} mod 341$只可能是$1$或$340$;当算得$2^{170} mod 341$确实等于$1$时,我们可以继续查看$2^{85}$除以$341$的结果。我们发现,$2^{85} mod 341=32$,这一结果摘掉了$341$头上的素数皇冠,面具后面真实的嘴脸显现了出来,想假扮素数和我的素MM交往的企图暴露了出来。
    ​ 这就是$Miller-Rabin$素性测试的方法。不断地提取指数$n-1$中的因子$2$,把$n-1$表示成$d*2^r$(其中$d$是一个奇数)。那么我们需要计算的东西就变成了$a$的$d*2^r$次方除以$n$的余数。于是,$a^{d * 2^{r-1}}$要么等于1,要么等于$n-1$。如果$a^{d * 2^{r-1}}$等于$1$,定理继续适用于$a^{d * 2^{r-2}}$,这样不断开方开下去,直到对于某个$i$满足$a^{d * 2^i} mod n = n-1$或者最后指数中的$2$用完了得到的$a^d mod n=1$或$n-1$。这样,$Fermat$小定理加强为如下形式:
    ​ 尽可能提取因子$2$, 把$n-1$表示成$d*2^r$,如果$n$是一个素数,那么或者$a^d mod n=1$,或者存在某个$i$使得$a^{d*2^i} mod n=n-1$ ($ 0<=i<r$ ) (注意$i$可以等于$0$,这就把$a^d mod n=n-1$的情况统一到后面去了)
    ​ $Miller-Rabin$ 素性测试同样是不确定算法,我们把可以通过以$a$为底的$Miller-Rabin$测试的合数称作以$a$为底的强伪素数($strong pseudoprime$)。第一个以$2$为底的强伪素数为$2047$。第一个以$2$和$3$为底的强伪素数则大到$1 373 653$。

    上述内容便于读者有个大概的认识

     

    下面是$Miller-Rabin$的算法流程

    1.先判断是不是$2$,是的话就返回$true$。

    2.判断是不是小于$2$的,或合数,是的话就返回$false$。

    3.令$n-1=u*2^t$,求出$u$,$t$,其中u是奇数。

    4.随机取一个$a$,且$1<a<n$

    (根据费马小定理,如果$a^{n-1}≡1(mod p)$那么$n$就极有可能是素数,如果等式不成立,那肯定不是素数了

    因为$n-1=u*2^t$,所以$a^{n-1}=a^{u*2^t}=(a^u)^{2^t}。)$

    5.所以我们令$x=a^umod n$

    6.然后是$t$次循环,每次循环都让$y=(x*x)mod n,x=y$,这样$t$次循环之后$x=a^{u*2^t}=a^{n-1}$了

    7.因为循环的时候$y=(x*x)mod n$,且$x$肯定是小于$n$的,正好可以用二次探测定理,

    如果$x^2mod n==1$,也就是$y$等于$1$的时候,假如$n$是素数,那么$x==1||x==n-1$,如果$x!=1$且$x!=n-1$,那么$n$肯定不是素数了,返回$false$。

    8.运行到这里的时候$x=a^{n-1}$,根据费马小定理,$x!=1$的话,肯定不是素数了,返回$false$

    9.因为$Miller-Rabin$得到的结果的正确率为75%,所以要多次循环步骤$4~8$来提高正确率

    10.循环多次之后还没返回,那么$n$几乎可以肯定是素数了,返回$true$

    代码

    #include<iostream>
    #include<cstdio>
    #include<algorithm>
    #include<time.h>
    #define ll long long
    using namespace std;
    
    ll qpow(ll x,ll y,ll mod)//快速幂取模 
    {
        ll r=1;
        for (;y;y>>=1,x=x*x%mod) if (y&1) r=r*x%mod;
        return r; 
    }
    bool is_prime(ll n)
    {
        if (n<2) return false;
        if (n==2) return true;
        if (!(n&1)) return false;
        ll u=n-1,t=0;
        while (!(u&1))//算出u,t 
        {
            t++;
            u>>=1;
        }
        for (int i=1;i<=20;i++)
        {
            ll a=rand()%(n-1)+1;
            ll x=qpow(a,u,n),y;
            for (int j=1;j<=t;j++)
            {
                y=x*x%n;
                if (y==1&&x!=1&&x!=n-1) return false;//二次探测定理 
                x=y;
            }
            if (x!=1) return false;//费马小定理 
        }
        return true;
    }
    int main()
    {
        srand(time(0));
        ll n;
        scanf("%lld",&n);
        if (is_prime(n)) puts("YES");
        else puts("NO");
        return 0;
    }
  • 相关阅读:
    MySQL的max()函数使用时遇到的小问题
    scp命令需要指定端口时要紧跟在scp后
    linux系统之间基于密钥对免输入密码登陆
    c++的引用用法
    预测模型
    mysql出现ERROR 1366 (HY000):的解决办法
    R语言可视化--颜色
    R语言可视化--ggplot函数
    R语言可视化--qplot函数
    R语言可视化二
  • 原文地址:https://www.cnblogs.com/xxzh/p/9328883.html
Copyright © 2011-2022 走看看