差分机

差分机是一种由 Charles Babbage 设计的，专门用于计算多项式的机器。

差分机一号设计于 19 世纪 20 年代，而更先进的差分机二号设计于 1846 年到 1849 年。但是，由于种种原因，它们并没有被成功地造出来。直到 1991 年，伦敦科学博物馆才根据 Babbage 的原始图纸造出差分机二号。

机械结构运算入门 🔗

差分机是怎么用机械结构计算多项式的？

含有单一未知数 $x$ 的多项式可以被表示为以下形式：

$$ a_nx^n+a_{n-1}x^{n-1}+\dotsb+a_2x^2+a_1x+a_0 $$

要计算多项式的值，涉及到加法和乘法。

首先要解决的问题是，怎么在机械结构里表示数字？

假设有一个 10 齿齿轮，我们可以沿顺时针方向，为每个齿分别标上数字 0~9。

上面两个齿轮正对着我的齿分别是 3 和 4，就可以认为它们分别表示 3 和 4。

如果齿轮逆时针转动一格，表示的数字就加 1，顺时针转动一格，数字就减 1。把两个齿轮拼在一起，顺时针转动右侧齿轮，当右侧齿轮的数字变为 0 时，就完成了加法。

可以看到，当右侧齿轮的数字变为 0 时，左侧齿轮的数字变成了 3+4=7。

利用差分代替乘法 🔗

乘法比加法复杂得多。Babbage 的思路是利用差分，化乘法为加法，这也是差分机名字的由来。

差分是离散化的微分。这里我们研究的是，对函数 $f(x)$，当 $\Delta x = 1$ 时的 $\Delta f(x)$。也就是说，$x$ 每增加 $1$，$f(x)$ 增加多少。

考虑多项式 $p(x) = 2x^2-3x+2$，其差分为

$$\begin{align} \Delta p(x) &= p(x+\Delta x) - p(x) \\ &= p(x+1) - p(x) \end{align} $$

可以发现

$$ \Delta p(0) = -1, \Delta p(1) = 3, \Delta p(2) = 7, \dotsb $$

如果我们对差分求差分

$$ \begin{align} \Delta(\Delta p(x)) &= \Delta p(x+1) - \Delta p(x) \end{align} $$

为了方便表示，我们将 $p(x)$ 的一阶差分记作 $\Delta_1 p(x)$，二阶差分记作 $\Delta_2 p(x)$。

可以发现

$$ \Delta_2 p(0) = 4, \Delta_2 p(1) = 4, \dotsb $$

$p(x)$ 的二阶差分 $\Delta_2p(x)$ 是常数 $4$。

这不难理解，如果我们对 $p(x)$ 连续求导两次，也能得到常数 4。

将上面的结果整理成表格

表格的最后一列为常数 4，从第 2 行开始，每一格的数字都是上一行的两数之和。比如，$ {\color{blue}7} = {\color{red}3} + {\color{red}4}$。

这个特性让我们只需要用加法，就可以从任意一行开始，不断地迭代出下一行。

因为每个数只与上方和右上方数有关，所以我们可以覆盖式地求解。只需要 3 个计数器（寄存器），分别存储 $p(x)$、$\Delta_1 p(x)$ 和 $\Delta_2 p(x)$，就可以不断地迭代出下一行。

同理，求解 $n$ 次多项式需要使用 $n+1$ 个计数器，分别存储 $p(x), \Delta_1 p(x), \Delta_2 p(x), \dotsb, \Delta_n p(x) $。

十进制加法器 🔗

如果把 $p(x)$ 记为 $\Delta_0 p(x)$，每迭代一个 $x$ 的过程需要做 $n$ 次加法： $$ \begin{align} \Delta_0 p(x) &\leftarrow \Delta_0 p(x) + \Delta_1 p(x) \\ \Delta_1 p(x) &\leftarrow \Delta_1 p(x) + \Delta_2 p(x) \\ & \dotsb \\ \Delta_{n-1} p(x) &\leftarrow \Delta_{n-1} p(x) + \Delta_{n} p(x) \end{align} $$ 其中，$ A \leftarrow A + B $ 表示计算 $A+B$ 的结果，放入 $A$。

这样的加法器无法满足需求，虽然左侧齿轮上的数字变成了两数之和 7，但是右侧齿轮上的数字变成了 0。我们期望加法完成后，右侧齿轮能保持原来的数字 4。

为此，我们需要在两个齿轮中间加上一个齿轮。这个齿轮通过升降来调整与两个齿轮是否相连。

1. 开始时，中间的齿轮与左右两个齿轮都相连，右侧齿轮带动左侧齿轮，将左侧齿轮的数字变为两数之和 7；
2. 然后，中间的齿轮上升，断开与左侧齿轮的连接，但保持与右侧齿轮相连；
3. 最后，中间的齿轮带动右侧齿轮变回原数字 4。

这就是一位十进制加法器的大致工作原理。如果对实现细节感兴趣，比如如何存储多位十进制数、如何进位，可以搜索差分机的相关视频。

性能优化 🔗

根据上面的原理，同一时间最多只能有一个加法器在工作。 $$ \begin{align} \Delta_0 p(x) &\leftarrow \Delta_0 p(x) + {\color{blue}\Delta_1 p(x)} \\ {\color{blue}\Delta_1 p(x)} &\leftarrow \Delta_1 p(x) + \Delta_2 p(x) \end{align} $$ 如果这两个加法同时进行，前者需要读 $\Delta_1 p(x)$，后者需要写 $\Delta_1 p(x)$，显然是有冲突的。

差分机二号使用了一种技巧，同时计算多轮迭代中的部分加法。考虑多项式 $q(x) = x^3+1$，可以这样求值：

把 $n$ 个加法操作一个隔一个，分成两组 $$ \begin{align} \Delta_0 p(x) &\leftarrow \Delta_0 p(x) + \Delta_1 p(x) \\ \Delta_2 p(x) &\leftarrow \Delta_2 p(x) + \Delta_3 p(x) \\ \Delta_4 p(x) &\leftarrow \Delta_4 p(x) + \Delta_5 p(x) \\ & \dotsb \end{align} $$ 和 $$ \begin{align} \Delta_1 p(x) &\leftarrow \Delta_1 p(x) + \Delta_2 p(x) \\ \Delta_3 p(x) &\leftarrow \Delta_3 p(x) + \Delta_4 p(x) \\ \Delta_5 p(x) &\leftarrow \Delta_5 p(x) + \Delta_6 p(x) \\ & \dotsb \end{align} $$

因为同一组中的加法不涉及读写同一个计数器，所以它们可以并行。

经过优化，差分机二号计算多项式的时间与多项式的次数无关，只跟要迭代的次数有关。

关于并行性的讨论 🔗

差分机二号的优化可以给我们一些利用并行优化算法的思路。类似的优化对现代 CPU 而言也是有效的。

考虑对数组 nums 求和：

int nums[1<<20];
int num_of_nums = sizeof nums / sizeof nums[0];

最简单的办法是使用循环，把每个数字加到累加器上：

int sum() {
  int s = 0;
  for (int i = 0; i < num_of_nums; i++) {
    s += nums[i];
  }
  return s;
}

手动进行循环展开理论上可以减少分支带来的开销：原先每次加法之前都需要检查循环条件是否满足，现在只需要每 4 次加法检查一次。

int unrolled_sum() {
    int s = 0;
    int i = 0;
    for (i = 0; i + 3 < num_of_nums; i += 4) {
      s += nums[i];
      s += nums[i + 1];
      s += nums[i + 2];
      s += nums[i + 3];
    }
    /* 别忘了其余元素 */
    for (; i < num_of_nums; i++) {
      s += nums[i];
    }
    return s;
}

但是实际上 CPU 的分支预测器会进行预测，如果条件大概率是满足的，就会预执行条件满足时的指令。

上面的代码也没有提高并行性，因为循环中的 4 条语句都需要读写同一个变量 s。

要进一步榨干 CPU 的性能，就需要去除对公共变量 s 的依赖，充分利用 CPU 的流水线。

一个简单的思路是让循环中的 4 条语句单独求和，最终再汇总起来：

int parallel_sum() {
    int s[4] = {0};
    int i;
    for (i = 0; i + 3 < num_of_nums; i += 4) {
      s[0] += nums[i];
      s[1] += nums[i + 1];
      s[2] += nums[i + 2];
      s[3] += nums[i + 3];
    }
    for (; i < num_of_nums; i++) {
      s[0] += nums[i];
    }
    return s[0] + s[1] + s[2] + s[3];
}

此时循环中的 4 条指令可以并行，而这个优化本质上跟差分机二号的优化一样，都是通过去除公共依赖，利用并行提升效率。

以下是使用 Clang 11.0，开启 O1 级别优化时的 benchmark 结果（越小越快）：

测试结果符合预期，手动循环展开的版本 unrolled_sum 相比 sum 快一些，而循环内可并行的版本 parallel_sum 则更快，效率是 sum 的 2.1 倍。

是不是突然觉得自己掌握了代码优化的精髓，想把以前的代码一顿优化？

别急着动手，我们把优化级别调整到 O2（GCC 11.1 需要开启 O3），看下现代编译器的功力：

什么情况？sum 反而是最快的，效率逼近其余两个版本的 5 倍。直接看产生的汇编，可以看出端倪。

paddd 指令是一条 SIMD（Single Instruction Multiple Data） 指令，一次可以操作 4 个 32 位整数（16 字节，也就是图中 nums+16、nums+32……的原因），简单理解就是向量加法。

xmm0 和 xmm1 是两个 128 位寄存器，类似 parallel_sum 中的 int s[4]。产生的汇编相当于把数组分成了 8 个部分并行求和，还顺便帮我们做了循环展开。

不能说是快，只能说是飞一般的速度。

得到 8 个和 $S_0, S_1, \dotsb, S_7$ 之后，xmm0 和 xmm1 寄存器是这样的：

$$ \begin{align}\text{xmm0} &= (S_0, S_1, S_2, S_3) \\ \text{xmm1} &= (S_4, S_5, S_6, S_7)\end{align} $$

在计算 $\sum_{i=0}^7{S_i}$ 时，生成的汇编代码用到了归并的思路。

paddd   xmm0, xmm1

先把 xmm1 加到 xmm0 上，此时

$$ \begin{align}\text{xmm0} = (&S_0 + {\color{blue}S_4}, \\ &S_1 + {\color{blue}S_5}, \\ &S_2 + {\color{blue}S_6}, \\ &S_3 + {\color{blue}S_7})\end{align} $$

然后目标就变成了求 xmm0 四个元素之和。

pshufd  xmm1, xmm0, 78   # xmm1 = xmm0[2,3,0,1]

再将 xxm0 的前半部分赋值给 xmm1 的后半部分，后半部分赋值给 xmm1 的前半部分，使前后交错。

$$ \begin{align}\text{xmm0} = (&{\color{blue}S_0 + S_4}, \\ &{\color{blue}S_1 + S_5}, \\ &{\color{red}S_2 + S_6}, \\ &{\color{red}S_3 + S_7}) \\ \text{xmm1} = (&{\color{red}S_2 + S_6}, \\ &{\color{red}S_3 + S_7}, \\ &{\color{blue}S_0 + S_4}, \\ &{\color{blue}S_1 + S_5})\end{align} $$

paddd   xmm1, xmm0

再次求和，此时 xmm1 的前半部分和即为结果：

$$ \begin{align}\text{xmm1} = (& {\color{blue}S_0} + S_2 + {\color{blue}S_4} + S_6, \\ & {\color{blue}S_1} + S_3 + {\color{blue}S_5} + S_7, \\ & \dotsb)\end{align} $$

pshufd  xmm0, xmm1, 229  # xmm0 = xmm1[1,1,2,3]

使用同样的套路，把 xmm1 的第二项赋值给 xmm0 的第一项：

$$ \begin{align}\text{xmm0} = (& S_1 + S_3 + S_5 + S_7, \\ & \ldots ) \\ \text{xmm1} = (& S_0 + S_2 + S_4 + S_6, \\ & \ldots ) \end{align} $$

paddd   xmm0, xmm1
movd    eax, xmm0

最后相加得到结果，放到函数返回值使用的寄存器里去：

$$ \begin{align}\text{xmm0} = (& {\color{blue}\sum_{i=0}^{7}{S_i}}, \ldots ) \end{align} $$

可以看出，现代编译器已经能够帮我们做很多优化了，甚至很多时候比我们自己优化的结果要好得多。

编译器通过模式匹配和重写来做优化，我们手动进行循环展开，或者企图榨干 CPU 性能的代码，对编译器而言反而不是常见的模式，最终都没有优化到与 sum 相同的性能水平。这也是「不要过早优化」的一个原因。如果基于某些场景做了优化，一定要做 benchmark。

两百年的时间，把齿轮变成了晶体管，把十进制变成了二进制，把多项式计算变成了通用计算。但人们对计算力的追求却像一台永无止境的差分机，构建了属于这个时代的摩尔定律。