我正在做一些需要我生成高达 10^12 的所有素数的事情。
由于我以前从未需要过这么多素数,所以我通常只是在此网页上实现算法
当然,这里的问题是 10^12 大于整数的最大值,因此我无法创建该大小的数组。
我不熟悉用于有效生成这么多素数的方法,并且想知道是否有人可以阐明这种情况。
您将需要使用分段筛。
分段筛的基本思想是选择小于 n 的平方根的筛分素数,选择一个相当大但仍适合内存的片段大小,然后依次筛选每个片段,从最小的片段开始。在第一段,计算段内每个筛素的最小倍数,然后以正常方式将筛素的倍数标记为合数;当所有的筛选素数都用完后,该段中剩余的未标记数是素数。然后,对于下一个片段,对于每个筛分素数,您已经知道当前片段中的第一个倍数(它是结束前片段中那个素数的筛分的倍数),所以您在每个筛分素数上进行筛分,依此类推直到你完成。
考虑以 20 为一组从 100 到 200 筛分的例子;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 < hi,并且 lo 必须大于 hi 的平方根。段大小是 delta 的两倍。长度为 m 的数组 ps 包含小于 hi 的平方根的筛选素数,由于偶数被忽略,因此删除了 2,这是通过正常的 Eratosthenes 筛法计算的。数组 qs 包含在相应筛分素数的当前段中最小倍数的筛位数组中的偏移量。在每个段之后,lo 前进两倍 delta,因此筛位数组的索引 i 对应的数字是 lo + 2 i + 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。
delta 的值很关键;您应该使 delta 尽可能大,只要它适合高速缓存内存,以提高速度。将您的语言库用于位数组,这样您只需为每个筛子位置取一个位。如果你需要一个简单的 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。在伪代码显示输出的地方,您可以打印素数,或将素数收集到数组中,无论您想对它们做什么。
我在我的博客上做了很多关于素数的工作,包括用素数编程的文章,其中包括最后一页的分段筛。
真正的解决方案是找到另一种方法来解决潜在问题,而无需生成完整的素数集。根据素数定理,素数之间的平均间隔为 ln(1e12),约为 27.6。这给出了超过 39e9 个小于 1e12 的素数的估计。
您可能不需要所有这些。考虑研究生成可能素数和/或素数测试的方法。当然,如果不知道您要解决的潜在问题,就不可能确切地知道要做什么。
这是我用于计算素数段的 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,基本筛子是一样的。