数学研发论坛

 找回密码
 欢迎注册
楼主: medie2005

[讨论] 阶乘和素数个数

[复制链接]
发表于 2008-11-17 19:22:53 | 显示全部楼层


关键是计算到那个程度
需要多少时间?
我想算10^16的就很长
否则不会在70年代算不出
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-17 20:30:12 | 显示全部楼层
原帖由 无心人 于 2008-11-17 15:28 发表
你还是想办法把
$O(n^{2/3})$
的$pi(n)$算法
先做好了吧


这篇论文我曾上传到论坛的;在 Donald E.Knuth 著的 TAOCP 中某道习题中提及。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-11-17 20:53:28 | 显示全部楼层
呵呵,gxq说的那篇文章很好找,google一下就能找到免费的原文。
我找了一下,发现还有时间复杂度为$O(n^(3/5))$的算法来计算$pi(n)$,
不过,好像是没什么人用那个算法来计算$pi(n)$,看样子只是时间复杂度好看点而已。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-11-17 21:03:51 | 显示全部楼层
另外,有一个网站:PrimeDatabase。里面可以查到第1到第$10^12$个素数。
我对网页不熟悉,是不是能写一个脚本,每隔一段时间就查一个数,然后记录返回的素数?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 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秒。

代码如下:
  1. #include   <iostream.h>   
  2. #include   <math.h>   
  3. #include   <windows.h>   

  4. const   unsigned   default_n=4000000000;   
  5. const   int   maxlen=1000;   

  6. int   primes[maxlen];   
  7. int   len;   
  8. int   sum;   
  9. unsigned   int   n;   
  10. int   m=7;   
  11. int   Q=1;   
  12. int   phiQ=1;   
  13. int*   v;   

  14. bool   string2int(unsigned   int&   n,char*   s);   
  15. void   init(void);   
  16. unsigned   int   phi(unsigned   int   x,int   a);   

  17. void   main(void){   
  18.         char   number_string[80];   
  19.         unsigned   int   time;   
  20.    
  21.         cout<<"calc   pi(n)\n";   
  22.         cout<<"n(default   =   "<<default_n<<")   =   ";   
  23.         cin.get(number_string,80);   
  24.         if(!string2int(n,number_string))n=default_n;   
  25.    
  26.         time=GetTickCount();   
  27.         init();   
  28.         int   num=(int)phi(n,len)-sum+len-1;   
  29.         time=GetTickCount()-time;   
  30.         GlobalFree(v);   
  31.         cout<<"pi(n)="<<num<<"("<<time<<"ms)\n";   
  32.         cout<<"press   any   key   to   continue...";   
  33.         cin.get();   
  34.         cin.get();   
  35. }   

  36. void   init(void){   
  37.         int   max;   
  38.         int   sqrt_max;   
  39.         bool*   mark;   
  40.         int   i,j;   
  41.         int   len2,len3;   
  42.         int   s;   
  43.    
  44.         max=(int)pow(1.0*n,2.0/3.0);   
  45.         sqrt_max=(int)sqrt(max);   
  46.         mark=(bool*)GlobalAlloc(GMEM_FIXED|GMEM_ZEROINIT,max);   
  47.    
  48.         for(i=2;i<sqrt_max;++i){   
  49.                 if(!mark[i]){   
  50.                         j=i*i;   
  51.                         while(j<max){   
  52.                                 mark[j]=true;   
  53.                                 j+=i;   
  54.                         }   
  55.                 }   
  56.         }   
  57.    
  58.         for(len=0,i=2;(unsigned)i*i*i<=n;++i){   
  59.                 if(!mark[i])primes[++len]=i;   
  60.         }   
  61.    
  62.         for(j=max-1,sum=0,s=0,len2=len;(unsigned)i*i<=n;++i){   
  63.                 if(!mark[i]){   
  64.                         ++len2;   
  65.                         while((unsigned)i*j>n){   
  66.                                 if(!mark[j])++s;   
  67.                                 --j;   
  68.                         }   
  69.                         sum+=s;   
  70.                 }   
  71.         }   
  72.    
  73.         for(len3=len2;i<max;++i){   
  74.                 if(!mark[i])++len3;   
  75.         }   
  76.    
  77.         sum=(len2-len)*len3-sum;   
  78.         sum+=len*(len-1)/2-len2*(len2-1)/2;   
  79.    
  80.         if(m>len)m=len;   
  81.         for(i=1;i<=m;++i){   
  82.                 Q*=primes[i];   
  83.                 phiQ*=primes[i]-1;   
  84.         }   
  85.    
  86.         v=(int*)GlobalAlloc(GMEM_FIXED,Q*sizeof(int));   
  87.         for(i=0;i<Q;++i)v[i]=i;   
  88.         for(i=1;i<=m;++i){   
  89.                 for(j=Q-1;j>=0;--j){   
  90.                         v[j]-=v[j/primes[i]];   
  91.                 }   
  92.         }   
  93.    
  94.         GlobalFree(mark);   
  95. }   

  96. bool   string2int(unsigned   int&   n,char*   s){   
  97.         unsigned   int   val=0;   
  98.         for(int   i=0;s[i];++i){   
  99.                 if(s[i]>'9'||s[i]<'0')return   false;   
  100.                 val*=10;   
  101.                 val+=s[i]-'0';   
  102.                 if(val>default_n)return   false;   
  103.         }   
  104.         n=val;   
  105.         return   true;   
  106. }   

  107. unsigned   int   phi(unsigned   int   x,int   a){   
  108.         if(a==m){   
  109.                 return   (x/Q)*phiQ+v[x%Q];   
  110.         }   
  111.         if(x<(unsigned)primes[a]){   
  112.                 return   1;   
  113.         }   
  114.         return   phi(x,a-1)-phi(x/primes[a],a-1);   
  115.   }   
复制代码
另外,哪位能写个可以自动到PrimeDatabase上查询,并能记录结果的脚本?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-19 13:34:59 | 显示全部楼层
原帖由 medie2005 于 2008-11-19 12:40 发表
关于计算$pi(n)$,目前有时间复杂度为$O(n^((1/2)+varepsilon))$,空间复杂度为$O(n^((1/4)+varepsilon))$的算法。

...


可否提供参考资料?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-19 14:35:41 | 显示全部楼层
回15楼, 张一飞我知道,张一飞是3届(2000,2001,2002)IOI国家集训队的成员,第14届国际信息学奥林匹克竞赛金牌获得者。intfree的这个代码我也看过(可惜我我看不懂), 但我不知道intfree就是张一飞。楼主怎么知道的?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-19 14:44:54 | 显示全部楼层
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-19 15:32:17 | 显示全部楼层
这个不难,在Linux下面直接用脚本调用wget就可以了。然后将返回的结果找关键字就可以了。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-19 15:37:20 | 显示全部楼层
如果上面intfree的代码是$O(n^{2/3})$,既然计算到$4*10^9$只需要0.5秒,那么计算到$4*10^12$也就需要100倍,50秒左右,还是可以容忍的
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2019-7-17 04:16 , Processed in 0.055177 second(s), 15 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表