洛谷 NOIP2003 麦森数

[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\sim 11$ 行：十进制高精度数 $2^P-1$ 的最后 $500$ 位数字。（每行输出 $50$ 位，共输出 $10$ 行，不足 $500$ 位时高位补 $0$）

不必验证 $2^P-1$ 与 $P$ 是否为素数。

样例 #1

样例输入 #1

1279

样例输出 #1

386
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000104079321946643990819252403273640855
38615262247266704805319112350403608059673360298012
23944173232418484242161395428100779138356624832346
49081399066056773207629241295093892203457731833496
61583550472959420547689811211693677147548478866962
50138443826029173234888531116082853841658502825560
46662248318909188018470682222031405210266984354887
32958028878050869736186900714720710555703168729087

解题

分析

很明显，这道题需要用到高精度运算，但是简单的是，我们只需要运算500位即可。

计算位数

这一个问题应该才是这道题里面最难的地方（看了很多资料才堪堪找到一个非常简洁的版本）
由于我们需要计算2^n - 1的位数，而通过列举可以知道，2^n最后一位不可能为0，所以，显而易见，2^n - 1的位数就是2^n的位数。
计算k^n的位数可以使用一个万能公式，即ceil(n * log10(k))（证明就免了，懒 (我) 得 (不) 写 (会)）
所以，可以直接cout << ceil(n * log10(2)) << endl;

这里有一个小坑，就是如果直接用printf("%d\n", ...);会无法正确输出（全部WA），换成%ld或%lld都不能，需要额外定义一个有类型的变量输出
而且，两个32位2进制的数相乘可能得到64位，所以最多30位，而且要用long long

计算后500位

如果，我只是说如果，你的第一反应是用两个高精度数来做，那一定是高精度乘法已经镌刻在脑海里了，挥之不去。但实际上，我们根本不需要这么麻烦，用高精度乘一个32位整形不好吗？
而且只需要后五百位，输出需要补齐前置0，连判断位数都不要了，这么还需要两个高精度数？

可是n可以达到3100000，总不能乘这么多次吧

观察数据，每次需要乘2，而32位有符号整数的数据范围是-2^31 ~ 2^31，所以，每次乘2^30不过分吧。

优化

当乘的代码出来之后，经过观察后可得，乘上一个2^n其实就是把数字向左位移n位
如果是位移的话，在multi函数中，每一位数又不会超过100 （有进位的情况），所以我们可以大胆的位移64 - 7位。但是经过实际测试，最多可以位移59位（数据的问题吧），所以，这时间就是这么节约的（减少了接近一半的时间）
**- - - - -
其实，还有另外一种利用快速幂的算法，这就需要高精度相乘了。我很懒（这次是真的），所以就不写具体的思考方法了。需要的话可以参看这位大佬的博客：快速幂算法详解
(～￣▽￣)～

---

代码

普通乘法

#include <iostream>
#include <cstring>
#include <cstdio>
// 包括了cmath才有log10这个函数！ 
#include <cmath>

using namespace std;

typedef long long Int;

// 为了不出现命名冲突，在全局的定义的用g作为前缀标志 
Int g_r[500] = {0};
// 简化memcpy时的代码 
#define numByte sizeof(Int)*500

void multi(Int* r, int x) {
	// 利用static减少栈空间的反复创建与释放 
	static Int tmp[501] = {0};
	
	// 清零，使用的是static！ 
	tmp[0] = 0;
	for (int i = 0; i < 500; i++) {
		tmp[i] += r[i] * x;
		tmp[i + 1] = tmp[i] / 10;
		tmp[i] %= 10;
	}
	
	memcpy(r, tmp, numByte);
}

// 从后往前输出，每50位要记得换行 
void printLargeNum(Int* n) { 
	int cnt = 0;
	for (int i = 500 - 1; i >= 0; i--) {
		putchar(n[i] + '0');
		if (++cnt % 50 == 0) putchar('\n');
	}
}

int main() {
	int p = 3021377;
	cin >> p;
	cout << ceil(p * log10(2)) << endl;
	// 如果用printf,代码如注释
	// Int width = ceil(p * log10(2));
	// printf("%lld\n", width);

	// 一开始的时候要设置为乘法元，1	
	++g_r[0];
	
	// 每次乘2^30次方 
	int m30 = 1 << 30;
	while (p >= 30) {
		multi(g_r, 1 << 30);
		p -= 30;
	}
	
	// 处理剩下的位数 
	if (p) multi(g_r, 1 << p); 
	// 最后记得是求 2^n - 1 
	--g_r[0];
	printLargeNum(g_r);
}

位移优化

#include <iostream>
#include <cstring>
#include <cstdio>
// 包括了cmath才有log10这个函数！ 
#include <cmath>

using namespace std;

typedef long long Int;

// 为了不出现命名冲突，在全局的定义的用g作为前缀标志 
Int g_r[500] = {0};
// 简化memcpy时的代码 
#define numByte sizeof(Int)*500

void multi(Int* r, int x) {
	// 利用static减少栈空间的反复创建与释放 
	static Int tmp[501] = {0};
	
	// 清零，使用的是static！ 
	tmp[0] = 0;
	for (int i = 0; i < 500; i++) {
		// tmp[i] += r[i] * x;
		tmp[i] += r[i] << x;
		tmp[i + 1] = tmp[i] / 10;
		tmp[i] %= 10;
	}
	
	memcpy(r, tmp, numByte);
}

// 从后往前输出，每50位要记得换行 
void printLargeNum(Int* n) { 
	int cnt = 0;
	for (int i = 500 - 1; i >= 0; i--) {
		putchar(n[i] + '0');
		if (++cnt % 50 == 0) putchar('\n');
	}
}

int main() {
	int p = 3021377;
	cin >> p;
	cout << ceil(p * log10(2)) << endl;
	// 如果用printf,代码如注释
	// Int width = ceil(p * log10(2));
	// printf("%lld\n", width);

	// 一开始的时候要设置为乘法元，1	
	++g_r[0];
	
	int shiftBits = 54;
	while (p >= shiftBits) {
		multi(g_r, shiftBits); 
		p -= shiftBits;
	}
	
	// 处理剩下的位数 
	if (p) multi(g_r, p); 
	// 最后记得是求 2^n - 1 
	--g_r[0];
	printLargeNum(g_r);
}

快速幂法

#include <iostream>
#include <cstring>
#include <cmath>

using namespace std;

#define MAXN 1001
#define MAXI 500
#define NUMSIZE sizeof(base)
int base[MAXN], nums[MAXN], tmp[MAXN], p;

void multi() {
	// nums乘以base 
	memset(tmp, 0, NUMSIZE);
	
	for (int i = 0; i < MAXI; i++) {
		for (int j = 0; j < MAXI; j++)
			tmp[i + j] += nums[i] * base[j];
	}

	for (int i = 0; i < MAXI; i++) {
		tmp[i + 1] += tmp[i] / 10;
		tmp[i] %= 10;
	}
	
	memcpy(nums, tmp, NUMSIZE);
}

void updatebase() {
	// base * base 
	memset(tmp, 0, NUMSIZE);

	for (int i = 0; i < MAXI; i++) {
		for (int j = 0; j < MAXI; j++)
			tmp[i + j] += base[i] * base[j];
	}

	for (int i = 0; i < MAXI; i++) {
		tmp[i + 1] += tmp[i] / 10;
		tmp[i] %= 10;
	}
	
	memcpy(base, tmp, NUMSIZE);
}

void printLargeNum(int * n) { 
	int cnt = 0;
	for (int i = 500 - 1; i >= 0; i--) {
		putchar(n[i] + '0');
		if (++cnt % 50 == 0) putchar('\n');
	}
}

int main() {
	cin >> p;
	cout << ceil(p * log10(2)) << endl;
	
	// 设为乘法元，1 
	nums[0] = 1;
	// 设2的1次方为初始值 
	base[0] = 2;
	while (p > 0) {
		if (p % 2 == 1) // 需要乘一次base 
			multi();
		p /= 2;
		updatebase();
	}
	
	// 最后记得是求 2^n - 1 
	nums[0] -= 1;
	printLargeNum(nums);
	return 0;
}