更严格的界限:
static const unsigned short primes_small[] = {0,2,3,5,7,11};
static unsigned long nth_prime_upper(unsigned long n) {
double fn = (double) n;
double flogn, flog2n, upper;
if (n < 6) return primes_small[n];
flogn = log(n);
flog2n = log(flogn);
if (n >= 688383) /* Dusart 2010 page 2 */
upper = fn * (flogn + flog2n - 1.0 + ((flog2n-2.00)/flogn));
else if (n >= 178974) /* Dusart 2010 page 7 */
upper = fn * (flogn + flog2n - 1.0 + ((flog2n-1.95)/flogn));
else if (n >= 39017) /* Dusart 1999 page 14 */
upper = fn * (flogn + flog2n - 0.9484);
else /* Modified from Robin 1983 for 6-39016 _only_ */
upper = fn * ( flogn + 0.6000 * flog2n );
if (upper >= (double) ULONG_MAX) {
/* Adjust this as needed for your type and exception method */
if (n <= 425656284035217743UL) return 18446744073709551557UL;
fprintf(stderr, "nth_prime_upper overflow\n"; exit(-1);
}
return (unsigned long) ceil(upper);
}
这些不应该小于实际的 nth_prime,应该适用于任何 64 位输入,并且比之前给出的 Robin 的公式(或 Wimblik 的复杂范围限制公式)一个数量级或更接近。对于我的使用,我有一个稍大的小素数表,因此可以进一步加强最后的估计。从技术上讲,我们可以使用 floor() 而不是 ceil() 的公式,但我担心精度。
编辑:改进这一点的另一个选择是实现良好的素数边界(例如 Axler 2014)并对它们进行二进制搜索。我的这种方法的代码比上面的代码要长约 10 倍(仍然运行在一毫秒以下),但可以将错误百分比降低一个数量级。
如果你想估计第 n 个素数,你可以这样做:
- Cipolla 1902(参见Dusart 1999第 12 页或本文。我发现三个项 (m=2) 加上一个三阶校正因子很有用,但如果项越多,它的振荡就太大了。维基百科链接中显示的公式就是这个公式(其中 m=2). 使用下面的两项逆 li 或逆黎曼 R 将给出更好的结果。
- 计算 Dusart 2010 的上限和下限并对结果进行平均。还不错,尽管我怀疑使用加权平均值会更好,因为界限并不同样严格。
- li^{-1}(n) 由于 li(n) 是素数的近似近似,因此逆是近似的 nth_prime 近似。这个,以及所有其他的,都可以作为函数的二分搜索相当快地完成。
- li^{-1}(n) + li^{-1}(sqrt(n))/4 更接近,因为这越来越接近 R(n)
- R^{-1} 逆黎曼 R 函数是我所知道的最接近的平均近似值,这是合理的。
最后,如果您有一个非常快速的素数计数方法,例如 LMO 实现之一(现在有三个开源实现),您可以编写一个快速精确的 nth_prime 方法。计算第 10^10 个素数可以在几毫秒内完成,第 10^13 个素数可以在几秒钟内完成(在现代快速机器上)。近似值在所有尺寸下都非常快,并且适用于更大的数字,但每个人对“大”的含义都有不同的想法。