筛选素数直到10^12

Question

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

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

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

我不熟悉人们会用来有效地生成这么多素数的方法，我想知道是否有人可以对这种情况有所了解。

Answer 1

您将需要使用分段筛子。

分段筛选的基本思想是选择小于n的平方根的筛选素数，选择一个合理的较大的段大小，但仍适合内存，然后依次筛选每个段，从最小的开始。在第一个片段中，计算片段中每个筛选素数的最小倍数，然后以正常方式将筛选素数的倍数标记为复合；当所有筛选素数都已使用时，片段中剩余的未标记数字为素数。然后，对于下一段，对于每个筛选质数，您已经知道当前段中的第一个倍数(正是该倍数结束了对上一段中该质数的筛选)，因此您筛选每个筛选质数，依此类推，直到完成。

考虑在20个片段中从100到200筛选的例子；5个筛选素数是3、5、7、11和13。在第一个从100到120的片段中，位数组有10个时隙，其中时隙0对应于101，时隙k对应于100 + 2k + 1，时隙9对应于119。片段中3的最小倍数为105，对应于插槽2；插槽2+3=5和5+3=8也是3的倍数；5的最小倍数在插槽2为105，插槽2+5=7也是5的倍数；7的最小倍数在插槽2为105，插槽2+7=9也是7的倍数。以此类推。

函数质数采用自变量lo，hi和δ；lo和hi必须是偶数，其中lo < hi，并且lo必须大于hi的平方根。段大小是增量的两倍。长度为m的数组ps包含小于hi的平方根的筛选素数，其中删除了2，因为忽略了偶数，由Eratosthenes的正常筛子计算。数组qs包含到相应筛选质数的当前段中的最小倍数的筛位数组中的偏移量。在每个段之后，L0前进两次δ，因此对应于筛位数组的索引i的数字是L0+2i+ 1。

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

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

增量的值是至关重要的；为了提高速度，您应该使增量尽可能大，因为它适合高速缓存。将您的语言库用于位数组，这样您就可以对每个筛子位置只取一个位。如果你需要一个简单的Eratosthenes筛子来计算筛选质数，这是我最喜欢的：

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，并使用适当的整数数据类型。伪代码表示输出，你可以打印素数，或者在数组中收集素数，无论你想用它们做什么。

我在我的博客上用素数做了很多工作，包括文章，它在最后一页包含了一个分段的筛子。

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，基本筛子是相同的。