擂台：筛法求质数Sieve of Eratosthenes

筛法求质数Sieve of Eratosthenes描述：

http://primes.utm.edu/glossary/page.php?sort=SieveOfEratosthenes

相关帖子:
擂台：求N!的最后九个非零尾数
http://expert.csdn.net/Expert/topic/2986/2986867.xml?temp=.511532

＃＃谁愿意帮我做个测试，最好是P4的CPU！＃＃ 
http://expert.csdn.net/Expert/topic/3027/3027142.xml?temp=.850918

程序下载（最新版，暂不提供源码，仍在不断改进中）：
http://www.77studio.net/files/PrimeLabPro.zip

参考源代码：
http://www.77studio.net/files/PrimeLab_Src.zip
http://www.77studio.net/files/PrimeLab_Src_New.zip （比前一个新）

另一种更快的方法：
Sieve of Atkins
介绍：http://primes.utm.edu/links/programs/sieves/binary_quadratic/index.html
程序：http://www.77studio.net/files/SieveAtkins.zip

特别感谢xstring(麻雀), liangbch(宝宝), languagec(各有所求)（是从你的帖子开始研究质数的^_^）等，只为交流，适当的时候会和大家分享代码。

注：
以上未嵌入汇编，未使用ICC(Intel C++ Compiler)优化，VS 6.0 + SP5编译

43 个解决方案

#1

提供一个测试数据，以供参考，我的CPU是赛扬II 900MHz，PC133 384M内存
L2 Cache为128KB全速Cache（和CPU主频一样），更详细的CPU以及Cache细节请下载
CPU-Z软件: http://www.skycn.com/soft/3943.html

测试数据：
PrimeLab Pro - Windows(R) CPU Identifier
-----------------------------------------

    Microsoft Windows 2000 Advanced Server Service Pack 4(Build 2195)
    -----------------------------------------------------------------
    TotalPhys: (383MB)392692KB, AvailPhys: 90700KB, FreePercent: 24%

    Pentium III (0.18 um) With 256 KB On-Die L2 Cache  - [CPU #0]
    -------------------------------------------------------------
        | Clock Frequency: 896MHz
        |
        + On-Chip Hardware Features
            | L1 Cache: 32KB
            | L2 Cache: 128KB

-----------------------------------------

Input a number [2-4294967295]: 4294967295

PrimeLab_Src.zip   约17.1秒
SieveAtkins.zip    约15.5秒

#2

上面打错了，更正：

PrimeLabPro.zip    约17.1秒
SieveAtkins.zip    约15.5秒

以上2个程序都是只统计质数个数，不保存质数，只要内存够大的话，保持标记位也不会慢很多（多几秒而已），鉴于多数程序都只统计个数，由此比较效率，所以不讨论保存标记位，如有需要，修改一下就是了

所提供的源码保存了质数的标记位，所以计算4294967295以内的质数需要分配约256M的内存，所以内存小于256M的机器可能会很慢，请自行修改一下代码，改为只统计质数个数

#3

我曾下载过一篇文章：computing-pi-x-the.pdf

COMPUTING p(x) x): THE MEISSEL-LEHMER METHOD

             J. C. Lagarias
          AT&T Bell Laboratories
          Murray Hill, New Jersey

             V. S. Miller
          I BM Watson Research Center
          Yorktown Heights, New York

             A. M. Odlyzko
          AT&T Bell Laboratories
          Murray Hill, New Jersey

里面介绍了一些方法，并列举了一些运行时间：

pi(  10^12)   2min    37,607,912,018
pi(2x10^12)   3min
pi(3x10^12)   4min
pi(4x10^12)   5min
pi(5x10^12)   6min
pi(6x10^12)   7min
pi(7x10^12)   7min
pi(8x10^12)   8min
pi(9x10^12)   8min

pi(  10^13)   9min    346,065,536,839

pi(  10^14)  39min    3,204,941,750,802

pi(  10^15) 165min    29,844,570,422,669

pi(  10^16) 718min    279,238,341,033,925

#4

顶。。目前只学了 the Sieve of Eratosthenes  找质数的方法。。。其他的学习中……

#5

机器不达标。关注之。

#6

to NowCan(((((( ★ ))))))：

PrimeLabPro.zip
SieveAtkins.zip

对内存没什么要求的，CPU也不会要求很高，但是最好有128KB的L2 Cache(这是唯一的要求，否则将会使程序慢2倍，即所花时间是正常的3倍）

#7

我按  shines(郭靖的郭，英雄的雄，光辉的辉)  给的网址上的程序改了一个 算10的10次方,能够输出的 要将近6分钟. shines 把你的代码公布一下吧,让大家研究研究.

#8

你算10的10次方，大约需要596M的内存，如果你内存太小，将会对速度有很大的影响；
而且要计算2^32次方以上的必须使用__int64，也将会变慢不少，最好的做法是把
2^32和2^32以后的筛选分为2个部分来计算，这样前面的就会快很多。

另外，我的网站今天被人利用动网的漏洞删了文件（本来以为那个门派的论坛不会有人来理的，就懒得升级了，没想到他还真无聊：），我明天再上传。过阵子，我将整理一下写个文档，一块把代码公布出来，目前我在计算超高精度大数的阶乘，已经有些接近擂主的速度了，参看：

http://community.csdn.net/Expert/topic/3049/3049738.xml?temp=.0685541

#9

To gxqcn(GxQ):
  以前你说过我的算法像是Sieve of Atkins，其实不是，说句实话，我还没看懂Sieve of Atkins呢，呵呵，不过粗看是有点像，它使用了60K+1, 60K+3....60K+59，我用的是30K+1, 30K+7,...30K+29，可以看到是不太一样的，我没有30K+3, 30K+5，因为这些数已经包含了3或5的倍数，是已知的合数，关于Sieve of Atkins，以后有空再研究，它的确很快。我原本的意思是先把传统的Sieve of Eratosthenes提高到接近或超过Sieve of Atkins的效率，不过现在看来是不太可能了，只能有些逼近。不过那个Sieve of Atkins代码过于依赖Cache，如果Cache Line不够（被其他程序所占）的时候，程序的效率会大打折扣的。

  其实我的方法还是传统的筛选法，只是对于3,5,7,11等筛选（2早就用数据存储的方法筛选掉了）做了特殊处理，即整理出来循环节，什么叫循环节，就是对2...N筛选的时候，你会发现筛选3的倍数的时候是会有循环节出现的，那么同时筛掉3,5的循环节，同时筛掉3,5,7的循环节，3,5,7,11。。。，都是可以计算出来的，因为筛选的数越小，则筛选的次数就越多，所以这个简单的优化是可以提高不少速度的，由于L2 Cache的限制，现在证实，3,5,7,11,13的循环节效率还没有3,5,7,11的循环节好。

  3,5,7,11的循环节需要3x5x7x11x4=4620字节存储，3,5,7,11,13的循环节则需要3x5x7x11x13x4=60060字节，至于筛选质数的质数倍30K+1,30K+7,30K+11,...30K+29, (K>=1)，原理也很简单，当初是受宝宝的启发，其实想想也很容易理解的，

liangbch(宝宝)：

to  NowCan 和 xstring(麻雀)
  在筛质数的倍数时，筛去质数的所有倍数自然是最慢的，筛去所有质数的质数倍是最快的，但这知识理论上的，象xstring(麻雀)说的， 实际上无法操作，典型的筛法是筛去质数的部分倍数，这样就有好多种筛法， xstring(麻雀)给出两种算法，我这里再总结和补充以下。
   方法1，筛去质数的 k 倍 (k>1)   (所有倍数)
   方法2，筛去质数的 2*k+1 倍(k>=1)（奇数)，xstring(麻雀)所说的方法1，
   方法3, 筛去质数的 6k+1,6K+5（ 非2和3 的倍数 ,k>=1 ) ，xstring(麻雀)所说的方法2

   方法4, 筛去质数的 30k+1,30K+7,30K+11,30K+13,30K+17,30K+19,30K+23,30K+29 （ 非2、3和5 的倍数  )
   方法5：筛去 非2，3，5，7的倍数

如果第一种算法耗费的时间为1，在方法2，3，4耗费的时间为1/2=0.5,1/3=0.333,8/30=0.267,48/210=0.229， 速度分别方法一的 2倍，3倍，3.75倍，4.375倍.

很容易就可以扩展到210K+1, 210K+11,210K+13,....210K+199，即筛去 非2，3，5，7的倍数，但由于编程的角度，选用30K+1,30K+7.....是比较合理的，某些时候可以考虑210K+1...，越往后提升的可能性就越低，同时的程序的复杂度就越高

#10

嘿嘿，楼上的兄弟都忽略了

你们的程序要求内存太多了

实际上，即使计算到10 ^ 16，也不用超过10 M内存

嘿嘿，分配 2 * 3 * 5 * 7 * 11 * 13 * 17 * 19 = 9699690的块就可以了，而且可以预先消掉2、3、5、7、11、13、17、19的所有倍数（这个过程不消耗时间，增加一个块复制时间，但是节约掉的时间更多），然后分块计算！！！

至于压缩表示，则节约更多内存

#11

1/2+1/3+1/5+1/7+1/11+1/13+1/17+1/19 = 1.45

谁计算一下到10000的素数倒数和

#12

To yaos(夜夜秋雨孤灯下·无心人):

你试试新的，已经不保存所有的标记位了，只记录分块的标记味，质数总个数在分块中统计出来并累加得到的，所以对内存没有什么要求，

http://www.77studio.net/files/PrimeLabPro.zip

另外，“分配 2 * 3 * 5 * 7 * 11 * 13 * 17 * 19 = 9699690的块就可以了”，9699690约等于9.6Mb，请注意一点，要让你的程序更快，必须让块的大小适合于L2 Cache，所以我为什么选择了3x5x7x11x13x4=60060的循环块，而我的CPU的L2 Cache是128KB，所以选择3x5x7x11x4 × 16 = 4620 × 16 = 73920最佳（即16个3x5x7x11的分块，没超过128Kb，但又比60060大，分块越大，分块计算的次数越少，速度就越快），如果是P4的CPU，L2 Cache大小是512KB，倒是可以3x5x7x11x13x4 × N = 60060 × N，只要60060 × N不倒于L2 Cache的大小就可以了，由于还有其他数据会保存在Cache中，所以经验得出60060 × N最好不要大于L2 Cache的2/3大小。

由于程序的角度，分块的大小应该是这么计算的，你可以自己实践一下就知道为什么是这样的了，比如你上面说的分块大小就应该是 3 * 5 * 7 * 11 * 13 * 17 * 19  *  4 = 19399380，为什么后面都会有个4呢，因为在现在的32Bit的CPU上，处理32Bit的数据是最快的，即一个DWORD，如果你是按照BYTE大小来保存循环块，那循环节是要小一点，但不利用数据的Copy和计算，因为CPU多数情况都要在数据4字节对齐的时候才最快（尤其是memcpy等）。

还有一点，现在的CPU的Cache都是32字节（256Bit）对齐，和4-way的，这也导致了，第一次开始筛选的倍数不能太大，否则就难以一开始就把分块的数据都load到Cache中了，如果Cache的命中率很低的话，程序的效率会是利用了Cache的1/3左右，即大约多花2倍的时间

#13

还有一个问题

如果分块太小，节约的时间不能抵消块复制的时间

这个没分析，请你看看

#14

够快的,佩服

我的算法被csdn贪污了,自己没保留,唉

我是计算最大素数间隔的时候做的程序,后来发现人家已经有10^16以下的结果了,就放弃了

#15

另外

知道n^(1/3)的所有素数，是可以计算出来小于n的素数的个数的

速度很快，计算10^20之类的素数个数就源于此算法

可惜我没参透

有了解的兄弟请指教

#16

实践证明，分块较小的时候（最适合L2 Cache的时候，当然也可以适合到L1 Cache（比L2快不少），可是由于L1 Cache实在太小，不适合），复制块的时间基本上可以忽略，大概是由于都在Cache中的原因，不必去RAM中去载入，增加的时间幅度不太记得了，大概在5-10%左右，相对于分块大于L2 Cache的时候增加2倍的时候是很值得的（也就是前者的300%左右）。

程序最好动手写起来，我的程序很乱，也没有什么技巧，就是为了给大家动手的机会
当然我新的程序比它考虑周全一点，其实给你们的源码里，只要把所有的
if(0 == _getisprime2(p, cursor))
注释掉就可以提高近一倍的速度，就是在写入标记位的时候先检查这个数是否是质数
本来这个技巧在所求的质数不是很大的时候是有很大的效率提高的，可是当N变大以后，
而且使用了6K+1, 6K+5等筛法以后，它却变位很费时间的一个操作，因为6K+1, 30K+1等方法
已经尽量是过滤的倍数是质数倍了，所以这个检查是没有太大意义的，去掉之后试试看，是不是快了很多，当然可以改进的地方还很多。。。

#17

考虑直接算的算法吧

另外，公开源代码吧

#18

对不起，已经看到题目的代码了

有点繁琐，嘿嘿

暂时看不懂，新程序的代码有么？

那个临时块可以程序生成的，写到程序中，似乎不好吧
不过节省0.01秒左右的时间 ：）

#19

to yaos(夜夜秋雨孤灯下·无心人)：

你说的对，知道n^(1/3)的所有素数，是可以计算出来小于n的素数的个数的

这个程序我有，也就是数学家所研究的PI(X)函数，PI是圆周率的π，但是这个只适合计算质数的个数，所以不拿来这里讨论，其他它也必须使用n^(1/3)的筛法

#20

另外，现代的CPU是读优化的，就是说，读不浪费时钟，即使Cache不命中
写不行，有Cache读入的延迟

#21

临时块可以程序生成的，我另外写了一个程序来生成，放在程序里只是当初就放了，嘿嘿：）
其实花不了多少程序的时间，基本可以忽略不计：）

等我写好了文档，我再公布源码吧，我最近在忙那个阶乘的程序，昨天才告一段落

#22

嘿嘿，我按照我的思路，你的分析做做，看结果如何
Cache按128kB计算

#23

好的，这样最好了

#24

#include <iostream.h>
#define  CacheSize 4 * 16 * 3 * 5 * 7 * 11

unsigned char PCache[CacheSize];
unsigned char PBuffer[CacheSize];
unsigned long n;
unsigned long sum = 0;
unsigned long P[]; //素数表

//分块筛法计算n --> n + p - 1的素数
void BlockSieve(unsigned long n, unsigned long m)
{
unsigned long i, j, k, l, r, t;

for (i = 0; i < CacheSize; i ++)
  PBuffer[i] = PCache[i];

k = sqrt(n);
l = sqrt(N + p - 1);

for (i = 5; P[i] <= k; i ++)
  {
  //标记P[i]的倍数, 在PBuffer中
  r = (n / P[i]) * P[i];
  if (r < n)
    r = r + P[i];
  if (r % 2 == 0)
    r = r + P[i];
  for (j = r - n; j <= m - 1; j + = 2 * P[i])
    PBuffer[j] = 0;  
  }

for (; P[i] <= l; i ++)
  {
  //标记P[i]的倍数
  for (j = P[i] * P[i] - n; j < m - 1; j += 2 * P[i])
    PBuffer[j] = 0;
  }

//输出
for (i = 0; i < CacheSize; i ++)
  if (PBuffer[i])
    {
    //输出
    sum ++;
    }  
}

void initPrime(unsigned long n)
{
unsigned long ps = n / ln(n) * 1.2; //大概素数个数
if (ps < 10) ps = 10;
P = new unsigned long[ps];
//同本程序的筛法
...

}

void initCache(void);
{
unsigned long i;
for (i = 0; i < CacheSize; i ++)
  { 
  if ((i % 2) == 0)
     PCache[i] = 0;
  else
  if ((i % 3) == 0)
     PCache[i] = 0;
  else
  if ((i % 5) == 0)
     PCache[i] = 0;
  else
  if ((i % 7) == 0)
     PCache[i] = 0;
  else 
  if ((i % 11) == 0)
     PCache[i] = 0;
  else
     PCache[i] = 1; 
  }
}

void main(void)
{
unsigned i, j, m, k;
initCache();

cin >> n;

//初始化素数表
initPrime(sqrt(n));

if (n < 1)
  cout << "0";
else
  if (n < 13)
    for (i = 0; i < n; i ++)
      {
      if (PCache[i])
        { 
        //输出
        sum ++;  
        } 
      }
else
  {
  m = n / (cacheSize);
  k = n % (CacheSize);
  if (m > 0)
    for (i = 0; i < m; i ++)
      BlockSieve(i * CacheSize, CacheSize);
  BlockSieve(m * CacheSize, k);
  }

//输出sum
cout >> endl >> sum;
}

先写提纲,后完善 :)

#25

嘿嘿，刚放上去，就发现了一个错误 :)
n < 13的计算是错的，似乎也不能根据n < 13划分情况

#26

我只能说很遗憾，你是用一个char来表示一个标记位（就是一个BYTE），而为了压缩数据，可以使用bit来做标记位，这样不是省了7/8的存储空间了吗？再者，偶数中只有2是质数，我们完全可以一开始就把偶数给去掉，标记位中第一个bit存储的是1的标记位，依次是3,5,7,9,......，偶数根本没有必要处理，这样就又省了1/2，为什么要这么省呢？虽然这样读写标记位处理是复杂了，而节省的空间却可以使用较大的循环块，因而总体的速度还是快了很多

#27

如果你上面的代码只是大体的提供的话，还可以，不过你的Cache应该是这样定义（按照你的分块的逻辑）才对：

#define  CacheSize (2 * 3 * 5 * 7 * 11 * 32)

还有，如果你要把第一个分块的位置放在从0开始的话，要注意第一个分块是特殊的，你注意到了吗？因为你把2,3,5,7,11都筛掉了（可是他们是质数），所以一般第一个分块不从0开始，或者做特别处理

#28

更正，大体的提纲

#29

多谢指教

确实是Bit压缩效率高，不过要牺牲点运算时间啊，除非用386 Bit运算指令

另外，昨天没想仔细，第一个块要做特殊处理，最好是直接数 ：）
素数表的生成也要特殊处理，打算利用筛的结果，需要一个初始素数表，大概
小于sqrt(CacheSize)就可以了，但是总的素数表至少要包含小于cacheSize的
所有素数！！

程序中，已经忽略了偶数了 ：）
j + = 2 * P[i]
至于你再次强调的Bit标志问题，还是用汇编做比较合适，性能提高很大的哦

#30

嘿嘿，我还不打算上汇编，调试比较麻烦，暂时用Byte吧 ：）

#31

WoW Beta2的私服我玩过了，打怪和不能打怪（看风景的），Beta3的私服也出来了，可惜硬盘最近有点紧张，暂时没有下载

其实你没有忽略偶数，你只是跳过，真正忽略应该都不保存它的标记位，这样你的CacheSize要大一倍，还有for中i++, i+=2的效率是不一样的

质数表是这样的，要求N内的质数，至少需要求出sqrt(N)内的所有质数，这个好办，对于2^32（即DWORD），只要求出2^16(65536)内的质数就可以了，所以我的循环块(3*5*7*4)×K=420×K（K为系数，使420*K刚好大于sqrt(N)，1<=K<=10），轻易就可以覆盖了这个65536内的质数，因为我的一个Byte覆盖了16个数的标记位（为什么看前面的bit压缩），所以420×K的循环块其实相当于420*16*K=6720*K，K=10时，＝67200，还比65536略多，所以除了前面一个420×K是特殊处理以外（这是一个小的循环节，重新看了下我的代码里是INIT_BLOCK_SIZE ((3*5*7 * 4) * 1)），真正开始循环我是从420*16×K+1开始的，循环块的大小是CACHE_BLOCK_SIZE (60060) = (3*5*7*11*4)*13，其实后来我改为73920稍快一些

至于Bit标志的问题，不用汇编也可以的，因为编译器基本上可以给你一个比较好（只能说还行，否则编译器也不是吃素的）的效率，当然我也还没有用汇编试试，多数的情况下没有很高的技巧或者明显汇编可以优化的地方，是不会有很大的效率优化

#32

386汇编就可以根据偏移来存取位信息，才1个时钟，严重优化阿

#33

今天查了一下，的确 BT 和 BTs (s为任一字母) 指令可以对位检测和赋值，应该有得优化

#34

你有时间做做吧

我最近没时间调试程序了

另外，我们能不能从16进制考虑问题 ：）

#35

16进制？什么意思？不明白？

你先看下我的程序，我的循环块都是16进制的啊

最近2天生活有点乱，下了beta3，可服务器就上去了一下就当了，就看见了个开场动画

#36

也就是有4种情况
位数组，小cache
位数组，大cache
字节数组，小cache
字节数组，大cache

#37

我使用字节数组，分段筛，每段为 65520 ，筛4*10^9以内的所有素数，筛过程用时22秒，统计用时6秒，共计28.5秒。楼主的程序SieveAtkins 为6.32秒。

#38

好像intfree的程序求10G内素数在 C2 535上求个数只要<10秒,求出素数只要50秒吧,咔咔

#39

如果秋初素数只要50秒的话，那求个数绝对不用10秒。

#40

这是 intfree 的程序 :

1 : 可以求出具体素数的 ;

{$APPTYPE CONSOLE}
{$R-,S-,Q-,O-}

program prime;

uses sysutils, windows;

const
  maxn = int64(10) shl 30;
  maxlen = 10000;
  m = 6;
  dots = 80;

type
  Pls = ^Tls;
  Tls = array[byte] of byte;

var
  mask: pointer;
  primes, mprimes, ans: array[1 .. maxlen] of integer;
  bits: array[byte] of integer;
  len, r: integer;
  start: TDateTime;
  
procedure solve(const a, b, c: integer; var x, y: integer);
begin
  if a = 0 then begin
    x := 0; y := c div b;
  end else begin
    solve(b mod a, a, c, y, x);
    x := x - b div a * y;
  end;
end;

procedure init;
var
  p, k: integer;
begin
  k := 2; len := 0;
  repeat
    len := len + 1; primes[len] := k;
    repeat
      k := k + 1; p := 1;
      while (sqr(primes[p]) < k) and (k mod primes[p] <> 0) do p := p + 1;
    until k mod primes[p] <> 0;
  until 1.0 * k * k > maxn;
  bits[0] := 0; for k := 1 to $FF do bits[k] := bits[k shr 1] + k and 1;
end;

procedure makeprimes;
var
  i, k, n, size, kmax, min, last, p, x, dot: integer;
begin
  n := 1; for i := 1 to m do n := n * primes[i];
  for i := m + 1 to len do begin
    p := primes[i]; solve(p, - n, 1, x, mprimes[i]);
    mprimes[i] := (mprimes[i] mod p + p) mod p;
  end;
  r := len; size := (maxn div n + 7) shr 3; last := 0; dot := 0;
  fillchar(ans, sizeof(ans), 0); getmem(mask, size);
  for i := 1 to n - 1 do begin
    k := 1; while (k < m) and (i mod primes[k] <> 0) do k := k + 1;
    if i mod primes[k] = 0 then continue;
    p := i - last; last := i; kmax := (maxn - i) div n;
    for k := m + 1 to len do ans[k] := (ans[k] + mprimes[k] * p) mod primes[k];
    fillchar(mask^, size, 0);
    for k := m + 1 to len do begin
      p := primes[k];
      if ans[k] = 0 then min := p - 1 else min := ans[k] - 1;
      asm
        mov EAX, mask
        mov EBX, kmax
        mov ECX, min
        mov EDX, p
@001:
        BTS [EAX], ECX
        add ECX, EDX
        cmp ECX, EBX
        jl @001
      end;
    end;
    for k := 0 to size - 1 do kmax := kmax - bits[Pls(mask)^[k]];
    r := r + kmax;
    if i / n > dot / dots then begin
      dot := dot + 1; write('.');
    end;
  end;
  freemem(mask, size);
  writeln;
end;

begin
  SetPriorityClass(GetCurrentProcess, REALTIME_PRIORITY_CLASS);
  start := now;  
  init;
  makeprimes;
  messageboxA(0,
              PChar(Format('sum = %d, Time = %.2f',
                           [r, (now - start) * SecsPerDay])),
              'prime', 0); 
end.

#41

2: 求素数个数的,咔咔,好像以前的大虾如 starfish , intfree, mathe 都不来了啊

{$R-,S-,Q-,O+}
{$APPTYPE CONSOLE}

program prime;

uses
  sysutils;

const
  maxlen3 = 1000;
  dn = int64(10) shl 30;
  dm = 7;

var
  primes: array[1 .. maxlen3] of integer;
  v: array of integer;
  mark: array of boolean;
  m, ans, Q, phiQ: integer;
  n, len2, len3, sum: int64;
  start: TDateTime;

procedure init;
var
  l, k, s, max, sqmax: integer;
  p: int64;
begin

  max := round(exp(2 / 3 * ln(1.0 * n))) + 1;
  sqmax := round(sqrt(max));
  setlength(mark, max);
  for s := 2 to sqmax do
    if not mark[s] then
    begin
      k := s * s;
      while k < max do
      begin
        mark[k] := true;
        k := k + s;
      end;
    end;

  l := 0; p := 2;
  while p * p * p <= n do
  begin
    if not mark[p] then
    begin
      l := l + 1;
      primes[l] := p;
    end;
    p := p + 1;
  end;
  len3 := l;

  k := max - 1; sum := 0; s := 0;
  while p * p <= n do
  begin
    if not mark[p] then
    begin
      l := l + 1;
      while p * k > n do
      begin
        if not mark[k] then
          s := s + 1;
        k := k - 1;
      end;
      sum := sum + s;
    end;
    p := p + 1;
  end;
  len2 := l;

  while p < max do
  begin
    if not mark[p] then
      l := l + 1;
    p := p + 1;
  end;

  sum := (len2 - len3) * l - sum;
  sum := len3 * (len3 - 1) div 2 - len2 * (len2 - 1) div 2 + sum;

  Q := 1;
  phiQ := 1;
  if len3 < dm then
    m := len3
  else
    m := dm;
  for s := 1 to m do
  begin
    Q := Q * primes[s];
    phiQ := phiQ * (primes[s] - 1);
  end;

  setlength(v, Q);
  for s := 0 to Q - 1 do
    v[s] := s;
  for s := 1 to m do
    for k := Q - 1 downto 0 do
      v[k] := v[k] - v[k div primes[s]];

end;

function phi(x: int64; a: integer): int64;
begin

  if a = m then
  begin
    result := (x div Q) * phiQ + v[x mod Q];
    exit;
  end;

  if x < primes[a] then
  begin
    result := 1;
    exit;
  end;

  result := phi(x, a - 1) - phi(x div primes[a], a - 1);

end;

var
  s: string;
begin
  writeln('calc pi(n)');
  write(format('n (default = %d) = ', [dn]));
  try
    readln(s);
    n := strtoint64(s);
  except
    n := dn;
  end;

  start := now;
  init;
  ans := phi(n, len3) - sum + len3 - 1;

  writeln(format('pi(%d) = %d (%.2fs)',
                 [n, ans, (now - start) * SecsPerDay]));
  writeln('Press enter');
  readln;
end.

#42

我改进了我的程序，但仍然使用字节数组，分段筛，每段为 262140，筛4*10^9以内的所有素数，筛过程用时15.38秒，统计用时4.26秒，共计20.10秒(有一部分计算时间不在前2者之内). 
    若改为每段120K,总时间为22.60秒.

    楼主的程序SieveAtkins 为在我的电脑上运行速度为6.32秒。

  字节数组最大的问题在于读写内存的数据量较大，特别是统计素数的个数部分，很难提高速度。

#43

按照CSDN的精神，此贴最近将结