【梅森素数】——nefu120——lucas-Lehmer判定法 & MIller-Robin素数测试法——***模板***

                                               梅森素数

Problem:120

Time Limit:1000ms

Memory Limit:65536K

Description

由于梅森学识渊博，才华横溢，为人热情以及最早系统而深入地研究2p－1 型的数（其中p为素数），为了纪念他，数学界就把这种数称为“梅森数”；并以Mp 记之（其中M为梅森姓名的首字母），即Mp=2p-1 。如果梅森数为素数，则称之为“梅森素数”。 比如p＝2，3，5，7时，Mp都是素数，但211-1 不是素数 。现在请你求出前N个梅森素数。

Input

有多组测试数据。
第一行是一个正整数T，表示测试数据的组数。接下来每组1个数p的值，这里2<= p <= 62。

Output

对于每组测试数据，判断Mp 是不是梅森素数，是就输出“yes ”，否就输出“no”，输出后要换行。

Sample Input

2
2
7

Sample Output

yes
yes

lucas-Lehmer判定法：

解题思路：1.梅森素数的判定方法之一——Lucas—Lehmer判别法：构造rk序列判断Mp是否是素数，同时应用大数乘积取模的算法，将乘法变成加法
附表：梅森素数表（1-47个）
序号 p Mp=（2^p）-1 Mp的位数 发现时间 发现者
1 2 3 1 古代 古人
2 3 7 1 古代 古人
3 5 31 2 古代 古人
4 7 127 3 古代 古人
5 13 8191 4 1456年 无名氏
6 17 131071 6 1588年 Cataldi
7 19 524287 6 1588年 Cataldi
8 31 2147483647 10 1772年 欧拉
9 61 2305843009213693951 19 1883年 Pervushin
10 89 618970019642690137449562111 27 1911年 Powers
11 107 162259276829213363391578010288127 33 1914年 Powers
12 127 170141183460469231731687303715884105727 39 1876年 卢卡斯
13 521 6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151 157 1952年1月30日 Robinson
14 607 531137992816767098689588206552468627329593117727031923199444138200403559860852242739162502265229285668889329486246501015346579337652707239409519978766587351943831270835393219031728127 183 1952年1月30日 Robinson
15 1279 10407932194664399081925240327364085538615262247266704805319112350403608059673360298012239441732324184842421613954281007791383566248323464908139906605677320762924129509389220345773183349661583550472959420547689811211693677147548478866962501384438260291732348885311160828538416585028255604666224831890918801847068222203140521026698435488732958028878050869736186900714720710555703168729087 386 1952年6月25日 Robinson
16 2203 1475979915214180235084898622737381736312066145333169775147771216478570297878078949377407337049389289382748507531496480477281264838760259191814463365330269540496961201113430156902396093989090226259326935025281409614983499388222831448598601834318536230923772641390209490231836446899608210795482963763094236630945410832793769905399982457186322944729636418890623372171723742105636440368218459649632948538696905872650486914434637457507280441823676813517852099348660847172579408422316678097670224011990280170474894487426924742108823536808485072502240519452587542875349976558572670229633962575212637477897785501552646522609988869914013540483809865681250419497686697771007 664 1952年10月7日 Robinson
17 2281 446087557183758429571151706402101809886208632412859901111991219963404685792820473369112545269003989026153245931124316702395758705693679364790903497461147071065254193353938124978226307947312410798874869040070279328428810311754844108094878252494866760969586998128982645877596028979171536962503068429617331702184750324583009171832104916050157628886606372145501702225925125224076829605427173573964812995250569412480720738476855293681666712844831190877620606786663862190240118570736831901886479225810414714078935386562497968178729127629594924411960961386713946279899275006954917139758796061223803393537381034666494402951052059047968693255388647930440925104186817009640171764133172418132836351 687 1952年10月9日 Robinson
18 3217 259117086013202627776246767922441530941818887553125427303974923161874019266586362086201209516800483406550695241733194177441689509238807017410377709597512042313066624082916353517952311186154862265604547691127595848775610568757931191017711408826252153849035830401185072116424747461823031471398340229288074545677907941037288235820705892351068433882986888616658650280927692080339605869308790500409503709875902119018371991620994002568935113136548829739112656797303241986517250116412703509705427773477972349821676443446668383119322540099648994051790241624056519054483690809616061625743042361721863339415852426431208737266591962061753535748892894599629195183082621860853400937932839420261866586142503251450773096274235376822938649407127700846077124211823080804139298087057504713825264571448379371125032081826126566649084251699453951887789613650248405739378594599444335231188280123660406262468609212150349937584782292237144339628858485938215738821232393687046160677362909315071 969 1957年9月8日 Riesel
19 4253 190797007524439073807468042969529173669356994749940177394741882673528979787005053706368049835514900244303495954950709725762186311224148828811920216904542206960744666169364221195289538436845390250168663932838805192055137154390912666527533007309292687539092257043362517857366624699975402375462954490293259233303137330643531556539739921926201438606439020075174723029056838272505051571967594608350063404495977660656269020823960825567012344189908927956646011998057988548630107637380993519826582389781888135705408653045219655801758081251164080554609057468028203308718724654081055323215860189611391296030471108443146745671967766308925858547271507311563765171008318248647110097614890313562856541784154881743146033909602737947385055355960331855614540900081456378659068370317267696980001187750995491090350108417050917991562167972281070161305972518044872048331306383715094854938415738549894606070722584737978176686422134354526989443028353644037187375385397838259511833166416134323695660367676897722287918773420968982326089026150031515424165462111337527431154890666327374921446276833564519776797633875503548665093914556482031482248883127023777039667707976559857333357013727342079099064400455741830654320379350833236245819348824064783585692924881021978332974949906122664421376034687815350484991 1,281 1961年11月3日 Hurwitz
20 4423 285542542228279613901563566102164008326164238644702889199247456602284400390600653875954571505539843239754513915896150297878399377056071435169747221107988791198200988477531339214282772016059009904586686254989084815735422480409022344297588352526004383890632616124076317387416881148592486188361873904175783145696016919574390765598280188599035578448591077683677175520434074287726578006266759615970759521327828555662781678385691581844436444812511562428136742490459363212810180276096088111401003377570363545725120924073646921576797146199387619296560302680261790118132925012323046444438622308877924609373773012481681672424493674474488537770155783006880852648161513067144814790288366664062257274665275787127374649231096375001170901890786263324619578795731425693805073056119677580338084333381987500902968831935913095269821311141322393356490178488728982288156282600813831296143663845945431144043753821542871277745606447858564159213328443580206422714694913091762716447041689678070096773590429808909616750452927258000843500344831628297089902728649981994387647234574276263729694848304750917174186181130688518792748622612293341368928056634384466646326572476167275660839105650528975713899320211121495795311427946254553305387067821067601768750977866100460014602138408448021225053689054793742003095722096732954750721718115531871310231057902608580607 1,332 1961年11月3日 Hurwitz
21 9689 478220278…225754111 2,917 1963年5月11日 Gillies
22 9941 346088282…789463551 2,993 1963年5月16日 Gillies
23 11213 281411201…696392191 3,376 1963年6月2日 Gillies
24 19937 431542479…968041471 6,002 1971年3月4日 布莱恩特·塔克曼
25 21701 448679166…511882751 6,533 1978年10月30日 Noll& Nickel
26 23209 402874115…779264511 6,987 1979年2月9日 Noll
27 44497 854509824…011228671 13,395 1979年4月8日 Nelson& Slowinski
28 86243 536927995…433438207 25,962 1982年9月25日 Slowinski
29 110503 521928313…465515007 33,265 1988年1月28日 Colquitt& Welsh
30 132049 512740276…730061311 39,751 1983年9月20日 Slowinski
31 216091 746093103…815528447 65,050 1985年9月6日 Slowinski
32 756839 174135906…544677887 227,832 1992年2月19日 Slowinski& Gage
33 859433 129498125…500142591 258,716 1994年1月10日 Slowinski& Gage
34 1257787 412245773…089366527 378,632 1996年9月3日 Slowinski& Gage
35 1398269 814717564…451315711 420,921 1996年11月13日 GIMPS/Joel Armengaud
36 2976221 623340076…729201151 895,932 1997年8月24日 GIMPS/Gordon Spence
37 3021377 127411683…024694271 909,526 1998年1月27日 GIMPS/Roland Clarkson
38 6972593 437075744…924193791 2,098,960 1999年6月1日 GIMPS/Nayan Hajratwala
39 13466917 924947738…256259071 4,053,946 2001年11月14日 GIMPS/Michael Cameron
40 20996011 125976895…855682047 6,320,430 2003年11月17日 GIMPS/Michael Shafer
41 24036583 299410429…733969407 7,235,733 2004年5月15日 GIMPS/Josh Findley
42* 25964951 122164630…577077247 7,816,230 2005年2月18日 GIMPS/Martin Nowak
43* 30402457 315416475…652943871 9,152,052 2005年12月15日 GIMPS/Curtis Cooper及Steven Boone
44* 32582657 124575026…053967871 9,808,358 2006年9月4日 GIMPS/Curtis Cooper及Steven Boone
45* 37156667 202254406…308220927 11,185,272 2008年9月6日 GIMPS/Hans-Michael Elvenich
46* 42643801 169873516…562314751 12,837,064 2009年4月12日 GIMPS/Odd M. Strindmo
47* 43112609 316470269…697152511 

代码如下：
#include<iostream>
using namespace std;

//快速乘
long long multi(long long a,long long b,long long mod)
{
    long long ans=0;

    while(b>0)
    {
        if(b&1)
            ans=(ans+a)%mod;
        b>>=1;//*****************移位运算符需要赋值*************
        a=(a<<1)%mod;
    }

    return ans;
}

long long data[100];//LUCAS序列

int main()
{

    long long sum;
    long long tmp;
    int p;
    int t;

    cin>>t;
    while(t--)
    {
        sum=1;

        cin>>p;

        if(p==2)
        {
            cout<<"yes"<<endl;
            continue;
        }

        data[1]=4;
        sum<<=p;//*****************移位运算符需要赋值*************
        sum-=1;

        //打Lucas序列的表
        for(int i=2;i<=p-1;i++)
        {
            tmp=multi(data[i-1],data[i-1],sum);
            data[i]=(tmp-2)%sum;
        }

        if(data[p-1]==0)
            cout<<"yes"<<endl;
        else
            cout<<"no"<<endl;
    }
    return 0;
}

 

 MIller-Robin素数测试法：

      费马小定理：对于素数p和任意整数a，有ap ≡ a(mod p)（同余）。反过来，满足ap ≡ a(mod p)，p也几乎一定是素数。

　　伪素数：如果n是一个正整数，如果存在和n互素的正整数a满足 an-1 ≡ 1(mod n)，我们说n是基于a的伪素数。如果一个数是伪素数，那么它几乎肯定是素数。

　　Miller-Rabin测试：不断选取不超过n-1的基b(s次)，计算是否每次都有bn-1 ≡ 1(mod n)，若每次都成立则n是素数，否则为合数。　

注意：Iller-Rabin测试是概率型的，不是确定型的，不过由于多次运行后出错的概率非常小，所以实际应用还是可行的。（一次Miller-Rabin测试其成功的概率为3/4） 

 

算法原理：

 

输入：一个大于3的奇整数n和一个大于等于1的安全参数t（用于确定测试轮数）。

 

输出：返回n是否是素数（概率意义上的，一般误判概率小于((1/2))80)即可)。

 

 

 

  1. 将n-1表示成2sr

  2. 对i从1到t做循环做以下操作：

    1. 选择一个随机整数a（2 ≤ a ≤ n−2）

    2. 计算y ← ar bmod n

    3. 如果y≠1并且y≠ n−1循环做下面的操作，否则转3：

      1. j ← 1

      2. 当j ≤ s−1并且y ≠ n−1循环做下面操作，否则跳到（iv.）

      3. 计算y ← y2 bmod n，如果y = 1返回“合数”，否则j ← j + 1

      4. 如果y ≠ n−1则返回“合数。

  3. 返回“素数”

 

    伪代码：

Function Miller-Rabin (n : longint) :boolean;
begin
    for i := 1 to s do
    begin
        a := random(n - 2) + 2;
        if mod_exp(a, n-1, n) <> 1 then return false;
    end;
    return true;
end;

 

二次探测定理：

　　如果p是奇素数，则 x2 ≡ 1(mod p)的解为 x = 1 || x = p - 1(mod p);

 

可以利用二次探测定理在实现Miller-Rabin上添加一些细节，具体实现如下：

bool miller_rabin(ll n) {
    if(n == 2 || n == 3 || n == 5 || n == 7 || n == 11)    return true;
    if(n == 1 || !(n%2) || !(n%3) || !(n%5) || !(n%7) || !(n%11))    return false;

    ll x, pre, u;
    int i, j, k = 0;
    u = n - 1;    //要求x^u % n

    while(!(u&1)) {    //如果u为偶数则u右移，用k记录移位数
        k++; u >>= 1;
    }

    srand((ll)time(0));
    for(i = 0; i < S; ++i) {    //进行S次测试
        x = rand()%(n-2) + 2;    //在[2, n)中取随机数
        if((x%n) == 0)    continue;

        x = mod_exp(x, u, n);    //先计算(x^u) % n，
        pre = x;
        for(j = 0; j < k; ++j) {    //把移位减掉的量补上，并在这地方加上二次探测
            x = mod_mul(x, x, n);
            if(x == 1 && pre != 1 && pre != n-1)    return false;    //二次探测定理，这里如果x = 1则pre 必须等于 1，或则 n-1否则可以判断不是素数
            pre = x;
        }
        if(x != 1)    return false;    //费马小定理
    }
    return true;
}

 

代码如下：

#include<iostream>
#include<time.h>
#include<cstdlib>
using namespace std;

#define TIMES 12

long long random(long long n)
{
    /*#include<cstdlib> 里 #define RAND_MAX 0x7FFF*/
    return (long long)((double)rand()/RAND_MAX*n+0.5);
}

long long multi(long long a,long long b,long long m)
{
    long long ans=0;

    while(b>0)
    {
        if(b&1)
            ans=(ans+a)%m;
        b>>=1;
        a=(a<<1)%m;
    }
    return ans;
}

long long quick_mod(long long a,long long b,long long m)
{
    long long ans=1;
    a%=m;

    while(b>0)
    {
        if(b&1)
        {
            ans=multi(ans,a,m);
            b--;
        }
        b>>=1;
        a=multi(a,a,m);
    }
    return ans;
}

bool witness(long long a,long long n)
{
    long long m=n-1;
    int j=0;

    while((m&1)==0)
    {
        j++;
        m>>=1;
    }

    long long x=quick_mod(a,m,n);

    if(x==1||x==n-1)
    {
        return false;
    }

    while(j--)
    {
        x=x*x%n;
        if(x==n-1)
            return false;
    }

    return true;
}

bool miller_rabin(long long n)
{
    if(n<2)return false;
    if(n==2)return true;
    if((n&1)==0)return false;

    for(int i=1;i<=TIMES;i++)
    {
        long long a=random(n-2)+1;
        if(witness(a,n))
            return false;
    }
    return true;
}

int main()
{
    int t,p;
    long long n;

    cin>>t;
    while(t--)
    {
        cin>>p;
        n=((long long)1<<p)-1;
        if(miller_rabin(n)==true)
            cout<<"yes"<<endl;
        else
            cout<<"no"<<endl;
    }
    return 0;
}

 

posted @ 2016-06-28 10:30  琥珀川||雨露晨曦  阅读(1319)  评论(0)    收藏  举报