梅森数是指形如2^p-1的正整数，其中指数p为素数，常记作Mp。并且当Mp也是素数时，则称Mp为梅森素数。

例如：2^2-1=3，2^3-1=7，2^5-1=31，M4=2^7-1=127，2^13-1=8191，2^17-1=131071，……等结果是素数；而2^11-1=2047=23*89，2^23-1=8388607=47*138481，2^29-1=536870911=233*1103*2089，……等结果却不是素数，所以一个梅森数是否为梅森素数，可以通过试分解检验，但是随着指数p的增大，梅森数2^p-1也将越来越大，甚至超过了计算机的处理范围。

因此，要分解梅森数必须解决好三个问题：一是能进行超大数运算（包括加、减、乘、除、乘方及分解质因数运算）；二是找出简便方法；三是尽量增加每单元存储整数的长度。

关于第一个问题，前面几篇文章中：超大数的加、减、乘、除、乘方运算已经解决了。

第二个问题请看下面：

2^11-1=2047=23*89,其中23=11*2+1,89=11*2*4+1;

2^23-1=8388607=47*178481,其中47=23*2+1，178481=23*2*3880 +1;

2^29-1=536870911=233*1103*2089，其中233=29*2*4+1,1103=29*2*19+1,2089=29*2*36+1;……

所以梅森数MP=2^p-1如果是合数，其因子一定是2p*k+1的形式（k为自然数），因此为了提高分解效率，就可以直接用（2p*k+1）去试除Mp；

第三个问题：关于增加每单元存储整数的长度，考虑到除法中被除数包含上次的余数，对32位操作系统最大能存储的无符号完整整数是9位，所以最大只能取到4位（因4*2=8，而5*2=10已经大于9了）；对64位操作系统最大存储的无符号完整整数是19位，最大只能取到9位（因9*2=18，而10*2=20已经大于19了）。

分解结果展示：

每单元存储9位整数分解梅森数程序：

//9位为1个单元，分解 Mp=2^p-1 (其中p为素数)

#include <stdio.h>

#include <math.h>

#include <time.h>

#define N 20

main ()

{ void shuchu(int,unsigned long long [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,积的单元数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 long long b[N]={0,1},c[N]={0,1},s[N]={},y[N*2]={},xj;

printf("请输入素数指数p：");scanf("%u",&p);

float t0=clock();

for(i=1;i<=p;i++)

{ jw=0;

for(k=1;k<=lb;k++)

{ b[k]=b[k]*2+jw;

if(b[k]<=999999999)jw=0; else {jw=1;b[k]-=1000000000;}

}

if(jw>=1){lb++;b[lb]=1;}

}

b[1]--;

printf("2^%u-1=",p); shuchu(lb,b); //调用输出函数：输出被分解数b的各单元

//开始分解：

printf(" = 1");

p*=2;c[1]+=p; c3=c[1];

while (lcz<=lbz)

{ lc2=lcz*2;

//做除法运算：

if(lc2<lbz||(lc2==lbz&&c3<=q))

{ lbz=lb*9; lb1=lb-1;

if(b[lb]>=100000000) {b5=b[lb]/10000;}

else if(b[lb]>=10000000) {lbz--;b5=b[lb]/1000;}

else if(b[lb]>=1000000) {lbz-=2;b5=b[lb]/100;}

else if(b[lb]>=100000) {lbz-=3;b5=b[lb]/10;}

else if(b[lb]>=10000) {lbz-=4;b5=b[lb];}

else if(b[lb]>=1000) {lbz-=5;b5=b[lb]*10+b[lb1]/100000000;}

else if(b[lb]>=100) {lbz-=6;b5=b[lb]*100+b[lb1]/10000000;}

else if(b[lb]>=10) {lbz-=7;b5=b[lb]*1000+b[lb1]/1000000;}

else {lbz-=8;b5=b[lb]*10000+b[lb1]/100000;}

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]*1000000000;y[i+1]=0; s[i]=0;

while(y[i]>c[lc])

{ if(y[i]>=18446744073) ss=y[i]/(c[lc]+1);

else ss=(y[i]*1000000000+y[i-1])/(c[lc]*1000000000+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<=999999999)jw=0; else {jw=xj/1000000000;xj%=1000000000;}

l=k+i-lc;

if(y[l]<xj) {y[l]+=1000000000;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++)

{ if(y[x]<c[x]) {y[x]+=1000000000;y[x+1]--;}

y[x]-=c[x];

}

y[lc]-=c[lc];

}

tc: //判断余数是否为0：

for(x=lc;x>=1;x--)

{ if(y[x]!=0) break;}

if(x!=0) // 1.2.1 余数 != 0,求新的除数：

{ c[1]+=p; g++;

if(jr!=0)

{ if(g%646969323!=0) //646969323=3*7*11*13*17*19*23*29

{ while((g%3==0||c[1]%5==0||g%7==0||g%11==0||g%13==0||g%17==0||g%19==0||g%23==0||g%29==0)==1)

{ g++;c[1]+=p; }

}

else {g=1;c[1]+=p;}

}

else {if(c[1]%646969323==0) {jr=1;g=1;c[1]+=p;} }

if(c[1]>=1000000000)

{ c[2]++;c[1]-=1000000000;

for(x=2;x<=lc;x++)

{ if(c[x]>=1000000000){c[x+1]++;c[x]-=1000000000;}

}

if(c[lc+1]>=1) lc++;

lcz=lc*9; lc1=lc-1;

if(c[lc]>=100000000) {c3=c[lc]/1000000;}

else if(c[lc]>=10000000) {lcz--;c3=c[lc]/100000;}

else if(c[lc]>=1000000) {lcz-=2;c3=c[lc]/10000;}

else if(c[lc]>=100000) {lcz-=3;c3=c[lc]/1000;}

else if(c[lc]>=10000) {lcz-=4;c3=c[lc]/100;}

else if(c[lc]>=1000) {lcz-=5;c3=c[lc]/10;}

else if(c[lc]>=100) {lcz-=6;c3=c[lc];}

else if(c[lc]>=10) {lcz-=7;c3=c[lc]*10+c[lc1]/100000000;}

else {lcz-=8;c3=c[lc]*100+c[lc1]/10000000;}

} }

else //输出：

{ printf(" * "); shuchu(lc,c); //调用输出函数：输出因数c的各单元

if(s[lb]==0) lb--;

for(x=lc;x<=lb-1;x++) //p=73时商有大于9位数的，所以必须处理商进位

{ if(s[x]>=1000000000) {s[x+1]++;s[x]-=1000000000;}

}

for(x=lc;x<=lb;x++) {b[x-lc1]=s[x];}

lb-=lc1;

}

}

else

{ printf(" * "); shuchu(lb,b); //调用输出函数：输出最后一个因数b的各单元

break;

}

}

printf("\n用时%.3f秒",(clock()-t0)/1000);

return 0;

}

//输出函数：

void shuchu(int lb,unsigned long long b[N])

{ int x;

printf("%u",b[lb]); lb--;

for(x=lb;x>=1;x--) printf("%9u",b[x]);

}