你的范围只有一千万,这对于这种事情来说很小.我有两个建议:
1)以方便的间隔创建pi(n)表,然后使用分段的Eratosthenes筛来计算包含所需值的两个表条目之间的素数.间隔的大小决定了所需表的大小和计算结果的速度.
2)使用勒让德的phi(x,a)函数和Lehmer的计数公式直接计算结果.phi功能需要一些存储,我不确定究竟有多少.
在这两个中,我可能会根据您的问题大小选择第一个替代方案.我的博客上提供了Eratosthenes的分段筛和Lehmer的计数功能的实现.
编辑1:
经过反思,我有第三种选择:
3)使用对数积分来估计pi(n).它是单调增加的,并且在你需要的时间间隔内总是大于pi(n).但差异很小,从不超过200左右.因此,您可以预先计算所有小于千万的值的差异,制作200个变化点的表格,然后在请求时计算对数积分并查找校正因子表.或者你可以用Riemann的R函数做类似的事情.
第三种方案占用的空间最少,但我怀疑第一种方案所需的空间不会太大,筛分可能比计算对数积分更快.所以我会坚持原来的建议.我的博客上有对数积分和Riemann 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
Eratosthenes的筛子分为两部分.首先,从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并在全局ps变量中存储446个筛选素数:
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执行筛分的循环,依次取每个筛分素数.当筛分素数的平方大于范围的极限时,循环终止,因为所有质数将在该点被识别(我们需要比平方根更大的下一个素数的原因是为了有筛选素数停止循环).神秘变量q是在lo到hi范围内p的最小倍数的筛子的偏移(注意它不是范围中p的最小倍数,而是p的最小倍数的偏移的索引)范围,这可能令人困惑).当if语句引用一个完美正方形的数字时,该语句会增加q.然后j上的循环从筛子上击中p的倍数.
我们count以两种方式使用该功能.第一次使用建立一个pi(n)值的表,其值为1000的倍数; 第二次使用在表格内插入.我们将表存储在全局变量piTable中:
piTable = [0] * 10000
我们根据原始请求选择参数1000和10000,以将内存使用量保持在50千字节以内.(是的,我知道原始的海报放宽了这个要求.但我们无论如何都可以尊重它.)一万个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上进行这些计算时,它们花了大约五秒钟,但是这个时间不计算,因为当程序员第一次编写代码时,它可以一次完成.作为一般规则,您应该寻找将代码从运行时间移动到编译时间的机会,以使您的程序快速运行.
剩下的唯一事情就是编写实际计算小于或等于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,并且return语句将piTable中的值添加到表值和请求值之间的素数计数中; hi限制是n + 1,因为否则在n是素数的情况下它将不被计数.举个例子,说:
print pi(6543223)
将导致号码447519显示在终端上.
该pi功能是非常快的.在ideone.com上,在大约半秒钟内计算了一千次随机调用pi(n),因此大约半个毫秒; 其中包括生成素数的时间和结果的总和,因此实际计算pi函数的时间甚至不到半毫秒.这对我们建造桌子的投资来说是一个相当不错的回报.
如果你对使用素数进行编程感兴趣,我在博客上做了很多工作.请过来参观.