分治 麦森数解题报告（转）

描述 Description        

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

任务：从文件中输入P（1000<P<3100000），计算2^P-1的位数和最后500位数字（用十进制高精度数表示）

输入格式 Input Format      

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

输出格式 Output Format    

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

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

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

样例输入 Sample Input      

1279

样例输出 Sample Output    

386

00000000000000000000000000000000000000000000000000

00000000000000000000000000000000000000000000000000

00000000000000104079321946643990819252403273640855

38615262247266704805319112350403608059673360298012

23944173232418484242161395428100779138356624832346

49081399066056773207629241295093892203457731833496

61583550472959420547689811211693677147548478866962

50138443826029173234888531116082853841658502825560

46662248318909188018470682222031405210266984354887

32958028878050869736186900714720710555703168729087

 

题目分析：

第一问是很简单的，只需要求一个对数而已，数学原理：十进制正整数n的位数为int(log10(n))+1。所以2^P-1的位数int(log10(2)*p)+1 。

第二问的关键是高精度乘法和指数幂的运算，而且由于题目要求最后500位数字，所以在计算乘法的时候我们只要求计算乘数的低500位就好了。

指数幂的运算不能硬乘，而要采用分治算法，否则就超时了。分治递归算法求指数幂是非常经典的，其数学原理是a^n = a^(n/2)*a^(n/2)*f(a)，其中f(a) = 1(n%2==0)或f(a) = a(n%2==1)。

另外我们也可以创建一个栈，记录每次执行(n /= 2)前n的值是奇数还是偶数，然后根据上面的数学原理，模仿递归的思路，从n=1或n=0开始逆向计算a^n。

采用递归算法的时候，由于存储高精度整数数组的大小是预置MAX = 1000，所以在调用递归函数的时候要按引用传递参数，否则到了后面空间就不够分配了。

为了满足“每行输出50位”的条件，我把存储高精度整数数组的元素设置成5位数，这样输出的时候只需每行输出10个元素就行了。

 

说明：

算法思想：递归分治和高精度。

数据结构：结构数组。

时间复杂度：O(LEN*LEN*lonN)，其中LEN = 500 / WIDTH + 1;

空间复杂度：O(MAX);

程序语言：分别用c++和pascal实现。

附注：分别提供了递归和非递归算法求指数幂的子函数。

关于高精度整数的四则运算请参考拙作《高精度整数运算改进版》：

http://blog.csdn.net/goal00001111/archive/2008/12/15/3522737.aspx

//PASCAL代码：PROGRAM EXAMBigInt(INPUT, OUTPUT);CONST    MAX = 1000; {整型数组的最大长度}    WIDTHMAX = 100000; {整型数组val[MAX]的元素上限}    WIDTH = 5; {输出整型数组val[MAX]的元素时的格式宽度，为方便输出，设成500的因数}    MAXLENGTH = 500; {输出数字的位数}type    BigInt = RECORD        val  : array [1..MAX] of int64; {用来存储高精度整数}        size : LongWord;  {整型数组的实际长度}    end;VAR    a, b : BigInt;    p : LongWord;PROCEDURE PrintBigInt(a : BigInt);    var        i : LongWord;    begin        while (a.size*WIDTH) < MAXLENGTH do {少则补足0}        begin            inc(a.size);            a.val[a.size] := 0;        end; {while}        while (a.size*WIDTH) > MAXLENGTH do {多则只取低500}            dec(a.size);        for i:=a.size downto 1 do        begin            if a.val[i] < 10 then                write('0000')            else if a.val[i] < 100 then                write('000')            else if a.val[i] < 1000 then                write('00')            else if a.val[i] < 10000 then                write('0');            write(a.val[i]);            if (i mod 10) = 1 then                writeln;        end; {for}     end; {PrintBigInt} FUNCTION MulBigInt(a, b : BigInt): BigInt;    const        LEN = MAXLENGTH div WIDTH;    var        c : BigInt;        i, j : LongWord;    begin        if (a.size = 1) and (a.val[1] = 0) then            MulBigInt := a        else if (b.size = 1) and (b.val[1] = 0) then            MulBigInt := b        else        begin            for i:=1 to MAX do {全部赋初值为0}                c.val[i] := 0;            i := 1;            while (i <= b.size) and (i <= LEN) do            begin                j := 1;                while (j <= a.size) and (j <= LEN) do                begin                    c.val[i+j-1] := c.val[i+j-1] + a.val[j] * b.val[i];                    c.val[i+j] := c.val[i+j] + c.val[i+j-1] div WIDTHMAX;                    c.val[i+j-1] := c.val[i+j-1] mod WIDTHMAX;                    inc(j);                end; {while}                c.size := i + j - 1;                if c.val[c.size] <> 0 then {最高位有进位}                    inc(c.size);                inc(i);            end; {while}            MulBigInt := c;        end; {else}    end; {MulBigInt}PROCEDURE PowBigInt(var c : BigInt; n : LongWord);    begin        if (n = 0) or (n = 1) then {指数为0，则幂等于1；指数为1，则幂等于底数2}        begin            c.size := 1;            c.val[1] := n + 1;            exit;        end; {if}        PowBigInt(c, n div 2); {递归求高精度整数幂}        c := MulBigInt(c, c); {a^n = a^(n/2)*a^(n/2)*f(a)}        if (n mod 2) = 1 then {其中f(a) = 1(n%2==0)或f(a) = a(n%2==1)}            c := MulBigInt(b, c);    end; {PowBigInt}{PROCEDURE PowBigInt(var c : BigInt; n : LongWord);    var        stack : array[1..MAX] of integer;        i, top : LongWord;    begin        top := 0;        while n > 0 do        begin            inc(top);            stack[top] := n mod 2;            n := n div 2;        end; {while}        c.size := 1;        c.val[1] := 1;        for i:=top downto 1 do        begin            c := MulBigInt(c, c);            if stack[i] = 1 then                c := MulBigInt(c, b);        end; {for}    end; {PowBigInt}}BEGIN    read(p);    writeln(trunc(p * ln(2) / ln(10)) + 1); {输出十进制高精度数2^P-1的位数}    b.size := 1;    b.val[1] := 2;    PowBigInt(a, p);    dec(a.val[1]);    PrintBigInt(a);    writeln;END.

//c++代码：#include<iostream>#include<math.h> using namespace std; const unsigned int MAX = 1000; //整型数组的最大长度const long long WIDTHMAX = 100000;//整型数组val[MAX]的元素上限const unsigned int WIDTH = 5; //输出整型数组val[MAX]的元素时的格式宽度，为方便输出，设成500的因数 typedef struct node{    long long val[MAX];//用来存储高精度整数    unsigned int size; //整型数组的实际长度} BigInt; void PrintBigInt(BigInt a);BigInt MulBigInt(const BigInt & a, const BigInt & b);void PowBigInt(BigInt & c, unsigned int n); int main(){    unsigned int p;    cin >> p;       cout << int(log10(2)*p) + 1 << endl;  //输出十进制高精度数2^P-1的位数       BigInt a;    PowBigInt(a, p);    a.val[0] -= 1;    PrintBigInt(a);       system("pause");       return 0;}/*函数名称：PowBigInt函数功能：递归高效算法求高精度整数幂，底数默认为2输入参数：BigInt & c：存储高精度整数幂的整型数组          unsigned int n：  指数输出参数：BigInt & c：存储高精度整数幂的整型数组 */void PowBigInt(BigInt & c, unsigned int n){    if (n == 0 || n == 1)//指数为0，则幂等于1；指数为1，则幂等于底数2    {        c.size = 1;        c.val[0] = n + 1;        return ;     }       PowBigInt(c, n/2); //递归求高精度整数幂       c = MulBigInt(c, c); //a^n = a^(n/2)*a^(n/2)*f(a)    if (n % 2 == 1)     //其中f(a) = 1(n%2==0)或f(a) = a(n%2==1)    {        BigInt b;        b.size = 1;        b.val[0] = 2; //底数默认为2        c = MulBigInt(b, c);    }}/**/ /*函数名称：PowBigInt函数功能：非递归高效算法求高精度整数幂，底数默认为2输入参数：BigInt & c：存储高精度整数幂的整型数组          unsigned int n：  指数输出参数：BigInt & c：存储高精度整数幂的整型数组 *void PowBigInt(BigInt & c, unsigned int n){    int stack[MAX] = {0};    int top = 0;    while (n > 0) //利用一个栈来存储n的状态：奇数还是偶数    {        stack[top++] = n % 2;        n /= 2;    }    BigInt b;    b.size = 1;    b.val[0] = 2; //底数默认为2    c.size = 1;    c.val[0] = 1;    for (int i=top-1; i>=0; i--)    {        c = MulBigInt(c, c);  //a^n = a^(n/2)*a^(n/2)*f(a)        if (stack[i] == 1)   //其中f(a) = 1(n%2==0)或f(a) = a(n%2==1)            c = MulBigInt(b, c);    }}/**//*函数名称：PrintBigInt函数功能：十进制高精度数2^P-1的最后500位数字，每行输出50位，共输出10行，不足500位时高位补0输入参数：BigInt a：存储高精度整数幂的整型数组输出参数：无*/void PrintBigInt(BigInt a){    while (a.size*WIDTH < 500)//少则补足0    {        a.val[a.size++] = 0;    }    while (a.size*WIDTH > 500)//多则只取低500    {        a.size--;    }    for (int i=a.size-1; i>=0; i--)    {        if (a.val[i] < 10)            cout << "0000";        else if (a.val[i] < 100)            cout << "000";         else if (a.val[i] < 1000)            cout << "00";         else if (a.val[i] < 10000)            cout << "0";              cout << a.val[i];        if (i % 10 == 0) //每行输出10个元素            cout << endl;    }}/*函数名称：MulBigInt函数功能：高精度整数乘法输入参数：const BigInt & a：用整型数组表示的高精度整数被乘数          const BigInt & b：用整型数组表示的高精度整数乘数输出参数：BigInt：返回用整型数组表示的高精度整数乘积*/BigInt MulBigInt(const BigInt & a, const BigInt & b){    if (a.size == 1 && a.val[0] == 0)        return a;    if (b.size == 1 && b.val[0] == 0)        return b;       const int LEN = 500 / WIDTH + 1;//只需取低500位相乘就好了    BigInt c;    for (int i=0; i<MAX; i++) //全部赋初值为0        c.val[i] = 0;    for (int i=0, j=0; i<b.size && i<LEN; i++)    {        for (j=0; j<a.size && j<LEN; j++)        {            c.val[i+j] += a.val[j] * b.val[i];            c.val[i+j+1] += c.val[i+j] / WIDTHMAX;            c.val[i+j] %= WIDTHMAX;        }        c.size = i + j;        if (c.val[c.size] != 0)//最高位有进位            c.size++;    }    return c;}