梅森素数(mason)

梅森素数(mason)

题目描述

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

样例输入

1279

样例输出

386
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000104079321946643990819252403273640855
38615262247266704805319112350403608059673360298012
23944173232418484242161395428100779138356624832346
49081399066056773207629241295093892203457731833496
61583550472959420547689811211693677147548478866962
50138443826029173234888531116082853841658502825560
46662248318909188018470682222031405210266984354887
32958028878050869736186900714720710555703168729087

题解

第一小题： $2^p-1$ 是多少位数。

已知：$10^2=100$ 是最小的三位数，$10^3=1000$ 是最小的四位数，所以若 $10^x$ 是一个三位数，则 $2\le x＜3$。所以若求 $10^x$ 的位数，则求出 $x$ ，取整数部分再加一即可。所以求 $2^p-1$ 的位数，可以先求 $2^p$ 的位数，$int(lo{g}_{10}{2}^p)+1=int(p\times lo{g}_{10}2)+1$，这样计算所得是 $2^p$ 的位数，还应该判断：在减 $1$ 之后对结果有无影响。如果最后一位是 $0$ 则减 $1$ 之后需要向前借位，如果一直向前借位则可能会出现比原来的数少一位的情况。但又恰巧在本题中，底数是 $2$，指数 $p>1000$ ，则末位只可能是 $2,4,6,8$ ，不会出现需要借位的情况。所以 $2^p-1$ 的位数和 $2^p$ 是相同的，所以 $2^p-1$ 的位数就是 $int(p\times lo{g}_{10}2)+1$。

第二小题：求 $2^p-1$ 的后 $500$位。

看到这么多位数基本可以确定这 $500$ 需要利用高精度算法执行。
直接使用高精度进行多次幂的计算，得60分，循环达到 $3100000\times 500$ 次，超时。

//TLE:60
#include <iostream>
#include <cmath>
using namespace std;
const int N=507;
int a[N];
int main()
{
    int p;
    cin>>p;
    cout<<int(p*log10(2))+1<<endl;
    a[0]=1;
    while(p--) //乘上p个2
    {
        for(int i=0; i<500; ++i)
            a[i]<<=1; //每个本数乘2
        for(int i=0; i<500; ++i)
            a[i+1]+=a[i]/10, a[i]%=10; //处理进位
    }
    a[0]--; //处理减一问题（不会出现负数情况2的p次幂结尾只可能是1、2、4、8）
    for(int i=0; i<10; ++i)
    {
        for(int j=0; j<50; ++j)
            cout<<a[499-i*50-j]; //注意前面是倒着存放，所以此处要倒着输出
        cout<<endl;
    }
    return 0;
}

对于 $a^b$ 类型的问题，且 $b$ 值很大时，想到使用快速幂算法：

方法一：二分快速幂

基于上述暴力做法，需要计算有：

$2^0,2^1,2^2,...,2^p$ 每次多乘一个 $2$，因此该处求幂的过程需要计算 $p$ 次。再加上需要求高精度，因此此处直接暴力求取 $2^p$ 是超时做法。可以如何优化呢？

$a^p$ 可以写成 $a^{\frac{p}{2}}\times a^{\frac{p}{2}}$ ；而$a^{\frac{p}{2}}$ 又可以写成 $a^{\frac{p}{4}}\times a^{\frac{p}{4}}$，$a^{\frac{p}{4}}$ 又可以写成 $a^{\frac{p}{8}}\times a^{\frac{p}{8}}$；同样$a^{\frac{p}{8}}$ 又可以写成 $a^{\frac{p}{16}}\times a^{\frac{p}{16}}$ …… 直到指数变为 $0$ 时，可以直接得到 $a^0=1$。

因为每次都会被除 $2$，所以自然要想到，如果当前指数是奇数该如何处理呢？

当指数部分为奇数时（例如：$a^p$ 其中 $p$ 为奇数），$a^p=a^{\frac{p}{2}}\times a^{\frac{p}{2}}\times a$，即如果当前指数为奇数时，只需要在指数除 $2$ 之后再在最后乘一个 $a$ 即可。

综合上述思路，就可以利用递归来解决，每次返回当前指数的一半幂的结果。代码如下：

二分快速幂模板

int power(int a, int p)
{
    if(p==0) return 1; //终止条件
    int res=power(a, p/2); //利用递归每次求取
    res=res*res;
    if(b&1) res*=a; //如果当前指数为奇数时需要再乘一个底数
    return res;
}

这样的做法时间复杂度有怎样的变化呢？

原本需要求的是：$a^0,a^1,a^2,...,a^p$ ，需要计算 $p$ 次，二分快速幂每次将需要求的指数砍掉一半，那么需要求的次幂就变为：$a^0,a^2,a^4,a^8,a^{16},...,a^p$ 最终约需要计算 $lo{g}_{2}p$ 次。（中间的 $2^{\frac{p}{2}+1},2^{\frac{p}{2}+2},...,2^{p-1}$ 都无需再计算）

代码1：

//二分快速幂:100
#include <iostream>
#include <cmath>
#include <array>
using namespace std;
const int N=507;
typedef array<int, N> Ar;
Ar mul1(Ar x, Ar y) //重新定义高精度与高精度相乘
{
    Ar res{0}; //初始化res数组的第0位为0
    for(int i=0; i<500; ++i)
        for(int j=0; j<500; ++j)
        {
            if(i+j>=500) continue; //结果只需要求后500位，所以多余500的部分不用计算，直接舍弃
            res[i+j]+=x[i]*y[j]; //x的第i位与y的第j位相乘的结果刚好存储在结果数组res的第i+j位
        }
    for(int i=0; i<500; ++i) //进位处理
        res[i+1]+=res[i]/10, res[i]%=10;
    return res;
}
Ar mul2(Ar x, int y) //重新定义高精度与整型变量相乘
{
    for(int i=0; i<500; ++i) //高精度x的每一位都要对应乘上整形变量y
        x[i]*=y;
    for(int i=0; i<500; ++i) //进位处理
        x[i+1]+=x[i]/10, x[i]%=10;
    return x;
}
Ar power(int a, int p) //二分快速幂
{
    if(p==0) return Ar{1}; //终止条件
    Ar res=power(a, p>>1); //每次递归求出其一半的结果
    res=mul1(res, res); //前面的res求出的是一半的结果，完整的还要再平方一下
    if(p&1) res=mul2(res, a); //如果是奇数，要在结果最后补上一个底数a
    return res;
}
int main()
{
    int p;
    cin>>p;
    cout<<int(p*log10(2))+1<<endl; //转换为int直接舍弃小数部分所以要加一
    Ar a=power(2, p); //利用二分快速幂求2^p
    a[0]--; //别忘了减一（2^p的末位数只可能是2, 4, 6, 8）
    for(int i=0; i<10; ++i) //按行输出
    {
        for(int j=0; j<50; ++j) //每行50个数
            cout<<a[499-i*50-j]; //倒着输出（第一个坐标是0，所以最后一个数的坐标为499）
        cout<<endl;
    }
    return 0;
}

优化后的时间复杂度分析：

高精度与高精度相乘的次数是 $500\times 500=2500$；又已知其中 $p<3100000$，再结合上面计算过的二分快速幂的次数大约是 $2500\times lo{g}_{2}3100000\approx 2500\times 21.564\approx 53910$，所以最终计算次数不会超过 $53910$ 次。

方法二：快速幂

二分快速幂使用了递归算法，而递归算法会大量占用栈空间，如果出现栈空间不足的情况，递归则无法进行下去。递归是从最终目标出发，将复杂问题化为简单问题，然后再通过回溯，将简单问题的解组成复杂问题的解，递归是一种逆向过程。而递推则是从简单问题出发，一步一步向前计算，最终求得问题的结果，是一种正向的过程。很多问题，当找到由简到繁的递推关系后，既可以逆向用递归求解，也可以正向用递推求解。所以快速幂，除了上一种方法利用递归解决的二分快速幂，也可以利用递推解决快速幂算法。

由上述的二分快速幂可知，$p$ 为奇数时的递推关系与 $p$ 为偶数时的递推关系并不相同：$p$ 为奇数时，$2^p=2^{\frac{p}{2}}\times 2^{\frac{p}{2}}\times 2$；$p$ 为偶数时，$2^p=2^{\frac{p}{2}}\times 2^{\frac{p}{2}}$。所以如果由此正向递推，需要从 $2^0$ 逐步推到 $2^p$，如果用这个递推式，需要考虑 $p$ 在不断的二分过程中是分成奇数还是分成偶数，会比较麻烦。换个思路，借鉴一下二进制的思想：任何一个整数，都可以转换成二进制的表达。

例如：$a^{11}=a^1\times a^2\times a^8$，这里把 $11$ 拆成了 $1,2,8$ ，都是 $2$ 的倍数。$11$ 用二进制为 $1011_{(2)}$，这就是其中第 $0,1,3$ 位上为 $1$，所以有 $a^{11}=a^{2^0}\times a^{2^1}\times a^{2^3}$。

快速幂模板

int fast_power(int a, int p)
{
    int res=1; //用res返回结果
    while(p) //将p看成二进制，逐个处理它的最后一位
    {
        if(p&1) res*=a; //p的最后一位为1，表示当前位需要乘
        a*=a; //递推：a^2 -> a^4 -> a^8 -> a^16 ...
        p>>=1; //p右移一位，把p的最后一位去掉
    }
    return res;
} 

将这个模板用在本题中，都是高精度与高精度相乘。

代码

#include <iostream>
#include <array>
#include <cmath>
using namespace std;
const int N=507;
typedef array<int, N> Ar;
Ar mul(Ar x, Ar y) //定义函数进行高精度与高精度的相乘
{
    Ar res{0}; //初始化res数组的第0位为0
    for(int i=0; i<500; ++i)
        for(int j=0; j<500; ++j)
        {
            if(i+j>=500) continue; //结果只需要求后500位，所以多余500的部分不用计算，直接舍弃
            res[i+j]+=x[i]*y[j]; //x的第i位与y的第j位相乘的结果刚好存储在结果数组res的第i+j位
        }
    for(int i=0; i<500; ++i) //进位处理
        res[i+1]+=res[i]/10, res[i]%=10;
    return res;
}
int main()
{
    int p;
    cin>>p;
    cout<<int(p*log10(2))+1<<endl; //强制类型转换之后小数部分直接舍弃，所以最后别忘了加一
    Ar ans{1}, base{2}; //将ans答案数组的第0位初始化为1，将base底数数组的第0位初始化为2
    while(p) //循环次数：p的二进制表示的位数（log2(p)）
    {
        if(p&1) ans=mul(ans, base); //当前二进制位为1乘上当前的基数
        base=mul(base, base); //基数继续进到下一位
        p>>=1; //每次右移一位
    }
    ans[0]--; //2^p-1的减一别忘了
    for(int i=0; i<10; ++i)
    {
        for(int j=0; j<50; ++j)
            cout<<ans[499-i*50-j]; //倒序输出
        cout<<endl;
    }
    return 0;
}

时间复杂度的分析：

快速幂执行的次数是由 $p$ 的二进制位数决定，所以执行次数为 $lo{g}_{2}p$ 与二分快速幂相同，但无需用到栈空间。同上一种方法的计算，计算次数也是小于 $53910$ 的。

参考视频：
快速幂：
蓝桥杯软件类备赛：快速幂和矩阵快速幂
麦森数：
信奥赛_NOIP/CSP_算法_高精快速幂_P1045_麦森数(上)
信奥赛_NOIP/CSP_算法_高精快速幂_P1045_麦森数(中)
信奥赛_NOIP/CSP_算法_高精快速幂_P1045_麦森数(下)