ACM-麦森数

来源：互联网发布：编程开发在哪里学编辑：程序博客网时间：2024/06/09 18:19

问题描述

形如 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)

输出要求

第 1 行：十进制高精度数 2^p-1 的位数。第 2-11 行：十进制高精度数 2^p-1 的最后 500位数字。（每行输出 50 位，共输出 10 行，不足 500 位时高位补 0）
不必验证 2^p-1 与 P 是否为素数。

输入样例

输出样例

38600000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010407932194664399081925240327364085538615262247266704805319112350403608059673360298012239441732324184842421613954281007791383566248323464908139906605677320762924129509389220345773183349661583550472959420547689811211693677147548478866962501384438260291732348885311160828538416585028255604666224831890918801847068222203140521026698435488732958028878050869736186900714720710555703168729087

解题思路

第一个问题，输出 2^p-1 有多少位。由于 2^p-1 的各位数只可能是 2, 4, 6, 8 所以 2^p-1 和 2^p的位数相同。使用 C/C++标准库中在 math.h 里声明的，求以 10 为底的对数的函数 doublelog10(double x) 函数，就能轻松求得 2^p-1 的位数。

        2^p 的值需要用一种高效率的办法来算。
        显然，对于任何 p>0，考虑 p 的二进制形式，则不难得到：
        p = a₀2⁰ +a₁2¹+ a₂2²……+ a_n-12^n-1 + 2ⁿ
        这里， ai 要么是 1，要么是 0。
        因而：

计算 2^p 的办法就是，先将结果的值设为 1，计算 2¹。如果 a₀ 值为 1，则结果乘以 2¹；计算 2²，如果 a₁ 为 1，则结果乘以 2²；计算 2⁴，如果 a₂ 为 1，则结果乘以 2⁴；……总之，第 i 步(i 从 0 到 n， a_n 是 1)就计算2^2ⁱ ，如果 a_i 为 1，则结果就乘以2^2ⁱ 。每次由2^2ⁱ×2^2ⁱ就能算出2^2ⁱ⁺¹ 。由于 p 可能很大，所以上面的乘法都应该使用高精度计算。由于题目只要求输出 500 位，所以这些乘法都是只须算出末尾的 500 即可。

在前面的高精度计算中，我们用数组来存放大整数，数组的一个元素对应于十进制大整数的一位。本题如果也这样做，就会超时。为了加快计算速度，可以用一个数组元素对应于大整数的 4 位，即将大整数表示为 10000 进制，而数组中的每一个元素就存放 10000 进制数的 1 位。例如，用 int 型数组 a 来存放整数 6373384，那么只需两个数组元素就可以了，a[0]=3384, a[1]=637。

由于只要求结果的最后 500 位数字，所以我们不需要计算完整的结果，只需算出最后500 位即可。因为用每个数组元素存放十进制大整数的 4 位，所以本题中的数组最多只需要 125 个元素。

参考程序

#include <iostream>#include <cstdio>#include <memory.h>#include <cmath>#define LEN 125using namespace std;void Multiply(int *a,int *b){int i,j;int nCarry;//存放进位int nTemp;int c[LEN];//存放结果的末 500 位memset(c,0,sizeof(int)*LEN);for(i=0;i<LEN;i++){nCarry=0;for(j=0;j<LEN-i;j++){nTemp = c[i+j]+a[j]*b[i]+nCarry;c[i+j]=nTemp%10000;nCarry=nTemp/10000;//算出进位}}memcpy(a,c,LEN*sizeof(int));}int main(int argc, char *argv[]){int i;int p;int anPow[LEN];//存放不断增长的 2 的次幂int aResult[LEN];//存放最终结果的末 500 位cin >> p;cout << (int)(p*log10(2))+1 <<endl;//下面将 2 的次幂初始化为 2^(2^0)(a^b 表示 a 的 b 次方)anPow[0]=2;aResult[0]=1;for(i=1;i<LEN;i++){anPow[i]=0;aResult[i]=0;}//下面计算 2 的 p 次方while(p > 0){// p = 0 则说明 p 中的有效位都用过了，不需再算下去if(p & 1){//判断此时 p 中最低位是否为 1Multiply(aResult,anPow);}p>>=1;Multiply(anPow,anPow);}aResult[0]--;//2 的 p 次方算出后减 1//输出结果for(i=LEN-1;i>=0;i--){if(i % 25 == 12){printf("%02d\n%02d",aResult[i]/100,aResult[i]%100);}else{printf("%04d",aResult[i]);if(i%25 == 0){cout<<endl;}}}return 0;}

第13之20行：j 只要算到 LEN - i - 1，是因为 b[i]×a[j]的结果总是加到 c[i+j]上， i+j 大于等于 LEN 时， c[i+j]是不需要的，也不能要，否则 c 数组就越界了。b[i]×a[j]的结果总是要加到 c[i+j]上，此外还要再加上上次更新 c[i+j-1]时产生的进位。
第39之45行：每次执行循环都判断 a_i(i 从 0 开始)的值是否为 1, 如果是，则将最终结果乘以2^2ⁱ 。接下来再由2^2ⁱ算出2^2ⁱ⁺¹ 。

第48之57行：输出从万进制数的第 124 位开始，万进制数的每一位输出为十进制数的 4 位，每行只能输出 50 个十进制位，所以发现当 i%25 等于 12 时，第 i 个万进制位会被折行输出，其对应的后两个十进制位会跑到下一行。

0 0