在古希臘時代人們就發現了4個完全數:6,28,496,8128,它們之間相差不大,但是直到15世紀才發現了第五個完全數33550336,和前面相差竟然如此遙遠,有沒有簡單法則尋找完全數呢?其實公元前300年歐幾裡得已經得出了一條定理:若2^p-1是素數(也叫質數),則2^(p-1)*(2^p-1)是完全數。因此,對于一個新的p,隻要能判斷2^p-1 是素數,就相當于找到了一個新的完全數。
梅森是17世紀歐洲科學界一位獨特的中心人物,是法蘭西學院的奠基人,被選為“100位在世界科學史上有重要地位的科學家”之一。他在歐幾裡得、費爾馬等人研究的基礎上對2^p-1型的數做了大量的計算和驗證,于1644年在他的《物理數學随感》一書中斷言:在不大于257的素數中,當p=2、3、5、7、13、17、19、31、67、127、257時,2^p-1是素數,其它都是合數。前面的7個數(即p=2、3、5、7、13、17、19)已經被前人所證實,而後面的4個數(p=31、67、127、257)則是梅森自己的推斷,由于梅森在科學界有崇高的學術地位,當時的人們對其斷言都深信不疑(後來人們發現其斷言包含若幹錯漏)。
由于梅森的工作極大地激發了人們研究2^p-1型素數的熱情,使其擺脫了作為“完全數”的附庸地位,是一個轉折點和裡程碑,為了紀念他,數學界就把2^p-1,其中p是素數的數稱為“梅森數”,并以Mp記之,即Mp=2^p-1,若p是素數,則2^p-1叫做梅森數,例如:2^2-1=3,2^3-1=7,2^5-1=31,2^7-1=127,前面4個都是素數;而2^11-1=2047=23*89是合數,2^13-1=8191,2^17-1=131071,2^19-1=52428,是素數;2^23-1=8388607=47*178481,又是合數……,如果2^p-1也是素數,則稱其為梅森素數。
分解2^p-1 的結果展示:
C語言分解梅森數程序如下:
//分解梅森數Mp=2^p-1,其中p為素數
#include <stdio.h>
#include <math.h>
#include <time.h>
#define N 20
main () //注:如果(2^p-1)是合數,一定能被(2p*k 1)整除,其中k 為自然數中的某一個
{ void shuchu(int,unsigned [N]); //輸出函數原型聲明
int i,k,x,lbz,lb=1,lcz=1,lc=1; //循環變量i,k,x;被除數總位數lbz,單元數lb,除數總位數lcz,單元數lc
int p,ss,l,jw,g=0,jr=0; //指數p,試商ss,積的單元數l,進位jw,個數g,進入标志jr
int lb1,lc1,lc2,b5,q,c3; //lb1=lb-1,lc1=lc-1,lc2=lcz*2-1,被除數前5位b5及其算術根q,除數前3位c3
unsigned b[N]={0,1},c[N]={0,1},s[N]={},y[N*2]={},xj; //被除數b,除數c,商s,餘數y,新積xj
printf("請輸入素數p:");scanf("%u",&p);
float t0=clock();
//求2^p-1:
for(i=1;i<=p;i )
{ jw=0;
for(k=1;k<=lb;k )
{ b[k]=b[k]*2 jw;
if(b[k]<=9999)jw=0; else{jw=1;b[k]-=10000;} //每4位為1個單元
}
if(jw>=1){lb ;b[lb]=1;}
}
b[1]--;
printf("2^%u-1=",p);
shuchu(lb,b); //調用輸出函數:輸出被分解數b的各單元(lb為單元數,b為數組名)
//開始分解:
printf(" = 1");
p*=2;c[1] =p; c3=c[1];//求(2*p 1)
while (lcz<=lbz)
{ lc2=lcz*2;
if(lc2<lbz||(lc2==lbz&&c3<=q))
{ lbz=lb*4;lb1=lb-1;
if(b[lb]>=1000) {b5=b[lb]*10 b[lb1]/1000;}
else if(b[lb]>=100) {lbz--;b5=b[lb]*100 b[lb1]/100;}
else if(b[lb]>=10) {lbz-=2;b5=b[lb]*1000 b[lb1]/10;}
else {lbz-=3;b5=b[lb]*10000 b[lb1];}
q=sqrt(b5 1); lc1=lc-1;
for(x=1;x<=lb;x ) {y[x]=b[x];} //做除法:
for(i=lb;i>=lc;i--)
{ y[i] =y[i 1]*10000;y[i 1]=0; s[i]=0;
while(y[i]>c[lc])
{ if(y[i]>=429496) ss=y[i]/(c[lc] 1); //2^32=42 9496 7296
else ss=(y[i]*10000 y[i-1])/(c[lc]*10000 c[lc1] 1);
if(ss==0) ss=1;
jw=0;s[i] =ss;
for(k=1;k<=lc1;k ) //lc1=lc-1
{ xj=c[k]*ss jw;
if(xj<=9999)jw=0; else{jw=xj/10000;xj%=10000;}
l=k i-lc;
if(y[l]<xj) {y[l] =10000;y[l 1]--;}
y[l]-=xj;
}
xj=c[lc]*ss jw;
y[i]-=xj;
}
}
while(y[lc]>=c[lc])
{ for(x=lc;x>=1;x--)
{ if(y[x]>c[x]) break;
if(y[x]<c[x]) goto tc;
}
s[lc] ;
for(x=1;x<=lc1;x ) //lc1=lc-1
{ if(y[x]<c[x]) {y[x] =10000;y[x 1]--;}
y[x]-=c[x];
}
y[lc]-=c[lc];
}
tc:
for(x=lc;x>=1;x--)
{ if(y[x]!=0) break;}
if(x!=0) //餘數 !=0,求新的除數:
{ c[1] =p; g ;
if(jr!=0)
{ if(gQ051!=0) //因51051=3*7*11*13*17
{ while((g%3==0||c[1]%5==0||g%7==0||g==0||g==0||g==0)==1)
{ g ;c[1] =p; }
}
else {g=1;c[1] =p;}
}
else {if((c[2]*10000 c[1])Q051==0) {jr=1;g=1;c[1] =p;}}
if(c[1]>=10000)
{ c[2] ;c[1]-=10000;
for(x=2;x<=lc;x ) { if(c[x]>=10000){c[x 1] ;c[x]-=10000;} }
if(c[lc 1]>=1) lc ;
lcz=lc*4; lc1=lc-1;
if(c[lc]>=1000) {c3=c[lc]/10;}
else if(c[lc]>=100) {lcz--;c3=c[lc];}
else if(c[lc]>=10) {lcz-=2;c3=c[lc]*10 c[lc1]/1000;}
else {lcz-=3;c3=c[lc]*100 c[lc1]/100;}
}
}
else // 餘數 =0,輸出因數:
{ printf(" * "); shuchu(lc,c); //調用輸出函數:輸出因數c的各單元
if(s[lb]==0) lb--;lb1=lb-1;
for(x=lc;x<=lb1;x )
{ if(s[x]>=10000) {s[x 1] ;s[x]-=10000;}
}
for(x=lc;x<=lb;x ) {b[x-lc1]=s[x];} //求新的被除數:
lb-=lc1; //lc1=lc-1
}
}
else //分解完了輸出最後一個因數:
{ printf(" * "); shuchu(lb,b); //調用輸出函數:輸出最後一個因數b的各單元
break;
}
}
printf("\n用時%.3f秒",(clock()-t0)/1000);
}
//輸出函數:
void shuchu(int lb,unsigned b[N])
{ int x;
printf("%u",b[lb]); lb--;
for(x=lb;x>=1;x--)
{ if(b[x]>=1000) printf(" %u",b[x]);
else if(b[x]>=100) printf(" 0%u",b[x]);
else if(b[x]>=10) printf(" 00%u",b[x]);
else printf(" 000%u",b[x]);
}
}
,更多精彩资讯请关注tft每日頭條,我们将持续为您更新最新资讯!