【轉】[NOIP2003普及組]麥森數

來源：http://vivid.name/tech/mason.html

 

不得不紀念一下這道題，因為我今天一整天的時間都花到這道題上了。因為這道題，我學會了快速冪，學會了高精度乘高精度，學會了靜態查錯，學會了一個小小的變量的使用可能會導致整個程序掛掉。。

 

Description

形如2^P-1的素數稱為麥森數，這時P一定也是個素數。但反過來不一定，即如果P是個素數，2^P-1不一定也是素數。到1998年底，人們已找到了37個麥森數。最大的一個是P=3021377，它有909526位。麥森數有許多重要應用，它與完全數密切相關。
任務：從文件中輸入P（1000 < P < 3100000），計算2^P-1的位數和最后500位數字（用十進制高精度數表示）

Input

只包含一個整數P（1000 < P < 3100000）

Output

第一行：十進制高精度數2^P-1的位數。
第2-11行：十進制高精度數2^P-1的最后500位數字。（每行輸出50位，共輸出10行，不足500位時高位補0）
不必驗證2^P-1與P是否為素數。

Sample Input

1279

Sample Output

386
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000104079321946643990819252403273640855
38615262247266704805319112350403608059673360298012
23944173232418484242161395428100779138356624832346
49081399066056773207629241295093892203457731833496
61583550472959420547689811211693677147548478866962
50138443826029173234888531116082853841658502825560
46662248318909188018470682222031405210266984354887
32958028878050869736186900714720710555703168729087

看到這題第一感覺是簡單，很容易就看懂了。然后，就沒有然后了。不知道怎么做了，直接不斷乘2肯定不僅超時而且超出范圍。。明顯用高精度。然后不知道用高精度乘這么多次會不會超時，於是想到快速冪。

想歸想，這倆算法我都不會。於是問百度，看題解，看了好久，終於在上午似懂非懂地編出了快速冪（求a^b%m)。下午繼續研究，終於又編出了高精度乘高精度的程序。這可完全是我自己蒙出來的，我沒找到高精度乘高精度的例程。所以毫無懸念地在2^p，p上萬時掛掉了。
就是因為自己想的高*高用了t、tt以及計算時處理進位，程序在次數較高的時候才掛掉了。這三個全部改掉才可以，只改掉任何一個仍然會掛。

本題分兩問，第一問求位數，可以證明：當x有n位時，必有10^(n-1)<=x<10^n（如x有3位時必有100=10^2<=x<1000=10^3），取常用對數，n-1<=lgx<n，即lgx的整數部分是n-1，也就是說數x的位數是lg(x)的整數部分+1。故欲求x的位數只需求floor(log10(x)+1).

  

---

PS：注意：原文這里描述是“故欲求x的位數只需求floor(log10(x)+1+0.5)（floor(x+0.5)是計算x的整數部分，這樣可以避免浮點誤差）。”

其實這個浮點數x加上0.5的最主要功能是四舍五入。這里多加0.5的話，會使得有些情況下的計算結果比正確結果多1. 

---

 

然后第二問是去掉了取模運算、用上高精度的快速冪：

while(p>0)
{
    if(p==1)
    {
　　　　mul(ans,a);
　　　　break;
　　}
    else
    {
        if(p%2) mul(ans,a);
        p/=2;
        mul(a,a);
    }
}

話說我看了很多題解，都發現里面有一句“末位不要忘了減1”，一直百思不得其解。直到最后我才明白，原來是因為題目讓算2^p-1。。還有，寫這個程序的時候除了很多小錯，比如遞減的for循環寫成了i++，x[i]*y[j]寫成了x[i]*y[i]……這些都是編譯器查不出來的，只有靜態查錯，也就是傳說中的用眼睛自己看，才可以查到。從此以后我再也不敢完全依靠編譯器了。

廢話不多說，直接貼AC程序：

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>

void mul(int x[],int y[])
{
    /*   x*y->x   */
    int tmp[520]={0},lx=500,ly=500,i,j,len;
    memset(tmp,0,sizeof(tmp));
    while(x[lx]==0&&lx>0) lx--; //計算x首位位置
    while(y[ly]==0&&ly>0) ly--; //計算y首位位置
    len=lx+ly;
    for(i=1;i<=ly;i++)
       for(j=1;j<=lx;j++)
           if(i+j-1<=500) tmp[i+j-1]+=y[i]*x[j];
    for(i=1;i<=500;i++) {tmp[i+1]+=tmp[i]/10; tmp[i]%=10; if(i<500&&tmp[i+1]==0) {len=i; break;}}
    for(i=500;i>0;i--) x[i]=tmp[i];
}
int main()
{
    long p;
    int ans[501]={0},a[501]={0},i;
    scanf("%ld",&p);
    printf("%ld\n",(long)floor(p*log10(2)+1));
    /*快速冪求2^p,同時用高精度乘法*/
    ans[1]=1; a[1]=2;
    while(p>0){
    if(p==1) {mul(ans,a); break;}
    else{
    if(p%2) mul(ans,a);
    p/=2;
    mul(a,a);}
    }
    ans[1]-=1;
    for(i=500;i>0;i--) {printf("%d",ans[i]); if((i-1)%50==0) printf("\n");}
    //system("pause");
    return 0;
}

 

 

參考原文的代碼，簡化了一些地方：

#include<stdio.h>
#include<math.h>
//#include<string.h>
void mul(int x[],int y[]);
int main()
{
    long p;
    int num=0;
    int ans[501]={0},a[501]={0},i;
    freopen("Mason.in","r",stdin);
    freopen("Mason.out","w",stdout);
    scanf("%ld",&p);
    num=(int)floor(p*log10(2)+1);
    printf("%d\n",num);
    
    /*快速冪求2^p,同時用高精度乘法*/
    ans[1]=1; a[1]=2;
    while(p>0)
    {
        if(p&1)      //if(p%2==1)
            mul(ans,a);
        p=p>>1;      //p=p/2;
        mul(a,a);    //a*a -> a
    }
    ans[1]-=1;     //這個地方其實直接減1會有bug。當ans[1]為0的時候是錯誤的結果。所以可以采用下面的方法減1。但是對這個題目而言，2^p-1必然是奇數，也即個位是不可能為0。(ans[1]不會為0.) 
    /*for(i=1;i<=500;i++)
    {
        if(ans[i]>0) {ans[i]--;break;}
    }
    for(;i>=1;i--) ans[i]=9;*/
    for(i=500;i>0;i--) {printf("%d",ans[i]); if((i-1)%50==0) printf("\n");}
    printf("\n");
    return 0;
}
void mul(int x[],int y[])
{
    /*   x*y->x   */
    int tmp[520]={0},lx=500,ly=500,i,j,len;  //tem[]一定要清零。
    //memset(tmp,0,sizeof(tmp));
    while(x[lx]==0&&lx>0) lx--; //計算x首位位置
    while(y[ly]==0&&ly>0) ly--; //計算y首位位置
    len=lx+ly;
    for(i=1;i<=ly;i++)
       for(j=1;j<=lx;j++)
           if(i+j-1<=500) tmp[i+j-1]+=y[i]*x[j];   //if語句是保證只保留500位 
    for(i=1;i<=500;i++) 
    {
        tmp[i+1]+=tmp[i]/10; //把進位的值加到高位 
        tmp[i]%=10;
        /*if(i<500&&tmp[i+1]==0)
        {len=i; break;}*/  //這個地方沒理解其用意，似乎只是優化循環次數。但len沒有使用的機會，所以干脆刪除掉比較好。
    }  
    for(i=500;i>0;i--) x[i]=tmp[i];  //把結果復制到x數組 
}