一、题目规定
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)
注意里面有两个汇编过程调用,仔细观察下名字发现是源码中取余和除法所造成。由于t是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编译出来的速度更快些。算法优化,生生不息…
膜拜之
没想到用C语言数组之前的知识就可以解决这么难的问题[e01]