【noip2003】 麥森數

題目描述

形如2^P-1的素數稱為麥森數，這時P一定也是個素數。但反過來不一定，即如果P是個素數，2^P-1不一定也是素數。到1998年底，人們已找到了37個麥森數。最大的一個是P=3021377，它有909526位。麥森數有許多重要應用，它與完全數密切相關。

任務：從文件中輸入P（1000<P<3100000），計算2^P-1的位數和最後500位數字（用十進位高精度數表示）

輸入

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

輸出

第一行：十進位高精度數2^P-1的位數。

第2-11行：十進位高精度數2^P-1的最後500位數字。（每行輸出50位，共輸出10行，不足500位時高位補0）

不必驗證2^P-1與P是否為素數。

樣例輸入

樣例輸出

題解

首先這道題得從數論的角度入手，平時我們估計一個十進位數x(x>0)的位數，往往潛在地使用了放縮法：

若10^k<=x<10^(k+1)，則x有k+1位。而當x=2^p-1時，由於x的末位不為0，因此不存在借位現象，所以x的位數等於2^p的位數。

又假設10^k<=2^p<10^(k+1)，此中的k，便等於[log10(2^p)]也就是int(log10(2^p))。

最後簡單變換一下，第一問的答案便是int(p*log10(2))+1。

第二問若用高精加法，時間複雜度為O(n)。但常數過大，超時在所難免。

於是想到，要用快速冪高精。

只有一點不同，此處的快速冪是遞歸形式的。

當由遞歸得到2^(p/2)的值時，便可以通過2^p=2^(p/2)*2^(p/2)得到2^p的值，若p是奇數，2^p=2^(p/2)*2^(p/2)*2即可。

根據2^p=2^(p/2)*2^(p/2)，由高精得到a=b*b; if(n%2==1)

根據2^p=2^(p/2)*2^(p/2)*2，由高精得到a=b*b*2;

處理a數組進位，更新b數組；

清空a數組；

 1 #include<cstdio>
 2 #include<iostream>
 3 #include<cstring>
 4 #include<cmath>
 5 #include<algorithm>
 6 using namespace std;
 7 long long p,a[10005]={0},b[10005]={0};
 8 void bigpow(long long n){
 9     if(n==0) return;
10     bigpow(n/2);
11     if(n%2==0){
12         for(int i=1;i<=500;i++)
13             for(int j=1;j<=500;j++)
14                 a[i+j-1]=a[i+j-1]+b[i]*b[j];
15     }
16     if(n%2==1){
17         for(int i=1;i<=500;i++)
18             for(int j=1;j<=500;j++)
19                 a[i+j-1]=a[i+j-1]+b[i]*b[j]*2;
20     }
21     for(int i=1;i<=500;i++){
22         b[i]=a[i]%10;
23         a[i+1]=a[i+1]+a[i]/10;
24     }
25     memset(a,0,sizeof(a));
26 }
27 int main(){
28     ios::sync_with_stdio(false);
29     cin>>p;
30     cout<<(int)(p*log10(2)+1)<<endl;
31     b[1]=1;
32     bigpow(p);
33     for(int i=500;i>1;i--){
34         cout<<b[i];
35         if(i%50==1) cout<<endl;
36     }
37     cout<<b[1]-1<<endl;
38     return 0;
39 }