免费注册 查看新帖 |

Chinaunix

  平台 论坛 博客 文库
12下一页
最近访问板块 发新帖
查看: 4754 | 回复: 19

计算 1000! 源代码 [复制链接]

论坛徽章:
0
发表于 2007-04-29 11:18 |显示全部楼层
看到论坛里好多讨论大数阶乘的话题的,我也来凑个热闹,贴一个自己前段时间刚写的。
文件内的注释已经很清楚了,懒得再用中文写一遍了,请大家看英文注释吧。
文件已包含测试代码,main函数计算了1000!。

测试结果:

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


放上计算10000!(恩,不是1000!, 是10000!)的运行时间, 供大家参考:


  1. [ykuang@qdcvs 0_UnSorted]$ time a.out
  2. indexEnd: 8914
  3. ======================== 10000! =======================
  4. 2846259680917054518906413212119868890148....
  5. /*结果太大,论坛发帖限制,这里省略,但运行时间是硬梆梆的!*/
  6. ======================== end ===========================

  7. real    0m2.196s
  8. user    0m1.730s
  9. sys     0m0.000s
复制代码


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

  1. #include <iostream>
  2. #include <string>
  3. #include <sstream>
  4. #include "LargeNumMultiply.h"

  5. using namespace std;

  6. /*
  7. * The factorial of big number: e.g. 20000!
  8. *
  9. * Parameter:
  10. *      N     -   The big number;
  11. *  factorial  -   The result of factorial.
  12. *  indexEnd  -   The valid max index of factResult array.
  13. *  baseSize  -   How many bits of factResult[i] should be saved.
  14. *                e.g. If N = 8, baseSize = 2,
  15. *                     then factResult[0] = 20
  16. *                          factResult[1] = 3, should be taken as 03 because baseSize if 2
  17. *                          factResult[2] = 4, should be taken as 04, but since indexEnd is 2, only 4 is good.
  18. *                          result = 40302
  19. *                     
  20. *                     If N = 8, baseSize = 3,
  21. *                     then factResult[0] = 320,
  22. *                          factResult[1] = 40, should be taken as 040, but since indexEnd is 1, only 40 is good.
  23. *                     
  24. *                     If N = 8, baseSize = 4,
  25. *                     then factResult[0] = 320, should be taken as 0320 because baseSize if 4.
  26. *                          factResult[1] = 4,   should be taken as 0004, but since indexEnd is 1, only 4 is good.
  27. *
  28. * Note:
  29. *    Make sure N is <= MAX_FACT.
  30. *    Because every multiply result is stored at an unsigned int var, so we need to make sure
  31. *    the every step of mutiply result should be less than UINT_MAX.
  32. * */

  33. int Factorial (unsigned int N, unsigned int *factResult, int *indexEnd, unsigned int baseSize)
  34. {
  35.     if (N > MAX_FACT)
  36.     {
  37.         cout << "N should not be larger than: " <<MAX_FACT <<endl;
  38.         return 0;
  39.     }

  40.     int valIndex = 0;
  41.     factResult[0] = 1;

  42.     unsigned int baseValue = 1;
  43.     for( unsigned int k = 0; k < baseSize; k++ )
  44.         baseValue *= 10;

  45.     /* If baseValue * N is larger than UINT_MAX, we will overflow at factResult[j] *= i;*/
  46.     if ( baseValue > (UINT_MAX/N))
  47.     {   
  48.         cout <<"baseSize: "<< baseSize<<" is too large, will make factResult[i] * N larger than UINT_MAX." <<endl;
  49.         cout <<"Please choose smaller baseSize"<<endl;
  50.         return 0;
  51.     }
  52.    
  53.     for( unsigned int i = 1; i <= N; i++ )
  54.     {   
  55.         unsigned int add = 0;
  56.         // Every time we multiply i, we should multily every section of factResult.
  57.         // valIndex stands for the lagest index of elem that is valid in factResult;
  58.         for( int j = 0; j <= valIndex ; j++ )
  59.         {   
  60.             // Multiply of result section j.
  61.             factResult[j] *= i;                 
  62.             // Add the 'add' from lower section.
  63.             factResult[j] += add;
  64.             // Caculate 'add' of section j.
  65.             add = factResult[j]/baseValue;
  66.             // Caculate value of section j.
  67.             factResult[j] %= baseValue;
  68.         }
  69.         if( add > 0 )
  70.             // if 'add' is not 0, then store 'add' to higher section.
  71.             factResult[++valIndex] = add;
  72.     }
  73.     // Store the valid index of array factResult;
  74.     *indexEnd = valIndex;
  75.     return 1;
  76. }

  77. string IntToStr (unsigned int value, unsigned int baseSize )
  78. {
  79.     ostringstream o;
  80.     string s="";
  81.     if( o<<value )
  82.         s = o.str();
  83.     else return "";
  84.     int j;
  85.     int l = baseSize - s.length();
  86.     for( j = 0; j < l; j++ )
  87.         s = "0" + s;
  88.     return s;
  89. }

  90. int main()
  91. {
  92.     unsigned int factResult[100000];
  93.     int indexEnd;
  94.     int baseSize = 4;
  95.     Factorial (1000, factResult, &indexEnd, baseSize);

  96.     cout <<"indexEnd: "<< indexEnd << endl;
  97.     cout <<"======================== start ======================="<<endl;
  98.     /*
  99.      * 8! = 40320, if baseSize is 2, then factResult:
  100.      * factResult[0] = 4
  101.      * factResult[1] = 3
  102.      * factResult[2] = 20
  103.      * Result: 40320, not 040320
  104.      * */
  105.     cout << factResult[indexEnd];
  106.     for (int i = indexEnd-1; i >= 0; i-- )
  107.         cout << IntToStr (factResult[i], baseSize);
  108.     cout<<"\n======================== end ==========================="<<endl;
  109. }

复制代码


/*LargeNumMultiply.h*/

  1. #ifndef __MYUTILSTUFF_H
  2. #define __MYUTILSTUFF_H

  3. #include <climits>
  4. #include <cmath>
  5. #include <iostream>
  6. using namespace std;

  7. const unsigned int MAX_FACT = (unsigned int)sqrt((double)UINT_MAX);
  8. int Factorial (unsigned int N, unsigned int *factResult, int *indexEnd , unsigned int baseSize = 5 );
  9. string IntToStr (unsigned int value, unsigned int baseSize);
  10. void BigNumMutiply (unsigned int N1, unsigned int N2, unsigned int *factResult, unsigned int n);

  11. #endif
复制代码

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

评分

参与人数 1可用积分 +1 收起 理由
win_hate + 1 原创内容

查看全部评分

论坛徽章:
38
2017金鸡报晓
日期:2017-02-08 10:39:4215-16赛季CBA联赛之深圳
日期:2023-02-16 14:39:0220周年集字徽章-年
日期:2022-08-31 14:25:28黑曼巴
日期:2022-08-17 18:57:0919周年集字徽章-年
日期:2022-04-25 13:02:5920周年集字徽章-20	
日期:2022-03-29 11:10:4620周年集字徽章-年
日期:2022-03-14 22:35:1820周年集字徽章-周	
日期:2022-03-09 12:51:3220周年集字徽章-年
日期:2022-02-10 13:13:4420周年集字徽章-周	
日期:2022-02-03 12:09:4420周年集字徽章-20	
日期:2022-01-25 20:14:2720周年集字徽章-周	
日期:2022-01-13 15:12:33
发表于 2007-04-29 11:28 |显示全部楼层
好贴,比GMP速度如何?

论坛徽章:
0
发表于 2007-04-29 11:33 |显示全部楼层
不知道GMP是什么,呵呵。
速度么,计算1000!几乎不需要等待。
可以自己编译测试一下。

原帖由 醉卧水云间 于 2007-4-29 11:28 发表
好贴,比GMP速度如何?

论坛徽章:
38
2017金鸡报晓
日期:2017-02-08 10:39:4215-16赛季CBA联赛之深圳
日期:2023-02-16 14:39:0220周年集字徽章-年
日期:2022-08-31 14:25:28黑曼巴
日期:2022-08-17 18:57:0919周年集字徽章-年
日期:2022-04-25 13:02:5920周年集字徽章-20	
日期:2022-03-29 11:10:4620周年集字徽章-年
日期:2022-03-14 22:35:1820周年集字徽章-周	
日期:2022-03-09 12:51:3220周年集字徽章-年
日期:2022-02-10 13:13:4420周年集字徽章-周	
日期:2022-02-03 12:09:4420周年集字徽章-20	
日期:2022-01-25 20:14:2720周年集字徽章-周	
日期:2022-01-13 15:12:33
发表于 2007-04-29 18:31 |显示全部楼层
不会吧,不知道GMP是什么?
那么你的代码计算10000!效果如何,时间多久?抗的住么?
我一般用GMP, windows和linux都好用.功能也全.

论坛徽章:
0
发表于 2007-04-29 18:35 |显示全部楼层
用 OpenSSL 的 BN 类型试一试。

论坛徽章:
0
发表于 2007-04-29 21:11 |显示全部楼层
醉卧水云间 兄说笑了,计算10000!很轻松啊:

下面这个有说服力了吧,呵呵

[ykuang@qdcvs 0_UnSorted]$ time a.out
indexEnd: 8914
======================== 10000! =======================
2846259680917054518906413212119868890148....
/*结果太大,论坛发帖限制,这里省略,但运行时间是硬梆梆的!*/
======================== end ===========================

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

原帖由 醉卧水云间 于 2007-4-29 18:31 发表
不会吧,不知道GMP是什么?
那么你的代码计算10000!效果如何,时间多久?抗的住么?
我一般用GMP, windows和linux都好用.功能也全.

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

论坛徽章:
0
发表于 2007-04-29 21:52 |显示全部楼层
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 时分组,使中间结果长度相近,然后当数比较大时,可以采用经典的快速乘法。

评分

参与人数 1可用积分 +2 收起 理由
langue + 2

查看全部评分

论坛徽章:
0
发表于 2007-04-29 21:55 |显示全部楼层
原帖由 win_hate 于 2007-4-29 21:52 发表
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 时分组,使中间结果长度相近,然后当数比较大时,可以采用经典的快速乘法。


divide-and-conquer?

.

论坛徽章:
38
2017金鸡报晓
日期:2017-02-08 10:39:4215-16赛季CBA联赛之深圳
日期:2023-02-16 14:39:0220周年集字徽章-年
日期:2022-08-31 14:25:28黑曼巴
日期:2022-08-17 18:57:0919周年集字徽章-年
日期:2022-04-25 13:02:5920周年集字徽章-20	
日期:2022-03-29 11:10:4620周年集字徽章-年
日期:2022-03-14 22:35:1820周年集字徽章-周	
日期:2022-03-09 12:51:3220周年集字徽章-年
日期:2022-02-10 13:13:4420周年集字徽章-周	
日期:2022-02-03 12:09:4420周年集字徽章-20	
日期:2022-01-25 20:14:2720周年集字徽章-周	
日期:2022-01-13 15:12:33
发表于 2007-04-29 22:02 |显示全部楼层
原帖由 softsongs 于 2007-4-29 21:11 发表
醉卧水云间 兄说笑了,计算10000!很轻松啊:

下面这个有说服力了吧,呵呵

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


回头我抽时间比较一下GMP和你的代码谁更快,谢谢!

论坛徽章:
0
发表于 2007-04-29 22:05 |显示全部楼层
原帖由 langue 于 2007-4-29 21:55 发表


divide-and-conquer?

.


嗯,那个父亲临终之前应该举计算阶乘为例子,而不应该举竹筷........
您需要登录后才可以回帖 登录 | 注册

本版积分规则 发表回复

  

北京盛拓优讯信息技术有限公司. 版权所有 京ICP备16024965号-6 北京市公安局海淀分局网监中心备案编号:11010802020122 niuxiaotong@pcpop.com 17352615567
未成年举报专区
中国互联网协会会员  联系我们:huangweiwei@itpub.net
感谢所有关心和支持过ChinaUnix的朋友们 转载本站内容请注明原作者名及出处

清除 Cookies - ChinaUnix - Archiver - WAP - TOP