【转】[NOIP2003普及组]麦森数

  1. 云栖社区>
  2. 博客>
  3. 正文

【转】[NOIP2003普及组]麦森数

华山青竹 2014-10-16 00:53:00 浏览658
展开阅读全文

来源: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程序:

 1 #include<stdio.h>
 2 #include<stdlib.h>
 3 #include<string.h>
 4 #include<math.h>
 5 
 6 void mul(int x[],int y[])
 7 {
 8     /*   x*y->x   */
 9     int tmp[520]={0},lx=500,ly=500,i,j,len;
10     memset(tmp,0,sizeof(tmp));
11     while(x[lx]==0&&lx>0) lx--; //计算x首位位置
12     while(y[ly]==0&&ly>0) ly--; //计算y首位位置
13     len=lx+ly;
14     for(i=1;i<=ly;i++)
15        for(j=1;j<=lx;j++)
16            if(i+j-1<=500) tmp[i+j-1]+=y[i]*x[j];
17     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;}}
18     for(i=500;i>0;i--) x[i]=tmp[i];
19 }
20 int main()
21 {
22     long p;
23     int ans[501]={0},a[501]={0},i;
24     scanf("%ld",&p);
25     printf("%ld\n",(long)floor(p*log10(2)+1));
26     /*快速幂求2^p,同时用高精度乘法*/
27     ans[1]=1; a[1]=2;
28     while(p>0){
29     if(p==1) {mul(ans,a); break;}
30     else{
31     if(p%2) mul(ans,a);
32     p/=2;
33     mul(a,a);}
34     }
35     ans[1]-=1;
36     for(i=500;i>0;i--) {printf("%d",ans[i]); if((i-1)%50==0) printf("\n");}
37     //system("pause");
38     return 0;
39 }
View Code

 

 

参考原文的代码,简化了一些地方:

 1 #include<stdio.h>
 2 #include<math.h>
 3 //#include<string.h>
 4 void mul(int x[],int y[]);
 5 int main()
 6 {
 7     long p;
 8     int num=0;
 9     int ans[501]={0},a[501]={0},i;
10     freopen("Mason.in","r",stdin);
11     freopen("Mason.out","w",stdout);
12     scanf("%ld",&p);
13     num=(int)floor(p*log10(2)+1);
14     printf("%d\n",num);
15     
16     /*快速幂求2^p,同时用高精度乘法*/
17     ans[1]=1; a[1]=2;
18     while(p>0)
19     {
20         if(p&1)      //if(p%2==1)
21             mul(ans,a);
22         p=p>>1;      //p=p/2;
23         mul(a,a);    //a*a -> a
24     }
25     ans[1]-=1;     //这个地方其实直接减1会有bug。当ans[1]为0的时候是错误的结果。所以可以采用下面的方法减1。但是对这个题目而言,2^p-1必然是奇数,也即个位是不可能为0。(ans[1]不会为0.) 
26     /*for(i=1;i<=500;i++)
27     {
28         if(ans[i]>0) {ans[i]--;break;}
29     }
30     for(;i>=1;i--) ans[i]=9;*/
31     for(i=500;i>0;i--) {printf("%d",ans[i]); if((i-1)%50==0) printf("\n");}
32     printf("\n");
33     return 0;
34 }
35 void mul(int x[],int y[])
36 {
37     /*   x*y->x   */
38     int tmp[520]={0},lx=500,ly=500,i,j,len;  //tem[]一定要清零。
39     //memset(tmp,0,sizeof(tmp));
40     while(x[lx]==0&&lx>0) lx--; //计算x首位位置
41     while(y[ly]==0&&ly>0) ly--; //计算y首位位置
42     len=lx+ly;
43     for(i=1;i<=ly;i++)
44        for(j=1;j<=lx;j++)
45            if(i+j-1<=500) tmp[i+j-1]+=y[i]*x[j];   //if语句是保证只保留500位 
46     for(i=1;i<=500;i++) 
47     {
48         tmp[i+1]+=tmp[i]/10; //把进位的值加到高位 
49         tmp[i]%=10;
50         /*if(i<500&&tmp[i+1]==0)
51         {len=i; break;}*/  //这个地方没理解其用意,似乎只是优化循环次数。但len没有使用的机会,所以干脆删除掉比较好。
52     }  
53     for(i=500;i>0;i--) x[i]=tmp[i];  //把结果复制到x数组 
54 }

 

网友评论

登录后评论
0/500
评论
华山青竹
+ 关注