ybt1925 麦森数（NOIP2003普及组第4题）（高×高）

ybt1925  麦森数（NOIP2003普及组第4题）

时空限制    1000ms/64MB

题目描述

        形如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是否为素数。

输入样例

1279

输出样例

386
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000104079321946643990819252403273640855
38615262247266704805319112350403608059673360298012
23944173232418484242161395428100779138356624832346
49081399066056773207629241295093892203457731833496
61583550472959420547689811211693677147548478866962
50138443826029173234888531116082853841658502825560
46662248318909188018470682222031405210266984354887
32958028878050869736186900714720710555703168729087

分析

        第一问：比较简单，求对数。数学原理：十进制正整数n的位数为int(log10(n))+1。所以2^p-1的位数为int(long10(2^p))+1 ==> int(long10(2)*p)+1。

        证明：当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。得证。

        第二问：要求最后500位，一定是高精度啦。如果p取1000000（1百万），2^p大约 8^(p/3)<10^(p/3)，约三十万位。如果依次乘2，在p很大的时候超时。

        这里要用到快速幂。原理：x^n=x^(n/2)*x^(n/2)*f(x),其中f(x)=1(x%2==0),或者f(x)=x(x%2==1)。模仿递归求快速幂的方法，从n=0和n=1开始计算x^n。（题目里x=2）。

 

代码

法一：数组模拟

#include<iostream>
#include<cmath>
#include<cstring>
using namespace std;
const int N = 510;
int a[N],b[N],c[N];

void mul(int a[],int b[]){
	memset(c,0,sizeof(c));
	for (int i=1; i<=500; ++i)
		for (int j=1; j<=500; ++j)
			if (i+j-1<=500) c[i+j-1]+=a[i]*b[j];
	for (int j=1; j<=500; j++)
		c[j+1]+=c[j]/10,c[j]%=10;
	copy(c,c+N,a);
}

void output(int a[]){
	int cnt=0;
	for (int i=500; i>=1; --i){
		cnt++;
		cout<<a[i];
		if (cnt%50==0) cout<<'\n';
	}
}

void mss(int n){
	memset(a,0,sizeof(a));
	memset(b,0,sizeof(b));
	a[1]=1; b[1]=2;
	while (n){
		if (n%2) mul(a,b);
		mul(b,b);
		n /= 2;
	}
	a[1]--;
	output(a);
}

int main(){
	int p;
	cin>>p;
	cout<<int(p*log10(2))+1<<endl;
	mss(p);
	return 0;
}

法二：重载运算符

#include<iostream>
#include<cmath>
#include<cstring>
using namespace std;
const int N = 505;
struct bign{
	int len,d[N];
	bign(){ memset(d,0,sizeof(d)); len=1; }	//构造函数，默认初始化
	bign(int x){							//构造函数，int初始化
		memset(d,0,sizeof(d));
		len = 0;
		do {
			d[++len] = x%10;
			x /= 10;
		} while (x);
	}
	bign operator * (const bign &b){		//重载* 高*高
		bign c;
		for (int i=1; i<=500; ++i)
			for (int j=1; j<=500; ++j)
				if (i+j-1<=500) c.d[i+j-1]+=d[i]*b.d[j];
		for (int i=1; i<=500; ++i)
			c.d[i+1]+=c.d[i]/10,c.d[i]%=10;
		return c;
	}
	bign operator *= (const bign &b){		//重载*=
		*this = *this*b;
		return *this;
	}
	friend ostream& operator << (ostream& output,const bign &x){//重载输出流
		int cnt=0;
		for (int i=500; i>1; --i){
			cnt++;
			cout<<x.d[i];
			if (cnt%50==0) cout<<'\n';
		}
		cout<<x.d[1]-1<<'\n';	//单独处理-1问题
		return output;
	}
};

int main(){
	int p;
	cin>>p;
	cout<<int(p*log10(2))+1<<endl;
	bign a=1,b=2;		//int初始化
	while (p){
		if (p%2) a*=b;	//需重载高*高、*=
		b *= b;
		p /= 2;
	}
	cout<<a<<'\n';		//需重载输出流
	return 0;
}