阶乘和素数个数

无心人

关键是计算到那个程度 需要多少时间？ 我想算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>
5. const unsigned default_n=4000000000;
6. const int maxlen=1000;
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;
17. bool string2int(unsigned int& n,char* s);
18. void init(void);
19. unsigned int phi(unsigned int x,int a);
21. void main(void){
22. char number_string[80];
23. unsigned int time;
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;
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. }
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;
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);
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. }
63. for(len=0,i=2;(unsigned)i*i*i<=n;++i){
64. if(!mark[i])primes[++len]=i;
65. }
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. }
78. for(len3=len2;i<max;++i){
79. if(!mark[i])++len3;
80. }
82. sum=(len2-len)*len3-sum;
83. sum+=len*(len-1)/2-len2*(len2-1)/2;
85. if(m>len)m=len;
86. for(i=1;i<=m;++i){
87. Q*=primes[i];
88. phiQ*=primes[i]-1;
89. }
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. }
99. GlobalFree(mark);
100. }
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. }
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秒左右，还是可以容忍的