c - 求素数的位置

Question

我需要做与找到第 N 个素数相反的事情，即给定一个素数，我需要找到它在

2, 3, 5, 7...

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

我有一个可以进行二分搜索的预先计算的素数索引，但我也有 50k 的空间限制!可以做筛子吗？或者有其他快捷方式吗？

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

Answer 1

你的范围只有一千万，这对于这种事情来说很小。我有两个建议:

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

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

在这两者中，鉴于您的问题规模，我可能会选择第一个选项。 segmented Sieve of Eratosthenes 的实现和 Lehmer's素数计数功能可在我的博客上找到。

编辑 1:

经过反射(reflection)，我有第三种选择:

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

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

编辑2:

正如评论所表明的那样，这并不奏效。请忽略我的第三个建议。

为了弥补我在提出一个不起作用的解决方案时所犯的错误，我编写了一个程序，该程序使用一个 pi(n) 值表和一个分段的 Eratosthenes 筛来计算 n < 10000000 时的 pi(n) 值。我将使用 Python，而不是原始海报要求的 C，因为 Python 更简单，更易于阅读以用于教学目的。

我们首先计算小于一千万平方根的筛选素数；这些素数将用于构建 pi(n) 的值表和执行计算最终答案的筛选。千万的平方根是3162.3。我们不想使用 2 作为筛选素数——我们将只筛选奇数，并将 2 视为特例——但我们确实希望下一个素数大于平方根，这样列表筛选素数永远不会用尽(这会导致错误)。所以我们使用这个非常简单的 Eratosthenes 筛法来计算筛分素数:

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 是 lo 到 hi 范围内 p 的最小倍数进入筛子的偏移量(注意不是范围内 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 个字节存储量会非常快。您可能想尝试其他参数以查看它们如何影响程序的空间和时间使用。 build piTable是通过拨打 count 来完成的功能一万次:

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计算为小于请求n的piTable的最大索引，return语句将piTable中的值与表值和请求值之间的素数计数相加； hi 限制是 n+1，因为否则在 n 是素数的情况下它不会被计算在内。举个例子，说:

print pi(6543223)

将导致数字 447519 显示在终端上。
pi功能非常快。在 ideone.com , 大约半秒内计算了一千个对 pi(n) 的随机调用，所以每个调用大约半毫秒；这包括生成素数和对结果求和的时间，因此实际计算 pi 函数的时间甚至不到半毫秒。这是我们 build table 的投资的相当不错的返回。

如果您对使用质数编程感兴趣，我已经在 blog 上做了很多工作。 .请来参观。