1.2 函数与它们产生的计算

至此我们已经考虑了程序设计中的一些要素。我们用过各种基本算术操作，知道如何组合使用这些运算，还考虑了通过把组合运算声明为复合函数，实现对复合运算的抽象。但是，即使知道了这些，我们还不能说自己已经理解了如何去编写程序。我们现在的状况很像初学象棋的人，他们已经知道了移动棋子的规则，但还不知道典型的开局、战术和策略。与之类似，我们现在还不知道本领域中有用的常见模式，缺少有关棋步价值的知识（哪些函数值得声明），缺少对所走棋步（执行一个函数）的后果做出预期的经验。

能看清所考虑的动作的后果，这种能力对于成为程序设计专家是至关重要的，就像类似能力在所有综合性的创造性活动中的作用一样。举个例子，要成为专业摄影师，就必须学会如何去观察景象，知道在各种可能的曝光和显影选择下，景象中的各个区域在影像中的明暗情况。只有在具有了这些能力后，我们才能去做反向推理，对得到所需效果应该做的取景、亮度、曝光和显影等做出规划。同样，在程序设计领域，我们需要对计算过程中各种动作的进行情况做出规划，用一个程序去控制这种过程的进展。要想成为专家，我们必须学习如何看清各种不同类型的函数产生的计算过程。只有掌握了这种技能，我们才能学会如何可靠地构造出程序，使之能表现出我们期望的行为。

一个函数就是一个计算过程的 局部演化 的模式，它描述了这个过程中的每个步骤如何在前面步骤的基础上进行。在用函数描述了一个计算过程的局部演化过程之后，我们当然希望能做出一些有关该计算过程的整体或 全局 行为的论断。一般而言这是非常困难的，但我们至少还是可以试着去描述这种过程演化的一些典型模式。

在这一节，我们要考察由一些简单函数产生的计算过程的“形状”，还要研究这些计算过程消耗各种重要计算资源（时间和空间）的速率。这里考察的函数都很简单，它们的用途就像摄影术中的测试模式，都是极度简化的，并不是很实际的例子。

1.2.1 线性递归和迭代

我们从阶乘函数开始，这种函数的定义是：

n != n ·（ n -1）·（ n -2）…3·2·1

存在许多计算阶乘的方法，一种最简单的方法就是利用下面的认识，对任意正整数 n ， n !就等于 n 乘以（ n -1）!：

n != n ·[（ n -1）·（ n -2）…3·2·1]= n ·（ n -1）！

这样，我们就能通过计算（ n -1）!，再将其结果乘以 n 的方式计算 n !。如果再注意到1!就是1，这些认识就可以直接翻译成一个计算机函数了：

我们可以利用1.1.5节介绍的代换模型，观察这个函数在计算6!时的行为，如图1.3所示。

现在我们从另一个不同的角度考虑阶乘计算。计算阶乘 n !的规则可以改述为：首先乘起1和2，然后把结果乘以3，然后再乘以4，这样下去直至达到 n 。更形式化地说，我们维持一个变动的乘积product和一个从1到 n 的计数器counter，这个计算过程可以用counter和product的如下变化规则描述。在计算中的每一步，它们同时按下面的规则改变：

可以看到， n !也就是计数器counter超过 n 时乘积product的值。

同样，我们又可以根据这个描述构造出一个计算阶乘的函数 ^[3] ：

与前面一样，我们也可以用代换模型查看6!的计算过程，如图1.4所示。

比较上面的两个计算过程。从一个角度看它们没什么不同：两者计算的是同一个定义域上的同一个数学函数，算出 n !都需要使用与 n 成正比的步骤数。实际上，这两个计算过程甚至采用了同样的乘法运算序列，得到同样的部分乘积序列。但在另一方面，如果我们考虑这两个计算过程的“形状”，就可以看到它们的进展情况大相径庭。

考虑第一个计算过程，代换模型揭示出一种先逐步展开而后收缩的形状，如图1.3中的箭头所示。在展开阶段，该计算过程构造起一个 推迟进行的运算 形成的链条（在这里是一个乘法的链条），收缩阶段则表现为这些运算的实际执行。这种由推迟执行的操作链条刻画的计算过程，称为 递归计算过程 。为了执行这种计算过程，解释器需要维护好以后将要执行的那些操作的轨迹。在阶乘 n !的计算中，推迟执行的乘法链条的长度就是记录其轨迹需要保存的信息量。这个长度随着 n 的值而线性增长（正比于 n ），与计算中的步骤数目一样。这种计算过程称为 线性的递归计算过程 。

与之相对，第二个计算过程里没有任何增长或收缩。对任意的 n ，在计算过程中的每一步，在我们需要保存的轨迹里，所有的东西就是名字product、counter和max_count的当前值。我们称这种过程为 迭代计算过程 。一般而言，迭代计算过程就是那种其中的状态可以用固定数目的 状态变量 描述的计算过程；而与此同时，又存在一套固定的规则，它们描述了该计算过程从一个状态到下一状态转换时这些变量的更新方法；还有一个（可能有的）结束检测，它说明了计算过程应该终止的条件。另外，在 n !的计算中，所需计算步骤随着 n 线性增长，因此，这种过程称为 线性的迭代计算过程 。

我们还可以从另一个角度对比这两个计算过程。在迭代的情况里，在计算过程中的任何一点，上述的几个程序变量都提供了有关计算状态的一个完整描述。如果我们令这一计算在某两个步骤之间停下来，再想重新唤醒这一计算，只需要为解释器提供这三个变量的值。而对于递归计算过程，还存在一些额外的“隐含”信息，它们并没有保存在程序变量里，而是由解释器维持。这些信息表明了在所推迟的操作形成的链条里漫游的过程中“计算过程正处于何处”。这个链条越长，需要保存的信息也就越多 。

在比较迭代与递归时，我们必须当心，不要搞混了递归 计算过程 的概念和递归 函数 的概念。当我们说一个函数是递归的时，讨论的是一个语法形式上的事实，说明在这个函数的声明中（直接或间接地）引用了该函数本身。而在说某一计算过程具有某种模式时（例如，线性递归），我们说的是这一计算过程的进展方式，而不是相应函数书写上的语法形式。我们说某个递归函数（例如fact_iter）将产生迭代计算过程时，可能会有人感到不舒服。然而，这一计算过程确实是迭代的，因为其状态由它的三个状态变量完全刻画，在执行这一计算过程时，解释器只需要维持这三个名字的轨迹就足够了。

区分计算过程和函数可能使人感到困惑，一个原因来自各种常见语言（包括C、Java和Python）的情况。这些语言的大多数实现在解释任何递归函数时，消耗的存储量总是与函数调用的次数成正比，即使从原理上看该函数描述的计算过程是迭代的。作为这一事实的后果，要想在这些语言里描述迭代计算过程，就必须借助某些特殊的“循环结构”，如do、repeat、until、for和while等。我们将在第5章考察的JavaScript实现则没有这一缺陷，它总能在常量空间中执行迭代计算过程，即使这一计算是用一个递归函数描述的。具有这种特性的实现称为 尾递归 。有了尾递归的实现，我们就可以采用常规的函数调用机制描述迭代，这也使得各种复杂的专用迭代结构变成不过是一些语法糖了 。

练习1.9 下面两个函数都是基于函数inc（它得到参数加1）和dec（它得到参数减1）声明的，它们各定义了一种得到两个正整数之和的方法。

请用代换模型描绘这两个函数在求值plus（4, 5）时产生的计算过程。这些计算过程是递归的或者迭代的吗？

练习1.10 下面的函数计算一个称为Ackermann函数的数学函数：

下面各语句的值是什么：

请考虑下面的函数，其中的A就是上面声明的函数：

请给出函数f、g和h对给定整数值 n 计算的函数的简洁数学定义。例如， k （ n ）计算5 n ² 。

1.2.2 树形递归

另一种常见计算模式称为 树形递归 。作为例子，现在考虑斐波那契（Fibonacci）数序列的计算，这一序列中的每个数都是前面两个数之和：

0, 1, 1, 2, 3, 5, 8, 13, 21,…

一般而言，斐波那契数由下面的规则定义：

我们立刻就能把这个定义翻译为一个计算斐波那契数的递归函数：

现在考虑这一计算的模式。为了计算fib（5），我们需要计算fib（4）和fib（3）。而为了计算fib（4），又需要计算fib（3）和fib（2）。一般而言，这一演化过程看起来像一棵树，如图1.5所示。请注意，这里的分支在每一层分裂为两个分支（除了最下面的），反映出对fib函数的每个调用里两次递归调用自身的事实。

上面这个函数作为典型的树形递归很有教育意义，但对计算斐波那契数而言，它却是一种很糟糕的方法，因为做了太多的冗余计算。从图1.5中可以看到，对fib（3）的全部计算重复做了两次，差不多占到所有工作的一半。事实上，不难证明，在整个计算过程中，计算fib（1）和fib（0）的次数（一般来说，也就是上面这种树里的树叶数）正好是Fib（ n +1）。要感受一下这种情况有多糟，我们可以证明Fib（ n ）值的增长相对于 n 是指数的。更准确地说（见练习1.13），Fib（ n ）等于最接近 的那个整数，其中：

这也就是 黄金分割 的值，它满足方程：

ϕ ² = ϕ +1

这样，这个计算过程中所用的步骤数将随着输入的增长而指数增长。在另一方面，其空间需求只是随着输入的增长而线性增长，因为，在计算中的每一点，我们都只需要保存树中在此结点之上的结点的轨迹。一般而言，树形递归计算过程里所需要的步骤数正比于树中的结点数，其空间需求正比于树的最大深度。

我们也可以规划出一种计算斐波那契数的迭代计算过程，其基本想法就是用一对整数 a 和 b ，把它们分别初始化为Fib（1）=1和Fib（0）=0，然后反复地同时应用下面的变换规则：

a ← a + b

b ← a

不难证明，在 n 次应用这些变换后， a 和 b 将分别等于Fib（ n +1）和Fib（ n ）。因此，我们可以定义下面这两个函数，以迭代方式计算斐波那契数：

这种计算Fib（ n ）的方法是线性迭代。这两种方法在所需步数上差异巨大——后一个相对 n 为线性的，前一个增长得如Fib（ n ）一样快——即使不大的输入也可能造成巨大差异。

然而，我们也不应该做出结论说树形递归计算过程根本没用。如果需要考虑操作具有层次结构的数据（而不是操作数值），我们可能发现树形递归计算过程是一种很自然的、威力强大的工具 。即使对于数的计算，树形递归计算过程也可能帮助我们理解和设计程序。以计算斐波那契数的程序为例，虽然第一个fib函数远比第二个低效，但它却更直截了当，几乎就是把斐波那契序列的定义直接翻译为JavaScript。而要规划出对应的迭代算法，就需要注意到该计算过程可以重塑为一个采用三个状态变量的迭代。

实例：换零钱方式的统计

要找到迭代的斐波那契算法，只需要不多的智慧。下面这个问题的情况与此不同：我们有一些50美分、25美分、10美分、5美分和1美分的硬币，要把1美元换成零钱，一共有多少种不同的方式？更一般的问题是，给了任意数量的美元现金，我们能写出一个程序，计算出所有换零钱方式的数目吗？

如果采用递归函数，这个问题有一种很简单的解法。假定我们把考虑的可用硬币类按某种顺序排好，于是就有下面的关系。

把总数为 a 的现金换成 n 种硬币的不同方式的数目等于

现金 a 换成除去第一种硬币之外的所有其他硬币的不同方式数目

加上现金 a-d 换成所有种类的硬币的不同方式数目，其中 d 是第一种硬币的币值

要弄清为什么这一说法正确，请注意这里把换零钱的方式分成两组，第一组的换法都没用第一种硬币，而第二组都用了第一种硬币。显然，换零钱的全部换法的数目，就等于完全不用第一种硬币的换法数，加上用了第一种硬币的换零钱的换法数。而后一个数目也就等于去掉一个第一种硬币值后，将剩下的现金换零钱的换法数。

这样，我们就把某个给定现金数的换零钱方式的问题，递归地归约为对更少现金数或者更少硬币种类的同一个问题。请仔细考虑上面的归约规则，设法使自己确信，如果采用下面的方法处理各种退化情况，我们就能利用上述规则写出一个算法 ：

●如果 a 就是0，应该算作有1种换零钱的方式。

●如果 a 小于0，应该算作有0种换零钱的方式。

●如果 n 是0，应该算作有0种换零钱的方式。

我们很容易把这些规则翻译成一个递归函数：

（函数first_denomination以可用硬币的种类数作为输入，返回第一种硬币的币值。这里认为硬币已经从最小到最大排好了，其实采用任何顺序都可以。）我们现在就能回答开始的问题了，下面是1美元换硬币的不同方法的数目：

函数count_change产生树形递归的计算过程，其中的冗余计算与前面fib的第一种实现类似。但另一方面，想设计一个能算出同样结果的更好的算法，就不那么明显了。我们把这个问题留给读者作为一个挑战。可以看到，树形递归的计算过程可能极其低效，但常常很容易描述和理解，这导致有人提出能否利用这两个世界里最好的东西设计一种“聪明编译器”，使之能把树形递归的函数翻译为能完成同样计算的更高效的函数 。

练习1.11 函数 f 由如下规则定义：如果 n <3，那么 f （ n ）= n ；如果 n ≥3，那么 f （ n ）= f （ n -1）+2 f （ n -2）+3 f （ n -3）。请写一个JavaScript函数，它通过一个递归计算过程计算 f 。再写一个函数，通过迭代计算过程计算 f 。

练习1.12 下面的数值模式称为 帕斯卡三角形 ：

三角形两个斜边上的数都是1，内部的每个数是位于它上面的两个数之和 ^[4] 。请写一个函数，它通过一个递归计算过程计算帕斯卡三角形。

练习1.13 证明Fib（ n ）是最接近 的整数，其中 。提示：利用归纳法和斐波那契数的定义（见1.2.2节），证明 ，其中 。

1.2.3 增长的阶

前面的例子说明，不同的计算过程，在消耗计算资源的速率上可能存在巨大差异。描述这种差异的一种方便方法是使用 增长的阶 记法，它能帮助我们理解在输入变大时，某一计算过程所需要的资源的粗略度量情况。

令 n 是一个参数，作为问题规模的一种度量，令 R （ n ）是一个计算过程在处理规模为 n 的问题时所需要的资源量。在前面的例子里，我们取 n 为给定函数需要计算的那个数值，当然也存在其他可能性。例如，如果我们的目标是计算数的平方根的近似值，那么就可以取 n 为所需精度的数字个数。对矩阵乘法，我们可以取 n 为矩阵的行数。一般而言，总存在某个与问题的特性有关的数值，使我们可以相对于它去分析给定的计算过程。类似的， R （ n ）也可以是所用内部存储寄存器的个数，也可能是需要执行的机器操作数目的度量值，或者其他类似的东西。对于每个单位时间只能执行固定数目的操作的计算机，所需的时间将正比于需要执行的基本机器指令条数。

我们称 R （ n ）具有 Θ （ f （ n ））的增长阶，记为 R （ n ）= Θ （ f （ n ））（读作“ f （ n ）的theta”），如果存在与 n 无关的整数 k ₁ 和 k ₂ ，使得：

k ₁ f （ n ）≤ R （ n ）≤ k ₂ f （ n ）

对任意足够大的 n 值成立（换句话说，对足够大的 n ， R （ n ）总位于 k ₁ f （ n ）和 k ₂ f （ n ）之间）。

举例说，在1.2.1节讨论的计算阶乘的线性递归计算过程里，操作步数的增长正比于输入 n 。也就是说，这一计算过程所需步骤的增长为 Θ （ n ），其空间需求的增长也是 Θ （ n ）。对迭代的阶乘，其步数还是 Θ （ n ）而空间是 Θ （1），即为常量 。树形递归的斐波那契计算需要 Θ （ ϕ ⁿ ）步和 Θ （ n ）空间，这里的ϕ就是1.2.2节里说的黄金分割率。

增长的阶为我们提供了对计算过程的行为的一种粗略描述。例如，某个计算过程需要 n ² 步，另一计算过程需要1000 n ² 步，还有一个计算过程需要3 n ² +10 n +17步，它们的增长的阶都是 Θ （ n ² ）。但另一方面，增长的阶也为我们在问题的规模改变时，预计一个计算过程的行为变化提供了有用的线索。对一个 Θ （ n ）（线性）的计算过程，规模增大一倍将使它所用的资源增加大约一倍。而对一个指数的计算过程，问题规模每增加1都将导致所用资源按某个常数倍增长。在1.2节的剩下部分，我们将考察两个算法，其增长的阶都是对数的，因此，当问题规模增大一倍时，所需资源量只增加一个常数。

练习1.14 请画出一棵树，展示1.2.2节的函数count_change把11美分换成硬币时产生的计算过程。相对于被换现金量的增加，这一计算过程的空间和步数增长的阶各是什么？

练习1.15 当角度（用弧度描述） x 足够小的时候，其正弦值可以用sin x ≈ x 计算，而三角恒等式

可以减小sin的参数的大小（为完成这一练习，如果一个角不大于0.1弧度，我们就认为它“足够小”）。这些想法都体现在下面的函数里：

a.在求值sine（12.15）时p将被调用多少次？

b.在求值sine（a）时，由函数sine产生的计算过程使用的空间和步数（作为 a 的函数）增长的阶是什么？

1.2.4 求幂

现在考虑对给定的数计算乘幂的问题，我们希望函数的参数是一个基数 b 和一个正整数的指数 n ，函数计算 b ⁿ 。做这件事的一种方法是根据下面的递归定义：

b ⁿ = b · b ^{n -1}

b ⁰ =1

它可以直接翻译为如下的JavaScript函数：

这个函数产生一个线性递归计算过程，需要 Θ （ n ）步和 Θ （ n ）空间。像阶乘一样，我们很容易将其形式化为一个等价的线性迭代：

这一版本需要 Θ （ n ）步和 Θ （1）空间。

我们可以通过反复求平方，以更少的步骤完成乘幂的计算。例如，不是采用下面这样的方式计算 b ⁸ ：

b ·（ b ·（ b ·（ b ·（ b ·（ b ·（ b · b ））））））

而是通过三次乘法算出它来：

b ² = b · b

b ⁴ = b ² · b ²

b ⁸ = b ⁴ · b ⁴

这种方法对指数为2的乘幂都适用。如果再利用下面规则，我们就可以借助连续求平方的方法，完成一般的乘幂计算了：

b ⁿ =（ b ^{n /2} ） ² 如果 n 是偶数

b ⁿ = b · b ^{n -1} 如果 n 是奇数

我们可以用如下的函数表述上面的方法：

其中检测一个整数是否偶数的谓词可以基于计算整数除法余数的运算符%定义：

由fast_expt产生的计算过程，在空间和步数上相对于 n 都是对数的。要看到这个情况，请注意，在用fast_expt计算 b ^{2 n} 时，只需要比计算 b ⁿ 多做一次乘法。每做一次新乘法，能使计算的指数值（大约）增大一倍。这样，计算指数 n 所需要的乘法次数的增长大约就是以2为底的 n 的对数值，也就是说，这一计算过程增长的阶为 Θ （log n ） ^[5] 。

随着 n 的变大， Θ （log n ）增长与 Θ （ n ）增长之间的差异也会变得非常明显。例如，对于 n =1000，fast_expt只需要做14次乘法 。我们也可能采用连续求平方的想法，设计一个具有对数步数的计算乘幂的迭代算法（见练习1.16）。但是，就像迭代算法的常见情况一样，写出这一算法就不像递归算法那样直截了当了 ^[6] 。

练习1.16 请设计一个函数，它使用一系列的求平方，产生一个迭代的求幂计算过程，但是就像fast_expt那样只需要对数的步数。（提示：请利用关系（ b ^{n /2} ） ² =（ b ² ） ^{n /2} ，除了指数 n 和基数 b 之外，还应该维持一个附加的状态变量 a ，并定义好状态变换，使得从一个状态转到另一状态时乘积 a · b ⁿ 不变。在计算过程开始时令 a 取值1，并用计算过程结束时 a 的值作为回答。一般而言，定义一个 不变量 ，要求它在状态之间保持不变，这一技术是思考迭代算法设计问题时的一种非常强有力的方法。）

练习1.17 本节中几个求幂算法的基础都是通过反复做乘法去求乘幂。与此类似，也可以通过反复做加法的方式求出乘积。下面的乘积函数与expt函数类似（在这里假定我们的语言里只有加法而没有乘法）：

这一算法所需的步骤数相对于b是线性的。现在假定除了加法外还有一个函数double，它求出一个整数的两倍；还有函数halve，它把一个（偶）数除以2。请用这些运算设计一个类似fast_expt的求乘积函数，使之只用对数的计算步数。

练习1.18 利用练习1.16和练习1.17的结果设计一个函数，它能产生一个基于加、加倍和折半运算的迭代计算过程，只用对数的步数就能求出两个整数的乘积 。

练习1.19 存在只需要对数步就能求出斐波那契数的巧妙算法。请回忆1.2.2节fib_iter产生的计算过程中状态变量 a 和 b 的变换规则 a ← a + b 和 b ← a ，我们把这种变换称为 T 变换。可以看到，从1和0开始把 T 反复应用 n 次，将产生一对数Fib（ n +1）和Fib（ n ）。换句话说，斐波那契数可以通过将 T ⁿ （变换 T 的 n 次方）应用于对偶（1,0）而得到。现在把 T 看作变换族 T _pq 中 p =0且 q =1的特殊情况，其中 T _pq 是对偶（ a , b ）按 a ← bq + aq + ap 和 b ← bp + aq 规则的变换。请证明，如果应用变换 T _pq 两次，其效果等同于应用同样形式的变换 T _p′q′ 一次，其中的 p′ 和 q′ 可以由 p 和 q 算出来。这就指明了一种计算这种变换的平方的路径，使我们可以通过连续求平方的方法计算 T ⁿ ，就像fast_expt函数里所做的那样。把所有这些放到一起，就得到了下面的函数，其运行只需要对数的步数 ：

1.2.5 最大公约数

两个整数 a 和 b 的最大公约数（GCD），就是能除尽这两个数的最大整数。例如，16和28的GCD是4。在第2章，当我们研究有理数算术的实现时，就需要计算GCD，以便把有理数约化到最简形式（要把有理数约化到最简形式，我们必须将其分母和分子同时除掉它们的GCD。例如，16/28将约简为4/7）。要得到两个整数的GCD，一种方法是对它们做因数分解，然后找出公共因子。然而，存在一个更高效的著名算法。

该算法的思想基于下面的观察：如果 r 是 a 除以 b 的余数，那么 a 和 b 的公约数正好也是 b 和 r 的公约数。因此我们可以使用等式：

GCD（ a , b ）=GCD（ b , r ）

这样就把一个GCD计算问题连续归约到越来越小的整数对的GCD计算问题。例如：

GCD（206, 40）=GCD（40, 6）

=GCD（6, 4）

=GCD（4, 2）

=GCD（2, 0）

=2

这里从GCD（206, 40）归约到GCD（2, 0），最终得到了2。可以证明，从任意两个正整数出发，反复执行这种归约，最终将产生一个整数对，其中第二个数是0，此时的GCD就是第一个数。这个计算GCD的方法称为 欧几里得算法 ^[7] 。

不难把欧几里得算法描述成一个函数：

这个函数将产生一个迭代计算过程，其步数按照所涉及的数的对数增长。

欧几里得算法所需步数是对数增长的，这一事实与斐波那契数有一种有趣的关系：

Lamé定理 ： 如果欧几里得算法计算出一对整数的GCD需要用 k 步，那么这对数中较小的那个数必然大于或等于第 k 个斐波那契数 ^[8] 。

我们可以利用这一定理，给出欧几里得算法的增长阶估计。令 n 作为函数输入的两个数中较小的那个，如果计算过程需要 k 步，那么我们就一定有 。这样，步数 k 的增长就是 n 的对数（对数的底是 ϕ ）。这样，算法增长的阶就是 Θ （log n ）。

练习1.20 一个函数产生的计算过程当然依赖解释器使用的规则。作为例子，请考虑上面给出的迭代式gcd函数，假定解释器采用1.1.5节介绍的正则序解释这个函数（对条件表达式的正则序求值规则在练习1.5中说明）。请采用（正则序的）代换方法展示求值表达式gcd（206, 40）时产生的计算过程，并标出实际执行的remainder运算。采用正则序求值gcd（206, 40），需要执行多少次remainder运算？如果采用应用序求值呢？

1.2.6 实例：素数检测

本节讨论检查整数 n 是否素数的两种方法，第一个具有 的增长阶，另一个“概率”算法具有 Θ （log n ）的增长阶。本节最后的练习中有几个基于这些算法的编程作业。

寻找因子

从古代开始，数学家就一直被有关素数的问题所吸引，许多人研究过确定整数是否素数的方法。检测一个数是否素数，一种方法就是找出它的因子。下面的程序能找出给定整数 n 的（大于1的）最小整数因子。这里采用一种直截了当的方法，从2开始一个个地检查整数，看它们能否整除 n 。

我们可以用如下方式检查一个数是否素数： n 是素数当且仅当它是自己的最小因子。

find_divisor的结束判断基于如下事实：如果 n 不是素数，它必然有一个小于等于 的因子 ^[9] 。这也意味着上面的算法中只需要在1和 之间检查因子。由此可知，这样确定是否素数需要的步数具有 的增长阶。

费马检查

Θ （log n ）的素数检查基于数论里著名的费马小定理的结果 。

费马小定理 ： 如果 n 是素数， a 是小于 n 的任意正整数，那么 a 的 n 次方与 a 模 n 同余。

（我们说两个数 模 n 同余 ，如果它们除以 n 的余数相同。数 a 除以 n 的余数也称为 a 取模 n 的余数 ，或者简称为 a 取模 n 。）

如果 n 不是素数，那么，一般而言，大部分的 a < n 都不满足上面的关系。由此就得到了下面的素数检查算法：对给定的整数 n ，随机取一个 a < n 并算出 a ⁿ 取模 n 的余数。如果得到的结果不等于 a ，那么 n 肯定不是素数。如果结果是 a ， n 就有很大的机会是素数。现在再随机取一个 a 并做同样检查。如果结果满足上面的关系，那么我们对 n 是素数就有了更大的信心。通过检查越来越多的 a 值，我们就可以不断增加对有关结果的信心。这一算法称为费马检查。

为了实现费马检查，我们需要一个函数来计算一个数的幂对另一个数取模的结果：

这个函数很像1.2.4节的函数fast_expt，其中采用了连续求平方的方法，因此，相对于计算中的指数，其执行步数具有对数的增长阶 ^[10] 。

执行费马检查时需要随机选取位于1和 n -1之间（包含这两者）的数 a ，然后检查 a 的 n 次幂取模 n 的余数是否等于 a 。随机数 a 的选取用基本函数math_random完成，该函数返回一个小于1的非负实数。这样，要得到一个1和 n -1之间的随机整数，我们把math_random的返回值乘以 n -1，再用基本函数math_floor四舍五入，最后再加一：

下面这个函数检查给定的数，它按第二个参数给定的次数运行上述检查。如果每次检查都成功，这一函数的值就是真，否则就是假：

概率方法

从特性上看，费马检查与我们前面熟悉的算法都不一样。前面那些算法都保证计算结果一定正确，而这里得到的结果则只是可能正确。说得更准确些，如果数 n 不能通过费马检查，我们可以确信它不是素数。而 n 通过了检查的事实只是一个很强的证据，仍然不能保证 n 为素数。我们只能说，对任何正整数 n ，如果执行这一检查的次数足够多，而且看到 n 通过了所有检查，就能把这一素数检查出错的概率减小到所需要的任意程度。

不幸的是，这一断言并不完全正确。因为确实存在一些能骗过费马检查的整数。存在一些数 n ，它们不是素数但却具有如下性质：对任意的整数 a < n ，都有 a ⁿ 与 a 模 n 同余。由于这种数极其罕见，因此费马检查在实践中还是很可靠的 。存在一些不会受骗的费马检查变形，它们也像费马方法一样，在检查整数 n 是否素数时，随机选择整数 a < n 并去检查某些依赖于 n 和 a 的关系（练习1.28是这类检查的一个例子）。另一方面，与费马检查不同，可以证明，对任意数 n ，相应条件对整数 a < n 中的大部分都不成立，除非 n 是素数。这样，如果对某个随机选出的 a ， n 能通过检查， n 是素数的机会就大于一半。如果 n 对两个随机选择的 a 能通过检查， n 是素数的机会就大于四分之三。通过用更多随机选择的 a 值运行这种检查，我们可以使出错的概率减小到所需要的任意程度。

存在这样的检查，能证明其出错的机会可以达到任意小，这一情况也引起了人们对这类算法的极大兴趣，目前已形成了称为 概率算法 的领域。这一领域已经有了大量研究工作，概率算法也被成功地应用于许多重要领域 。

练习1.21 请用smallest_divisor函数找出下面各数的最小因子：199，1 999，19 999。

练习1.22 假设有一个无参数的基本函数get_time，它返回一个整数，表示从1970年1月1日的00:00:00 起到现在已经过去的微秒数。如果对整数 n 调用下面的timed_prime_test函数，它将打印出 n ，然后检查 n 是否素数。如果 n 是素数，函数将打印出三个星号 ，随后是执行这一检查所用的时间量。

请利用这个函数写一个search_for_primes函数，它检查给定范围内连续的各个奇数的素性。请用你的函数找出大于1 000、大于10 000、大于100 000和大于1 000 000的最小的三个素数。请注意检查每个素数所需的时间。因为这一检查算法具有 的增长阶，你可以期望在10 000附近的素数检查耗时大约为在1 000附近的素数检查的 倍。你得到的数据确实如此吗？由100 000和1 000 000得到的数据，对这一 预测的支持情况如何？概念上说，在你的机器上运行的时间正比于计算所需的步数，你的结果符合这一说法吗？

练习1.23 本节开始时给出的smallest_divisor函数做了许多无用检查：在检查了一个数能否被2整除后，完全没必要再检查它是否能被任何偶数整除。这说明test_divisor用的值不该是2, 3, 4, 5, 6, …，而应该是2, 3, 5, 7, 9, …。请实现这种修改。其中声明一个函数next，用2调用时它返回3，否则返回其输入值加2。修改smallest_divisor函数，让它用next（test_divisor）而不是test_divisor+1。请用结合了这一smallest_divisor版本的timed_prime_test运行练习1.22里那个找12个素数的测试。这样修改使检查的次数减半，你可能期望其运行速度快一倍。实际情况符合这一预期吗？如果不符，你观察到的两个算法速度的比值是什么？你如何解释比值不是2的事实？

练习1.24 请修改练习1.22的timed_prime_test函数，让它使用fast_is_prime（费马方法），并检查你在该练习中找出的12个素数。因为费马检查具有 Θ （log n ）的增长速度，对于检查接近1 000 000的素数与检查接近1000的素数，你预期两个时间之间的比较应该怎样？你的数据符合这一预期吗？你能解释所发现的任何不符吗？

练习1.25 Alyssa P.Hacker提出，expmod做了过多的额外工作。她说，毕竟我们已经知道怎样计算乘幂，因此只需要简单地写：

她说的对吗？这个函数能很好地服务于我们的快速素数检查程序吗？请解释这些问题。

练习1.26 Louis Reasoner在做练习1.24时遇到很大困难，他的fast_is_prime检查看起来运行得比他写的is_prime还慢。Louis请他的朋友Eva Lu Ator过来帮忙。在检查Louis的代码时，两人发现他重写了expmod函数，其中采用显式的乘法而没有调用square：

“我看不出来这会造成什么不同，”Louis说。“我能看出，”Eva说，“用这种方式写这个函数，你就把一个 Θ （log n ）的计算过程变成 Θ （ n ）的了。”请解释这个问题。

练习1.27 请证明脚注47中列出的那些Carmichael数确实能骗倒费马检查。也就是说，写一个函数，它以整数 n 为参数，对每个 a < n 检查 a ⁿ 是否与 a 模 n 同余。然后用你的函数去检查前面给出的那些Carmichael数。

练习1.28 费马检查的一种不会被骗的变形称为Miller-Rabin 检查 （Miller 1976, Rabin 1980），它源于费马小定理的一个变形。该变形断言，如果 n 是素数， a 是任何小于 n 的正整数，则 a 的（ n -1）次幂与1模 n 同余。用Miller-Rabin检查考察数 n 的素数性时，我们随机选取数 a < n 并用函数expmod求 a 的（ n -1）次幂对 n 的模。然而，在执行expmod中的平方步骤时，我们要查看是否发现了一个“1取模 n 的非平凡平方根”，也就是说，是否存在不等于1或 n -1的数，其平方取模 n 等于1。可以证明，如果1的这种非平凡平方根存在，那么 n 就不是素数。还可以证明，如果 n 是非素数的奇数，那么，至少存在一半的 a < n ，按这种方式计算 a ^{n -1} ，会遇到1取模 n 的非平凡平方根。这也是Miller-Rabin检查不会受骗的原因。请修改expmod函数，让它在发现1的非平凡平方根时报告失败。利用它实现一个类似于fermat_test的函数完成Miller-Rabin检查。通过检查一些已知的素数和非素数的方式考验你的函数。提示：让expmod送出失败信号的一种方便方法是让它返回0。 PvUl7KJ4Z8ydkXn4MPFZnNNPqJOoCpiiJb0jTNCsXxqQrRSsBVs6C3GfNGqyLhq+