无心人
发表于 2008-11-17 19:22:53
:b:
关键是计算到那个程度
需要多少时间?
我想算10^16的就很长
否则不会在70年代算不出
gxqcn
发表于 2008-11-17 20:30:12
原帖由 无心人 于 2008-11-17 15:28 发表 http://bbs.emath.ac.cn/images/common/back.gif
你还是想办法把
$O(n^{2/3})$
的$pi(n)$算法
先做好了吧
这篇论文我曾上传到论坛的;在 Donald E.Knuth 著的 TAOCP 中某道习题中提及。
medie2005
发表于 2008-11-17 20:53:28
呵呵,gxq说的那篇文章很好找,google一下就能找到免费的原文。
我找了一下,发现还有时间复杂度为O(n^(3/5))的算法来计算pi(n),
不过,好像是没什么人用那个算法来计算pi(n),看样子只是时间复杂度好看点而已。
medie2005
发表于 2008-11-17 21:03:51
另外,有一个网站:PrimeDatabase。里面可以查到第1到第10^12个素数。
我对网页不熟悉,是不是能写一个脚本,每隔一段时间就查一个数,然后记录返回的素数?
medie2005
发表于 2008-11-19 12:40:53
关于计算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秒。
代码如下:#include <iostream.h>
#include <math.h>
#include <windows.h>
const unsigned default_n=4000000000;
const int maxlen=1000;
int primes;
int len;
int sum;
unsigned int n;
int m=7;
int Q=1;
int phiQ=1;
int* v;
bool string2int(unsigned int& n,char* s);
void init(void);
unsigned int phi(unsigned int x,int a);
void main(void){
char number_string;
unsigned int time;
cout<<"calc pi(n)\n";
cout<<"n(default = "<<default_n<<") = ";
cin.get(number_string,80);
if(!string2int(n,number_string))n=default_n;
time=GetTickCount();
init();
int num=(int)phi(n,len)-sum+len-1;
time=GetTickCount()-time;
GlobalFree(v);
cout<<"pi(n)="<<num<<"("<<time<<"ms)\n";
cout<<"press any key to continue...";
cin.get();
cin.get();
}
void init(void){
int max;
int sqrt_max;
bool* mark;
int i,j;
int len2,len3;
int s;
max=(int)pow(1.0*n,2.0/3.0);
sqrt_max=(int)sqrt(max);
mark=(bool*)GlobalAlloc(GMEM_FIXED|GMEM_ZEROINIT,max);
for(i=2;i<sqrt_max;++i){
if(!mark){
j=i*i;
while(j<max){
mark=true;
j+=i;
}
}
}
for(len=0,i=2;(unsigned)i*i*i<=n;++i){
if(!mark)primes[++len]=i;
}
for(j=max-1,sum=0,s=0,len2=len;(unsigned)i*i<=n;++i){
if(!mark){
++len2;
while((unsigned)i*j>n){
if(!mark)++s;
--j;
}
sum+=s;
}
}
for(len3=len2;i<max;++i){
if(!mark)++len3;
}
sum=(len2-len)*len3-sum;
sum+=len*(len-1)/2-len2*(len2-1)/2;
if(m>len)m=len;
for(i=1;i<=m;++i){
Q*=primes;
phiQ*=primes-1;
}
v=(int*)GlobalAlloc(GMEM_FIXED,Q*sizeof(int));
for(i=0;i<Q;++i)v=i;
for(i=1;i<=m;++i){
for(j=Q-1;j>=0;--j){
v-=v];
}
}
GlobalFree(mark);
}
bool string2int(unsigned int& n,char* s){
unsigned int val=0;
for(int i=0;s;++i){
if(s>'9'||s<'0')return false;
val*=10;
val+=s-'0';
if(val>default_n)return false;
}
n=val;
return true;
}
unsigned int phi(unsigned int x,int a){
if(a==m){
return (x/Q)*phiQ+v;
}
if(x<(unsigned)primes){
return 1;
}
return phi(x,a-1)-phi(x/primes,a-1);
} 另外,哪位能写个可以自动到PrimeDatabase上查询,并能记录结果的脚本?
gxqcn
发表于 2008-11-19 13:34:59
原帖由 medie2005 于 2008-11-19 12:40 发表 http://bbs.emath.ac.cn/images/common/back.gif
关于计算$pi(n)$,目前有时间复杂度为$O(n^((1/2)+varepsilon))$,空间复杂度为$O(n^((1/4)+varepsilon))$的算法。
...
可否提供参考资料?
liangbch
发表于 2008-11-19 14:35:41
回15楼, 张一飞我知道,张一飞是3届(2000,2001,2002)IOI国家集训队的成员,第14届国际信息学奥林匹克竞赛金牌获得者。intfree的这个代码我也看过(可惜我我看不懂), 但我不知道intfree就是张一飞。楼主怎么知道的?
gxqcn
发表于 2008-11-19 14:44:54
Google 了一下:http://topic.csdn.net/t/20020827/09/972666.html
mathe
发表于 2008-11-19 15:32:17
这个不难,在Linux下面直接用脚本调用wget就可以了。然后将返回的结果找关键字就可以了。
mathe
发表于 2008-11-19 15:37:20
如果上面intfree的代码是$O(n^{2/3})$,既然计算到$4*10^9$只需要0.5秒,那么计算到$4*10^12$也就需要100倍,50秒左右,还是可以容忍的