阶乘和素数个数

无心人

关键是计算到那个程度
需要多少时间？
我想算10^16的就很长
否则不会在70年代算不出

gxqcn

原帖由 无心人 于 2008-11-17 15:28 发表
你还是想办法把
$O(n^{2/3})$
的$pi(n)$算法
先做好了吧

这篇论文我曾上传到论坛的；在 Donald E.Knuth 著的 TAOCP 中某道习题中提及。

medie2005

呵呵，gxq说的那篇文章很好找，google一下就能找到免费的原文。
我找了一下，发现还有时间复杂度为$O(n^(3/5))$的算法来计算$pi(n)$，
不过，好像是没什么人用那个算法来计算$pi(n)$，看样子只是时间复杂度好看点而已。

medie2005

另外，有一个网站：PrimeDatabase。里面可以查到第1到第$10^12$个素数。
我对网页不熟悉，是不是能写一个脚本，每隔一段时间就查一个数，然后记录返回的素数？

medie2005

关于计算$pi(n)$，目前有时间复杂度为$O(n^((1/2)+varepsilon))$，空间复杂度为$O(n^((1/4)+varepsilon))$的算法。

至于$O(n^((2/3)+varepsilon))$的算法，我找到了intfree（张一飞，mathe应该还记得^_^）写的一个dephi版本的C++版，在我机器上，计算$pi(4*10^9)$需要0.5秒。

代码如下：

1. #include   <iostream.h>   
   
2. #include   <math.h>   
   
3. #include   <windows.h>   
   
4. 
   
5. const   unsigned   default_n=4000000000;   
   
6. const   int   maxlen=1000;   
   
7. 
   
8. int   primes[maxlen];   
   
9. int   len;   
   
10. int   sum;   
    
11. unsigned   int   n;   
    
12. int   m=7;   
    
13. int   Q=1;   
    
14. int   phiQ=1;   
    
15. int*   v;   
    
16. 
    
17. bool   string2int(unsigned   int&   n,char*   s);   
    
18. void   init(void);   
    
19. unsigned   int   phi(unsigned   int   x,int   a);   
    
20. 
    
21. void   main(void){   
    
22.         char   number_string[80];   
    
23.         unsigned   int   time;   
    
24.    
    
25.         cout<<"calc   pi(n)\n";   
    
26.         cout<<"n(default   =   "<<default_n<<")   =   ";   
    
27.         cin.get(number_string,80);   
    
28.         if(!string2int(n,number_string))n=default_n;   
    
29.    
    
30.         time=GetTickCount();   
    
31.         init();   
    
32.         int   num=(int)phi(n,len)-sum+len-1;   
    
33.         time=GetTickCount()-time;   
    
34.         GlobalFree(v);   
    
35.         cout<<"pi(n)="<<num<<"("<<time<<"ms)\n";   
    
36.         cout<<"press   any   key   to   continue...";   
    
37.         cin.get();   
    
38.         cin.get();   
    
39. }   
    
40. 
    
41. void   init(void){   
    
42.         int   max;   
    
43.         int   sqrt_max;   
    
44.         bool*   mark;   
    
45.         int   i,j;   
    
46.         int   len2,len3;   
    
47.         int   s;   
    
48.    
    
49.         max=(int)pow(1.0*n,2.0/3.0);   
    
50.         sqrt_max=(int)sqrt(max);   
    
51.         mark=(bool*)GlobalAlloc(GMEM_FIXED|GMEM_ZEROINIT,max);   
    
52.    
    
53.         for(i=2;i<sqrt_max;++i){   
    
54.                 if(!mark[i]){   
    
55.                         j=i*i;   
    
56.                         while(j<max){   
    
57.                                 mark[j]=true;   
    
58.                                 j+=i;   
    
59.                         }   
    
60.                 }   
    
61.         }   
    
62.    
    
63.         for(len=0,i=2;(unsigned)i*i*i<=n;++i){   
    
64.                 if(!mark[i])primes[++len]=i;   
    
65.         }   
    
66.    
    
67.         for(j=max-1,sum=0,s=0,len2=len;(unsigned)i*i<=n;++i){   
    
68.                 if(!mark[i]){   
    
69.                         ++len2;   
    
70.                         while((unsigned)i*j>n){   
    
71.                                 if(!mark[j])++s;   
    
72.                                 --j;   
    
73.                         }   
    
74.                         sum+=s;   
    
75.                 }   
    
76.         }   
    
77.    
    
78.         for(len3=len2;i<max;++i){   
    
79.                 if(!mark[i])++len3;   
    
80.         }   
    
81.    
    
82.         sum=(len2-len)*len3-sum;   
    
83.         sum+=len*(len-1)/2-len2*(len2-1)/2;   
    
84.    
    
85.         if(m>len)m=len;   
    
86.         for(i=1;i<=m;++i){   
    
87.                 Q*=primes[i];   
    
88.                 phiQ*=primes[i]-1;   
    
89.         }   
    
90.    
    
91.         v=(int*)GlobalAlloc(GMEM_FIXED,Q*sizeof(int));   
    
92.         for(i=0;i<Q;++i)v[i]=i;   
    
93.         for(i=1;i<=m;++i){   
    
94.                 for(j=Q-1;j>=0;--j){   
    
95.                         v[j]-=v[j/primes[i]];   
    
96.                 }   
    
97.         }   
    
98.    
    
99.         GlobalFree(mark);   
    
100. }   
     
101. 
     
102. bool   string2int(unsigned   int&   n,char*   s){   
     
103.         unsigned   int   val=0;   
     
104.         for(int   i=0;s[i];++i){   
     
105.                 if(s[i]>'9'||s[i]<'0')return   false;   
     
106.                 val*=10;   
     
107.                 val+=s[i]-'0';   
     
108.                 if(val>default_n)return   false;   
     
109.         }   
     
110.         n=val;   
     
111.         return   true;   
     
112. }   
     
113. 
     
114. unsigned   int   phi(unsigned   int   x,int   a){   
     
115.         if(a==m){   
     
116.                 return   (x/Q)*phiQ+v[x%Q];   
     
117.         }   
     
118.         if(x<(unsigned)primes[a]){   
     
119.                 return   1;   
     
120.         }   
     
121.         return   phi(x,a-1)-phi(x/primes[a],a-1);   
     
122.   }   

复制代码

另外，哪位能写个可以自动到PrimeDatabase上查询，并能记录结果的脚本？

gxqcn

原帖由 medie2005 于 2008-11-19 12:40 发表
关于计算$pi(n)$，目前有时间复杂度为$O(n^((1/2)+varepsilon))$，空间复杂度为$O(n^((1/4)+varepsilon))$的算法。

...

可否提供参考资料？

liangbch

回15楼， 张一飞我知道，张一飞是3届(2000,2001,2002)IOI国家集训队的成员，第14届国际信息学奥林匹克竞赛金牌获得者。intfree的这个代码我也看过(可惜我我看不懂), 但我不知道intfree就是张一飞。楼主怎么知道的？

gxqcn

Google 了一下：http://topic.csdn.net/t/20020827/09/972666.html

mathe

这个不难，在Linux下面直接用脚本调用wget就可以了。然后将返回的结果找关键字就可以了。

mathe

如果上面intfree的代码是$O(n^{2/3})$,既然计算到$4*10^9$只需要0.5秒，那么计算到$4*10^12$也就需要100倍，50秒左右，还是可以容忍的