筛选质数高达 10^12

Question

我正在做一些事情，需要我生成最多 10^12 的所有素数。

因为我以前从来没有需要这么多素数，所以我通常只是在这个网页上实现算法here

当然，这里的问题是 10^12 大于整数的最大值，因此我无法创建该大小的数组。

我不熟悉有效生成这么多素数的方法，想知道是否有人可以阐明这种情况。

Answer 1

您需要使用分段筛。

分段筛的基本思想是选择小于 n 的平方根的筛分素数，选择一个相当大但适合内存的片段大小，然后依次筛选每个片段，从最小的片段开始.在第一段，计算段内每个筛素的最小倍数，然后以正常方式将筛素的倍数标记为合数；当所有的筛选素数都用完后，该段中剩余的未标记数是素数。然后，对于下一个片段，对于每个筛分素数，您已经知道当前片段中的第一个倍数（它是结束前片段中那个素数的筛分的倍数），所以您在每个筛分素数上进行筛分，依此类推直到你完成。

考虑将 100 到 200 分成 20 段的示例； 5个筛选素数分别是3、5、7、11和13。在从100到120的第一段中，位数组有10个槽，槽0对应101，槽k对应100 + 2k + 1，槽9对应119。段中3的最小倍数是105，对应slot 2；槽 2+3=5 和 5+3=8 也是 3 的倍数。槽 2 处 5 的最小倍数是 105，槽 2+5=7 也是 5 的倍数。7 的最小倍数是 105在槽 2 处，槽 2+7=9 也是 7 的倍数。以此类推。

函数 primes 接受参数 lo、hi 和 delta； lo 和 hi 必须是偶数，lo

function primes(lo, hi, delta)
    sieve := makeArray(0..delta-1)
    ps := tail(primes(sqrt(hi)))
    m := length(ps)
    qs := makeArray(0..m-1)
    for i from 0 to m-1
        qs[i] := (-1/2 * (lo + ps[i] + 1)) % ps[i]
    while lo < hi
        for i from 0 to delta-1
            sieve[i] := True
        for i from 0 to m-1
            for j from qs[i] to delta step ps[i]
                sieve[j] := False
            qs[i] := (qs[i] - delta) % ps[i]
        for i from 0 to delta-1
            t := lo + 2*i + 1
            if sieve[i] and t < hi
                output t
        lo := lo + 2*delta

对于上面给出的示例，这称为 primes(100, 200, 10)。在上面给出的示例中，qs 最初是 [2,2,2,10,8]，对应于最小倍数 105、105、105、121 和 117，并且对于第二段重置为 [1,2,6, 0,11]，对应最小倍数123、125、133、121和143。

delta 的值很关键；您应该使 delta 尽可能大，只要它适合高速缓存内存，以提高速度。将您的语言库用于位数组，这样您只需为每个筛子位置取一个位。如果你需要一个简单的埃拉托色尼筛来计算筛分质数，这是我最喜欢的：

function primes(n)
    sieve := makeArray(2..n, True)
    for p from 2 to n step 1
        if sieve(p)
            output p
            for i from p * p to n step p
                sieve[i] := False

这些函数都是伪代码；您必须使用适当的整数数据类型转换为 Java。在伪代码显示输出的地方，您可以打印素数，或将素数收集到一个数组中，无论您想对它们做什么。

我在博客上对素数做了很多工作，包括文章 Programming with Prime Numbers，其中包括最后一页上的分段筛。

Answer 2

真正的解决方案是找到另一种方法来解决潜在问题，而无需生成完整的素数集。根据Prime Number Theorem，素数之间的平均差距为 ln(1e12)，约为 27.6。这给出了超过 39e9 个小于 1e12 的素数的估计。

您可能不需要所有这些。考虑研究生成可能素数和/或素数测试的方法。当然，如果不了解您要解决的潜在问题，就不可能确切地知道该怎么做。

Answer 3

这是我用于计算素数段的 Java 代码：

/**
 * Computes the primes in a range using the sieve of Eratosthenes.
 * The size of the range must not exceed Integer.MAX_VALUE.
 *
 * @param start  The start index of the prime sieve.
 * @param limit  Primes will be sieved up to but not including this limit.
 *
 * @return  A bit set representing the integer range from start to limit.
 *     Each bit in this set is set to true if and only if
 *     the corresponding integer is prime.
 */
public static BitSet computePrimes(long start, long limit)
{
    if (limit - start > Integer.MAX_VALUE)
    {
        throw new IllegalArgumentException();
    }

    final long sqrtLimit = sqrtCeil(limit);
    final BitSet primes = computePrimes((int) sqrtLimit);

    final BitSet segment = new BitSet();
    if (0 - start >= 0)
    {
        segment.set((int) (0 - start), false);
    }
    if (1 - start >= 0)
    {
        segment.set((int) (1 - start), false);
    }
    segment.set((int) (Math.max(0, 2 - start)), (int) (limit - start), true);
    for (int d = 2; d < sqrtLimit; d++)
    {
        if (primes.get(d))
        {
            final int remainder = (int) (start % d);
            final long mStart = start - remainder + (remainder == 0 ? 0 : d);
            for (long m = Math.max(mStart, d * d); m < limit; m += d)
            {
                segment.clear((int) (m - start));
            }
        }
    }
    return segment;
}

它需要一个标准筛子来计算筛分片段的素数（它会为每个片段重新计算它，你应该改变它）：

/**
 * Computes the primes using the sieve of Eratosthenes.
 *
 * @param limit  Primes will be sieved up to but not including this limit.
 *
 * @return  A bit set where exactly the elements with prime index
 *     are set to true.
 */
public static BitSet computePrimes(int limit)
{
    final BitSet primes = new BitSet();
    primes.set(0, false);
    primes.set(1, false);
    primes.set(2, limit, true);
    for (int d = 2; d < sqrtCeil(limit); d++)
    {
        if (primes.get(d))
        {
            for (int m = d * d; m < limit; m += d)
            {
                primes.clear(m);
            }
        }
    }
    return primes;
}

请注意，车轮分解可以将速度提高三倍。另见this answer，基本筛子是一样的。