为了计算第n个素数,我知道两个主要变体。
直截了当的方式
也就是说,当你找到它们时,从2开始计算所有素数,直到你达到所需的nth。
这可以通过不同程度的复杂程度和效率来完成,并且有两种概念上不同的方法可以做到这一点。首先是
测试序列中所有数字的素数
这将通过驱动程序功能实现,例如
public static int nthPrime(int n) {
int candidate, count;
for(candidate = 2, count = 0; count < n; ++candidate) {
if (isPrime(candidate)) {
++count;
}
}
// The candidate has been incremented once after the count reached n
return candidate-1;
}
决定效率的有趣部分是功能。isPrime
鉴于素数的定义是一个大于1的数字,并且我们在学校¹中学到的,因此,素数检查的明显方法是
审判庭
将定义直接转换为代码是
private static boolean isPrime(int n) {
for(int i = 2; i < n; ++i) {
if (n % i == 0) {
// We are naive, but not stupid, if
// the number has a divisor other
// than 1 or itself, we return immediately.
return false;
}
}
return true;
}
但是,如果您尝试一下,您很快就会发现,它的简单性伴随着缓慢。通过这个素数检验,你可以在几毫秒内找到第1000个素数,7919(在我的计算机上大约20个),但是找到第10000个素数,104729,需要几秒钟(~2.4s),第10000个素数,1299709,几分钟(约5),第百万个素数,15485863,大约需要八个半小时,第一千万个素数,179424673, 周,依此类推。运行时复杂度比二次 Θ(n² * log n) 差。
因此,我们希望稍微加快一些素数测试的速度。许多人采取的一个步骤是意识到(除了它本身)的除数最多可以是 。如果我们使用这个事实,让试验除法循环只运行到 而不是 ,算法的运行时间如何变化?对于复合数,循环下限不会更改任何内容。对于素数,试验区的数量减半,因此总体而言,运行时间应减少一个略小于2的系数。如果你尝试一下,你会发现运行时间几乎完全减半,所以几乎所有的时间都花在验证素数的素数上,尽管复合物比素数多得多。n
n
n/2
n/2
n-1
现在,如果我们想找到第一亿个素数,这并没有多大帮助,所以我们必须做得更好。为了进一步减少循环限制,让我们看看实际需要哪些数字的上限。如果 是 的除数,则是整数,换句话说,可被 2 整除。但是,循环不会超过 2,因此它永远不会(除了 )到达 。很高兴,那么的下一个最大可能的除数是什么?当然,为什么。但只能是整数的除数,换句话说,if 可以被 3 整除。然后循环将在 3 处退出(或更早,在 2 处),并且永远不会到达(除了 )。下一个可能的除数...n/2
n/2
n
n/2
n
n = 4
n/2
n
n/3
n/3
n
n
n/3
n = 9
等一下!我们有 和 .n 的除数成对出现。2 <-> n/2
3 <-> n/3
如果我们考虑 的对应除数对,则 ,即 ,或其中一个,比如说,小于另一个。但是然后和.每对相应的除数包含(至少)一个不超过 的除数。(d, n/d)
n
d = n/d
d = √n
d
d*d < d*(n/d) = n
d < √n
n
√n
如果 为复合,则其最小的非平凡除数不超过 。n
√n
因此,我们可以将循环限制减小到 ,从而降低算法的运行时复杂性。它现在应该是Θ(n1.5 * √(log n)),但从经验上讲,它似乎扩展得更好一点 - 但是,没有足够的数据从经验结果中得出可靠的结论。√n
这将在大约16秒内找到第100万个素数,在不到9分钟的时间内找到第1000万个素数,它将在大约四个半小时内找到第1亿个素数。这仍然很慢,但与十年左右的天真审判部门相去甚远。
由于存在素数的平方和两个接近素数的乘积,例如 323 = 17*19,因此我们无法将试验除法循环的极限降低到下面。因此,在坚持试验分工的同时,我们现在必须寻找其他方法来改进算法。√n
一个容易看到的事情是,除了2之外,没有素数是偶数,所以我们只需要在处理完2之后检查奇数。不过,这并没有太大的区别,因为偶数是最便宜的,可以找到复合数 - 而且大部分时间仍然花在验证素数的素数上。但是,如果我们将偶数视为候选除数,我们会看到,如果可以被偶数整除,则本身必须是偶数,因此(2除外)在尝试除以任何大于2的偶数之前,它将被识别为复合数。因此,算法中出现的所有大于 2 的偶数除法都必须留下非零余数。因此,我们可以省略这些除法,并仅检查2的可整除性以及从3到的奇数。这将确定一个数为素数或复合数所需的除法数减半(不完全正确),从而缩短运行时间。这是一个良好的开端,但我们能做得更好吗?n
n
√n
另一个大的数字家族是3的倍数。我们执行的每三个除法都是3的倍数,但如果可以被其中一个除以,它也可以被3整除,因此不能除以9,15,21,...我们在算法中执行的将留下剩余的0。那么,我们怎样才能跳过这些分歧呢?好吧,可以被2和3整除的数字恰恰是形式的数字。从5开始(因为我们只对大于1的数字感兴趣),它们是5,7,11,13,17,19,...,从一个到下一个的步长在2和4之间交替,这很容易,所以我们可以使用n
6*k ± 1
private static boolean isPrime(int n) {
if (n % 2 == 0) return n == 2;
if (n % 3 == 0) return n == 3;
int step = 4, m = (int)Math.sqrt(n) + 1;
for(int i = 5; i < m; step = 6-step, i += step) {
if (n % i == 0) {
return false;
}
}
return true;
}
这给我们带来了另一个(接近)1.5倍的加速,所以我们需要大约一个半小时才能达到第一亿个素数。
如果我们继续这条路线,下一步是消除5的倍数。共素数为 2、3 和 5 的数字是形式的数字
30*k + 1, 30*k + 7, 30*k + 11, 30*k + 13, 30*k + 17, 30*k + 19, 30*k + 23, 30*k + 29
因此,我们只需要除以每三十个数字中的八个(加上三个最小的素数)。从一个步骤到下一个,从 7 开始,循环到 4、2、4、2、4、6、2、6。这仍然很容易实现,并且又将速度提高了1.25倍(对于更复杂的代码,减去一点)。更进一步,7的倍数将被消除,每210个数字中有48个被除以,然后是11(480/2310),13(5760/30030)依此类推。每个被消除倍数的素数都会产生(几乎)的加速,因此回报会降低,而成本(代码复杂性,步骤的查找表的空间)会随着每个素数的增加而增加。p
p/(p-1)
一般来说,在消除六到七个素数(甚至更少)的倍数之后,人们很快就会停下来。然而,在这里,我们可以一直坚持到最后,当所有素数的倍数都被消除并且只有素数被保留为候选除数时。由于我们按顺序查找所有素数,因此在需要将其作为候选除数之前找到每个素数,然后可以存储以供将来使用。这将算法复杂性降低到 - 如果我没有计算错误 - O(n1.5 / √(log n))。以存储素数的空间使用为代价。
使用试验除法,这是尽可能好的,你必须尝试除以所有素数或第一除数以确定 的素数。这在这里大约半小时内找到了第一个亿个素数。√n
n
n
那么怎么样
快速素数测试
素数具有其他数论性质,而不是没有非平凡除数,而复合数通常没有。如果这些属性可以快速检查,则可以构成概率或确定性素数检验的基础。这种性质的原型与皮埃尔·德·费马的名字有关,他在17世纪初发现
如果 是素数,则是所有 的 (p-a) 的除数。p
p
a
这个 - 费马所谓的“小定理” - 在等价公式中
设为素数,不可被 整除。然后除以p-1 - 1。p
a
p
p
大多数广泛传播的快速素数检验(例如Miller-Rabin)的基础以及其变体或类似物出现在更多(例如Lucas-Selfridge)中。
因此,如果我们想知道一个不太小的奇数是否是素数(偶数和小数被试验除法有效地处理),我们可以选择任何不是,例如2的倍数的数字(>1),并检查是否除以n-1 - 1。由于n-1 变得很大,因此最有效的方法是检查是否,即通过模幂。如果这种一致性不成立,我们知道这是复合的。然而,如果它成立,我们就不能断定它是素数,例如,但它是复合的。这样的复合数称为基数的费马伪素数。n
a
n
n
a^(n-1) ≡ 1 (mod n)
n
n
2^340 ≡ 1 (mod 341)
341 = 11 * 31
n
a^(n-1) ≡ 1 (mod n)
a
但这种情况很少见。给定任何基数,尽管有无限数量的费马伪素数到基数,但它们比实际的素数稀有得多。例如,只有78个基数-2的费马伪素数和76个基数-3的费马伪素数低于100000,但有9592个素数。因此,如果一个人选择一个任意的奇数和一个任意的基并找到,那么很有可能实际上是素数。a > 1
a
n > 1
a > 1
a^(n-1) ≡ 1 (mod n)
n
但是,我们处于一个稍微不同的情况,我们被给予并且只能选择 。那么,对于一个奇数复合,对于多少,可以成立?不幸的是,有合数 - 卡迈克尔数 - 使得同余对于每个共素数都成立。这意味着,要用费马检验将卡迈克尔数识别为复合数,我们必须选择一个基数,该基数是素数之一的倍数 - 可能没有很多这样的倍数。n
a
n
a
1 < a < n-1
a^(n-1) ≡ 1 (mod n)
a
n
n
但是我们可以加强费马测试,以便更可靠地检测复合材料。如果 是奇数素数,请写 .然后,如果,p
p-1 = 2*m
0 < a < p
a^(p-1) - 1 = (a^m + 1) * (a^m - 1)
并正好除以两个因子中的一个(两个因子相差 2,因此它们的最大公约数是 1 或 2)。如果是均匀的,我们可以以同样的方式分裂。继续,如果带有奇数,则写入p
m
a^m - 1
p-1 = 2^s * k
k
a^(p-1) - 1 = (a^(2^(s-1)*k) + 1) * (a^(2^(s-2)*k) + 1) * ... * (a^k + 1) * (a^k - 1)
然后精确地除以其中一个因素。这就产生了强烈的费马检验,p
设为奇数。用奇数书写。给定任何带有 ,如果n > 2
n-1 = 2^s * k
k
a
1 < a < n-1
-
a^k ≡ 1 (mod n)
或
-
a^((2^j)*k) ≡ -1 (mod n)
对于任何与j
0 <= j < s
那么是一个强(费马)可能素数的碱。复合强碱(费马)可能素数称为基的强(费马)伪素数。强费马伪素数甚至比普通费马伪素数更罕见,低于10000000,有78498个素数,245个基数-2个费马伪素数,只有46个基数-2基数强费马伪素数。更重要的是,对于任何奇数复合,大多数基都是强费马伪素数。n
a
a
a
n
(n-9)/4
1 < a < n-1
n
因此,如果 是一个奇数复合,则通过强费马检验且随机选择的基数在 1 和(排除边界)之间的概率小于 。n
n
k
n-1
1/4^k
强费马检验需要O(log n)步,每个步骤涉及一个或两个带有O(log n)位的数字的乘法,因此复杂度是O((log n)^3)和朴素乘法[对于巨大,更复杂的乘法算法可能是值得的]。n
米勒-拉宾检验是随机选择碱基的k倍强费马检验。这是一个概率检验,但对于足够小的边界,已知基数的短组合,从而给出确定性结果。
强费马检验是确定性APLC试验的一部分。
建议在此类测试之前先对前几个小素数进行试验除法,因为除法相对便宜,并且会清除大多数复合材料。
对于找到第三个素数的问题,在测试所有数的素数是可行的范围内,有已知的碱基组合使多个强费马检验正确,因此这将给出一个更快的 - O(n*(log n)4) - 算法。n
因为,基数 2、7 和 61 足以验证素数。使用这个,在大约六分钟内找到第一亿个素数。n < 2^32
消除素除数的复合物,埃拉托斯特尼筛
与其按顺序研究这些数字并从头开始检查每个数字是否是素数,不如将整组相关数字视为一个整体,并一次性消除给定素数的倍数。这被称为埃拉托斯特尼的筛子:
找到不超过的素数N
- 列出从 2 到以下的所有数字
N
- 对于从 2 到 : 如果尚未划掉,则为素数;划掉作为复合材料的所有倍数
k
N
k
k
素数是列表中未划掉的数字。
该算法与试验除法有根本的不同,尽管两者都直接使用素数的可除性表征,与费马检验和使用素数其他性质的类似检验相反。
在试验除法中,每个数与所有素数配对,其数不超过 的较小素数和最小素数除数。由于大多数复合材料具有非常小的素除数,因此平均而言,检测复合材料在这里是便宜的。但是测试素数是昂贵的,因为下面有相对多的素数。尽管复合数比素数多得多,但测试素数的成本非常高,以至于它完全主导了整体运行时间,并使试验除法成为一种相对较慢的算法。所有小于以下数字的试验除法需要 O(N1.5 / (log N)²) 步长。n
√n
n
√n
N
在筛子中,每个复合物都与其所有素除数配对,但仅与这些除数配对。因此,素数是廉价的数字,它们只被看过一次,而复合材料更昂贵,它们被划掉了多次。有人可能会认为,由于筛子包含比“便宜”数字更多“昂贵”的数字,因此它总体上是一个糟糕的算法。然而,一个合数没有许多不同的素数除数 - 的不同素数除数的数目以 为界,但通常它要小得多,这些数的不同素数除数的平均值是 - 所以即使是筛子中的“昂贵”数字平均也不会比试验除法的“便宜”数字更昂贵(或几乎不更贵)。n
n
log n
<= n
log log n
筛分到,对于每个素数,都有要划掉的倍数,所以划线的总数是。这产生了比试验除法或具有更快素数检验的顺序测试更快的算法来查找素数。N
p
Θ(N/p)
Θ(∑ (N/p)) = Θ(N * log (log N))
N
然而,筛子有一个缺点,它使用内存。(但是使用分段筛,可以将其减少到不增加时间复杂性。O(N)
O(√N)
为了找到第th个素数,而不是直到的素数,还有一个问题是事先不知道筛子应该达到多远。n
N
后者可以使用素数定理求解。PNT说
π(x) ~ x/log x (equivalently: lim π(x)*log x/x = 1),
其中 是不超过素数的数目(这里和下面,必须是自然对数,对于算法复杂性,为对数选择哪个基数并不重要)。由此可以推断,哪里是第素数,并且从更深入的分析中知道有很好的上限,特别是π(x)
x
log
p(n) ~ n*log n
p(n)
n
p(n)
n*(log n + log (log n) - 1) < p(n) < n*(log n + log (log n)), for n >= 6.
因此,可以将其用作筛分限制,它不会超过目标。
通过使用分段筛可以克服空间要求。然后,可以记录下面的素数以消耗内存,并使用长度增加的段(当筛子接近N时,使用O(√N))。O(N)
√N
O(√N / log N)
如上所述,对算法有一些简单的改进:
- 开始划掉仅处 的倍数,而不是
p
p²
2*p
- 从筛子中消除偶数
- 从筛子中消除更多小素数的倍数
这些都没有降低算法的复杂性,但它们都显着减少了常数因子(与试验除法一样,消除倍数的产生较小的速度,而增加代码复杂性的幅度大于较小的速度)。p
p
p
使用前两个改进可产生效果
// Entry k in the array represents the number 2*k+3, so we have to do
// a bit of arithmetic to get the indices right.
public static int nthPrime(int n) {
if (n < 2) return 2;
if (n == 2) return 3;
int limit, root, count = 1;
limit = (int)(n*(Math.log(n) + Math.log(Math.log(n)))) + 3;
root = (int)Math.sqrt(limit) + 1;
limit = (limit-1)/2;
root = root/2 - 1;
boolean[] sieve = new boolean[limit];
for(int i = 0; i < root; ++i) {
if (!sieve[i]) {
++count;
for(int j = 2*i*(i+3)+3, p = 2*i+3; j < limit; j += p) {
sieve[j] = true;
}
}
}
int p;
for(p = root; count < n; ++p) {
if (!sieve[p]) {
++count;
}
}
return 2*p+1;
}
它在大约18秒内找到第一亿个素数,2038074743。通过存储打包的标志(每个标志一位,而不是作为 s),可以将此时间减少到大约 15 秒(此处为 YMMV),因为减少的内存使用量提供了更好的缓存局部性。boolean
包装标志,也消除了3的倍数,并使用位微调以更快,更快地计数,
// Count number of set bits in an int
public static int popCount(int n) {
n -= (n >>> 1) & 0x55555555;
n = ((n >>> 2) & 0x33333333) + (n & 0x33333333);
n = ((n >> 4) & 0x0F0F0F0F) + (n & 0x0F0F0F0F);
return (n * 0x01010101) >> 24;
}
// Speed up counting by counting the primes per
// array slot and not individually. This yields
// another factor of about 1.24 or so.
public static int nthPrime(int n) {
if (n < 2) return 2;
if (n == 2) return 3;
if (n == 3) return 5;
int limit, root, count = 2;
limit = (int)(n*(Math.log(n) + Math.log(Math.log(n)))) + 3;
root = (int)Math.sqrt(limit);
switch(limit%6) {
case 0:
limit = 2*(limit/6) - 1;
break;
case 5:
limit = 2*(limit/6) + 1;
break;
default:
limit = 2*(limit/6);
}
switch(root%6) {
case 0:
root = 2*(root/6) - 1;
break;
case 5:
root = 2*(root/6) + 1;
break;
default:
root = 2*(root/6);
}
int dim = (limit+31) >> 5;
int[] sieve = new int[dim];
for(int i = 0; i < root; ++i) {
if ((sieve[i >> 5] & (1 << (i&31))) == 0) {
int start, s1, s2;
if ((i & 1) == 1) {
start = i*(3*i+8)+4;
s1 = 4*i+5;
s2 = 2*i+3;
} else {
start = i*(3*i+10)+7;
s1 = 2*i+3;
s2 = 4*i+7;
}
for(int j = start; j < limit; j += s2) {
sieve[j >> 5] |= 1 << (j&31);
j += s1;
if (j >= limit) break;
sieve[j >> 5] |= 1 << (j&31);
}
}
}
int i;
for(i = 0; count < n; ++i) {
count += popCount(~sieve[i]);
}
--i;
int mask = ~sieve[i];
int p;
for(p = 31; count >= n; --p) {
count -= (mask >> p) & 1;
}
return 3*(p+(i<<5))+7+(p&1);
}
在大约9秒内找到第一亿个素数,这不是令人难以忍受的长。
还有其他类型的素筛,特别有趣的是阿特金筛,它利用了这样一个事实,即(有理)素数的某些同余类是Q的某些二次扩展的代数整数环中的复合。这里不是扩展数学理论的地方,足以说阿特金筛子的算法复杂性低于埃拉托斯特尼筛,因此对于大极限更可取(对于小限制,没有过度优化的阿特金筛具有更高的开销,因此可能比相对优化的埃拉托斯特尼筛慢)。D. J. Bernstein的素数库(用C语言编写)针对低于232的数字进行了很好的优化,并在大约1.1秒内找到了第一亿个素数(这里)。
快速的方式
如果我们只想找到第 个素数,那么也找到所有较小的素数就没有内在价值。如果我们能跳过其中的大部分,我们可以节省大量的时间和工作。给定第 th 个素数的良好近似值,如果我们有一种快速的方法来计算不超过的素数,那么我们可以筛选高于或低于的一个小范围,以识别 和 之间缺失或多余的几个素数。n
a(n)
n
p(n)
π(a(n))
a(n)
a(n)
a(n)
p(n)
我们已经看到了一个容易计算的相当好的近似值,我们可以采取p(n)
a(n) = n*(log n + log (log n))
例如。
一个好的计算方法是Meissel-Lehmer方法,它大致在一段时间内进行计算(确切的复杂性取决于实现,Lagarias,Miller,Odlyzko,Deléglise和Rivat的改进使人们以O(x2/3 / log² x)时间计算)。π(x)
π(x)
O(x^0.7)
π(x)
从简单的近似开始,我们计算。根据素数定理,附近的素数密度大约是 ,所以我们期望接近,我们会筛出一个略大于 的范围。为了在筛分范围内获得更大的信心,可以将范围增加2倍,例如,这几乎肯定会足够大。如果范围看起来太大,可以用更好的近似值代替 、 计算 和 进行迭代。通常,将比 小得多。如果 是近似值,则将是 更好的近似值。只有在非常接近(并且不是非常接近0)的极不可能的情况下,找到一个足够好的近似值,即最终的筛分阶段可以在与计算相当的时间内完成,这成为一个问题。a(n)
e(n) = π(a(n)) - n
a(n)
1/log a(n)
p(n)
b(n) = a(n) - log a(n)*e(n)
log a(n)*e(n)
p(n)
b(n)
a(n)
π(b(n))
f(n) = π((b(n)) - n
|f(n)|
|e(n)|
f(n)
-e(n)
c(n) = (a(n) + b(n)) / 2
p(n)
f(n)
e(n)
p(n)
π(a(n))
通常,在对初始近似进行一两次改进后,要筛分的范围足够小,以使筛分阶段的复杂度为O(n^0.75)或更高。
这种方法在大约40毫秒内找到第一亿个素数,在八秒内找到第10个12个素数,29996224275833。
tl;dr:找到第三个素数可以有效地完成,但是你越想要它,涉及的数学就越多。n
我在这里准备了大多数讨论的算法的Java代码,以防有人想玩它们。
¹ 对过度感兴趣的灵魂的旁注:现代数学中使用的素数的定义是不同的,适用于更一般的情况。如果我们调整学校的定义以包括负数 - 所以一个数字是素数,如果它既不是1也不是-1,并且只能被1,-1,本身及其负数整除 - 这定义了(对于整数)现在称为Z的不可约元素,但是,对于整数,素数和不可约元素的定义是一致的。