計算 1000! 源代碼

softsongs

看到論壇里好多討論大數(shù)階乘的話題的，我也來湊個熱鬧，貼一個自己前段時間剛寫的。
文件內(nèi)的注釋已經(jīng)很清楚了，懶得再用中文寫一遍了，請大家看英文注釋吧。
文件已包含測試代碼，main函數(shù)計算了1000!。

測試結(jié)果：

1. 
   
2. [ykuang@qdcvs 0_UnSorted]$ a.out
   
3. indexEnd: 641
   
4. ======================== 1000! =======================
   
5. 402387260077093773543702433923003985719374864210714632543799910429938512398629020592044208486969404800479988610197196058631666872994808558901323829669944590997424504087073759918823627727188732519779505950995276120874975462497043601418278094646496291056393887437886487337119181045825783647849977012476632889835955735432513185323958463075557409114262417474349347553428646576611667797396668820291207379143853719588249808126867838374559731746136085379534524221586593201928090878297308431392844403281231558611036976801357304216168747609675871348312025478589320767169132448426236131412508780208000261683151027341827977704784635868170164365024153691398281264810213092761244896359928705114964975419909342221566832572080821333186116811553615836546984046708975602900950537616475847728421889679646244945160765353408198901385442487984959953319101723355556602139450399736280750137837615307127761926849034352625200015888535147331611702103968175921510907788019393178114194545257223865541461062892187960223838971476088506276862967146674697562911234082439208160153780889893964518263243671616762179168909779911903754031274622289988005195444414282012187361745992642956581746628302955570299024324153181617210465832036786906117260158783520751516284225540265170483304226143974286933061690897968482590125458327168226458066526769958652682272807075781391858178889652208164348344825993266043367660176999612831860788386150279465955131156552036093988180612138558600301435694527224206344631797460594682573103790084024432438465657245014402821885252470935190620929023136493273497565513958720559654228749774011413346962715422845862377387538230483865688976461927383814900140767310446640259899490222221765904339901886018566526485061799702356193897017860040811889729918311021171229845901641921068884387121855646124960798722908519296819372388642614839657382291123125024186649353143970137428531926649875337218940694281434118520158014123344828015051399694290153483077644569099073152433278288269864602789864321139083506217095002597389863554277196742822248757586765752344220207573630569498825087968928162753848863396909959826280956121450994871701244516461260379029309120889086942028510640182154399457156805941872748998094254742173582401063677404595741785160829230135358081840096996372524230560855903700624271243416909004153690105933983835777939410970027753472000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
   
6. ======================== end ===========================
   

復(fù)制代碼

放上計算10000�。ǘ�，不是1000！， 是10000�。┑倪\行時間， 供大家參考：

1. 
   
2. [ykuang@qdcvs 0_UnSorted]$ time a.out
   
3. indexEnd: 8914
   
4. ======================== 10000！ =======================
   
5. 2846259680917054518906413212119868890148....
   
6. /*結(jié)果太大，論壇發(fā)帖限制，這里省略，但運行時間是硬梆梆的！*/
   
7. ======================== end ===========================
   
8. 
   
9. real    0m2.196s
   
10. user    0m1.730s
    
11. sys     0m0.000s
    

復(fù)制代碼

/*                 LargeNumFactorial.cpp
  * #g++ LargeNumFactorial.cpp
  * #a.out
  */

1. 
   
2. #include <iostream>
   
3. #include <string>
   
4. #include <sstream>
   
5. #include "LargeNumMultiply.h"
   
6. 
   
7. using namespace std;
   
8. 
   
9. /*
   
10. * The factorial of big number: e.g. 20000!
    
11. *
    
12. * Parameter:
    
13. *      N     -   The big number;
    
14. *  factorial  -   The result of factorial.
    
15. *  indexEnd  -   The valid max index of factResult array.
    
16. *  baseSize  -   How many bits of factResult[i] should be saved.
    
17. *                e.g. If N = 8, baseSize = 2,
    
18. *                     then factResult[0] = 20
    
19. *                          factResult[1] = 3, should be taken as 03 because baseSize if 2
    
20. *                          factResult[2] = 4, should be taken as 04, but since indexEnd is 2, only 4 is good.
    
21. *                          result = 40302
    
22. *                     
    
23. *                     If N = 8, baseSize = 3,
    
24. *                     then factResult[0] = 320,
    
25. *                          factResult[1] = 40, should be taken as 040, but since indexEnd is 1, only 40 is good.
    
26. *                     
    
27. *                     If N = 8, baseSize = 4,
    
28. *                     then factResult[0] = 320, should be taken as 0320 because baseSize if 4.
    
29. *                          factResult[1] = 4,   should be taken as 0004, but since indexEnd is 1, only 4 is good.
    
30. *
    
31. * Note:
    
32. *    Make sure N is <= MAX_FACT.
    
33. *    Because every multiply result is stored at an unsigned int var, so we need to make sure
    
34. *    the every step of mutiply result should be less than UINT_MAX.
    
35. * */
    
36. 
    
37. int Factorial (unsigned int N, unsigned int *factResult, int *indexEnd, unsigned int baseSize)
    
38. {
    
39.     if (N > MAX_FACT)
    
40.     {
    
41.         cout << "N should not be larger than: " <<MAX_FACT <<endl;
    
42.         return 0;
    
43.     }
    
44. 
    
45.     int valIndex = 0;
    
46.     factResult[0] = 1;
    
47. 
    
48.     unsigned int baseValue = 1;
    
49.     for( unsigned int k = 0; k < baseSize; k++ )
    
50.         baseValue *= 10;
    
51. 
    
52.     /* If baseValue * N is larger than UINT_MAX, we will overflow at factResult[j] *= i;*/
    
53.     if ( baseValue > (UINT_MAX/N))
    
54.     {   
    
55.         cout <<"baseSize: "<< baseSize<<" is too large, will make factResult[i] * N larger than UINT_MAX." <<endl;
    
56.         cout <<"Please choose smaller baseSize"<<endl;
    
57.         return 0;
    
58.     }
    
59.    
    
60.     for( unsigned int i = 1; i <= N; i++ )
    
61.     {   
    
62.         unsigned int add = 0;
    
63.         // Every time we multiply i, we should multily every section of factResult.
    
64.         // valIndex stands for the lagest index of elem that is valid in factResult;
    
65.         for( int j = 0; j <= valIndex ; j++ )
    
66.         {   
    
67.             // Multiply of result section j.
    
68.             factResult[j] *= i;                 
    
69.             // Add the 'add' from lower section.
    
70.             factResult[j] += add;
    
71.             // Caculate 'add' of section j.
    
72.             add = factResult[j]/baseValue;
    
73.             // Caculate value of section j.
    
74.             factResult[j] %= baseValue;
    
75.         }
    
76.         if( add > 0 )
    
77.             // if 'add' is not 0, then store 'add' to higher section.
    
78.             factResult[++valIndex] = add;
    
79.     }
    
80.     // Store the valid index of array factResult;
    
81.     *indexEnd = valIndex;
    
82.     return 1;
    
83. }
    
84. 
    
85. string IntToStr (unsigned int value, unsigned int baseSize )
    
86. {
    
87.     ostringstream o;
    
88.     string s="";
    
89.     if( o<<value )
    
90.         s = o.str();
    
91.     else return "";
    
92.     int j;
    
93.     int l = baseSize - s.length();
    
94.     for( j = 0; j < l; j++ )
    
95.         s = "0" + s;
    
96.     return s;
    
97. }
    
98. 
    
99. int main()
    
100. {
     
101.     unsigned int factResult[100000];
     
102.     int indexEnd;
     
103.     int baseSize = 4;
     
104.     Factorial (1000, factResult, &indexEnd, baseSize);
     
105. 
     
106.     cout <<"indexEnd: "<< indexEnd << endl;
     
107.     cout <<"======================== start ======================="<<endl;
     
108.     /*
     
109.      * 8! = 40320, if baseSize is 2, then factResult:
     
110.      * factResult[0] = 4
     
111.      * factResult[1] = 3
     
112.      * factResult[2] = 20
     
113.      * Result: 40320, not 040320
     
114.      * */
     
115.     cout << factResult[indexEnd];
     
116.     for (int i = indexEnd-1; i >= 0; i-- )
     
117.         cout << IntToStr (factResult[i], baseSize);
     
118.     cout<<"\n======================== end ==========================="<<endl;
     
119. }
     
120. 
     

復(fù)制代碼

/*LargeNumMultiply.h*/

1. 
   
2. #ifndef __MYUTILSTUFF_H
   
3. #define __MYUTILSTUFF_H
   
4. 
   
5. #include <climits>
   
6. #include <cmath>
   
7. #include <iostream>
   
8. using namespace std;
   
9. 
   
10. const unsigned int MAX_FACT = (unsigned int)sqrt((double)UINT_MAX);
    
11. int Factorial (unsigned int N, unsigned int *factResult, int *indexEnd , unsigned int baseSize = 5 );
    
12. string IntToStr (unsigned int value, unsigned int baseSize);
    
13. void BigNumMutiply (unsigned int N1, unsigned int N2, unsigned int *factResult, unsigned int n);
    
14. 
    
15. #endif
    

復(fù)制代碼

[ 本帖最后由 softsongs 于 2007-4-29 21:14 編輯 ]

醉臥水云間

好貼，比GMP速度如何？

softsongs

不知道GMP是什么，呵呵。
速度么，計算1000！幾乎不需要等待。
可以自己編譯測試一下。

原帖由 醉臥水云間 于 2007-4-29 11:28 發(fā)表
好貼，比GMP速度如何？

醉臥水云間

不會吧,不知道GMP是什么?
那么你的代碼計算10000!效果如何,時間多久?抗的住么?
我一般用GMP, windows和linux都好用.功能也全.

langue

用 OpenSSL 的 BN 類型試一試。

softsongs

醉臥水云間 兄說笑了，計算10000！很輕松�。�

下面這個有說服力了吧，呵呵

[ykuang@qdcvs 0_UnSorted]$ time a.out
indexEnd: 8914
======================== 10000！ =======================
2846259680917054518906413212119868890148....
/*結(jié)果太大，論壇發(fā)帖限制，這里省略，但運行時間是硬梆梆的！*/
======================== end ===========================

real    0m2.196s
user    0m1.730s
sys     0m0.000s

原帖由 醉臥水云間 于 2007-4-29 18:31 發(fā)表
不會吧,不知道GMP是什么?
那么你的代碼計算10000!效果如何,時間多久?抗的住么?
我一般用GMP, windows和linux都好用.功能也全.

[ 本帖最后由 softsongs 于 2007-4-29 21:13 編輯 ]

win_hate

gmp 是這樣處理的：

1） n! 分解成 2^m p_1^a_1 p_2^a_2 ... p_r^a_r

2）2^m 只不過是移位

3）其余素因子按方冪分組，即同次的放在一起，變成 (q_1 q_2 ... q_t)^a 的形式。由于我們能快速計算方冪，所以先求出  q_1 .. q_r 的值就可以了

4）求 q_1 ... q_r 時分組，使中間結(jié)果長度相近，然后當(dāng)數(shù)比較大時，可以采用經(jīng)典的快速乘法。

langue

原帖由 win_hate 于 2007-4-29 21:52 發(fā)表
gmp 是這樣處理的：

1） n! 分解成 2^m p_1^a_1 p_2^a_2 ... p_r^a_r

2）2^m 只不過是移位

3）其余素因子按方冪分組，即同次的放在一起，變成 (q_1 q_2 ... q_t)^a 的形式。由于我們能快速計算方冪，所以先求出  q_1 .. q_r 的值就可以了

4）求 q_1 ... q_r 時分組，使中間結(jié)果長度相近，然后當(dāng)數(shù)比較大時，可以采用經(jīng)典的快速乘法。

divide-and-conquer？

.

醉臥水云間

原帖由 softsongs 于 2007-4-29 21:11 發(fā)表
醉臥水云間 兄說笑了，計算10000！很輕松啊：

下面這個有說服力了吧，呵呵

[ykuang@qdcvs 0_UnSorted]$ time a.out
indexEnd: 8914
======================== 10000！ =======================
2846259 ...

回頭我抽時間比較一下GMP和你的代碼誰更快，謝謝！

win_hate

原帖由 langue 于 2007-4-29 21:55 發(fā)表


divide-and-conquer？

.

嗯，那個父親臨終之前應(yīng)該舉計算階乘為例子，而不應(yīng)該舉竹筷........