阶乘和圆周率

无心人

我重读这个帖子，期待对另外一个问题程序的改进

wayne

我期待这个题有一个新的突破

云梦

如果允许小数的话，问题就变得很简单了。
比如 n!=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446 E1000000000000
那么n=94851898540.1396977006676482198750221184274755718919148371815155857816611145105279808242712122064124653944534478504904850688162913

clockman

1# medie2005

做新手任务，凑字数

郭先抢

这个题目真的很有意思!

clockman

1# medie2005


难度确实挺大

gua64

我曾发过帖子说2的幂能以任何数字开始
这个问题和那帖子上的问题类似
但更难做了
无心人 发表于 2008-12-4 11:38

Matlab?

liangbch

62# 无心人

今天仔细考虑了这个题目。我认为无需使用大数库。无心人的程序中调用了大数库，而且用到对数，故速度较慢。其实这个题目虽然需要做高精度计算，但只用到小整数乘以大数的运算，完全可以自己实现。我刚刚写了一个程序，未用到任何大数库。这个程序从1开始，一直查找到2^32-1,在我的电脑上只需要3分钟左右。

说明：
    1. 我在程序中设定的精度是4个DWORD,每个DWORD表示9位10进制数，故最大精度不超过36位，若需更高精度，则只需改动一处即可。
    2. 目前的代码最多可计算到2^32-1，若需要扩大搜索范围，需要实现一个64bit数和大数相乘的算法。
    3. 我的程序依次计算每个32位整数n,求n!的科学计数法表示和pi的差，若其差小于当前的最小值，则输出n和n!.

下面是我的程序的输出结果：

1. 
   
2. 20!=2.432902008176640000*10^18
   
3. 23!=2.5852016738884976640000*10^22
   
4. 28!=3.04888344611713860501504000000*10^29
   
5. 154!=3.0897696138473508879585646636*10^271
   
6. 332!=3.1034430734498676144426322732*10^694
   
7. 612!=3.1345284031431087019726722580*10^1441
   
8. 2703!=3.1398094010536875630936179148846*10^8104
   
9. 8945!=3.139993064736104560624185498*10^31464
   
10. 17254!=3.14084642621024495446191444792*10^65612
    
11. 28726!=3.14106463595708329690285209393642766*10^115595
    
12. 50583!=3.1415379577347882738422653077918*10^215977
    
13. 527600!=3.14155838836242816120923043886*10^2789957
    
14. 780018!=3.1415890794236473559982729593497*10^4257193
    
15. 1812538!=3.141589535267260763834817102767*10^10556211
    
16. 5385973!=3.141591234607732809306850985*10^33915312
    
17. 5573998!=3.14159238184880918674222498012222406*10^35182367
    
18. 52604972!=3.1415925665166406146782547740*10^383318353
    
19. 65630447!=3.141592602756036723755404377401952*10^484537182
    
20. 187024036!=3.141592617297502030148476882134*10^1465820139
    
21. 410390476!=3.141592622630656164649753699858562*10^3356543814
    
22. 464736214!=3.14159264915921444226889681595*10^3826132373
    
23. 688395641!=3.141592651962886290907283002629*10^5784962808
    
24. It took 188 secs

复制代码

liangbch

1. #include "stdlib.h"
   
2. #include "stdio.h"
   
3. #include "string.h"
   
4. #include "windows.h"
   
5. 
   
6. typedef __int64  INT64;
   
7. 
   
8. typedef unsigned long DWORD;
   
9. 
   
10. #define DEC_RADIX 1000000000
    
11. 
    
12. #define PRECISION_DIV9 4                //10进制数表示的数字个数除以9
    
13. 
    
14. 
    
15. const DWORD PI_VALUE[][2]={
    
16.         {3,141592653},
    
17.         {31,415926535},
    
18.         {314,159265358},
    
19.         {3141,592653589},
    
20.         {31415,926535897},
    
21.         {314159,265358979},
    
22.         {3141592,653589793},
    
23.         {31415926,535897932},
    
24.         {314159265,358979323},
    
25. };
    
26. 
    
27. double scales[]=
    
28. {
    
29.         1.00,
    
30.         10.00,
    
31.         100.00,
    
32.         1000.00,
    
33.         10000.00,
    
34.         100000.00,
    
35.         1000000.00,
    
36.         10000000.00,
    
37.         100000000.00,
    
38. };
    
39. 
    
40. 
    
41. int Uint32BitCount(DWORD n)
    
42. {
    
43.         _asm
    
44.         {
    
45.                 bsr eax,n
    
46.                 inc eax
    
47.         }
    
48. }
    
49. 
    
50. //******************************************************
    
51. DWORD DEC_Array_Mul_DWORD_ALU(DWORD *a, size_t len, DWORD num)
    
52. //******************************************************
    
53. {
    
54.         static DWORD _s_TEN9=DEC_RADIX;
    
55.         //---------------------------       
    
56.         _asm
    
57.         {
    
58.                 push esi
    
59. 
    
60.                 mov ecx,len
    
61.                 mov esi,a
    
62.                 lea esi,[esi+ecx*4-4]
    
63.                 xor ecx,ecx                        //carry
    
64. 
    
65.                 test esi,3                        //the a must be 4byte align
    
66.                 jnz  thisExit
    
67. 
    
68.                 jmp cmp000
    
69. 
    
70. loop_start:
    
71.                 mov eax,[esi]
    
72.                 mul dword ptr num
    
73.                 add eax,ecx
    
74.                 adc edx,0
    
75.                 div dword ptr _s_TEN9
    
76.                 mov [esi],edx
    
77.                 mov ecx,eax
    
78.                 sub esi,4
    
79. cmp000:
    
80.                 cmp esi,a
    
81.                 jae loop_start
    
82. thisExit:
    
83.                 mov eax,ecx
    
84.                 pop esi
    
85.         }
    
86. }
    
87. 
    
88. double check_diff(DWORD num[],int shift)
    
89. {
    
90.         int t,carry;
    
91.         DWORD tArr[2];
    
92. 
    
93.         tArr[0]=PI_VALUE[shift][0];
    
94.         tArr[1]=PI_VALUE[shift][1];
    
95.        
    
96.         t=tArr[1]-num[1];
    
97.         if (t<0)
    
98.         {
    
99.                 carry=1;
    
100.                 t+=DEC_RADIX;
     
101.         }
     
102.         else
     
103.                 carry=0;
     
104. 
     
105.         tArr[1]=t;
     
106.        
     
107.         t=tArr[0]-num[0]-carry;
     
108.         tArr[0]=t;
     
109. 
     
110.         if (t<0)
     
111.         {
     
112.                 return 999;
     
113.         }
     
114.         else
     
115.         {
     
116.                 double diff;
     
117.                 diff =tArr[0] + tArr[1] * 0.000000001;
     
118.                 diff /= scales[shift];
     
119.                 return diff;
     
120. 
     
121.         }
     
122. }
     
123. 
     
124. void print_number(DWORD n,DWORD res[],int len,double exp)
     
125. {
     
126.         char buff[PRECISION_DIV9*9+9];
     
127.         char *p;
     
128.         int i,s;
     
129.         INT64 e;
     
130. 
     
131.         p=buff;
     
132.         p+=sprintf(p,"%d",res[0]);
     
133.         s=strlen(buff);
     
134. 
     
135.         if (len>1)
     
136.         {
     
137.                 for (i=1;i<len;i++)
     
138.                         p+=sprintf(p,"%09d",res[ i ]);
     
139.         }
     
140.        
     
141.         printf("\n%u!=",n);
     
142.         for (i=0;i<strlen(buff);i++)
     
143.         {
     
144.                 printf("%c",buff[ i ]);
     
145.                 if (i==0)
     
146.                         printf(".");
     
147.         }
     
148. 
     
149.         e=(INT64)(exp+(len-1)*9+s-1);
     
150.         printf("*10^%I64d\n",e);
     
151. }
     
152. 
     
153. 
     
154. void search_n()
     
155. {
     
156.         DWORD i,j;
     
157.         DWORD carry;
     
158. 
     
159.         int shift;
     
160.         int bc,fac_len;                // fac(i)的结果的数组的长度
     
161.         double diff,min_diff;
     
162.         double res_exp;
     
163.        
     
164.         DWORD fac[PRECISION_DIV9+1]; //fac(i)的结果的数组
     
165.         DWORD time;
     
166. 
     
167.        
     
168.         time=GetTickCount();
     
169.         min_diff=1;
     
170.         res_exp=0;
     
171.         fac[0]=1;
     
172.         fac_len=1;
     
173.        
     
174.         i=1;
     
175.         //while (i<DEC_RADIX)
     
176.         while (i!=0)
     
177.         {
     
178.                 carry=DEC_Array_Mul_DWORD_ALU(fac,fac_len,i);
     
179.                 if (carry>=DEC_RADIX)
     
180.                 {
     
181.                         for (j=PRECISION_DIV9-2;j>=2;j--)
     
182.                                 fac[j]=fac[j-2];
     
183.                         fac[0]=carry / DEC_RADIX;
     
184.                         fac[1]=carry % DEC_RADIX;
     
185.                        
     
186.                         if ( fac_len+2 <= PRECISION_DIV9)
     
187.                                 fac_len+=2;
     
188.                         else if (fac_len+1 ==PRECISION_DIV9)
     
189.                         {
     
190.                                 fac_len=PRECISION_DIV9;
     
191.                                 res_exp+=9;
     
192.                         }
     
193.                         else
     
194.                                 res_exp+=18;
     
195.                 }
     
196.                 else if (carry>0)
     
197.                 {
     
198.                         for (j=PRECISION_DIV9-1;j>=1;j--)
     
199.                                 fac[j]=fac[j-1];
     
200.                         fac[0]=carry;
     
201.                         if (fac_len<PRECISION_DIV9)
     
202.                                 fac_len++;
     
203.                         else
     
204.                                 res_exp+=9;
     
205.                 }
     
206.                 bc=Uint32BitCount(fac[0]);
     
207.                
     
208.                 switch (bc)
     
209.                 {
     
210.                 case 2:
     
211.                          shift=0;
     
212.                          diff=check_diff(fac,shift);
     
213.                          if (diff<min_diff)
     
214.                          {
     
215.                                  min_diff=diff;
     
216.                                  if (fac_len>=2)
     
217.                                          print_number(i,fac,fac_len,res_exp);
     
218.                          }
     
219.                          break;
     
220.                 case 5:
     
221.                         shift=1;
     
222.                         diff=check_diff(fac,shift);
     
223.                         if (diff<min_diff)
     
224.                         {
     
225.                                  min_diff=diff;
     
226.                                  if (fac_len>=2)
     
227.                                          print_number(i,fac,fac_len,res_exp);
     
228.                          }
     
229.                          break;
     
230.                 case 9:
     
231.                          shift=2;
     
232.                         diff=check_diff(fac,shift);
     
233.                          if (diff<min_diff)
     
234.                          {
     
235.                                  min_diff=diff;
     
236.                                  if (fac_len>=2)
     
237.                                          print_number(i,fac,fac_len,res_exp);
     
238.                          }
     
239.                          break;
     
240.                 case 12:
     
241.                          shift=3;
     
242.                         diff=check_diff(fac,shift);
     
243.                          if (diff<min_diff)
     
244.                          {
     
245.                                  min_diff=diff;
     
246.                                  if (fac_len>=2)
     
247.                                          print_number(i,fac,fac_len,res_exp);
     
248.                          }
     
249.                          break;
     
250.                 case 15:
     
251.                          shift=4;
     
252.                         diff=check_diff(fac,shift);
     
253.                          if (diff<min_diff)
     
254.                          {
     
255.                                  min_diff=diff;
     
256.                                  if (fac_len>=2)
     
257.                                          print_number(i,fac,fac_len,res_exp);
     
258.                          }
     
259.                          break;
     
260.                 case 19:
     
261.                          shift=5;
     
262.                         diff=check_diff(fac,shift);
     
263.                          if (diff<min_diff)
     
264.                          {
     
265.                                  min_diff=diff;
     
266.                                  if (fac_len>=2)
     
267.                                          print_number(i,fac,fac_len,res_exp);
     
268.                          }
     
269.                          break;
     
270.                 case 22:
     
271.                          shift=6;
     
272.                         diff=check_diff(fac,shift);
     
273.                          if (diff<min_diff)
     
274.                          {
     
275.                                  min_diff=diff;
     
276.                                  if (fac_len>=2)
     
277.                                          print_number(i,fac,fac_len,res_exp);
     
278.                          }
     
279.                          break;
     
280.                 case 25:
     
281.                          shift=7;
     
282.                         diff=check_diff(fac,shift);
     
283.                          if (diff<min_diff)
     
284.                          {
     
285.                                  min_diff=diff;
     
286.                                  if (fac_len>=2)
     
287.                                          print_number(i,fac,fac_len,res_exp);
     
288.                          }
     
289.                          break;
     
290.                 case 29:
     
291.                         shift=8;
     
292.                         diff=check_diff(fac,shift);
     
293.                          if (diff<min_diff)
     
294.                          {
     
295.                                  min_diff=diff;
     
296.                                  if (fac_len>=2)
     
297.                                          print_number(i,fac,fac_len,res_exp);
     
298.                          }
     
299.                          break;
     
300.                 default:
     
301.                         break;
     
302.                 }
     
303.                 i++;
     
304.         }
     
305.        
     
306.         time=(GetTickCount()-time);
     
307.         printf("It took %d secs\n",time/1000);
     
308. 
     
309. }
     
310. 
     
311. int main(int argc, char* argv[])
     
312. {
     
313.         search_n();
     
314.         return 0;
     
315. }

复制代码

liangbch

刚才仔细看了一下无心人的代码，明白了其所用的算法。其主要算法如下，
从start开始，逐个计算每个整数n的阶乘的常用对数的小数部分，记作log(fac(n)), 和PI的常用对数的小数部分比较，符合条件则输出。
    当n是1000的整数倍时，计算base=n, lf0=log10(fac(base))的小数部分，
    若n不是base的整数倍时，则计算从base+1到n之间所有数的累乘积的自然对数的小数部分lf1, 则 n!的对数的小数部分为 lf= lf0+lf1。

无心人的算法主要受制于 c语言库中(或者CPU浮点指令)中的双精度数的对数函数精度的制约，很难提高精度。


受无心人的启发，我打算编写一个不使用求双精度对数和多重精度对数函数的程序。
  主要特点有：
  1. 可计算到10^15或者更高。
  2. 可分段计算。
  3. 当n小于10亿时，使用自己写的多重精度大数乘以32bit整数的子程序,计算n的阶乘，fac(n),以10^9进制表示.
     当n>=10亿且n是100万的倍数时，使用斯特林公式的普通形式(而不是指数形式)计算 fac(n)，以10^9进制表示,记作fac_0
     当n>10亿且n不是100万的倍数时，令n=base+d(base是100万的整数倍，d>0 且 d<base)，使用自己写的多重精度大数乘以64bit整数的子程序,计算base+1到n的累乘积，记作fac_1, 则n的阶乘fac(n)=fac0*fac1

  我的算法和无心人算法的差别，我的记作B，无心人的记作A。
    1. B使用斯特林公式的普通形式，而A使用斯特林公式的对数形式。
    2. B不计算n的阶乘的常用对数，而是直接计算n的阶乘的值，以10亿进制的表示形式。
    3. A的核心部分为计算64bit整数的的常用对数，使用C语言函数库，A受制于double型数的精度，无法提高精度。 B的和核心部分计算一个多精度数与64bit整数的乘积，可通过增加多精度数数的长度来提高精度。
    4. 性能的PK. 主要看 多精度数与64bit整数的乘积 VS log10(double),那个算得更快。