3.3 实现该算法并优化其代码
一开始我们可能认为要用两个数组才能实现该算法,其中一个数组保存有待筛选的数字,另一个数组保存对应的Boolean标志,用以表示相应的数字有没有划掉。然而只要稍微思考一下就会发现,我们实际上不需要存储有待筛选的数字。因为其中的大多数值(也就是那些非素数)根本用不到。如果确实需要用到某个值,那我们可以根据它的位置来计算。由于首个值是3,而且后一个值比前一个值大2,因此第i个值就是2i + 3。
于是,只需要把筛选过程中用到的Boolean标志保存下来就可以了,我们用true表示素数,用false表示合数,并且把划掉非素数的操作称为给筛子做记号(mark the sieve)。下面这个函数可以根据给定的素因子(factor)来标注与之相关的所有非素数:
上面这段代码“声明”了一些带有要求的模板参数,而这些要求则称为概念(concept),
我们将在第10章详细讲解此内容,大家现在可以先参考附录C,来熟悉这个术语。(如果你对C++的模板不太熟悉,那么也请查阅附录C,该附录详细解释了模板机制。)
稍后我们将会看到,在调用上面这个函数时,first参数会指向首个还未划掉的factor倍数,而根据刚才的讨论,大家可以知道,这个值就是factor的平方。对于last参数来说,我们将遵循STL的约定,把刚刚超过表格中最后一个元素的那个迭代器传给它,这样一来,last-first就可以表示元素的数量了。
-
- *
在开始编写筛选所用的代码之前,我们先考虑下面几条引理(lemma):
- *
对于合数c来说,其最小的素因子的平方值,肯定小于等于c。
任何一个比p2小的合数,都会为某个比p小的素数所划掉(也就是说,它是p的倍数)。
如果某一轮是以p为素因子进行筛选的,那么该轮应该从p2开始执行。
如果要找的是小于等于m的所有素数,那么当p2>m时,整个筛选过程就应该停止。
下面这两条公式,用来在筛选的过程中进行计算:
第i个元素的值:value (i) = 3 + 2i = 2i + 3
值为v的元素所对应的下标:index (v) =
对于第i个元素来说,其数值的k倍与k+2倍之间所间隔的元素数量是:
现在,我们可以试着来实现这个素数筛选函数了:
有人可能认为,既然这个筛选函数必须从头开始完整地运用到某个序列上面,那么它就应该接受一个指向某种数据结构的引用,那种数据结构里面含有由Boolean值所构成的序列。但实际上,我们并没有那么做,而是令调用者把一个指向某段范围开始处的迭代器,以及该范围的长度传进来,以便能够应对更多的数据结构。这些数据既可以存放在STL容器里面,也可以存放在内存块里面,我们并不对此做出限定。请注意,该函数的第二个参数指的是表格的大小n,而不是筛选的上限m。
当前这一轮所要划掉的首个数值,就是筛选所用的那个素因子的平方,而该值在表格中的索引则用index_square变量表示。值得注意的是,在每一轮循环时,我们可能都要用
i + i + 3这个表达式,来计算当前这一轮筛选所使用的素因子,而且还要用另外一些表达式来计算其他的量(这些内容在代码中以斜体标出)。为此,我们可以把每轮循环都要用到的那些式子提取到循环的外面。修改之后的代码,用粗体标出了相关的改动:
敏锐的读者可能会发现,修改之后的代码要对factor变量进行运算,这样做的效率实际上比修改之前还要差一些。修改代码之前,我们只需要在if测试为true的时候才去计算i + i + 3的值,而修改完之后则需要在每次执行循环时,都把这个式子计算一遍。但是大家稍后就会看到,单独用一个factor变量来表示这个式子是很有意义的。更大的问题在于,对index_square变量所做的运算,其开销相对来说比较高,因为要执行两次乘法才行。针对这个问题,我们可以从编译器的优化技术中获得灵感。有一种优化技术叫做强度折减(strength reduction),也就是用开销相对较低的运算,来等效地取代开销较高的运算,例如用加法来实现乘法。既然编译器可以执行这样的优化,那我们同样可以采用手动的办法来实现它。
现在就来详细考虑这两个式子。假如我们把下面两种运算:
那么其中的δfactor和δindex_square,就分别表示factor和index_square变量在前后两轮之间的差值(也就是第i轮和第i+1轮之间的差值)。这两个差值可以这样来计算:
δfactor这一部分很好处理,由于其中的变量i可以消去,因此只需要用常量2来表示它就行了。接下来主要应该考虑的是,怎样对计算δindex_square所需的表达式进行简化。我们对表达式中的各项重新整理之后,可以看出,这条表达式相当于factor(i)与factor(i + 1)的和,其中的factor(i)是当前这一轮已经计算好的,而factor(i + 1)则正是我们紧接着就要计算的。(如果你要在代码中计算多个值,那就应该像本例这样,看看能不能用其中一个值来表示另外的那些值,做到了这一点,或许就可以减少一些计算量。)
用常量2来替换δfactor,并通过变量factor来表示δfactor之后,我们就得到了最终版本的sift函数。这次也用粗体表示其中的改进之处:
习题3.2 用二进制位(以std::vector来实现)、uint8_t、uint16_t、uint32_t及uint64_t等大小不同的数据来衡量素数筛选算法所耗的时间。
习题3.3 通过素数筛选算法来绘制下列函数的图像:
π (n) =小于n的素数个数
n一直取到107,并找出此函数的解析近似(analytic approximation)函数。
如果一个素数的数位按照相反顺序排列之后和原数一样,那么这种素数就叫做回文素数(palindromic prime)。下面列出2至1000之间的素数,并用方框来标注其中的回文素数:
习题3.4 有没有大于1000的回文素数?[1000, 2000]这个区间内,为什么不会出现回文素数?如果我们不用十进制,而是改用十六进制,那么情况是否有所变化?如果改用以任意数n为底的计数方式(也就是n进制)呢?