为了计算第n个质数,我知道两种主要方法。
直接的方法
从2开始,依次找到所有的质数,并计数,直到达到所需的第n个质数。
这可以通过不同程度的复杂性和效率来完成,有两种概念上不同的方式来实现。第一种是
顺序测试所有数字的素数性
这将通过类似以下驱动函数来实现
public static int nthPrime(int n) {
int candidate, count;
for(candidate = 2, count = 0; count < n; ++candidate) {
if (isPrime(candidate)) {
++count;
}
}
return candidate-1;
}
有趣的部分是决定效率的isPrime
函数。
根据我们在学校学到的质数定义,即大于1且只能被1和它本身整除的数字,进行素性检查的明显方法是
试除法
将定义直接翻译成代码:
private static boolean isPrime(int n) {
for(int i = 2; i < n; ++i) {
if (n % i == 0) {
return false;
}
}
return true;
}
但是,如果你尝试使用它,你很快就会发现它的简单性伴随着缓慢的速度。通过这个素性测试,你可以在几毫秒内(大约在我的电脑上是20毫秒)找到第1000个质数7919,但是要找到第10000个质数104729需要几秒钟(约2.4秒),第100000个质数1299709需要几分钟(约5分钟),第一百万个质数15485863需要大约八个半小时,第一千万个质数179424673需要数周等等。运行时间复杂度比二次函数更糟糕 - Θ(n² * log n)。
因此,我们希望加快素性测试的速度。许多人采取的一步是认识到除了n
本身之外的n
的除数最多只能为n/2
。
如果我们使用这个事实,并让试除法循环仅运行到n/2
而不是n-1
,该算法的运行时间将如何更改?
对于合成数,下限循环不会改变任何东西。对于质数,试除的数量减半,因此整体而言,运行时间应该缩短不到2倍。如果您尝试它,您会发现运行时间几乎减半,因此几乎所有时间都用于验证质数的素性,尽管比质数更多。
现在,如果我们想找到第1亿个质数,那么这并没有太大帮助,所以我们需要做得更好。试图进一步减少循环限制,让我们看看对于哪些数字实际上需要
n/2
的上限。如果
n/2
是
n
的除数,则
n/2
是一个整数,换句话说,
n
可被2整除。但是,循环不会超过2,因此它永远不会(除非
n = 4
)达到
n/2
。非常好,那么
n
的下一个最大可能的除数是什么呢?当然是
n/3
。但是,
n/3
只有在它是整数时才能成为
n
的除数,换句话说,如果
n
可被3整除。然后循环将在3处退出(或者在2之前),永远不会到达
n/3
(除非
n = 9
)。下一个最大可能的除数...
Hang on a minute! We have
2 <-> n/2
and
3 <-> n/3
. The divisors of n come in pairs.
If we consider the pair
(d, n/d)
of corresponding divisors of
n
, either
d = n/d
, i.e.
d = √n
, or one of them, say
d
, is smaller than the other. But then
d*d < d*(n/d) = n
and
d < √n
. Each pair of corresponding divisors of
n
contains (at least) one which does not exceed
√n
.
If
n
is composite, its smallest nontrivial divisor does not exceed
√n
.
所以我们可以将循环限制减小到
√n
,这降低了算法的运行时复杂度。现在应该是Θ(n
1.5 * √(log n)),但根据经验,它似乎比较好扩展 - 但是,没有足够的数据从经验结果中得出可靠的结论。
这将在大约16秒内找到第一百万个质数,在不到9分钟内找到第一千万个质数,并且它会在大约4个半小时内找到第一亿个质数。虽然仍然很慢,但与朴素的试除法需要大约10年左右相比,要快得多。
由于存在平方数和两个接近质数的乘积,例如323 = 17 * 19,因此我们不能将试除法循环的限制减小到
√n
以下。因此,在保持试除法的同时,我们必须寻找其他改进算法的方法。
一个容易看到的事情是,除了2以外没有其他质数是偶数,所以我们只需要在处理完2之后检查奇数。但这并没有太大的区别,因为偶数是最容易被判定为合数的 - 并且大部分时间仍然花费在验证质数的素性上。然而,如果我们把偶数看作候选因子,我们会发现,如果n是偶数,那么n本身必须是偶数,因此(除了2以外)在尝试除以任何大于2的偶数之前,它将被认为是合数。因此,算法中发生的所有大于2的偶数除法都必须留下非零余数。因此,我们可以省略这些除法,并仅检查2和从3到√n的奇数是否可整除。这将减少(不完全准确地)确定数字为质数或合数所需的除法次数,从而缩短运行时间。这是一个好的开始,但我们能做得更好吗?
另一个庞大的数字族群是3的倍数。我们每执行三次除法就会涉及到3的倍数,但如果n可被其中之一整除,则它也可被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倍,因此我们需要大约一个半小时来找到第一亿个质数。
如果我们继续这个方法,下一步就是消除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(n
1.5 / √(log n))。代价是用于存储质数的空间使用。
使用试除法,这是最好的方法,您必须尝试并除以所有小于等于根号n的质数或者首先除以n的第一个除数来确定n是否为质数。这样可以在大约半小时内找到第一亿个质数。
那么快速素性测试呢?
质数除了通常的无非平凡因子之外,还有其他与数论相关的性质,而合数通常没有这些性质。如果这些性质可以快速检查,就可以成为概率或确定性素性测试的基础。原型性质与Pierre de Fermat的名字相关,他在17世纪早期发现了以下内容:如果p是一个质数,则对于所有a,p都是(a^p-a)的因子。这被称为Fermat的小定理,等价于以下表述:设p是一个质数,a不是p的倍数,则p整除a^(p-1)-1。这是大多数广泛使用的快速素性测试(例如Miller-Rabin)和更多变体或类似物的基础(例如Lucas-Selfridge)。
如果我们想知道一个不太小的奇数n是否为素数(偶数和小数字可以通过试除法高效处理),我们可以选择任何不是n的倍数的数字a(> 1),例如2,并检查n是否能整除a^(n-1) - 1。由于a^(n-1)变得非常大,因此最有效的方法是通过检查a^(n-1) ≡ 1 (mod n)来进行模幂运算。如果这个同余式不成立,我们就知道n是合数。然而,如果它成立,我们不能得出n是质数的结论,例如2^340 ≡ 1 (mod 341),但是341 = 11 * 31是合数。对于基数a而言,使得a^(n-1) ≡ 1 (mod n)的合数n被称为费马伪素数。
但这种情况很少见。对于任意基数
a > 1
,虽然存在无限多个以费马伪素数为基数
a
,但它们比实际质数要稀少得多。例如,在小于100000的数中,只有78个以2为基数的费马伪素数和76个以3为基数的费马伪素数,但有9592个质数。因此,如果选择任意奇数
n > 1
和任意基数
a > 1
并找到
a^(n-1) ≡ 1 (mod n)
,那么
n
实际上很可能是质数。
然而,我们面临着稍微不同的情况,我们被给定了n
,只能选择a
。所以,对于一个奇合数n
,有多少个a
,使得a^(n-1) ≡ 1 (mod n)
,其中1 < a < n-1
?
不幸的是,存在复合数 - 卡迈克尔数 - 使得对于每一个与n
互质的a
,同余式都成立。这意味着要用费马测试来识别卡迈克尔数为复合数,我们必须选择一个基数,它是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)
对于任意满足0 <= j < s
的j
如果
n
是以
a
为基数的
强费马概率素数,则称复合数的强基数
a
(费马)概率素数为基数
a
的
强费马伪素数。强费马伪素数比普通费马伪素数更加稀有,在小于1000000的范围内,有78498个质数、245个2进制费马伪素数和仅有46个2进制强费马伪素数。更重要的是,对于任何奇合数
n
,存在最多
(n-9)/4
个基数
1 < a < n-1
,使得
n
是一个强费马伪素数。
因此,如果n
是一个奇合数,则n
通过在1到n-1
之间随机选择基数进行k
次强费马测试的概率小于1/4^k
。
强费马测试需要O(log n)步,每一步涉及一个或两个具有O(log n)位数的数字乘法,因此复杂度是O((log n)^3),使用朴素乘法[对于巨大的n,更复杂的乘法算法可能是值得的]。
米勒-拉宾测试是随机选择基数的k倍强费马测试。这是一种概率性测试,但对于足够小的界限,已知短的基数组合可以给出确定性结果。
强费马测试是确定性APRCL测试的一部分。
建议在前几个小质数的试除之前进行这样的测试,因为除法相对便宜,并且可以淘汰大多数复合物。
对于在可行范围内测试所有数字是否为素数的问题,已知基数组合可以使多重强费马测试正确,因此可以得到更快的-O(n*(log n)^4)-算法来找到第n个素数。
对于<2^32的n,基数2、7和61足以验证素性。使用这个方法,第一亿个素数可以在大约六分钟内找到。
使用质因数消除法,埃拉托斯特尼筛法
与其按顺序逐个检查每个数字是否为素数,不如将所有相关数字视为一个整体,一次性消除给定质数的所有倍数。这就是所谓的埃拉托斯特尼筛法:
要找出不大于N
的素数
- 列出从2到
N
的所有数字列表
- 对于从2到
N
的每个k
:如果k
还没有被划掉,那么它是素数;划掉k
的所有倍数作为合数
素数是未被划掉的列表中的数字。
这种算法与试除法根本不同,尽管两者都直接使用素数的可除性特征,但与费马测试和类似测试使用素数的其他属性形成鲜明对比。
在试除法中,每个数字n都与不超过较小的√n和n的最小素数因子的所有质数配对。由于大多数合数具有非常小的素数因子,因此在平均情况下检测合数是便宜的。但是测试质数很昂贵,因为在√n以下有相对较多的质数。尽管合数比质数多得多,但测试质数的成本非常高,完全支配了整个运行时间,并使试除法成为相对较慢的算法。对于小于N的所有数字的试除法需要O(N^1.5 / (log N)²)步骤。
在筛法中,每个合数n都与其所有质因子配对,但仅限于那些。因此,质数是便宜的数字,它们只被查看一次,而合数更昂贵,它们被多次划掉。有人可能认为,由于筛法包含比便宜数字更多的昂贵数字,它总体上将是一个糟糕的算法。然而,合数没有许多不同的质因子 - n的不同质因子数量受log n的限制,但通常要小得多,小于等于n的数字的不同质因子数量的平均值是log log n - 因此,即使在筛法中,昂贵的数字的平均成本也不会比试除法中的便宜数字更高(或几乎相等)。
筛法可以筛选到
N
,对于每个质数
p
,有
Θ(N/p)
个倍数需要被排除,因此总共需要排除的数量是
Θ(∑ (N/p)) = Θ(N * log (log N))
。这比试除法或使用更快速的素性测试进行顺序测试要快得多。
然而,筛法有一个缺点,它需要O(N)
的内存。(但是,使用分段筛法可以将其降低到O(√N)
,而不增加时间复杂度。)
对于查找第n
个素数,与查找素数N
不同,问题在于事先不知道筛法应该筛选到哪里。
后者可以通过使用素数定理来解决。素数定理说
π(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.
因此,可以将其用作筛选限制,它不会远超目标。
使用分段筛法可以克服O(N)的空间要求。然后可以记录小于√N的质数,以O(√N / log N)的内存消耗为代价,并使用递增长度的片段(当筛子接近N时为O(√N))。
上述算法有一些简单的改进:
- 仅从p²开始划去p的倍数,而不是从2*p开始
- 从筛子中删除偶数
- 从筛子中删除更小的其他质数的倍数
这些都不会减少算法复杂度,但它们都会显着减少常数因子(与试除法一样,划去p的倍数对于较大的p产生较小的加速效果,同时使代码复杂性比较小的p更高)。
使用前两个改进可得到:
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。通过将标记以每个标记一个比特的方式打包存储,而不是使用布尔值,可以将此时间缩短至约15秒(这取决于您的计算机配置),因为减少的内存使用可以提供更好的缓存局部性。
通过打包标记、消除3的倍数并使用位运算进行更快速的计数,可以进一步优化。
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;
}
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亿个质数,这并不是难以忍受的长时间。
还有其他类型的质数筛法,特别有趣的是Atkin筛法,它利用了某些(有理)质数的同余类在某些二次扩张的代数整数环中是合成数的事实。这里不是展开数学理论的地方,简单来说,Atkin筛法的算法复杂度比Eratosthenes筛法低,因此对于大限制而言更可取(对于小限制,一个未过度优化的Atkin筛法具有更高的开销,因此可能比一个相当优化的Eratosthenes筛法慢)。
D.J. Bernstein的primegen库(用C编写)针对小于232的数字进行了良好的优化,在这里大约1.1秒内找到了第1亿个质数。
快速方法
如果我们只想找到第n个质数,那么也没有必要找到所有更小的质数。如果我们可以跳过大部分质数,那么就可以节省很多时间和工作量。给定一个良好的近似值a(n)来表示第n个质数p(n),如果我们有一种快速计算不超过a(n)的质数数量π(a(n))的方法,那么我们就可以筛选出在a(n)和p(n)之间少量缺失或过剩的质数。
我们已经看到了一个容易计算出较好的p(n)近似值,我们可以采用:
a(n) = n*(log n + log (log n))
举个例子。
计算素数计数函数
Meissel-Lehmer方法是一种很好的方法,它可以在大致上
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) 近似值以使最终筛选阶段所需时间与计算
π(a(n))
相当的问题。
一般情况下,在对初始近似值进行一到两次改进后,需要筛选的范围足够小,以便筛选阶段具有O(n^0.75)或更好的复杂度。
这种方法可以在大约40毫秒内找到第一亿个质数,以及在不到8秒钟内找到10^12-1=29996224275833。
简而言之:寻找第n
个质数可以高效地完成,但是你想要的越高效,就需要涉及更多的数学知识。
如果有人想尝试这些算法,我已经为大多数讨论的算法准备好了Java代码在这里。
对于过于热心的人的一则旁注:现代数学中使用的质数定义不同,适用于更广泛的情况。如果我们将学校定义调整为包括负数-因此,如果一个数字既不是1也不是-1,并且只能被1,-1,本身及其负数整除,则定义了ℤ的不可约元素(对于整数而言),然而,对于整数,质数和不可约元素的定义是相同的。