Eratosthenes 筛法的 Free-Pascal 实现

Question

我的老师给我布置了这样的作业：使用给定的数n，求最大质数p且p<=n 和 n<=10^9。我尝试使用以下函数执行此操作：

Const amax=1000000000
Var i,j,n:longint;
    a:array [1..amax] of boolean;
Function lp(n:longint):longint;
 Var max:longint;
 Begin
  For i:=1 to n do a[i]:=true;
  For i:=2 to round(sqrt(n)) do
   If (a[i]=true) then
    For j:=1 to n div i do
     If (i*i+(j-1)*i<=n) then
      a[i*i+(j-1)*i]:=false;
  max:=0;
  i:=n;
  While max=0 do
   Begin
    If a[i]=true then max:=i;
    i:=i-1;
   End;
  lp:=max;
 End;

此代码对于诸如 100 万之类的数字可以完美运行，但是当我尝试 n=10^9 时，程序需要很长时间才能打印输出。所以这是我的问题：有什么方法可以改进我的代码以降低延迟吗？或者可能是不同的代码？

Answer 1

这里最重要的方面是不大于 n 的最大素数必须非常接近 n。快速浏览 The Gaps Between Primes (at The Prime Pages - 总是值得一看与素数有关的一切）表明对于 32 位数字，素数之间的差距不能大于 335。这意味着最大的素数不大于 n 必须在 [n - 335, n] 范围内。换句话说，最多需要检查 336 个候选人——例如通过试分——这肯定比筛选十亿个数字要快得多。

对于这类任务，试分是一个合理的选择，因为扫描的范围太小了。在我对 Prime sieve implementation (using trial division) in C++ 的回答中，我分析了几种加快速度的方法。

Eratosthenes筛法也是一个不错的选择，只是需要修改为只筛选感兴趣的范围，而不是从1到n的所有数字。这称为 'windowed sieve'，因为它只筛选 window。由于 window 很可能 而不是 包含直到 n 的平方根的所有素数（即所有可能是范围内复合材料的潜在最小素数因子的素数被扫描）最好通过一个单独的、简单的埃拉托色尼筛法筛选因子素数。

首先，我展示了一个简单的普通（非windowed）筛法，作为比较windowed 代码的基准。我使用 C# 是为了比使用 Pascal 更清楚地展示算法。

List<uint> small_primes_up_to (uint n)
{
    if (n == uint.MaxValue)
        throw new ArgumentOutOfRangeException("n", "n must be less than UINT32_MAX");

    var eliminated = new bool[n + 1];  // +1 because indexed by numbers

    eliminated[0] = true;
    eliminated[1] = true;

    for (uint i = 2, sqrt_n = (uint)Math.Sqrt(n); i <= sqrt_n; ++i)
        if (!eliminated[i])
            for (uint j = i * i; j <= n; j += i)
                eliminated[j] = true;

    return remaining_unmarked_numbers(eliminated, 0);
}

该函数的名称中有 'small'，因为它不太适合筛选大范围；我使用类似的代码（带有一些花里胡哨的东西）仅用于筛选更高级的筛子所需的小因子素数。

提取筛选素数的代码同样简单：

List<uint> remaining_unmarked_numbers (bool[] eliminated, uint sieve_base)
{
    var result = new List<uint>();

    for (uint i = 0, e = (uint)eliminated.Length; i < e; ++i)
        if (!eliminated[i])
            result.Add(sieve_base + i);

    return result;
}

现在，windowed 版本。一个区别是潜在的最小因子素数需要单独筛选（通过刚刚显示的函数），如前所述。另一个区别是给定素数的交叉序列的起点可能位于要筛选的范围之外。如果起点位于 window 的开始之前，那么需要一点模数魔法才能找到落在 window 中的第一个 'hop'。从此一切如常。

List<uint> primes_between (uint m, uint n)
{
    m = Math.Max(m, 2);

    if (m > n)  
        return new List<uint>();  // empty range -> no primes

    // index overflow in the inner loop unless `(sieve_bits - 1) + stride <= UINT32_MAX`
    if (n - m > uint.MaxValue - 65521)  // highest prime not greater than sqrt(UINT32_MAX)
        throw new ArgumentOutOfRangeException("n", "(n - m) must be <= UINT32_MAX - 65521");

    uint sieve_bits = n - m + 1;
    var eliminated = new bool[sieve_bits];

    foreach (uint prime in small_primes_up_to((uint)Math.Sqrt(n)))
    {
        uint start = prime * prime, stride = prime;

        if (start >= m)
            start -= m;
        else
            start = (stride - 1) - (m - start - 1) % stride;

        for (uint j = start; j < sieve_bits; j += stride)
            eliminated[j] = true;
    }

    return remaining_unmarked_numbers(eliminated, m);
}

模计算中的两个“-1”项可能看起来很奇怪，但它们将逻辑向下偏置 1 以消除需要映射到 0 的不方便情况 stride - foo % stride == stride。

这样，不超过n的最大素数可以这样计算：

uint greatest_prime_not_exceeding (uint n)
{
    return primes_between(n - Math.Min(n, 335), n).Last();
}

这总共花费了不到一毫秒的时间，包括因子素数的筛选等等，即使代码不包含任何优化。在我对 prime number summing still slow after using sieve 的回答中可以找到对适用优化的一个很好的概述；使用此处显示的技术，可以在大约半秒内筛选高达 10^9 的整个范围。

Eratosthenes 筛法的 Free-Pascal 实现

Free-Pascal Implementation of the Sieve of Eratosthenes

freepascal

sieve-of-eratosthenes