13

我需要找到第 N 个素数的反向操作，即给定一个素数，我需要找到它在

2, 3, 5, 7...

素数可以很大，大约为10^7. 此外，还有很多。

我有一个预先计算的素数索引，可以进行二进制搜索，但我也有 50k 的空间限制！可以过筛吗？或者还有什么快速的方法？

编辑：非常感谢所有精彩的答案，我没想到他们！我希望它们对寻求相同的其他人有用。

4

7 回答 7

9

你的范围只有一千万，这对于这种东西来说是很小的。我有两个建议：

1) 以方便的间隔创建一个 pi(n) 表，然后使用分段的 Eratosthenes 筛来计算包含所需值的两个表条目之间的素数。间隔的大小决定了所需表的大小和计算结果的速度。

2) 使用勒让德的 phi(x,a) 函数和 Lehmer 的素数计数公式直接计算结果。phi 函数需要一些存储空间，我不确定到底需要多少。

在这两者中，考虑到您的问题规模，我可能会选择第一个替代方案。我的博客上提供了 Eratosthenes 分段筛法和Lehmer素数计数功能的实现。

编辑1：

经过反思，我有第三种选择：

3) 使用对数积分来估计 pi(n)。它是单调递增的，并且在您需要的时间间隔内始终大于 pi(n)。但差异很小，永远不会超过 200。因此，您可以预先计算所有小于 1000 万的值的差异，制作 200 个变化点的表格，然后在需要时计算对数积分并在桌子。或者你可以用黎曼的 R 函数做类似的事情。

第三种选择占用的空间最少，但我怀疑第一种选择所需的空间不会太大，而且筛分可能比计算对数积分更快。所以我会坚持我原来的建议。在我的博客中有对数积分和黎曼 R 函数的实现。

编辑2：

正如评论所表明的那样，这并没有很好地工作。请忽略我的第三个建议。

为了弥补我在提出一个不起作用的解决方案时的错误，我编写了一个程序，该程序使用一个 pi(n) 值表和一个分段的 Eratosthenes 筛来计算 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:]

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

sieve = [True] * 500

该count函数使用分段的 Eratosthenes 筛来计算从 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 是 p 在 lo 到 hi 范围内的最小倍数入筛的偏移量（注意不是 p 在范围内的最小倍数，而是 p 在范围内的最小倍数的偏移量的索引范围，这可能会造成混淆）。if当该语句引用一个完全平方的数字时，该语句递增 q。然后 j 上的循环从筛子中击出 p 的倍数。

我们count以两种方式使用该功能。第一次使用建立一个 pi(n) 值的表，该表是 1000 的倍数；第二种用途在表内插值。我们将表存储在一个全局变量 piTable 中：

piTable = [0] * 10000

我们根据原始请求选择参数 1000 和 10000 以将内存使用量保持在 50 KB 以内。（是的，我知道最初的发布者放宽了这个要求。但我们无论如何都可以兑现它。）一万个 32 位整数将占用 40,000 字节的存储空间，从 lo 到 hi 的 1000 范围内筛选只需要 500 个字节存储空间，速度非常快。您可能想尝试其他参数以查看它们如何影响程序的空间和时间使用。构建piTable是通过调用count函数一万次完成的：

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

到目前为止，所有计算都可以在编译时而不是运行时完成。当我在ideone.com上进行这些计算时，它们花费了大约 5 秒钟，但那段时间不算在内，因为当程序员第一次编写代码时，它可以永远完成一次。作为一般规则，您应该寻找机会将代码从运行时移动到编译时，以使您的程序运行得非常快。

唯一剩下的就是编写实际计算小于或等于 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 计算为小于请求 n 的 piTable 的最大索引，return 语句将 piTable 中的值与 table 值和请求值之间的素数计数相加；hi 限制是 n+1，否则在 n 是素数的情况下，它不会被计算在内。例如，说：

print pi(6543223)

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

该pi功能非常快。在ideone.com，对 pi(n) 的一千次随机调用在大约半秒内被计算出来，因此每个大约是半毫秒；这包括生成素数和求和结果的时间，因此实际计算 pi 函数的时间甚至不到半毫秒。这是我们在建造桌子上的投资的相当不错的回报。

如果您对使用素数编程感兴趣，我已经在我的博客上做了很多工作。请前来参观。

于 2013-01-02T18:32:44.470 回答
4

如果您先验地知道输入是素数，则可以使用近似值 pi(n) ≈ n / log n 和一个小的修正表来查找素数，其中四舍五入的结果不足以获得正确的值n. 我认为这是您在大小限制内的最佳选择，除了缓慢的蛮力方法。

于 2013-01-02T18:23:49.677 回答
4

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

于 2013-01-02T18:32:57.457 回答
2

这是一些有效的代码。您应该使用适用于您的输入范围的确定性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));
}
于 2013-01-02T19:44:02.440 回答
1

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

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

于 2013-01-02T20:13:39.707 回答
0

你的建议是最好的。预先计算（或下载）小于 10^7 的素数列表，然后对它们进行二进制搜索。

只有 664579 个小于 10^7 的素数，因此该列表将消耗约 2.7 MB 的空间。解决每个实例的二进制搜索将非常快速 - 只需约 20 次操作。

于 2013-01-02T18:35:28.443 回答
0

我就这样做过一次。写了一段代码，即给定n，可以很快找到第n个素数，最多n=203542528，所以约2e8。或者，它可以倒退，对于任何数字 n，可以知道有多少质数小于 n。

使用数据库。我将所有素数存储到某个点（我的上限的 sqrt）。在您的情况下，这意味着您将所有素数存储到 sqrt(1e7)。其中有 446 个，您可以以压缩形式存储该列表，因为到该点的最大差异仅为 34。超出该点，存储每个第 k 个素数（对于某个 k 值）。然后快速筛子就足够了在很短的时间间隔内生成所有素数。

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

nthprime(1e7)
ans =
   179424673

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

nthprime(1e7,1)
ans =
      664579

关键是，这样的数据库易于构建和搜索。如果你的数据库可以不超过50k，应该没有问题。

于 2013-01-03T00:26:53.980 回答