为了计算第 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 且只能被 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;
}
但是,如果您尝试一下,您很快就会发现,它的简单性伴随着缓慢。
通过素数测试,您可以在几毫秒内找到第 1000th 素数 7919(在我的计算机上大约为 20),但找到第 10000th 素数 104729,需要几秒钟(~2.4s),第 100000 个素数,1299709,几分钟(大约 5),第 100 万个素数,15485863,大约需要 8 个半小时,第 1000 万个素数, 179424673、周等。运行时复杂度比二次方差 - Θ(n² * log n)。
所以我们想稍微加快素数测试的速度。许多人采取的一个步骤是意识到n(除了n 本身)的除数最多可以是n/2。
如果我们利用这个事实,让试除循环只运行到n/2而不是n-1,那么算法的运行时间会如何变化?
对于合数,循环下限不会改变任何内容。对于素数,试除的次数减半,所以总体来说,运行时间应该比 2 小一点。如果你试一试,你会发现运行时间几乎完全减半,所以 几乎所有的时间都花在验证素数的素数上,尽管复合数比素数多得多。
现在,如果我们想找到第 1 亿个素数,这并没有多大帮助,所以我们必须做得更好。尝试进一步减少循环限制,让我们看看实际需要n/2 的上限的数字。如果n/2 是n 的除数,那么n/2 是一个整数,换句话说,n 可以被 2 整除。但是循环不会超过 2,所以它永远不会(@ 除外987654337@) 到达n/2。好极了,那么n 的下一个最大可能除数是什么?
为什么,n/3 当然。但是n/3只能是n的除数,如果它是一个整数,换句话说,如果n可以被3整除。那么循环将在3(或之前,在2)处退出并且永远不会到达@ 987654344@(n = 9 除外)。下一个可能的最大除数 ...
等一下!我们有2 <-> n/2 和3 <-> n/3。 n 的除数成对出现。
如果我们考虑n 的相应除数对(d, n/d),d = n/d,即d = √n,或其中一个,例如d,小于另一个。但是然后d*d < d*(n/d) = n 和d < √n。 n 的每一对对应除数包含(至少)一个不超过√n。
如果 n 是复合的,它的最小非平凡除数不超过 √n。
所以我们可以将循环限制减少到√n,从而降低算法的运行时复杂度。它现在应该是 Θ(n1.5 * √(log n)),但从经验上看,它的扩展性似乎更好一些 - 然而,没有足够的数据可以从经验结果中得出可靠的结论。
这将在大约 16 秒内找到第 100 万个质数,在不到 9 分钟内找到第 1000 万个质数,在大约 4 个半小时内找到第 1 亿个质数。这仍然很慢,但与幼稚的审判部门需要十年左右的时间相去甚远。
由于存在素数平方和两个相近素数的乘积,例如 323 = 17*19,我们不能将试除法循环的限制降低到 √n 以下。因此,我们在继续试用除法的同时,现在必须寻找其他方法来改进算法。
很容易看出,除了 2 之外没有素数是偶数,所以我们只需要在处理完 2 后检查奇数。不过,这并没有太大区别,因为偶数是最便宜的找到复合 - 大部分时间仍然花在验证素数的素数上。但是,如果我们将偶数视为候选除数,我们会看到如果n 可以被偶数整除,则n 本身必须是偶数,因此(除了 2)它在被除数之前会被识别为合数尝试任何大于 2 的偶数。因此,算法中发生的所有大于 2 的偶数除法都必须留下非零余数。因此,我们可以省略这些除法,只检查 2 和 3 到 √n 之间的奇数是否可整除。这将(不完全准确)确定一个数为素数或合数所需的除法数减半,因此减少了运行时间。这是一个好的开始,但我们能做得更好吗?
另一个大数字族是 3 的倍数。我们执行的每三次除法都是 3 的倍数,但如果 n 能被其中一个整除,它也能被 3 整除,因此不能被 3 整除我们在算法中执行的 9, 15, 21, ... 将永远留下余数 0。
那么,我们怎样才能跳过这些部门呢?嗯,既不能被 2 也不能被 3 整除的数字正是 6*k ± 1 形式的数字。从 5 开始(因为我们只对大于 1 的数字感兴趣),它们是 5, 7, 11, 13, 17, 19, ...,从一个到下一个的步骤在 2 和 4 之间交替,即很简单,所以我们可以使用
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 倍的加速,所以我们需要大约一个半小时才能到达第 1 亿个素数。
如果我们继续这条路线,下一步是消除 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
所以我们只需要除以每 30 个数字中的 8 个(加上三个最小的素数)。从一个到下一个的步骤,从 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 的素数。在这里大约半小时内找到第 1 亿个素数。
怎么样
快速素性测试
除了没有复合数通常没有的非平凡除数之外,素数还有其他数论属性。如果可以快速检查这些属性,则可以构成概率或确定性素性检验的基础。此类财产的原型与 Pierre de Fermat 的名字有关,他在 17th 世纪初发现
如果p 是素数,则p 是所有a 的(ap-a) 的除数。
这 - 费马所谓的“小定理” - 是,在等价公式中
设p 为质数且a 不能被p 整除。然后p 除以一个p-1 - 1。
大多数广泛使用的快速素性测试的基础(例如 Miller-Rabin)以及更多出现的变体或类似物(例如 Lucas-Selfridge)。
所以如果我们想知道一个不太小的奇数n是否是质数(偶数和小数可以通过试除法有效地处理),我们可以选择任何不是的数a (> 1) n 的倍数,例如 2,并检查 n 是否除以 an-1 - 1。由于 an-1 变得很大,所以这是最通过检查是否有效地完成
a^(n-1) ≡ 1 (mod n),即通过模幂运算。如果这种一致性不成立,我们就知道n 是复合的。但是,如果它成立,我们不能断定n 是质数,例如2^340 ≡ 1 (mod 341),但341 = 11 * 31 是合数。合数n 使得a^(n-1) ≡ 1 (mod n) 被称为基数a 的费马伪素数。
但这种情况很少见。给定任何基数a > 1,尽管基数a 的费马伪素数是无限的,但它们比实际素数要少得多。例如,100000 以下的基数为 2 的费马伪素数只有 78 个,基数为 3 的费马伪素数只有 76 个,但素数是 9592 个。因此,如果选择任意奇数n > 1 和任意基数a > 1 并找到a^(n-1) ≡ 1 (mod n),则很有可能n 实际上是素数。
但是,我们的情况略有不同,给我们n,只能选择a。那么,对于一个奇数复合n,a^(n-1) ≡ 1 (mod n) 可以保持多少个a、1 < a < n-1?
不幸的是,存在合数 - Carmichael 数 - 使得 每个 a 与 n 互质。这意味着,要将卡迈克尔数识别为与费马检验的复合数,我们必须选择一个基数,该基数是 n 的一个主要因数的倍数 - 这样的倍数可能并不多。
但我们可以加强费马检验,以便更可靠地检测复合材料。如果p 是奇数素数,则写成p-1 = 2*m。那么,如果0 < a < p,
a^(p-1) - 1 = (a^m + 1) * (a^m - 1)
和p 正好除两个因子之一(两个因子相差 2,因此它们的最大公约数是 1 或 2)。如果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
那么n 是基数a 的强(费马)可能素数。复合强基a(费马)可能素数称为基a的强(费马)伪素数。强费马赝素数比普通费马赝素数还要少,在1000000以下,有78498个素数,245个base-2费马赝素,只有46个base-2强费马赝素。更重要的是,对于任何奇数复合n,最多有(n-9)/4 基1 < a < n-1,其中n 是一个强费马赝素。
因此,如果n 是一个奇数组合,则n 通过k 强费马测试的概率会小于1/4^k 之间随机选择的碱基。
强费马测试需要 O(log n) 步,每一步都涉及一个或两个数字与 O(log n) 位的乘法,因此复杂度为 O((log n)^3) 与天真的乘法 [for巨大的n,更复杂的乘法算法可能值得]。
米勒-拉宾检验是随机选择碱基的 k 倍强费马检验。这是一种概率测试,但对于足够小的界限,已知的碱基组合短,可以得出确定性结果。
强费马检验是确定性 APRCL 检验的一部分。
建议在此类测试之前先用前几个小素数进行试除法,因为除法相对便宜,并且可以排除大多数复合物。
对于寻找nth素数的问题,在可以测试所有数的素数的范围内,有已知的碱基组合可以使多重强费马测试正确,所以这将提供更快的 - O(n*(log n)4) - 算法。
对于n < 2^32,基数 2、7 和 61 足以验证素数。使用它,大约 6 分钟就可以找到第 1 亿个素数。
通过素数除数消除复合材料,埃拉托色尼筛法
除了按顺序研究数字并从头开始检查每个数字是否为素数,还可以将整组相关数字视为一个整体,并一次性消除给定素数的倍数。这被称为埃拉托色尼筛:
求不超过N的素数
- 列出从 2 到
N 的所有数字
- 对于从 2 到
N 的每个 k:如果 k 尚未被划掉,则它是素数;将k 的所有倍数划掉为复合词
素数是列表中未被划掉的数字。
该算法与试除法有着根本的不同,尽管两者都直接使用素数的可除性表征,与使用素数其他性质的费马测试和类似测试形成对比。
在试除法中,每个数n 都与不超过√n 和n 的最小因数除数中较小者的所有素数配对。由于大多数复合材料的主要因数非常小,因此平均而言,检测复合材料的成本很低。但是测试素数是昂贵的,因为在√n 以下的素数相对较多。尽管复合比素数多得多,但测试素数的成本如此之高,以至于它完全支配了整个运行时间,并使试除法成为一种相对较慢的算法。小于N 的所有数字的试除法需要 O(N1.5 / (log N)²) 步。
在筛子中,每个复合 n 与其所有的主要除数配对,但仅与这些除数配对。因此,素数是便宜的数字,它们只被查看一次,而复合物更昂贵,它们被多次划掉。有人可能会认为,由于筛子包含的“昂贵”数字比“便宜”数字多得多,因此总体而言它是一种糟糕的算法。然而,一个合数没有很多不同的素因数 - n 的不同素因数的数量以log n 为界,但通常它要小得多,数字<= n 的不同素数除数是log log n - 所以即使是筛子中的“昂贵”数字平均也不比试除法的“便宜”数字贵(或几乎不贵)。
筛选到N,对于每个素数p,有Θ(N/p)的倍数要划掉,所以划掉的总数是Θ(∑ (N/p)) = Θ(N * log (log N))。与使用更快的素数测试进行试除法或顺序测试相比,这会产生更多更快的算法来找到高达N 的素数。
但是,筛子有一个缺点,它使用O(N) 内存。 (但使用分段筛,可以减少到O(√N),而不会增加时间复杂度。)
为了找到nth素数,而不是直到N的素数,还有一个问题是事先不知道筛子应该到达多远。
后者可以用素数定理解决。 PNT 说
π(x) ~ x/log x (equivalently: lim π(x)*log x/x = 1),
其中π(x) 是不超过x 的素数数(这里和下面,log 必须是自然对数,对于算法的复杂性,为对数选择哪个底并不重要)。由此得出p(n) ~ n*log n,其中p(n) 是nth 素数,从更深入的分析可知,p(n) 有很好的上限,特别是
n*(log n + log (log n) - 1) < p(n) < n*(log n + log (log n)), for n >= 6.
因此可以将其作为筛分极限,它不会超过目标。
O(N) 空间要求可以通过使用分段筛来克服。然后可以记录√N 下面的素数,用于O(√N / log N) 内存消耗,并使用增加长度的段(当筛子接近 N 时,O(√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 秒内找到第 1 亿个素数 2038074743。通过存储打包的标志,每个标志一位,而不是 booleans,这个时间可以减少到大约 15 秒(这里是 YMMV),因为减少的内存使用提供了更好的缓存局部性。
打包标志,同时消除 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 秒内找到第 1 亿个素数,这不是令人难以忍受的长。
还有其他类型的素数筛,特别有趣的是阿特金筛,它利用了这样一个事实,即(有理)素数的某些同余类是 ℚ 的某些二次扩展的代数整数环中的复合物。这里不是扩展数学理论的地方,只要说阿特金筛的算法复杂度低于埃拉托色尼筛,因此更适合大限制(对于小限制,没有过度优化的阿特金筛具有更高的开销,因此可能比同等优化的 Eratosthenes 筛子慢)。
D. J. Bernstein 的 primegen 库(用 C 编写)针对 232 以下的数字进行了很好的优化,并在大约 1.1 秒内找到了第 1 亿个素数(此处)。
快捷方式
如果我们只想找到nth 素数,那么同时找到所有较小的素数没有内在价值。如果我们可以跳过其中的大部分,我们可以节省大量的时间和工作。给定a(n) 与nth 素数p(n) 的良好近似,如果我们有一种快速的方法来计算素数π(a(n)) 不超过a(n),那么我们可以筛选高于或低于a(n) 的小范围,以识别a(n) 和p(n) 之间的少数缺失或多余素数。
我们已经看到了上面p(n) 的一个很容易计算的相当好的近似值,我们可以采用
a(n) = n*(log n + log (log n))
例如。
计算π(x) 的一个好方法是Meissel-Lehmer method,它在大约O(x^0.7) 时间内计算π(x)(确切的复杂性取决于实现,Lagarias、Miller、Odlyzko、Deléglise 和 Rivat 的改进让在 O(x2/3 / log² x) 时间内计算一次 π(x)。
从简单的近似 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) 在筛选范围内,可以将范围扩大 2 倍,例如,这几乎肯定会足够大。如果范围似乎太大,可以用更好的近似值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)(并且不是非常接近 0),才能找到与 p(n) 足够好的近似值,以便最终筛分阶段可以及时完成,与计算 @987654525 相当@ 成为一个问题。
一般来说,对初始近似值进行一两次改进后,要筛分的范围小到足以使筛分阶段的复杂度达到O(n^0.75)或更好。
此方法在大约 40 毫秒内找到第 1 亿个素数,并在不到 8 秒的时间内找到第 1012 个素数 29996224275833。
tl;dr: 找到nth 素数可以很有效地完成,但是您想要的效率越高,涉及的数学就越多。
我为大多数讨论的算法准备了 Java 代码here,以防有人想玩弄它们。
¹ 顺便说一句,对于过度感兴趣的灵魂:现代数学中使用的素数定义不同,适用于更一般的情况。如果我们调整学校的定义以包括负数 - 所以一个数是素数,如果它既不是 1 也不是 -1 并且只能被 1、-1、自身及其负数整除 - 这定义了(对于整数)现在所谓的 不可约元素,但是对于整数,素数和不可约元素的定义是一致的。