你的范围只有一千万,这对于这种东西来说是很小的。我有两个建议:
1) 以方便的间隔创建一个 pi(n) 表,然后使用分段的 Eratosthenes 筛来计算包含所需值的两个表条目之间的素数。区间的大小决定了所需表的大小和计算结果的速度。
2) 使用勒让德的 phi(x,a) 函数和 Lehmer 的素数计数公式直接计算结果。 phi 函数需要一些存储空间,我不确定具体需要多少。
考虑到您的问题规模,我可能会选择第一个替代方案。我的博客上提供了segmented Sieve of Eratosthenes 和Lehmer's 素数计数功能的实现。
编辑 1:
经过反思,我有第三种选择:
3) 使用对数积分来估计 pi(n)。它是单调递增的,并且在您需要的时间间隔内始终大于 pi(n)。但差异很小,永远不会超过 200。因此,您可以预先计算所有小于 1000 万的值的差异,制作 200 个变化点的表格,然后在需要时计算对数积分并在桌子。或者你可以用黎曼的 R 函数做类似的事情。
第三种选择占用的空间最少,但我怀疑第一种选择所需的空间不会太大,而且筛分可能比计算对数积分更快。所以我会坚持我原来的建议。 my blog 处有对数积分和黎曼 R 函数的实现。
编辑 2:
正如 cmets 所指出的那样,这并不是很好。请忽略我的第三条建议。
为了弥补我在提出一个不起作用的解决方案时的错误,我编写了一个程序,该程序使用 pi(n) 值表和 Eratosthenes 的分段筛来计算 n 的 pi(n) 值10000000。我将使用 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 个字节存储空间,速度非常快。您可能想尝试其他参数以查看它们如何影响程序的空间和时间使用。通过调用count函数一万次来构建piTable:
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 中的值与表值和请求值之间的素数计数相加; hi 限制是 n+1,否则在 n 是素数的情况下,它不会被计算在内。例如,说:
print pi(6543223)
将导致数字 447519 显示在终端上。
pi 函数非常快。在ideone.com,大约半秒内计算了对 pi(n) 的一千次随机调用,因此每个调用大约半毫秒;这包括生成素数和求和结果的时间,因此实际计算 pi 函数的时间甚至不到半毫秒。这对我们在建表方面的投资来说是一个相当不错的回报。
如果您对使用素数进行编程感兴趣,我在blog 上做了很多工作。请前来参观。