跳至内容
Topos Uranos
中文 (中国)

数值算法

摘要

本单元前面的几课教了如何编写循环，以及如何借助不变式和变式对循环进行推理；本课把这种推理应用于四个经典的数值算法，它们比计算机还古老，至今仍在使用。用试除法做素性检验：它在 n\sqrt{n} 处停止，其理由可以证明；如果写得不小心，它的停止条件会溢出。求最大公因数的辗转相除法：它的不变式、终止性、负数带来的陷阱，以及它与 std::gcd 的比较。用逐次除法进行的进制转换，以及秦九韶算法（霍纳法则）：它由各位数字重建一个数，并用同一格式计算多项式的值。还有二分法：它把一个区间一分为二来逼近一个根，由介值定理保证，所需轮数可以事先算出。每个算法都作为数学对象来介绍：证明它做到了它所说的，证明它会终止，并初步数一数它的步数。

学习目标

  1. 实现试除到 n\sqrt{n} 为止的素性检验，并说明为什么到 n\sqrt{n} 为止就够了。
  2. 实现辗转相除法，陈述它的不变式（最大公因数不变），并把它与 std::gcd 比较。
  3. 用逐次除法把一个数转换为另一种进制，并用秦九韶算法计算一个位置表示或一个多项式的值。
  4. 用二分法把一个根逼近到给定的容差，并计算所需的轮数。

因数与素数

试除法

整数 n≥2n \geq 2 若其正因数只有 1 和 nn，就称为素数；若还有别的因数，就称为合数；0 和 1 既不是素数也不是合数。这个定义立刻提示了一个算法：逐个试验每个候选数 dd（从 2 到 n−1n - 1），一旦有一个能整除，就宣布这个数是合数。这正是循环一课中求最小因数的循环，它是正确的。然而，对于一个十位数的素数，它要做将近一百亿次除法，而且几乎全都是无用的，原因如下。

命题只需试到 n\sqrt{n}

若 n≥2n \geq 2 是合数，则它有一个因数 dd，满足 2≤d≤n2 \leq d \leq \sqrt{n}。

证明

  1. n=a⋅b,2≤a≤bn = a \cdot b, \quad 2 \leq a \leq b

    若 nn 是合数，它可以写成两个因子的乘积，其中没有一个等于 1 或 nn。把两者中较小的那个记为 aa。

  2. a2≤a⋅b=na^2 \leq a \cdot b = n

    把 a≤ba \leq b 两边乘以正数 aa，得到 a2≤ab=na^2 \leq ab = n，由此 a≤na \leq \sqrt{n}：任何因式分解中较小的因子都位于 n\sqrt{n} 的左边。

  3. d≤n  ⟺  nd≥nd \leq \sqrt{n} \iff \frac{n}{d} \geq \sqrt{n}

    因此，因数成对出现，即 (d,n/d)(d, n/d)，关于 n\sqrt{n} 对称：当 n=36n = 36 时，这些数对是 (2,18)(2, 18)、(3,12)(3, 12)、(4,9)(4, 9) 和 (6,6)(6, 6)。如果到 n\sqrt{n} 为止没有一个候选数整除 nn，那么超过它的候选数也没有一个能整除，因为与之配对的数本应更早出现。

由此可知，除法次数从大约 nn 次降到大约 n\sqrt{n} 次：对于接近二十亿的 nn，从二十亿次降到四万六千次。程序使用设计正确的循环一课中的标志模式：isPrime 一开始为真，如果循环找到一个因数，就把它清除。

C++ 程序
#include <iostream>#include <print>​int main(){    int n{};    std::cin >> n;    bool isPrime{n >= 2};    for (int d{2}; isPrime && d <= n / d; ++d) {        if (n % d == 0) {            isPrime = false;        }    }    if (isPrime) {        std::println("{} 是素数", n);    } else {        std::println("{} 不是素数", n);    }}
程序输入
2147483647
输出
2147483647 是素数

有三个决定值得评说。标志用 n >= 2 初始化，因此 0、1 和负数不进入循环就被判定为非素数：忘记这种情形是典型错误，因为当 n=1n = 1 时循环一轮也不执行，程序就会把 1 判定为素数。for 的条件有两部分：标志，它让循环在第一个因数处停下；以及界限。而界限写成 d <= n / d，对于正整数，它等价于 d2≤nd^2 \leq n（整数除法会截断，但 dd 不超过截断后的商，恰好当它不超过 n/dn/d 时成立，因为 dd 是整数）。输入的数 231−12^{31} - 1 是最大的 int，而且是素数；程序用 46 33946\,339 次求余来判定这一点，即对从 2 到 46 34046\,340 的各候选数求余。

循环的不变式是：只要标志仍为真，从 2 到 d−1d - 1 的整数都不整除 nn。进入时它成立，此时从 2 到 1 的区间是空的；它得以保持，因为如果 dd 不整除 nn，区间就只增加了一个不整除的数；而带着为真的标志退出时，不变式与为假的条件 d2>nd^2 > n 合起来说明：到 n\sqrt{n} 为止没有一个候选数整除 nn，而这正是命题所要求的。变式是还剩下待试的候选数的个数，也就是介于 dd 与 n\sqrt{n} 之间的整数个数：它每轮减少一，并且不会低于零。

会溢出的条件

写界限的自然方式是 d * d <= n，而这是错的。对同一个输入，循环到达 d=46 341d = 46\,341，这时乘积 46 3412=2 147 488 28146\,341^2 = 2\,147\,488\,281 超过了 int 的最大值。

C++ 程序
#include <iostream>#include <print>​int main(){    int n{};    std::cin >> n;    bool isPrime{n >= 2};    for (int d{2}; isPrime && d * d <= n; ++d) {        if (n % d == 0) {            isPrime = false;        }    }    std::println("素数？{}", isPrime);}
程序输入
2147483647
输出
main.cpp:9:33: runtime error: signed integer overflow: 46341 * 46341 cannot be represented in type 'int'

程序以代码 1 结束。

未定义行为检测器在乘积处停止程序。没有它，有符号溢出就是未定义行为，标准并不规定会发生什么。在 Linux 上用 GCC 观察到的情况，无论是否开启优化，都很有启发：乘积发生回绕，变成负数或小于 nn 的数，而由于没有任何 int 超过最大值，条件 d * d <= n 就再也不会不成立；循环继续试验候选数，大约二十亿轮，直到到达 d=nd = n（它整除 nn），于是程序把一个素数判定为合数。溢出错误并不表现为溢出，而是表现为错误的答案。d <= n / d 这种形式不可能溢出，因为它做的是除法而不是乘法；位与数学库一课给出过另一条出路，即用 std::sqrt 一次性算出界限，但那样界限的精确性就取决于浮点数的舍入，而整数除法没有这个问题。

代数与几何 I 课程中的“素数与算术基本定理”一课证明了，每个大于 1 的整数都能以唯一的方式写成素数的乘积；试除法是求出这个乘积最直接的方式：用找到的每个因数去除 nn，再对商继续进行。