第08章 数学函数实现：除法、泰勒展开与 CORDIC

这一篇在干嘛？

FPGA 里没有”除法器”和”sin/cos”这种现成的好用资源——加、减、乘都有高效电路，唯独除法和三角函数又慢又费面积。本章教你三种硬件除法方案（乘移法、迭代除法、Goldschmidt）和两种超越函数方案（泰勒级数、CORDIC），并告诉你每种方案适合什么场景。

一个有趣的事实：绝大多数现实中的数学问题，都可以用”移位 + 加法”的组合来解决。本章围绕除法和三角函数展开，最后再推广到更一般的函数类。同一个问题往往有多种解法，需要针对应用做优化。

8.1 硬件除法：为什么这么头疼

除法一直是数字逻辑设计者的老大难。与加法、减法、乘法不同，没有一个简单的逻辑运算能直接产生商。更微妙的是：定点数的除法不保证产生”有限且可预测”的定点结果（比如 1÷3 是无限循环小数）。解决方案有好几种，简单程度完全取决于应用约束——我们从最简单、约束最多的方案讲起，逐步过渡到通用方案。

8.1.1 乘法 + 移位法（Multiply and Shift）

这是最简单的除法方案，本质是乘以除数的倒数。它利用了二进制数的基本性质：右移一位等于除以 2。

先看定点表示的背景。假设一个 8 位寄存器，小数点固定在第 4 位（即低 4 位是小数部分），里面存着整数 3：

这个寄存器的值等于 2¹ + 2⁰ = 3。如果除以 2，只需右移 1 位：

移位后寄存器的值变为 2⁰ + 2⁻¹ = 1.5。

乘移法的思路是：假定除数是一个固定常数，把”除以 D”变成”乘以某个数，再除以 2 的幂”。书中给的例子是除以 7：

• 1/7 ≈ 73/512（73/512 = 0.14257…，而 1/7 ≈ 0.14285…）
• 所以 x ÷ 7 ≈ (x × 73) >> 9（右移 9 位即除以 512）
• 实际效果是除以 7.013…，有一定误差

想要更高精度？增大 2 的幂，并相应增大乘法因子即可。比如用 x × 585 >> 12（585/4096 ≈ 0.14282）误差就更小，但乘法器也更大。

点评一下这个方法的适用条件：

• 优点：除数固定时速度极快，只有一个乘法器加一次移位，非常适合高速流水线。
• 限制：只对”特定形式”的除数有效——本质上是要先把除数求倒数、转换成一个定点常数。这个求倒数通常由外部设备（如微处理器）预先算好。如果除数是设计中其他逻辑动态产生的、没法预先取倒数，就得用下面的通用方法。

8.1.2 迭代除法（Iterative Division）

这一节的算法属于**逐位循环法（digit recurrence）**一族。它的原理和小学学的十进制长除法一模一样，只不过用二进制还有额外的优化空间。看一个二进制长除法的例子——2 除以 3：

定点除法器的硬件架构如图 8.4，核心是一个比较器 + 减法器：

工作流程是这样的：

1. 预规格化（prenormalize）：把被除数的小数点移到某个位置，使被除数严格小于除数的 2 倍。这样做的好处是：之后每一步移位产生的部分被除数，必然小于除数的 2 倍。
2. 于是每一轮迭代中，“除数进入部分被除数”的次数只可能是 1 次或 0 次——如果除数 ≤ 当前部分被除数，商寄存器移入逻辑 1 并做减法；否则移入逻辑 0。
3. 部分被除数每轮左移 1 位，重复迭代直到达到所需精度。
4. 后规格化（postnormalize）：把结果移回正确的定点位置（抵消第 1 步的预移位）。

为什么二进制这么爽？ 十进制长除法每轮要试商 0~9，而二进制每轮只有 0 或 1 两种可能，试商就退化成一次比较。

顺带一提：因为这种结构出现频率太高，Synplify Pro 这类高级综合工具会对定点除法自动生成这种结构；如果声明的是 integer，综合器会用 32 位字宽并自动优化掉未使用的位。

这套紧凑架构适合能忍受较大延迟（迭代多轮才能出一个结果）的定点除法。如果面积允许、又需要对任意数做快速除法，就需要下面更高级的技术了。

8.1.3 Goldschmidt 方法：流水线高速除法

当需要把除法流水线化以最大化吞吐率时，轮到 Goldschmidt 出场。它属于逐次逼近（successive approximation）类算法——每迭代一次就更接近真商。它最大的好处是：硬件可以构造得每个时钟沿完成一次除法，而且面积比把上一节迭代法暴力展开要小得多。

核心思想：要算 Q = N/D，先近似求出 1/D，然后乘上 N，再通过逐次逼近收敛到真商。这对 IEEE 754 浮点这类大数特别有用——64 位浮点尾数可能有 50 多位，直接建 2⁵⁰ 大小的查找表不可行，但 2¹⁰ 量级的缩小查找表（Goldschmidt 典型配置）就实用了。

算法步骤（求 Q = N/D）：

1. 规格化：移动 N 和 D 的小数点，使得 N ≥ 1 且 D < 2（浮点语境下叫对分子分母做 normalizing）。
2. 查表得初值：用查找表得到 1/D 的初始近似 L₁，通常 8~16 位精度就够。
3. 第一轮近似：算出 q₁ = L₁·N 和误差项 e₁ = L₁·D（随着迭代趋于无穷，e₁ 会逼近 1）。
4. 开始迭代：令 L₂ = −e₁（二进制补码，实际含义是 2 − e₁）。
5. 更新：e₂ = e₁·L₂，q₂ = q₁·L₂。
6. 令 L₃ = −e₂（同第 4 步），如此继续迭代。

每迭代一轮，eᵢ 就更接近 1（也就是”除数 × 1/D = 1”），qᵢ 就更接近真正的商 Q。实践中，4~5 轮迭代通常就能满足 64 位浮点（53 位定点）运算的精度要求。

一个手算数值例子

以 N = 3.0、D = 1.5（真商 Q = 2）为例，假设查找表只给 4 位精度，输出 L₁ = 0.625（1/1.5 ≈ 0.6667 的粗略近似）：

迭代轮次	Lᵢ = 2 − eᵢ₋₁	商近似 qᵢ = qᵢ₋₁·Lᵢ	误差项 eᵢ = eᵢ₋₁·Lᵢ
第 1 轮（查表）	L₁ = 0.625	q₁ = 3.0 × 0.625 = 1.875	e₁ = 1.5 × 0.625 = 0.9375
第 2 轮	L₂ = 2 − 0.9375 = 1.0625	q₂ = 1.875 × 1.0625 = 1.9921875	e₂ = 0.9375 × 1.0625 = 0.99609375
第 3 轮	L₃ = 2 − 0.99609 = 1.00390625	q₃ ≈ 1.99997	e₃ ≈ 0.99998

可以看到 qᵢ 一路奔向 2，eᵢ 一路奔向 1——这就是”逐次逼近”的含义。误差的严格上界与迭代轮数、位宽的关系，可以参考专门讲 Goldschmidt 的论文。

Verilog 实现：53 位全流水线除法器

下面的例子实现 53 位除法（IEEE 754 双精度浮点所需），采用全流水线架构——输入假定已规格化（规格化由层次中的其他模块完成），不做溢出检查：

module div53(
    output [105:0] o,          // quotient
    input clk,
    input [52:0] a, b);        // dividend and divisor
    reg [261:0] mq5;
    reg [65:0] k2;
    reg [130:0] k3;
    reg [130:0] k4;
    reg [130:0] k5;
    reg [52:0] areg, breg;
    reg [65:0] r1, q1;
    reg [130:0] r2, q2;
    reg [130:0] r3, q3;
    reg [130:0] q4;
    wire [13:0] LutOut;
    wire [13:0] k1;
    wire [66:0] mr1, mq1;
    wire [131:0] mr2, mq2;
    wire [261:0] mr3, mq3;
    wire [261:0] mr4, mq4;
 
    gslut gslut(.addr(b[51:39]),
        .clk(clk),
        .dout(LutOut));        // 用除数高 13 位查表得初始近似
 
    assign k1 = LutOut;
 
    assign o = mq5[261-1:261-1-105];
 
    assign mr1 = breg * k1;    // e1 = L1 * D
    assign mq1 = areg * k1;    // q1 = L1 * N
 
    assign mr2 = r1 * k2;      // 第二轮迭代
    assign mq2 = q1 * k2;
 
    assign mr3 = k3 * r2;      // 第三轮迭代
    assign mq3 = k3 * q2;
 
    assign mr4 = k4 * r3;      // 第四轮迭代
    assign mq4 = k4 * q3;
 
always @(posedge clk) begin
    areg <= a;
    breg <= b;
    r1 <= mr1[65:0];
    k2 <= ~mr1[65:0] + 1;      // L2 = 取补码，即 2 - e1
    q1 <= mq1[65:0];
    r2 <= mr2[130:0];
    k3 <= ~mr2[130:0] + 1;     // L3 = 2 - e2
    q2 <= mq2[130:0];
    r3 <= mr3[260:130];
    k4 <= ~mr3[260:130] + 1;   // L4 = 2 - e3
    q3 <= mq3[260:130];
    k5 <= ~mr4[260:130] + 1;   // L5 = 2 - e4
    q4 <= mq4[260:130];
    mq5 <= k5 * q4;            // 最后一轮乘出商
end
endmodule

代码点评几个关键行：

• gslut 查找表：用除数 b 的高位（b[51:39]，13 位）查 1/D 的初始近似，这就是前面说的”2¹⁰ 量级的缩小查找表”。
• k2 <= ~mr1[65:0] + 1：按位取反加一就是二进制补码（即取负）。在定点语境下它实现的正是 L = 2 − e 这一步，是每轮迭代的核心。
• mr*/mq* 组合逻辑乘法 + 寄存器打拍：每轮迭代夹一级乘法器和一级寄存器，形成流水线，每个时钟都能接受新输入。
• assign o = mq5[261-1:261-1-105]：从最终 262 位乘积中截取商所在的位段。

注意：这个实现是完全展开、追求最大速度（每时钟一次除法），综合后面积也相对大。面积优化手段包括：用单个乘法器加状态机轮流迭代各轮系数，或者用重复移位-加法的紧凑乘法器。速度和面积，永远是取舍。

8.2 泰勒与麦克劳林级数展开

泰勒（Taylor）和麦克劳林（Maclaurin）级数可以把指数、三角函数、对数这类”硬件不认识的函数”分解成乘法和加法——这才是硬件擅长的东西。

泰勒展开的一般形式（式 8.1）：

$$\mathrm{T}(\mathrm{x})=\sum _0^{\infty }\frac{{\mathrm{f}}^{(\mathrm{n})}(\mathrm{a})(\mathrm{x}-\mathrm{a}{)}^{\mathrm{n}}}{\mathrm{n}!}$$

其中 f⁽ⁿ⁾ 是 f 的 n 阶导数。实践中常取 a = 0，此时简化为麦克劳林级数（式 8.2）：

$$\mathrm{M}(\mathrm{x})=\sum _0^{\infty }\frac{{\mathrm{f}}^{(\mathrm{n})}(0)(\mathrm{x}{)}^{\mathrm{n}}}{\mathrm{n}!}$$

几个最常用的展开式：

$$\sin (x)=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\dots$$ $$\cos (x)=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\dots$$ $$\ln (1+x)=x-\frac{x^2}{2}+\frac{x^3}{3}-\dots$$ $${\mathrm{e}}^{\mathrm{x}}=1+\mathrm{x}+\frac{{\mathrm{x}}^2}{2!}+\frac{{\mathrm{x}}^3}{3!}+\dots$$

展开的阶数越高，近似越准，但覆盖的角度范围也要考虑。下图展示了正弦波近似随展开阶数增加而逐步贴合真值的过程：

对硬件设计者来说，这套级数最妙的地方在于：所有分母（3!、5!、2! 等）都是固定常数，可以预先求倒数，然后用上一节的”定点乘法”来实现除法。于是 sin、cos、ln、eˣ 全都变成了”乘加链条”。

泰勒级数的主要缺点有两个：

1. 乘法次数多——每一项都是一次乘法；
2. 迭代耗时长——紧凑架构（移位-加法乘法器）会把算法拖到几百个时钟周期。

这两个问题往往是绑在一起的。下一节的 CORDIC 用向量旋转做二进制近似，能大幅提速。

8.3 CORDIC 算法

CORDIC（Coordinate Rotation Digital Computer，坐标旋转数字计算机）是一种逐次逼近算法，用于非常高效地计算正弦和余弦。它的原理：用一连串的向量旋转来逼近目标角度。

先建立直观理解——用画图法求 sin/cos：

1. 在 x-y 平面上画一个模长为 1、相位为 0 的向量（图 8.7）；
2. 逆时针旋转这个向量直到达到目标角度，全程保持模长为 1（图 8.8）；
3. 读出目标角度处的 (x, y) 坐标：y 就是 sin，x 就是 cos（斜边为 1，所以 sin = y/1，cos = x/1）（图 8.9）。

硬件实现时，旋转的角度增量是 90° 除以 2 的逐级递增的幂（也就是越转步子越小）：第一次转 45°，再转 22.5°，再转 11.25°……每跳一步就更新一次 x-y 坐标。加还是减当前增量，取决于当前累计角度相对目标角度的位置——这样就以越来越小的步长逐次逼近目标角度。迭代方程如下：

$${\mathrm{x}}_{\mathrm{i}+1}={\mathrm{K}}_{\mathrm{i}}[{\mathrm{X}}_{\mathrm{i}}-{\mathrm{y}}_{\mathrm{i}}{\mathrm{d}}_{\mathrm{i}}{2}^{-\mathrm{i}}]$$ $${\mathrm{y}}_{\mathrm{i}+1}={\mathrm{K}}_{\mathrm{i}}[{\mathrm{y}}_{\mathrm{i}}+{\mathrm{x}}_{\mathrm{i}}{\mathrm{d}}_{\mathrm{i}}{2}^{-\mathrm{i}}]$$ $${\mathrm{K}}_{\mathrm{i}}=(1+2^{-2\mathrm{i}}{)}^{-1/2}$$

怎么理解这几个式子：

• 2⁻ⁱ 就是”乘以越来越小的角度步长”——注意乘以 2 的负幂在硬件里就是右移，不花资源！
• dᵢ 是决策项：目标角度大于累计角度时 dᵢ = +1（增大角度，y 增 x 减）；小于时 dᵢ = −1（减小角度，y 减 x 增）。
• Kᵢ 是模长校正因子，由迭代级数唯一决定，当 i → ∞ 时收敛到 0.60725…。每旋转一次模长会稍微变大（乘了 √(1+2⁻²ⁱ)），所以要用 Kᵢ 补偿。实践中设计者会根据迭代轮数提前算好 Kᵢ，在计算的最后一次性乘上这个常数因子——一次常数乘法，其余全程无乘法。

对比一下：CORDIC 的每轮迭代只有加/减法和比较，全部可以单周期完成；唯一的乘法只在最后出现一次，而且是常数乘法还能进一步优化。所以 CORDIC 要么跑得比泰勒展开快得多，要么同速度下门数少得多。

结论：计算 sin/cos，CORDIC 应优先于泰勒展开。

8.4 要点总结

• 乘移法是做除法最省事的办法，但只适用于除数可预先取倒数（固定常数）的场景。
• 紧凑的迭代除法架构适合能忍受较大延迟的定点除法。
• Goldschmidt 方法把除法流水线化，比暴力展开迭代法高效得多。
• 泰勒/麦克劳林级数能把复杂函数分解成硬件友好的乘加操作。
• 算 sin/cos 时，CORDIC 应当优先于泰勒展开。

常见坑：用乘移法做除法时忘了误差评估

乘移法的”除以 7”实际上是”除以 7.013…”。很多工程师直接抄了近似系数就上项目，结果在精度敏感的场合（如 DSP 滤波器累加路径）误差逐级放大。正确做法：先根据系统精度要求确定需要保留的有效位，再反推 2 的幂次和乘法因子；另外别忘了乘法结果的位宽膨胀——x × 73 会让中间结果比输入宽好几个比特，截位/舍入策略也要一并设计。

通关标准：

学完本篇你应该能做到：

• 说出三种硬件除法方案（乘移法、迭代除法、Goldschmidt）各自的适用条件：除数是否固定、延迟要求、吞吐率要求；
• 手工推导”除以常数 D”的乘移系数（乘法因子 + 右移位数）；
• 解释 Goldschmidt 中误差项 eᵢ 为什么趋近 1，并手算 2~3 轮迭代验证收敛；
• 写出 CORDIC 的迭代方程，说明 dᵢ 如何决策、Kᵢ 为什么可以最后统一补偿；
• 给一个具体需求（如 64 位浮点除法）选择合适的算法并说明理由。

自测：为什么二进制的迭代除法比十进制长除法在硬件上更简单？

因为二进制下每轮试商只有 0 或 1 两种可能：预规格化保证部分被除数 < 2×除数后，每轮只需一次比较（除数是否 ≤ 部分被除数），是则商 1 并相减，否则商 0。十进制则要试探 0~9，硬件上代价大得多。

自测：Goldschmidt 为什么需要查找表？查找表为什么可以做得这么小？

Goldschmidt 需要一个 1/D 的初始近似 L₁ 作为迭代起点，初值越准收敛越快。但 D 是 50 多位的浮点尾数，直接建 2⁵⁰ 的表不可行；而迭代算法对初值精度要求不高（816 位就够），所以只用除数的高几位（如高 13 位）做地址、存 2¹⁰ 量级的表即可，剩下的精度靠 45 轮迭代补回来。

自测：CORDIC 中 2⁻ⁱ 因子为什么"不花硬件资源"？

乘以 2⁻ⁱ 等价于右移 i 位，在硬件中只是连线重新接线（wiring），不需要任何逻辑资源。这正是 CORDIC 高效的关键：迭代只含加减、比较和移位，全部可单周期完成。

自测：Kᵢ 校正因子为什么不在每轮迭代中相乘，而是最后统一乘一次？

Kᵢ = (1+2⁻²ⁱ)⁻¹/² 只由迭代级数 i 决定，与输入无关，所有轮次的 Kᵢ 连乘是一个固定常数。把它合并成末尾的一次常数乘法，可以把迭代路径完全做成无乘法的加减移位链，既省资源又提速。

自测：泰勒级数方案在什么情况下仍比 CORDIC 合适？

当目标函数 CORDIC 覆盖不了时（CORDIC 天生擅长 sin/cos/atan 等旋转类函数），比如 ln、eˣ 这类函数用泰勒级数更自然；另外当输入范围很窄、展开两三项就够精度时，泰勒级数的乘法次数很少，实现也很简单。