HDU1042解题报告

by davelv Posted on 2010/02/12

一、题目规定
Time Limit:
10000/5000 MS (Java/Others)    Memory Limit: 65536/32768 K (Java/Others)
Problem Description
Given an integer N(0 ≤ N ≤ 10000), your task is to calculate N!

Input
One N in one line, process to the end of file.

Output
For each N, output N! in one line.

Sample Input
1
2
3

Sample Output
1
2
6

Author
JGShining（极光炫影）
URL
http://acm.hdu.edu.cn/showproblem.php?pid=1042

二、题目分析：
题目本身很好理解，也就是求所输出数字N的阶乘，0 ≤ N ≤ 10000 ,故 1 ≤ N! ≤ 10^35660. N!的值比较大，远远超过了普通整型的范围，虽然没有明确说，但题目也不可能弱智的让你计算出来一个浮点值就够了。这里是高精度的计算，必须精确到个位，所以我想到的就是利用普通大整数乘法来计算这个问题。
三、解题步骤：
解法：
既然确定了大整数乘法，那么就分析下精度。如果我们利用一个int32_t/uint32_t来存储整数中的一段部分的话，我们能最多能存储9位数(|UINT_MAX|>4E9，在此感谢52computer同学，我原先一直写成8位来着)，虽然前面有个2个bit空出来很可惜，但是也不能存储10位数，否则就溢出了。既然每个int32_t存储9位数字，结果最多35660位，那么我们需要约为4000个int32_t来存储数字(记为数组r)，数据模型建立OK。
结构模型如下：由于最小的阶乘结果为1，我们先r[0]=1;然后用i = 2 .. N (N>=2) 来依次跟r相乘，相乘的方法就是 j = 0 .. p-1 { i*r[j]} (p为当前r中使用中元素的个数) i*r[j]可能大于INT_MAX,所以我们需要一个int64_t (变量t)来存储中间结果，t加上上次的进位，然后除以10^9。余数就是当前r[j]的值，而商则是这次乘法的进位。 这样计算完毕后，输出很简单,从r[p-1]递减输出到r[0]即可，需要注意中间的数字如果不够9位的话就得补0。
题目代码：

/*CODE00*/
/*HDU1042-davelv-2010-02-10*/

#include <stdio.h>
#include <stdint.h>
#define    MOD    1000000000
#define LEN    4000
/*计算 n!*/
int mut(int r[], int n)
{
    int i, j, p=1;
    int64_t t = 0;
    r[0] = 1;
    for(i=2; i<=n; i++)
    {
        /*对r中每个可用元素依次和i相乘*/
        for(j=0; j<p; j++)
        {
            t = (int64_t)r[j] * i + t;/*防止r[j]*i结果被int32_t截断，先转化为int64再计算*/
            r[j] = t % MOD ;
            t = t / MOD;
        }
        /*一轮乘法过后，是否产生新的进位*/
        if (t > 0)
        {
            r[p++] = t;
            t = 0;
        }
    }
    return p;
}
/*输出*/
void print(int r[], int p)
{
    for(printf("%d",r[--p]); p--;)
    {
        printf("%09d", r[p]);/*中间数字补0*/
    }
    puts("");
}

int main()
{
    int r[LEN], n, p;/*大整数数组，当前输入的n，数组中已使用的元素个数*/

    while(scanf("%d",&n)!=EOF)
    {
        p = mut(r, n);
        print(r, p);
    }
    return 0;
}

结果：
状态：AC，典型时间：500MS， 内存：220K， 代码长度787B ，编译器GCC 。这个结果是超时时间的1/10,已是不错。 不过我们还有优化的余地。

优化1：
通过汇编源代码我CODE00们可以看到这样的情况：
C源码：

/*CODE01*/
t = (int64_t)r[j] * i + t;
r[j] = t % MOD ;
t = t / MOD;

 

对应AT&T汇编码如下：

/*CODE02*/
movl    8(%ebp), %edx
movl    -28(%ebp), %ecx
movl    -32(%ebp), %eax
movl    (%edx,%ebx,4), %esi
movl    $1000000000, 8(%esp)
movl    $0, 12(%esp)
movl    %esi, %edi
sarl    $31, %edi
imull    %edi, %eax
imull    %esi, %ecx
addl    %eax, %ecx
movl    -32(%ebp), %eax
mull    %esi
leal    (%ecx,%edx), %edi
movl    %eax, %esi
addl    -40(%ebp), %esi
adcl    -36(%ebp), %edi
movl    %esi, (%esp)
movl    %edi, 4(%esp)
call    __moddi3
movl    8(%ebp), %edx
movl    %eax, (%edx,%ebx,4)
movl    $1000000000, 8(%esp)
movl    $0, 12(%esp)
movl    %esi, (%esp)
movl    %edi, 4(%esp)
call    __divdi3
movl    %eax, -40(%ebp)
movl    %edx, -36(%ebp)

注意里面有两个汇编过程调用，仔细观察下名字发现是源码中取余和除法所造成。由于ｔ是int64_t，所以作除法前，MOD也转化为int64_t与之兼容(参考<ISO/IEC 9899:1999>6.3 Conversions小节)。所以我们的除法是64位除法，但是32位机没有64位指令，或者说是ISO规范不够智能也好，只能用汇编过程来模拟运算。而这块是时间复杂度最高的，所以汇编过程会很占用时间，我们想办法把它去掉。
如果我们都用int32_t的话就不会有如此过程了，但是上面说的int32_t可能会乘法溢出。所以我们把每个int32_t存储的位数变小，把整数数组变长，由于max(N)=1^4，所以我们的整数可以存储5位数，这样最大是10^9-1不会溢出。代码稍微修改如下:main()函数参见CODE00

/*CODE03*/
#include <stdio.h>
#include <stdint.h>
#define    MOD    100000
#define LEN    8000
/*计算 n!*/
int mut(int r[], int n)
{
    int i, j, p=1, t = 0;
    r[0] = 1;
    for(i=2; i<=n; i++)
    {
        /*对r中每个可用元素依次和i相乘*/
        for(j=0; j<p; j++)
        {
            t = r[j] * i + t;
            r[j] = t % MOD ;
            t = t / MOD;
        }
        /*一轮乘法过后，是否产生新的进位*/
        if (t > 0)
        {
            r[p++] = t;
            t = 0;
        }
    }
    return p;
}
/*输出*/
void print(int r[], int p)
{
    for(printf("%d",r[--p]); p--;)
    {
        printf("%05d", r[p]);/*中间数字补0*/
    }
    puts("");
}

结果：
状态：AC，典型时间：218MS， 内存：236K， 代码长度731B ，编译器GCC 。时间效率提高1倍，优化效果明显。

优化2:
上面的结果虽然比较好了，但是int32_t却浪费了近一倍的空间。我们看不能在最初的代码上想办法把那两个汇编过程去掉。intel指令系统对32位乘除法是把结果扩展为64的，但是由于ISO规定，或者别的原因编译器没有很好的利用这点，我们可以绕开C编译器，使用终极奥义:内嵌汇编。
intel的32乘法指令mul的大体规定是，被乘数在寄存器EAX,乘数是指定的32位其他寄存器或者内存空间，结果是64位，由高到底分别存储在EDX:EAX中。除法div则是，64位被除数在EDX:EAX中，除数是指定的32位其他寄存器或内存空间,商放到EAX，余数放到EDX。(参考 <Intel汇编语言程序设计> 附录B)。
这样我们就可以利用这两个指令对乘法和除法进行操作，也不怕溢出了。代码如下(其余部分参考CODE00，阅读代码可对照CODE00 mut())：

/*CODE04*/
int mut(int r[], int n)
{
    __asm__(" /n/t"
    "push %ebx /n/t"/*保存需要修改的寄存器,ebp已经被编译器保存过*/
    "push %ecx /n/t"
    "push %edx /n/t"
    "push %esi /n/t"
    "push %edi /n/t"/*寄存器变量对应关系ebx:i,ecx:j,esi:p,edi:t,ebp:r*/
    "movl $0, %edi /n/t"/*t=0*/
    "movl 28(%esp), %ebp /n/t"/*ebp==r*/
    "movl $1, (%ebp) /n/t"/*r[0]=1*/
    "movl $4, %esi /n/t"/*p=1,esi=4*p int32_t占4个字节，所以4为步长*/

    "movl $2, %ebx /n/t"/*i=2*/
        "CMP01:cmpl %ebx, 32(%esp) /n/t"/*i<=n?*/
        "jl RET /n/t"/*i>n 跳出外循环*/

        "movl $0, %ecx /n/t"/*j=0*/            
            "CMP02:cmpl %ecx, %esi /n/t"/*j<p*/
            "jle IF /n/t"/*j>=p 跳出内循环*/
            "movl (%ebp, %ecx), %eax /n/t" /*r[j]作为被乘数*/
            "mull %ebx /n/t"/*r[j]*i*/
            "addl %edi, %eax /n/t"/*r[j]*i+t*/
            "adcl $0, %edx /n/t"/*对32位加法产生的进位加到高位上*/
            "movl $1000000000, %edi /n/t"/*MOD*/
            "divl %edi /n/t"/*(r[j]*i+5)除以MOD*/
            "movl %edx, (%ebp, %ecx) /n/t"/*余数到r[j]*/
            "movl %eax, %edi /n/t"/*商防到t*/
            "addl $4, %ecx /n/t"/*p++*/
            "jmp CMP02 /n/t"/*跳至内循环判断条件*/

        "IF:cmpl $0, %edi /n/t"/*t==0?*/
        "je EQZERO /n/t"/*相等则跳过下面三条指令*/
        "movl %edi, (%ebp, %esi) /n/t"/*r[p]==t*/
        "addl $4, %esi /n/t"/*p++*/
        "movl $0, %edi /n/t"/*t=0*/
        "EQZERO:incl %ebx /n/t"/*i++*/
        "jmp CMP01 /n/t"/*跳至外循环判断条件*/

    "RET: shrl $2, %esi /n/t"/*esi==p*/
    "movl %esi, %eax /n/t"/*p放入eax准备返回*/
    "popl %edi /n/t"
    "popl %esi /n/t"
    "popl %edx /n/t"
    "popl %ecx /n/t"
    "popl %ebx /n/t"
    "popl %ebp /n/t"/*恢复保存的寄存器*/
    "ret");/*返回*/

    return 0;/*不可抵达语句，防止编译器报错*/
}

结果：
状态:AC,典型时间:187MS,内存：220K，代码长度:2019B 。这次优化效果不尽人意。估计是由于我ASM功力不行，不能达到时间减半的效果，对ASM感兴趣的同学可以把这个汇编码再优化。

优化3:
这种辗转相乘的方法似乎已经到达极限了，我又不会别的阶乘算法。怎么办？只能用无耻的查表法了，又称暴破法。就是把运算出来的结果作为静态数据保存到程序中，这样程序每次计算通过查找保存的数据即可得到结果，这样一次计算时间复杂度可以最小至O(1)，是不是很无耻。但是我们这里大整数数字太多不能可完全存储，于是我选了几个自认为的关键点保存，这样以后计算都在这些关键点基础上进行，虽然没有O(1)那么夸张，但是也能减少很多计算量。
关键点数据取了6个，n={1，100，500,1000,3000,6000}时候的n,p,r全部存储下来。由于我们需要根据以前的数据来计算，所以mut()，i和p不能每次都从头开始了，需要修改mut()把保存的n,p都传进去即可。代码如下(由于整数数组数据庞大，就不写在解题报告中):

/*CODE05*/
#include <stdio.h>
#include <stdint.h>

#define    LEN    4000
#define    MOD    1000000000

int mut(uint32_t r[], int n, int i ,int pos)
{
    __asm__ (" /n/t"
    "push %ebx /n/t"
    "push %ecx /n/t"
    "push %edx /n/t"
    "push %esi /n/t"
    "push %edi /n/t"
    "movl $0, %edi /n/t"
    "movl 28(%esp), %ebp /n/t"
    "movl 40(%esp), %esi /n/t"
    "shll $2, %esi /n/t"
    "movl 36(%esp), %ebx /n/t"
        "CMP01:cmpl %ebx, 32(%esp) /n/t"
        "jl RET /n/t"
        "movl $0, %ecx /n/t"
            "CMP02:cmpl %ecx, %esi /n/t"
            "jle IF /n/t"
            "movl (%ebp, %ecx), %eax /n/t"
            "mull %ebx /n/t"
            "addl %edi, %eax /n/t"
            "adcl $0, %edx /n/t"
            "movl $1000000000, %edi /n/t"
            "divl %edi /n/t"
            "movl %edx, (%ebp, %ecx) /n/t"
            "movl %eax, %edi /n/t"
            "addl $4, %ecx /n/t"
            "jmp CMP02 /n/t"
        "IF:cmpl $0, %edi /n/t"
        "je EQZERO /n/t"
        "movl %edi, (%ebp, %esi) /n/t"
        "addl $4, %esi /n/t"
        "movl $0, %edi /n/t"
        "EQZERO:incl %ebx /n/t"
        "jmp CMP01 /n/t"
        "RET: shrl $2, %esi /n/t"
        "movl %esi, %eax /n/t"
        "popl %edi /n/t"
        "popl %esi /n/t"
        "popl %edx /n/t"
        "popl %ecx /n/t"
        "popl %ebx /n/t"
        "popl %ebp /n/t"
        "ret"
        );
    return 0;
}

void print(uint32_t r[], int pos)
{
    int i = pos -1;
    printf("%d",r[i--]);
    for( ; i>=0; i--)
    {
        printf("%09d", r[i]);
    }
    puts("");
}
uint32_t ro[6][LEN]; /*6个关键点的r值，数据太长,就不给出了*/
int lo[7][2]={{0,1},{100,18},{500,127},{1000,286},{3000,1015},{6000,2230},{10001,0}};/*关键点n,p值*/
int main()
{
    uint32_t  n;
    int pos,i;
    while(scanf("%d", &n)!=EOF)
    {
        for(i=0; i<6 ; i++)/*循环判断当前n需要用的关键点*/
        {
            if(n>=lo[i][0] && n<lo[i+1][0]) //判断是否可以使用
            {    pos = mut(ro[i], n, lo[i][0]+1, lo[i][1]);
                lo[i][0] = n;
                lo[i][1] = pos;//保存计算后新关键点n,p值
                print(ro[i],pos);
                break;
            }
        }
        if(i == 6) /*没有可以使用的关键点就从1开始计算，并把结果覆盖至关键点0*/
        {
            ro[0][0]=1;
            pos = mut(ro[0], n, 2, 1);
            lo[0][0] = n;
            lo[0][1] = pos;
            print(ro[0],pos);
        }
    }
    return 0;
}

结果：
状态：AC，时间:109MS，内存：240K,代码长度：35957B。这些 HDU1042 NO.1就被我无耻的窃取了。。

杂技：这里有个HDU1042的最短代码，核心内容参考雨中飞燕同学的代码,dave修改。原理跟我们上面用的一样，就是代码太短，很难读懂，仅供娱乐。
/*CODE06*/

#define N 10000
int n;main(){l:while(scanf("%d",&n)!=-1){int s[N]={1},t=2,a=0,b=0,m=0;if(n<2){puts("1");goto l;}for(;a<=m||++t<=n&&(a=b=0,1);m==a++&&b&&m++)s[a]=(b+=s[a]*t)%N,b/=N;for(printf("%d",s[m]);m--;printf("%04d",s[m]));puts("");}}

结果：
状态：AC，时间：187MS，内存：284K，代码长度：239B。 这个代码长度是1042最短的，但是不晓得为什么系统排名没有显示。。。

至此HDU1042算是告一段落了，dave在解题的过程享受了很大的精神愉悦，在此感谢52同学，感谢CFAN编程班所有同学。上面的代码仍然有改进的余地，如i=1 to n依次相乘的方法和查表的关键点位置和个数，而且根据dave观察，HDU上的G++似乎要比GCC编译出来的速度更快些。算法优化，生生不息…

This entry was posted in 程序设计 and tagged . Bookmark the permalink.

2 Responses to HDU1042解题报告

  2. gigglesun says:
    2010/04/30 at 00:45

    没想到用C语言数组之前的知识就可以解决这么难的问题[e01]

    回复

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注