HDU1042解题报告

一、题目规定
　　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), %edxmovl -28(%ebp), %ecxmovl -32(%ebp), %eaxmovl (%edx,%ebx,4), %esimovl $1000000000, 8(%esp)movl $0, 12(%esp)movl %esi, %edisarl $31, %ediimull %edi, %eaximull %esi, %ecxaddl %eax, %ecxmovl -32(%ebp), %eaxmull %esileal (%ecx,%edx), %edimovl %eax, %esiaddl -40(%ebp), %esiadcl -36(%ebp), %edimovl %esi, (%esp)movl %edi, 4(%esp)call __moddi3movl 8(%ebp), %edxmovl %eax, (%edx,%ebx,4)movl $1000000000, 8(%esp)movl $0, 12(%esp)movl %esi, (%esp)movl %edi, 4(%esp)call __divdi3movl %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.c):

/*CODE05*/#include <stdio.h>#include <stdint.h>#define LEN 4000#define MOD 1000000000int 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值，数据参考附件CODE05.c*/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 10000int 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编译出来的速度更快些。算法优化，生生不息...