查找质数的位置

Question

我需要做的是反向寻找第N个质数，也就是说，给定一个质数，我需要找到它在质数序列中的位置。

2, 3, 5, 7...

质数可以很大，达到10^7的数量级。而且它们的数量很多。

我有一个预计算质数的索引，可以进行二进制搜索，但我还有50k的空间限制！能使用筛法吗？或者其他快速的方法？

编辑： 非常感谢所有精彩的回答，我没有期望这么好！我希望它们对寻找相同答案的其他人有用。

- max

---

首先，你确定了可能的答案集中有多少个吗？解决这个问题的方法有几种。空间/速度权衡、压缩和数学技巧（如孪生质数）。 - Clinton Pierce

当我再次检查时，小于10^7的素数有664570个。在预计算的素数表上进行二分搜索在这里不是一个选项。 - nhahtdh

参见我的stackoverflow答案 https://dev59.com/WG865IYBdhLWcg3wOcDy，关于存储质数表和快速筛法的思路。您应该能够每个字节获取 <s>一个</s> 两个质数（请参阅 http://en.wikipedia.org/wiki/Prime_gap）。二分搜索是不必要的，您可以根据P(n)≈nlogn估计在表中开始查找的位置，在那里进行线性搜索。 - Potatoswatter

我认为筛法在这里不太实用。你可以看看达里克·亨利·勒默的计算素数方法：http://en.wikipedia.org/wiki/Prime-counting_function 他在1959年就能够在IBM 701上进行10^10次计算。 - hexist

非常感谢大家。@ColonelPanic，我所说的50k实际上是指源文件大小限制，这是一个在线评测问题的子问题！ - max

Answer 1

你的范围只能达到一千万，对于这种情况来说太小了。我有两个建议：

1）在方便的间隔处创建一个pi(n)表，然后使用分段埃拉托色尼筛选法来计算两个表项之间的质数总数，这些表项包含所需值。区间的大小决定所需表格的大小和计算结果的速度。

2）使用Legendre的phi(x，a)函数和Lehmer的计算素数公式直接计算结果。phi函数需要一些存储空间，但我不确定具体需要多少。

在这两个选择中，考虑到你的问题规模，我可能会选择第一个方案。我的博客上提供了分段埃拉托色尼筛选法和Lehmer's计算素数函数的实现。

编辑1：

经过思考，我有第三个选择：

3) 使用对数积分来估算pi(n)。它是单调递增的，并且始终大于您需要的区间内的pi(n)。但是差异很小，从不超过200。因此，您可以预先计算所有小于一千万的值的差异，制作一个包含200个变化点的表格，然后在请求时计算对数积分并查找表中的校正因子。或者您可以使用黎曼R函数进行类似的操作。

第三种选择所需的空间最少，但我怀疑第一种选择所需的空间也不会太大，并且筛法可能比计算对数积分更快。因此，我将坚持我的最初建议。这里有对数积分和黎曼R函数的实现：我的博客。

编辑2：

那并没有很好地发挥作用，正如评论所指出的那样。请忽略我的第三个建议。

作为我因提出一个无法工作的解决方案而做出的赎罪，我编写了一个程序，使用pi(n)值表和分段埃拉托色尼筛法计算n < 10000000的pi(n)值。出于教学目的，我会使用Python而不是原帖中要求的C，因为Python更简单易读。

我们首先计算小于一千万平方根的筛选质数；这些质数将用于构建pi(n)值表和执行计算最终答案的筛子。一千万的平方根是3162.3。我们不想使用2作为筛选质数——我们只对奇数进行筛选，并将2作为特殊情况处理——但我们确实需要比平方根大的下一个质数，以便筛选质数列表永远不会耗尽（这会导致错误）。所以我们使用这个非常简单的埃拉托色尼筛法版本来计算筛选质数：

def primes(n):
    b, p, ps = [True] * (n+1), 2, []
    for p in xrange(2, n+1):
        if b[p]:
            ps.append(p)
            for i in xrange(p, n+1, p):
                b[i] = False
    return ps

埃拉托斯特尼筛法分为两个部分。首先，从2开始，列出小于目标数的数字列表。然后，重复遍历该列表，从第一个未被划掉的数字开始，并将该数字的所有倍数从列表中划掉。最初，2是第一个未被划掉的数字，所以划掉4、6、8、10等数字。然后3是下一个未被划掉的数字，所以划掉6、9、12、15等数字。然后4作为2的倍数被划掉，下一个未被划掉的数字是5，所以划掉10、15、20、25等数字。继续进行，直到考虑完所有未被划掉的数字；剩下未被划掉的数字是素数。循环p逐个考虑每个数字，如果它未被划掉，则循环i会划掉它的倍数。 primes函数返回一个447个素数的列表：2、3、5、7、11、13、...、3121、3137、3163。我们从列表中去除2，并将446个筛选素数存储在全局ps变量中。

ps = primes(3163)[1:]

我们需要的主要功能是在一个范围内计算质数。它使用筛法，我们将其存储在全局数组中，以便可以重复使用，而不是在每次调用计数函数时重新分配：

sieve = [True] * 500

count函数使用分段埃拉托色尼筛法来计算从lo到hi（lo和hi都包含在范围内）的质数数量。该函数有四个for循环：第一个清除筛子，最后一个计算质数，另外两个以类似于上面展示的简单筛法的方式进行筛选：

def count(lo, hi):
    for i in xrange(500):
        sieve[i] = True
    for p in ps:
        if p*p > hi: break
        q = (lo + p + 1) / -2 % p
        if lo+q+q+1 < p*p: q += p
        for j in xrange(q, 500, p):
            sieve[j] = False
    k = 0
    for i in xrange((hi - lo) // 2):
        if sieve[i]: k += 1
    return k

该函数的核心是循环for p in ps，它执行筛选操作，依次取每个筛选素数p。当筛选素数的平方大于范围限制时，循环终止，因为此时所有素数都将被识别出来（我们需要比平方根大的下一个素数，以便有筛选素数停止循环）。神秘变量q是偏移量，指的是在范围lo到hi内p的最小倍数的筛子中的位置（请注意，它不是范围内最小的p的倍数，而是范围内p的最小倍数的偏移量的索引，这可能会令人困惑）。当引用完全平方数时，if语句会增加q。然后，j的循环会从筛子中删除p的倍数。

我们使用count函数有两种方式。第一种用法是在1000的倍数处建立pi(n)值的表；第二种用法是在表中进行插值。我们将表存储在全局变量piTable中：

piTable = [0] * 10000

我们选择参数1000和10000，基于原始请求以保持内存使用量在50千字节以内。(是的，我知道原帖作者放宽了这个要求。但我们仍然可以遵守。)一万个32位整数将占用40000字节的存储空间，而仅从lo到hi筛选1000的范围只需要500字节的存储空间，并且速度非常快。您可能想尝试其他参数，以查看它们如何影响程序的空间和时间使用情况。通过调用count函数10,000次来构建piTable。

for i in xrange(1, 10000):
    piTable[i] = piTable[i-1] + \
        count(1000 * (i-1), 1000 * i)

到目前为止，所有的计算都可以在编译时而不是运行时完成。当我在ideone.com上进行这些计算时，它们大约需要五秒钟，但这个时间不计算在内，因为程序员在编写代码时可以一次性完成。通常，您应该寻找将代码从运行时移动到编译时的机会，以使您的程序运行非常快。

唯一剩下的就是编写实际计算小于或等于n的质数数量的函数：

def pi(n):
    if type(n) != int and type(n) != long:
        raise TypeError('must be integer')
    if n < 2: return 0
    if n == 2: return 1
    i = n // 1000
    return piTable[i] + count(1000 * i, n+1)

第一个if语句进行类型检查。第二个if语句针对荒谬的输入返回正确响应。第三个if语句特殊处理2；我们的筛法使1成为质数，2成为合数，两者都是不正确的，所以我们在这里进行修正。然后，i被计算为piTable中小于请求n的最大索引，而返回语句将piTable的值加到表值和请求值之间的质数计数中；上限hi为n+1，因为如果n是质数，则不会被计算。例如，说：

print pi(6543223)

这将导致数字447519在终端上显示。

pi函数非常快速。在ideone.com，对pi(n)的一千个随机调用在大约半秒钟内计算完成，因此每个调用大约需要半毫秒; 这包括生成质数和求和结果的时间，因此实际计算pi函数的时间甚至少于半毫秒。这是我们在构建表格方面投资的相当不错的回报。

如果您对使用质数编程感兴趣，我在我的博客上做了很多工作。请来访问。

- user448810

---

不，我没有。感谢您的纠正。无论如何，更改的数量很少，因此这只需要一点空间。然而，我仍然认为筛法是更可取的。 - user448810

那比我想象的要多得多。我同意这不会起作用。那我们就坚持使用分段筛法吧。 - user448810

我喜欢#1，它让我想起了用巨步-小步算法解决不同问题的方法。就像那个算法一样，我期望你可以获得大约sqrt(N)的存储和工作量。然而，有一个问题：如何找到两个表项来框定该值？我认为你可以使用素数定理来猜测，然后进行修正（也许只需线性探测即可找到框定条目）。 - President James K. Polk

@GregS：假设你有10^3个表项，间隔为10^4。要找到小于567890的质数数量，请筛选从560000到570000的区间，计算从560000到567890的质数数量，并加上表中pi（560000）的值。 - user448810

@GregS：你可以通过简单的除法来完成，就像user448810所暗示的那样。（请查看我的答案代码。） - tmyklebu

Answer 2

如果你事先知道输入是质数，你可以使用近似公式 pi(n) ≈ n / log n，加上一张针对质数的小修复表格，在结果四舍五入不足以得到正确的值n时使用。我认为这是在大小限制内最好的选择，除了缓慢的暴力方法。

- R.. GitHub STOP HELPING ICE

---

那个表的大小非常关键。质数定理主要适用于大数，它可能在这里不太适用，但只有尝试才能告诉我们。 - Potatoswatter

我已经有一段时间没玩质数定理的数字了，但我认为这值得一试。 - R.. GitHub STOP HELPING ICE

我认为他正在寻找一个确切的数字，并且已知有好的方法来计算它（Derrick Henry Lehmer的方法对我来说看起来最好，Legendre也有一个看起来在这些限制下可能可行）。 - hexist

这行不通。质数的分布太复杂了。素数定理告诉你，pi(n)渐近于n/log(n)。这意味着它们之间的乘法差异趋近于零。你需要找到一个与pi(n)相差一个加性常数的东西，才有希望让修复方法奏效。（而n/log(n)并不在pi(n)的加性常数范围内。） - tmyklebu

显然，我的方法对于计算pi（n）没有用处。然而，OP的问题是计算pi^-1(p)，似乎在假设输入p为质数的情况下。我仍然不清楚我的提议是否可以用于解决这个问题。 - R.. GitHub STOP HELPING ICE

Answer 3

我建议在这里使用一种启发式混合模型。存储每第n个质数，然后通过质性测试进行线性搜索。为了加速，您可以使用快速简单的质性测试（例如Fermat测试，其中a==2），并预先计算错误的结果。基于输入的最大大小和存储限制的微调应该很容易解决。

- user295691

---

如果按照您所说的实施，这两个数字之间大约会有60个质数。 - nhahtdh

还有大约800个小于10^8的通过测试的非质数。计算2^n (mod n)是O(log n)的；在实际应用中，这意味着测试约1000个数字以获取计数。这种方法更适合一次性测试而不是批量测试，但仍应该很快。 - user295691

这是一个好主意。您可以使用Miller-Rabin的确定性变体来使其更好，并避免出现错误和Carmichael数; 检查证人2、7和61对于任何小于几十亿的输入都有效。 - tmyklebu

这是一个时间与空间的权衡；确定性费马测试非常快速（而且实现简单），很难被超越，而存储2064个数字的成本非常小。（另外，结果表明2064是小于10^8的正确数字；其中850个小于10^7。如果我们测试2和3，那么它会降至492）。 - user295691

Miller-Rabin同样快速，而且您不需要“备用”表查找或试除法。使用Miller-Rabin而不是直接使用Fermat测试没有理由。 - tmyklebu

我的错误——我假设米勒-拉宾算法所做的额外工作会使其平均速度变慢，但一些实验表明，早期终止可能性加速了事情，足以弥补这一点。 - user295691

Answer 4

这里有一些可用的代码。你应该使用一个适用于你输入范围的确定性Miller-Rabin测试来替换基于试除法的素数测试。在适当的小范围内筛选素数比试除法更好，但这是朝着错误的方向迈出的一步。

#include <stdio.h>
#include <bitset>
using namespace std;

short smallprimes[549]; // about 1100 bytes
char in[19531]; // almost 20k

// Replace me with Miller-Rabin using 2, 7, and 61.
int isprime(int j) {
 if (j<3) return j==2;
 for (int i = 0; i < 549; i++) {
  int p = smallprimes[i];
  if (p*p > j) break;
  if (!(j%p)) return 0;
 }
 return 1;
}

void init() {
 bitset<4000> siv;
 for (int i = 2; i < 64; i++) if (!siv[i])
  for (int j = i+i; j < 4000; j+=i) siv[j] = 1;
 int k = 0;
 for (int i = 3; i < 4000; i+=2) if (!siv[i]) {
  smallprimes[k++] = i;
 }

 for (int a0 = 0; a0 < 10000000; a0 += 512) {
  in[a0/512] = !a0;
  for (int j = a0+1; j < a0+512; j+=2)
   in[a0/512] += isprime(j);
 }
}

int whichprime(int k) {
 if (k==2) return 1;
 int a = k/512;
 int ans = 1 + !a;
 for (int i = 0; i < a; i++) ans += in[i];
 for (int i = a*512+1; i<k; i+=2) ans += isprime(i);
 return ans;
}

int main() {
 int k;
 init();
 while (1 == scanf("%i", &k)) printf("%i\n", whichprime(k));
}

- tmyklebu

---

我喜欢这个——在固定间隔内存储小于x的质数数量表比我的方法存储质数本身更有效，因为查找是微不足道的。 - user295691

Answer 5

以下是您正在寻找的内容。http://www.geekviewpoint.com/java/numbers/index_of_prime。在那里，您将找到代码和单元测试。由于您的列表相对较小（即10^7），因此它应该可以处理。

基本上，您需要找到2和n之间的所有质数，然后计算小于n的所有质数以查找索引。此外，如果n不是质数，则函数返回-1。

- Konsol Labapen

Answer 6

我曾经做过这件事。编写了一段代码，可以快速找到第n个质数，直到n = 203542528，大约2e8。或者，它可以向后查找，对于任何数字n，可以告诉有多少个质数小于n。

使用数据库。我存储了所有的质数，直到某个点（我的上限的平方根）。在你的情况下，这意味着你将存储所有的质数，直到sqrt(1e7)。有446个质数，你可以以压缩形式存储该列表，因为到那个点为止的最大差异仅为34。超过那个点，存储每k个质数（对于某个k值）。然后，一个快速筛法就足以在短时间内生成所有质数。

所以在MATLAB中，要找到第1e7个质数：

nthprime(1e7)
ans =
   179424673

或者，它可以找到小于1e7的质数数量：

nthprime(1e7,1)
ans =
      664579

关键是，这样的数据库易于构建和搜索。如果你的数据库不超过50k，那么就没有问题了。

- user85109

Answer 7

您建议的是最好的。预先计算（或 下载）小于10^7的质数列表，然后进行二分搜索。

小于10^7的质数只有664579个，因此列表将占用约2.7 MB的空间。每个实例的二分搜索将非常快速-仅需大约20次操作。

- Colonel Panic

---

但他说他有一个50k的空间限制。 - im so confused

“OP写道”（https://dev59.com/JmYq5IYBdhLWcg3w30bc#G47jnYgBc1ULPQZFfrXo）：“我实际上是指源文件限制为50k”。 - Colonel Panic