ACM-麦森数

问题描述

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

输入样例

1279

输出样例

386 00000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000 00000000000000104079321946643990819252403273640855 38615262247266704805319112350403608059673360298012 23944173232418484242161395428100779138356624832346 49081399066056773207629241295093892203457731833496 61583550472959420547689811211693677147548478866962 50138443826029173234888531116082853841658502825560 46662248318909188018470682222031405210266984354887 32958028878050869736186900714720710555703168729087

解题思路

        第一个问题，输出 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²ⁱ ，如果 a_i 为 1，则结果就乘以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 125 using 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]=nTemp000; 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 中最低位是否为 1 Multiply(aResult,anPow); } p>>=1; Multiply(anPow,anPow); } aResult[0]--;//2 的 p 次方算出后减 1 //输出结果 for(i=LEN-1;i>=0;i--){ if(i % 25 == 12){ printf("d\nd",aResult[i]/100,aResult[i]0); }else{ printf("d",aResult[i]); if(i% == 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 ²ⁱ⁺¹ 。

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