【RQNOJ】PID54 / 麦森数

题目描述

形如2^P-1的素数称为麦森数，这时P一定也是个素数。但反过来不一定，即如果P是个素数，2^P-1不一定也是素数。到1998年底，人们已找到了37个麦森数。最大的一个是P=3021377，它有909526位。麦森数有许多重要应用，它与完全数密切相关。

任务：从文件中输入P（1000

输入格式

文件中只包含一个整数P（1000

输出格式

第一行：十进制高精度数2^P-1的位数。

第2-11行：十进制高精度数2^P-1的最后500位数字。（一行输出，不足500位时高位补0）

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

分析
1. 500位的数只能用数组或字符串表示，由于涉及到计算，用数组表示更方便；
2. 如果M=10^500, 那么对于任意大的整数N，N%M 即为N的后500位数；
3. 假如T1 = N1%M，其中N1是500位以内的数，那么T1也是500位以内的数；令N2是一个32bit可以表达的整数，那么（N1*N2)%M=(T1*N2)%M
4. 对于题目中给的N=2^P（P可能是一个较大的数），令K=2^27 (32bit只能最大表达 2^32 -1，那么实际上我们只能用到2^31, 考虑到K会乘一个10以内的数，需要预留4bit，那么K的指数只能取到27=31-4）， 那么N可以表示为：N=K*K*...*K*(2^(P%27)。

    那么：
    (2^P)%M = {K*K*...*K*(2^(P%27) }%M
    将上式分步递归计算：T=（T*K）%M
    最终的T就是我们期望的2^P的后500位数

经过上面分析，获取2^P的后500位就不难实现了。

另外，我还有一个问题，那就是计算2^P的位数。
2^P是一个很大的数，怎么用10进制数表达？最常用的就是科学计数法：
    2^P=r * 10^H ( 1.0 <= r < 10.0)
那么H+1就是我们所求的位数。
对上式两边取10的对数：
    H+log10(r) = log10(2^P) = P*log10(2)
    因为0<=log10(r) < 1, 那么 H <=  P*log10(2) < H+1
    那么:
    H =floor(P*log10(2))

实现代码：

#incldue <math.h>
//Mersenne number
static void bigNumberMutil_500limit(unsigned int *rcd, int len, unsigned int key)
{
    int i;
    unsigned long long carry = 0;
    unsigned long long sum;

    for (i = 0; i < len; i++)
    {
        sum = rcd[i] * key + carry;
        rcd[i] = sum % 10;
        carry = sum / 10;
    }
}

void printMersenneNumber(int exp)
{
    int len = (int)floor(log10(2) * exp) + 1;
    unsigned int rcds[501] = { 0 };
    unsigned int mt = (unsigned int)(exp2(27));
    int i;

    rcds[0] = 1;
    while (exp >= 27)
    {
        bigNumberMutil_500limit(rcds, 501, mt);
        exp -= 27;
    }
    if (exp > 0)
    {
        mt = (unsigned int)(exp2(exp));
        bigNumberMutil_500limit(rcds, 501, mt);
    }

    if (rcds[0] > 0)
    {
        rcds[0] -= 1;
    }
    else
    {
        int carry = -1;
        int sum;

        for (i = 0; i < 500 && carry < 0; i++)
        {
            sum = rcds[i] + carry + 10;
            rcds[i] = sum % 10;
            carry = sum / 10 - 1;
        }
    }

    printf("%d\n", len);
    for (i = 499; i >= 0; i--)
    {
        printf("%c", '0' + rcds[i]);
    }
    printf("\n");
}