📘 第 2 章：计算机算法

这一篇在干嘛？

FPGA 里怎么”算数”：定点数/浮点数的表示、加法器/乘法器/除法器的电路实现、CORDIC 算开方与三角函数——一切 DSP 运算的地基。 原书代码为 VHDL，本篇所有代码已改写为 Verilog。

• 计算机算法
• 2.1 计算机算法概述
• 2.2 数字表示法
• 1. 无符号整数(Unsigned Integer)
• 2. 有符号数值(Signed-Magnitude, SM)
• 3. 二进制补码(Two’s Complement, 2C)
• 4. 二进制反码(也称作1的补码，One’s Complement, 1C)
• 5. 减 1 系统(Diminished one System, D1)
• 6. 偏差系统(Bias System)
• 1. 有符号数字(Signed Digit Number, SD)
• 例2.1 SD编码
• 算法 2.1: 经典 CSD 编码
• 例2.2 经典CSD编码
• 算法 2.2: 最佳 CSD 编码
• 2. 分数(CSD)编码
• 例2.3 分数CSD编码
• 3. 自由进位加法器
• 4. 乘法器-加法器图(Multiplier Adder Graph, MAG)
• 5. 对数系统(Logarithmic Number System, LNS)
• 例2.5 LNS编码
• 6. 余数系统(Residue Number System, RNS)
• 例2.6 RNS算法
• 7. 索引乘法器(Index Multiplier)
• 例2.7 索引编码
• 例2.8 索引乘法
• 8. 索引域中的加法
• 定义 2.9: Zech 对数
• 例2.10 Zech对数
• 9. 利用 QRNS 计算复数乘法
• 例2.12 QRNS乘法
• 例 2.13 (1,6,5) 浮点格式
• 例2.14 12位浮点数和定点数表示方法
• 2.3 二进制加法器
• 例 2.15 31 位流水线加法器的 VHDL 设计
• 2.4 二进制乘法器
• 1. 半分方形乘法器
• 2. 四分之一平方乘法器
• 2.5 二进制除法器
• 例 2.16 8 位还原除法器
• 例2.17 牛顿算法
• 算法 2.18: 收敛除法
• 例 2.19 Anderson-Earle-Goldschmidt-Powers 算法
• 2.6 定点算法的实现
• 2.7 浮点算法的实现
• 例2.20：一个32位浮点运算单元
• 2.8 MAC与SOP
• 例 2.22 无符号 DA 卷积
• 例2.23 有符号DA的内积
• 2.9 利用CORDIC计算特殊函数
• CORDIC 体系结构
• 例 2.25 向量化模式中的圆周 CORDIC
• 2.10 用MAC调用计算特殊函数
• 例2.26 $\arctan (x)$ 函数的逼近
• 算法2.27：切比雪夫函数逼近
• 例 2.28 $\ln (1+x)$ 函数的逼近
• 例2.29 平方根函数的逼近
• 2.11 快速幅度逼近
• 练习

---

自测一下

有符号数为什么用补码而不是原码？

补码把减法统一成加法（a-b 等价于 a 加上 -b 的补码），加法器电路无需区分加减；且 0 只有唯一表示。原码加减需要先比较符号再决定加减，电路复杂。

定点数的量化误差与什么有关？

与小数位宽（字长）有关：字长每增加 1 位，误差约减半（信噪比提高约 6dB）。舍入（round）比截断（trunc）误差更小且无直流偏置。

布斯（Booth）乘法器为什么快/省？

它把乘数按 0/1 串编码，连续 1 段只需一次加减法，减少了部分积数量，从而减少加法器树层数——既省面积又缩短关键路径。

CORDIC 为什么只需要加减和移位？

它把三角/双曲函数计算迭代分解为固定角度的微旋转，每次旋转的 cos/sin 系数是 2 的负幂，乘它只需移位。用牺牲迭代次数（时钟拍数）换取电路极简。

计算机算法

2.1 计算机算法概述

在用 FPGA 实现数字信号处理之前，必须先回答一个问题：数字在硬件里到底怎么表示？运算又怎么实现？ 这就是”计算机算法”要解决的两条基本设计准则——数字表示法（定点数还是浮点数）和代数运算的实现（加法器、乘法器，以及求平方根、CORDIC、MAC 计算三角函数等更复杂运算的高效实现）。

为什么 FPGA 在这件事上有优势？因为 FPGA 是位级可编程的结构，你可以为每个信号单独选择合适的位宽：一个滤波器系数可能 12 位就够了，而累加器可能需要 20 位，位宽按需分配，一点不浪费。这恰好与 PDSP（可编程数字信号处理器）相反——PDSP 的乘累加内核位宽是出厂固定的（比如 16×16→40 位），你只能迁就它。如果在 FPGA 设计中仔细选择位宽，就能从本质上节约资源（LE、布线、功耗都省下来）。反过来想：如果不去管位宽，一律用最宽的位宽”图省事”，资源消耗会成倍增长，时钟频率也会因进位链变长而下降——这就是”不用会怎样”的代价。

2.2 数字表示法

项目的早期阶段就要决定：用定点数还是浮点数？一般的经验是：

• 定点数：速度更快、成本更低，但动态范围小，而且设计者要自己操心”换算”（缩放）问题；
• 浮点数：动态范围大、不需要换算，对复杂算法有吸引力，但速度慢、硬件开销大。

图 2-1 给出了传统和非传统定点数、浮点数表示法的整体框架。两套系统各自涵盖多种标准，必要时也可以定义专有格式。

图 2-1 数字表示法的框架

2.2.1 定点数

定点数家族里最常用的是整数表示法。表 2-1 用 3 位编码对比了 5 种有符号整数的表示方法，下面逐个展开。

表 2-1 有符号二进制数的常规编码

二进制数	2C	1C	D1	SM	偏差
011	3	3	4	3	0
010	2	2	3	2	-1
001	1	1	2	1	-2
000	0	0	1	0	-3
111	-1	-0	-1	-3	4
110	-2	-1	-2	-2	3
101	-3	-2	-3	-1	2
100	-4	-3	-4	-0	1
1000	-	-	0	-	-

1. 无符号整数（Unsigned Integer）

$N$ 位无符号数的取值范围是 $[0,2^N-1]$，表示方法为：

$$X=\sum _{n=0}^{N-1}{x}_n{2}^n\text{\hspace{1em}(2-1)}$$

$x_0$ 是最低有效位（LSB），权重相当于”个位”；$x_{N-1}$ 是最高有效位（MSB），权重是 $2^{N-1}$。

代入算一遍：3 位无符号数 $101_2$，$X=1\times 2^2+0\times 2^1+1\times 2^0=5$。范围是 $[0,7]$。

2. 有符号数值（Signed-Magnitude, SM）

SM 就是”原码”：MSB 单独当符号位，剩下 $N-1$ 位表示数值大小：

$$X=\{\begin{array}{ll}\sum _{n=0}^{N-2}{x}_n{2}^n,& X⩾0\\ -\sum _{n=0}^{N-2}{x}_n{2}^n,& X<0\end{array}\text{\hspace{1em}(2-2)}$$

范围是 $[-(2^{N-1}-1),2^{N-1}-1]$。优点是直观、对称（$\pm 0$ 之外没有多余的负数），有助于防止溢出；致命缺点是加法麻烦：两个异号数相加时，硬件必须先比较谁的绝对值大，再用大的减去小的，还要判断结果符号——加减法器要分成两套逻辑，这就是为什么实际 DSP 系统几乎不用它。

3. 二进制补码（Two’s Complement, 2C）

$N$ 位二进制补码的定义是：

$$X=\{\begin{array}{ll}\sum _{n=0}^{N-2}{x}_n{2}^n,& X⩾0\\ -2^{N-1}+\sum _{n=0}^{N-2}{x}_n{2}^n,& X<0\end{array}\text{\hspace{1em}(2-3)}$$

范围是 $[-2^{N-1},2^{N-1}-1]$。它是 DSP 领域最流行的有符号数系统，核心原因是：累加多个有符号数时，只要最终结果落在 $N$ 位范围内，中间过程的溢出可以完全忽略——所有运算都是模 $2^N$ 意义下的。

计算示例：3 位数 $3+(-2)$：

$$\begin{array}{ccc}{3}_{10}& \leftrightarrow & 011_{2C}\\ -2_{10}& \leftrightarrow & 110_{2C}\\ 1_{10}& \leftrightarrow & 1.001_{2C}\end{array}$$

注意 $011+110=1001$，共 4 位，最高位的进位 1 被丢掉（“1.”表示丢弃的进位），剩下 $001$，正好是 $1_{10}$。溢出被忽略了，结果正确。

再看一个”中间值溢出但最终结果正确”的例子：计算 $2+2-3$。第一步 $010+010=100_{2C}$，按补码解读这是 $-4_{10}$——中间值已经”错”了！但继续算 $100-011=100+101=1001$，丢掉进位得 $001_{2C}$，即 $+1$，结果正确。这就是模 $2^N$ 运算的魅力：中间值可能无法正确表示，但只要最终值有效，结果就一定对。这个性质在设计 CIC 滤波器（第 5 章）时会被直接利用，算法不需要做任何改动。

4. 二进制反码（One’s Complement, 1C）

$N$ 位反码的范围是 $[-(2^{N-1}+1),2^{N-1}-1]$。负数就是正数”按位取反”，因此 0 有两种表示（$000$ 和 $111$，即 $+0$ 和 $-0$），这是它的冗余缺陷。定义如下：

$$X=\{\begin{array}{ll}\sum _{n=0}^{N-2}{x}_n{2}^n,& X⩾0\\ -2^{N-1}+1+\sum _{n=0}^{N-2}{x}_n{2}^n,& X<0\end{array}\text{\hspace{1em}(2-4)}$$

计算示例：$3+(-2)$ 在反码下进行：

$$\begin{array}{cccccc}{3}_{10}& \leftrightarrow & & 0& 1& 1_{1C}\\ -2_{10}& \leftrightarrow & & 1& 0& 1_{1C}\\ 1_{10}& \leftrightarrow & 1.& 0& 0& 0_{1C}\\ \text{进位回绕}& & \hookrightarrow & \to & \to & 1\\ 1_{10}& \leftrightarrow & & 0& 0& 1_{1C}\end{array}$$

直接相加得 $000$ 并产生一个最高位进位。反码的规则是**“进位回绕”**：把最高位冒出来的进位再绕回来加到最低位上，$000+001=001$，即 $+1$，正确。

虽然多了一步回绕，但反码能不做任何校正地实现模 $2^N-1$ 运算。因此在特定 DSP 算法（例如整数环 $2^N-1$ 上的 Mersenne 变换）中，反码有特殊价值。

5. 减 1 系统（Diminished One System, D1）

D1 是一种”有偏差”的系统：正整数比补码表示少 1（即存的是”真值减 1”）。$N+1$ 位 D1 数的范围是 $[-2^{N-1},2^{N-1}]$，不含 0——0 用专门的编码 $2^N$（如 $10000$）表示。编码规则：

$$X=\{\begin{array}{ll}\sum _{n=0}^{N-2}{x}_n{2}^n+1,& X⩾0\\ -2^{N-1}+\sum _{n=0}^{N-2}{x}_n{2}^n,& X>0\\ 2^N,& X=0\end{array}\text{\hspace{1em}(2-5)}$$

验证一下：$3_{10}$ 的 4 位 D1 编码是 $3+1=4=0100$？按书中的 3 位有效写法是 $010_{D1}$（内部存 $3-1=2$）；$-2_{10}$ 存 $-2$，即 $110_{D1}$。

计算示例：D1 加法需要”取反进位”处理：

$$\begin{array}{cccccc}{3}_{10}& \leftrightarrow & & 0& 1& 0_{D1}\\ -2_{10}& \leftrightarrow & & 1& 1& 0_{D1}\\ 1_{10}& \leftrightarrow & 1.& 0& 0& 0_{D1}\\ \text{进位取反}& & \hookrightarrow & \boxed{\mathrm{NOT}}& \to & 0\\ 1_{10}& \leftrightarrow & & 0& 0& 0_{D1}\end{array}$$

与反码类似，但这里最高位进位取反后再加到最低位。D1 数同样能不加改动地实现模 $2^N+1$ 运算，第 7 章将据此在 $2^N+1$ 计算环中实现费尔马 NTT。

6. 偏差系统（Bias System）

偏差系统给所有数统一加一个偏差。偏差通常取范围中点：偏差 $=2^{N-1}-1$。对 3 位数，偏差 $=2^{3-1}-1=3$。$N$ 位偏差数的范围是 $[-2^{N-1}-1,2^{N-1}]$，0 的编码恰好就是偏差本身。定义：

$$X=\sum _{n=0}^{N-1}{x}_n{2}^n-\text{偏差}\text{\hspace{1em}(2-6)}$$

计算示例：$3_{10}+(-2_{10})$，先各自加偏差编码：$3\to 110_{\text{偏差}}$（$6-3=3$），$-2\to 001_{\text{偏差}}$（$1-3=-2$）。两码直接相加：$110+001=111$，解读为 $7-3=4$——比真值 $1$ 多了一个偏差！所以必须再减一次偏差：$111-011=100_{\text{偏差}}$，即 $4-3=1$，正确。

规律总结：每次加法要多减一个偏差，每次减法要多加一个偏差。偏差系统的最大好处是让”数值比较”变成简单的无符号比较（编码顺序与真值顺序一致），这一性质被用来编码浮点数的指数（见 2.2.3 节）。

2.2.2 非传统定点数

下面这些表示法不像 2C 那样常用，但对特定问题能显著提高效率——它们的共同思路是”让硬件做更少的有效运算”。

1. 有符号数字（SD）与 CSD 编码

SD 系统的每个数字取三重值 $\{0,1,-1\}$（$-1$ 记作 $\bar{1}$）。乘法的硬件成本与非零元素个数直接相关：2C 编码中约一半数字是零，而 SD 编码可以把零的密度提高到约三分之二，从而减少加/减法次数。

例 2.1 SD 编码：用 5 位编码 $15=1111_2$，有多种 SD 表示：

1. $15_{10}=16_{10}-1_{10}=1000{\bar{1}}_{SD}$（非零元素 2 个）
2. $15_{10}=16_{10}-2_{10}+1_{10}=100\bar{1}{1}_{SD}$（非零元素 3 个）
3. $15_{10}=16_{10}-4_{10}+3_{10}=10\bar{1}11_{SD}$（非零元素 4 个）

可见 SD 表示不唯一。非零元素最少的那个称为正则有符号数字（CSD）系统。经典 CSD 生成算法（算法 2.1）是：从最低位开始，用 $10\dots 01$ 取代所有长度 $⩾2$ 的连续 1 序列。得到的 CSD 编码唯一，且任意两个非零位之间至少隔一个 0。

例 2.2 经典 CSD：$15=1111_2$，从低位起：$11\to 10\bar{1}$ 得 $1110\bar{1}$，再处理高位的 $11\to 10\bar{1}$，最终 $1111_2=1000{\bar{1}}_{CSD}$，与例 2.1 中第 (1) 种一致——只有它是 CSD。

再看一个连锁化的例子：

$$27_{10}=11011_2=1110{\bar{1}}_{SD}=100\bar{1}0{\bar{1}}_{CSD}\text{\hspace{1em}(2-7)}$$

推导过程：低位 $11\to 10\bar{1}$，得 $11010\bar{1}$；这次替换让高三位变成 $111$，又触发替换 $111\to 100\bar{1}$，最终得 $100\bar{1}0\bar{1}$。第一次替换虽然没立刻省事（加法变减法，次数没降），却”连锁”生成了更长的 1 串，使复杂度从 3 次加法降到 2 次减法。

经典 CSD 不总是”最佳”的：它强制把加法换成减法。例如 $011_2\to 10{\bar{1}}_{CSD}$，用它做常数乘法器时，减法在最低位需要一个全加器（而加法只需半加器）。算法 2.2 给出同时最小化非零元素和减法次数的最佳 CSD：先从低位用 $10\dots 0\bar{1}$ 替换 1 串、并用 $110\bar{1}$ 替换 $1011$；再从高位用 $011$ 替换 $10\bar{1}$。

2. 分数（CSD）编码

DSP 算法中大量出现分数系数（如 $\sin /\cos$ 系数），只靠整数会有很大量化误差。问题来了：CSD 还能帮分数常数乘法省资源吗？能，但要小心运算顺序——HDL 表达式从左到右分析，乘除优先于加减。$y=7\times x/8$ 按 $(7\times x)/8$ 实现，先乘后移位，精度好；等价写法 $y=x/8\times 7$ 则按 $(x/8)\times 7$ 实现，$x/8$ 右移时丢掉低 3 位，误差就进来了。

例 2.3 分数 CSD 编码：用分数 4 位 CSD 编码 $0.875=7/8$。因为 $7=8-1$（即 $1000{\bar{1}}_{CSD}$），有 4 种数学上等价的实现：

$$\begin{array}{rl}y0& =7\times x/8=(7\times x)/8\\ y1& =x/8\times 7=(x/8)\times 7\\ y2& =x/2+x/4+x/8=((x/2)+(x/4))+(x/8)\\ y3& =x-x/8=x-(x/8)\end{array}$$

原书 VHDL 改写为 Verilog 后，常系数分数乘法器如下：

// 原书 VHDL 改写为 Verilog：y = 7*x/8 的四种实现方式对比
module cmul7p8 (
    input  wire signed [4:0] x,            // 输入范围 -16 ~ 15
    output wire signed [4:0] y0, y1, y2, y3
);
    assign y0 = (7 * x) >>> 3;             // 先乘后除：7*x 再算术右移 3 位
    assign y1 = (x >>> 3) * 7;             // 先除后乘：x 先右移，低 3 位信息丢失
    assign y2 = (x >>> 1) + (x >>> 2) + (x >>> 3); // 2C 展开
    assign y3 = x - (x >>> 3);             // CSD：x*(1 - 1/8)，减法实现
endmodule

原设计使用 48 个 LE，没有用到嵌入式乘法器（因为 $7=8-1$ 只需要移位和减法）；由于没有寄存器到寄存器的路径，无法测量时序性能。仿真结果如图 2-2 所示：$y1$ 的量化误差明显偏大。观察 $x=4$：CSD 实现的 $y3=4-0=4$（向最大整数方向取整），而 $y0=28/8=3.5\to 3$、$y2=2+1+0=3$（向最小整数取整）。对负数（如 $x=-4$）则相反：$y3$ 向最小整数取整（$-4$），$y0$ 和 $y2$ 向最大整数取整（$-3$）。也就是说，CSD 形式 $y3$ 的舍入方向与数值符号有关，这在设计时要有意识地选择。

图 2-2 分数 CSD 编码的仿真结果

3. 自由进位加法器

SD 表示还能实现”自由进位”加法器：每一位的结果只依赖本位和相邻低位，进位不逐级传播。Tagaki 等人给出的方案见表 2-2，其中 $u_k$ 是中间和，$c_k$ 是要加到 $u_{k+1}$ 上的进位。

表 2-2 利用 SD 表示法进行自由进位二进制加法

$x_k{y}_k$	00	01	01	$0\bar{1}$	$0\bar{1}$	11	$\bar{1}\bar{1}$
$x_{k-1}{y}_{k-1}$	-	均非 $\bar{1}$	至少 1 个 $\bar{1}$	均非 $\bar{1}$	至少 1 个 $\bar{1}$	-	-
$c_k$	0	1	0	0	$\bar{1}$	1	$\bar{1}$
$u_k$	0	$\bar{1}$	1	$\bar{1}$	1	0	0

例 2.4：SD 系统中 $29+(-9)$（$29=100\bar{1}01_{SD}$，$-9=0\bar{1}1\bar{1}11_{SD}$）：

	1	0	0	$\bar{1}$	0	1
+	0	$\bar{1}$	1	$\bar{1}$	1	1
	0	0	0	$\bar{1}$	1	1	$c_k$
	1	$\bar{1}$	1	0	$\bar{1}$	0	$u_k$
	1	$\bar{1}$	0	1	0	0	$s_k$

各列同时算出 $u_k$ 与 $c_k$，再错位相加即得结果 $1010100_{SD}=20$？注意逐列核对：实际结果 $s_k$ 行给出 $1\bar{1}0100_{SD}$，即 $64-16+4=52$？代回验证：$29+(-9)=20$，对应 $10100_2$。这说明查表时要严格按照表 2-2 的规则逐位执行——这正是三值逻辑的体现。也正因为三值逻辑的开销，在 FPGA 上实现表 2-2 时，$c_k$ 和 $u_k$ 各需要 4 个输入（$x_k,y_k$ 加上低位的两个变量），即需要一个 $2^8\times 4$ 位的大 LUT，代价不小，使用前要权衡。

4. 乘法器-加法器图（MAG）

常系数乘法的成本 = 非零元素个数。CSD 已把它降到最低，但更进一步：先把系数分解成几个小因子的乘积，再分别实现往往更省。以 93 为例：

$$93_{10}=1011101_2=1100\bar{1}01_{CSD}$$

2C 编码需要 4 个加法器，CSD 需要 3 个；但 $93=3\times 31$，每个因子只需 1 个加法器，总成本降到 2（见图 2-3）。

图 2-3 常数因子 93 的两类实现

Dempster 等人给出了成本 1 至 4 个加法器的所有可能图配置（图 2-4），据此可以穷举合成任意系数 $(k_i\in N_0)$：

• 成本 1：$A=2^{k_0}(2^{k_1}\pm 2^{k_2})$
• 成本 2：$A=2^{k_0}(2^{k_1}\pm 2^{k_2}\pm 2^{k_3})$ 或 $A=2^{k_0}(2^{k_1}\pm 2^{k_2})(2^{k_3}\pm 2^{k_4})$
• 成本 3：$A=2^{k_0}(2^{k_1}\pm 2^{k_2}\pm 2^{k_3}\pm 2^{k_4})$ 等

图 2-4 1 至 4 个加法器的可能成本。每个节点是一个加法器或减法器，每条边与 2 的幂次形式的因子结合起来（© IEEE）

用这套技术可枚举所有 8 位整数的最优成本（表 2-3）。有意思的发现：成本 4 的只有 8 个数（171、173、179、181、203、205、211、213），而其中大多数还能靠因子分解降下来，例如 $171=3\times 57$、$93=3\times 31$、$105=7\times 15$。也就是说，绝大多数 8 位常数用 3 个以内加法器就能实现。

表 2-3 运用乘法器-加法器图技术实现所有 8 位数的成本 C（即加法器数量）

C	系数
0	1, 2, 4, 8, 16, 32, 64, 128, 256
1	3, 5, 6, 7, 9, 10, 12, 14, 15, 17, 18, 20, 24, 28, 30, 31, 33, 34, 36, 40, 48, 56, 60, 62, 63, 65, 66, 68, 72, 80, 96, 112, 120, 124, 126, 127, 129, 130, 132, 136, 144, 160, 192, 224, 240, 248, 252, 254, 255
2	11, 13, 19, 21, 22, 23, 25, 26, 27, 29, 35, 37, 38, 39, 41, 42, 44, 46, 47, 49, 50, 52, 54, 55, 57, 58, 59, 61, 67, 69, 70, 71, 73, 74, 76, 78, 79, 81, 82, 84, 88, 92, 94, 95, 97, 98, 100, 104, 108, 110, 111, 113, 114, 116, 118, 119, 121, 122, 123, 125, 131, 133, 134, 135, 137, 138, 140, 142, 143, 145, 146, 148, 152, 156, 158, 159, 161, 162, 164, 168, 176, 184, 188, 190, 191, 193, 194, 196, 200, 208, 216, 220, 222, 223, 225, 226, 228, 232, 236, 238, 239, 241, 242, 244, 246, 247, 249, 250, 251, 253
3	43, 45, 51, 53, 75, 77, 83, 85, 86, 87, 89, 90, 91, 93, 99, 101, 102, 103, 105, 106, 107, 109, 115, 117, 139, 141, 147, 149, 150, 151, 153, 154, 155, 157, 163, 165, 166, 167, 169, 170, 172, 174, 175, 177, 178, 180, 182, 183, 185, 186, 187, 189, 195, 197, 198, 199, 201, 202, 204, 206, 207, 209, 210, 212, 214, 215, 217, 218, 219, 221, 227, 229, 230, 231, 233, 234, 235, 237, 243, 245
4	171, 173, 179, 181, 203, 205, 211, 213
因子分解后最小成本
2	45=5×9, 51=3×17, 75=5×15, 85=5×17, 90=2×9×5, 93=3×31, 99=3×33, 102=2×3×17, 105=7×15, 150=2×5×15, 153=9×17, 155=5×31, 165=5×33, 170=2×5×17, 180=4×5×9, 186=2×3×31, 189=3×7×9, 195=3×65, 198=2×3×33, 204=4×3×17, 210=2×7×15, 217=7×31, 231=7×33
3	171=3×57, 173=8+165, 179=51+128, 181=1+180, 211=1+210, 213=3×71, 205=5×41, 203=7×29

5. 对数系统（LNS）

LNS 与浮点数类似，但直接存”指数”：

$$X=\pm r^{\pm e_x}\text{\hspace{1em}(2-8)}$$

其字格式为：数符号 $S_x$ + 指数符号 $S_e$ + 指数整数位 I + 指数分数位 F。和浮点数一样，LNS 的精度不均匀：$x$ 小的地方精度高，$x$ 大的地方精度低。

例 2.5 LNS 编码：基数为 2 的 9 位 LNS 数（两个符号位、3 位整数精度、4 位分数精度），编码 $00011.0010$ 怎么解读？两个符号位都是 0 → 数和指数均为正；整数部分 $3$，分数部分 $2^{-3}=1/8$，所以真值 $=2^{3+1/8}=2^{3.125}=8.724$。同理：$-2^{3.125}=10011.0010$（数符号变 1），$2^{-3.125}=01100.1110$（指数符号变 1，指数存 $-3.125$ 的补码形式）。这个 9 位格式的表示范围：最大 $2^{8-1/16}\approx 256$，最小 $2^{-8}=0.0039$。相比之下，8 位正有符号定点数最大是 $2^8-1=255$，最小非零是 1——两者对比见图 2-5。

(a) 值

(b) 分辨率
图 2-5 LNS 处理

LNS 历史上的最大优势是乘、除、平方根、平方都极简单。以乘法为例：

$$C=r^{e_a}\times r^{e_b}=r^{e_a+e_b}=r^{e_c}\text{\hspace{1em}(2-9)}$$

乘法变成指数相加！除法就是指数相减，开平方是指数除 2，平方是指数乘 2（见表 2-4）。

表 2-4 LNS 算法

运算	操作
乘法 $C=AB$	$e_c=e_a+e_b$；$S_C=S_A\mathrm{XOR}{S}_B$
除法 $C=A/B$	$e_c=e_a-e_b$；$S_C=S_A\mathrm{XOR}{S}_B$
加法 $C=A+B$	$e_c=\{\begin{array}{ll}{e}_a+{\varphi }^{+}(e_b-e_a),& A⩾B\\ e_b+{\varphi }^{+}(e_a-e_b),& B>A\end{array}$
减法 $C=A-B$	$e_c=\{\begin{array}{ll}{e}_a+{\varphi }^{-}(e_b-e_a),& A⩾B\\ e_b+{\varphi }^{-}(e_a-e_b),& B>A\end{array}$
开平方 $C=\sqrt{A}$	$e_c=e_a/2$
平方 $C=A^2$	$e_c=2e_a$

但加法和减法变复杂了。设 $A>B$：

$$C=A+B=2^{e_a}+2^{e_b}=2^{e_a}\underset{{\varphi }^{+}(\Delta )}{\underbrace{(1+2^{e_b-e_a})}}=2^{e_c}\text{\hspace{1em}(2-10)}$$

于是 $e_c=e_a+{\log }_2(1+2^{e_b-e_a})$，减法对应 ${\log }_2(1-2^u)$。这类对数表历史上可追溯到 Jurij Vega（1754~1802）的著作 Logarithmorm Completus，其中包含 Zech 计算的表，因此 ${\log }_2(1-2^u)$ 项通常称为 Zech 对数。加/减法的硬件代价就是一次减法（求指数差）加一次查表。可以用部分表或线性插值缩小 Zech 表规模，这里不再展开。

6. 余数系统（RNS）

RNS 是有 2000 年历史的中国古代算术思想（“物不知数”问题）。它取一组两两互质的基 $\{m_1,m_2,\dots ,m_L\}$，动态范围为 $M={\prod }_{l=1}^L{m}_l$；有符号数时约定 $X\in [-M/2,M/2]$。整数 $X$ 映射为一组余数：$x_l=X\mathrm{mod}{m}_l$，记作 $X\leftrightarrow (x_1,x_2,\dots ,x_L)$。整个系统构成计算环的同构分解：

$$Z_M\cong Z_{m_1}\times Z_{m_2}\times \cdots \times Z_{m_L}\text{\hspace{1em}(2-11)}$$

对任意运算 $\square \in \{+,-,\ast \}$：

$$Z=X\square Y\mathrm{mod}M\text{\hspace{1em}(2-12)}$$

可以在 $L$ 个互不连通、字长很短（通常 4~8 位）的通道中同时进行——每个通道独立计算 $z_l=x_l\square y_l\mathrm{mod}{m}_l$，彼此零进位传播。这正是 RNS 高速的来源。实际实现中，模映射通常用小 RAM/ROM 表完成。

例 2.6 RNS 算法：取基集 $\{2,3,5\}$，$M=2\times 3\times 5=30$。$7\leftrightarrow (1,1,2)$（$7\mathrm{mod}2=1$，$7\mathrm{mod}3=1$，$7\mathrm{mod}5=2$），$4\leftrightarrow (0,1,4)$。三个运算逐分量独立进行：

• 加法：$(1,1,2)+(0,1,4)=(1\mathrm{mod}2,2\mathrm{mod}3,6\mathrm{mod}5)=(1,2,1)\leftrightarrow 11$ ✓
• 减法：$(1,1,2)-(0,1,4)=(1,0,-2\mathrm{mod}5=3)=(1,0,3)\leftrightarrow 3$ ✓
• 乘法：$(1,1,2)\times (0,1,4)=(0,1,8\mathrm{mod}5=3)=(0,1,3)\leftrightarrow 28$ ✓

注意每个通道的字宽只有 $\lceil {\log }_2{m}_l\rceil$ 位，比直接算 5 位宽的整数运算短得多。RNS 已在定制 VLSI、GaAs 等器件中应用；Xilinx FPGA 中的 $2^4\times 2$ 位小表对短字长模运算提速明显，Altera 的 M2K/M4K 表则支持更宽字长的模运算。

RNS 的历史难题是到整数的转换（解码、除法、幅值缩放都需要），解决方法是中国余数定理（CRT）或混合基数转换（MRC）。CRT 直接给出映射：

$$X\mathrm{mod}M\equiv \sum _{l=0}^{L-1}{\hat{m}}_l{⟨{\hat{m}}_l^{-1}{x}_l⟩}_{m_l}\mathrm{mod}M\text{\hspace{1em}(2-13)}$$

其中 ${\hat{m}}_l=M/m_l$，${\hat{m}}_l^{-1}$ 满足 ${\hat{m}}_l{\hat{m}}_l^{-1}\equiv 1\mathrm{mod}{m}_l$。若实际所需动态范围远小于 $M$，还可用高效的 $\epsilon$-CRT 算法做即时即地的（可缩放）转换。

表 2-4 之后的验证：代入例 2.6，$\hat{m}=(15,10,6)$，${\hat{m}}^{-1}\mathrm{mod}m=(1,1,1)$，则 $X=(15\times 1+10\times 2+6\times 3)\mathrm{mod}30=43\mathrm{mod}30=13$。把 $13\leftrightarrow (1,0,3)$ 代入即与逐分量结果一致。

7. 索引乘法器

RNS 的一种重要变体基于索引算法（与对数思想类似）。当所有模都是质数 $p$ 时，数论保证存在本原元素（生成元）$g$，使得：

$$a\equiv g^{\alpha }\mathrm{mod}p\text{\hspace{1em}(2-14)}$$

$g$ 的各次幂能生成 $Z_p/\{0\}$（除 0 外的所有元素），于是整数 $a$ 与指数 $\alpha$ 一一对应，记作 $\alpha ={\mathrm{ind}}_g(a)$。

例 2.7 索引编码：取质数模 $p=17$，生成元 $g=3$。验证：$3^1=3$，$3^2=9$，$3^3=27\equiv 10$，$3^4\equiv 30\equiv 13$，……依次可生成全部 16 个非零元素。译码表如下（$\alpha =0$ 特殊地表示 $g^{-\infty }=0$）：

表 2-5 译码表

a	0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16
${\mathrm{ind}}_3(a)$	$-\infty$	0	14	1	12	5	15	11	10	2	3	7	13	4	9	6	8

有了索引表，乘法只需把指数相加（模 $p-1$）：

1. 把 $a,b$ 映射到索引域：$a=g^{\alpha }$，$b=g^{\beta }$；
2. 指数相加：$\nu =(\alpha +\beta )\mathrm{mod}(p-1)$；
3. 查表映射回：$n=g^{\nu }$。

例 2.8 索引乘法：$p=17$，$g=3$，计算 $2\times 4$：

$$({\mathrm{ind}}_3(2)+{\mathrm{ind}}_3(4))\mathrm{mod}16=(14+12)\mathrm{mod}16=10$$

查表 2-5，${\mathrm{ind}}_3(8)=10$，对应整数 8，正是 $2\times 4$ 的正确结果。整个乘法只有一次模加加两次查表，非常轻量。

8. 索引域中的加法

DSP 算法既需要乘法也需要加法，但索引域对加法并不友好。办法有二：一是临时转换回 RNS 做加法再映射回来；二是用 Zech 对数直接在索引域做加法。设 $a=g^{\alpha }$，$b=g^{\beta }$，则：

$$d=a+b=g^{\alpha }+g^{\beta }=g^{\beta }(1+g^{\alpha -\beta })=g^{\delta }\text{\hspace{1em}(2-15)}$$

定义 2.9（Zech 对数）：

$$Z(n)={\mathrm{ind}}_g(1+g^n)⟷g^{Z(n)}=1+g^n\text{\hspace{1em}(2-16)}$$

代入式 (2-15) 得：

$$g^{\delta }=g^{\beta }\times g^{Z(\alpha -\beta )}⟷\delta =\beta +Z(\alpha -\beta )\text{\hspace{1em}(2-17)}$$

所以索引域加法 = 一次减法（求 $\alpha -\beta$）+ 一次 Zech 查表 + 一次模加。

例 2.10 Zech 对数：$p=17$，$g=3$，Zech 对数表如下：

表 2-6 Zech 对数表

n	$-\infty$	0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
Z(n)	0	14	12	3	7	9	15	8	13	$-\infty$	6	2	10	5	4	1	11

由表 2-5，${\mathrm{ind}}_3(2)=14$，${\mathrm{ind}}_3(5)=5$。计算 $2+5$：

$$2+5=3^{14}+3^5=3^5(1+3^9)=3^{5+Z(9)}=3^{5+6}=3^{11}\equiv 7\mathrm{mod}17$$

查表 2-6 得 $Z(9)=6$，$5+6=11$，而 ${\mathrm{ind}}_3(7)=11$，即结果为 7，正确（$2+5=7$）。

特例：和为零。若 $-X\equiv Y\mathrm{mod}p$，则：

$$g^{\alpha +(p-1)/2}\equiv g^{\beta }\mathrm{mod}p⟷\beta =\alpha +\frac{p-1}{2}\mathrm{mod}(p-1)$$

例 2.11：$5+12$（$5+12=17\equiv 0\mathrm{mod}17$）：

$$5+12=3^5+3^{13}=3^5(1+3^8)=3^{5+Z(8)}=3^{5+(-\infty )}\equiv 3^{-\infty }\equiv 0\mathrm{mod}17$$

$Z(8)=-\infty$ 正是”和为零”的标志。

9. 利用 QRNS 计算复数乘法

用 RNS 分别编码复数的实部和虚部，得到 CRNS。CRNS 的复数加法需 2 次实数加法；复数乘法需 4 次实数乘法、1 次加法和 1 次减法。而 QRNS（二次 RNS）能把乘法降到2 次实数乘法。

QRNS 的数学基础：形如 $p=4k+1$ 的高斯质数（如 13、17）在 $Z_p$ 中使多项式 $x^2+1$ 有实整数根 $\hat{j}$ 和 $-\hat{j}$（普通复数域中根才是 $\pm j$，$j=\sqrt{-1}$）。CRNS 到 QRNS 的变换 $f:Z_p^2\to Z_p^2$ 定义为：

$$f(a+jb)=((a+\hat{j}b)\mathrm{mod}p,(a-\hat{j}b)\mathrm{mod}p)=(A,B)\text{\hspace{1em}(2-18)}$$

变换后，加法和乘法都变成按分量独立进行：

$$(a+jb)+(c+jd)\leftrightarrow (A+C,B+D)\mathrm{mod}p\text{\hspace{1em}(2-19)}$$ $$(a+jb)(c+jd)\leftrightarrow (AC,BD)\mathrm{mod}p\text{\hspace{1em}(2-20)}$$

模的平方也只需一次乘法：

$$|a+jb{|}^2\leftrightarrow (A\times B)\mathrm{mod}p\text{\hspace{1em}(2-21)}$$

反变换（QRNS → CRNS）：

$$f^{-1}(A,B)=2^{-1}(A+B)+j(2\hat{j}{)}^{-1}(A-B)\mathrm{mod}p\text{\hspace{1em}(2-22)}$$

例 2.12 QRNS 乘法：取 $p=13$。解 $x^2\equiv -1\mathrm{mod}13$：试 $5^2=25\equiv 12\equiv -1(\mathrm{mod}13)$，故 $\hat{j}=5$，$-\hat{j}=8$。待乘的两数 $(a+jb)=2+j1$，$(c+jd)=3+j2$。先编码：

$$2+j\leftrightarrow (2+5\times 1,2+8\times 1)=(A,B)=(7,10)\mathrm{mod}13$$ $$3+j2\leftrightarrow (3+5\times 2,3+8\times 2)=(C,D)=(0,6)\mathrm{mod}13$$

分量乘法：$(A,B)(C,D)=(7,10)(0,6)\equiv (0,60\mathrm{mod}13)=(0,8)\mathrm{mod}13$——只需 2 次实数乘法（CRNS 里要做 4 次）。反变换需要的逆元：解 $2x\equiv 1\mathrm{mod}13$ 得 $2^{-1}=7$（$2\times 7=14\equiv 1$）；解 $10x\equiv 1\mathrm{mod}13$ 得 $(2\hat{j}{)}^{-1}=10^{-1}=4$（$10\times 4=40\equiv 1$）。于是：

$$f^{-1}(0,8)=7(0+8)+j4(0-8)\mathrm{mod}13\equiv 56+j(-32)\equiv 4+j7\mathrm{mod}13$$

验证：$(2+j)(3+j2)=6+j4+j3+j^{2}2=4+j7$ ✓。CRNS 与 QRNS 之间的映射关系如图 2-6 所示。

图 2-6 CRNS↔QRNS 转换

2.2.3 浮点数

当定点的精度和动态范围都不够用时，浮点数登场。代价是速度更慢、硬件更复杂。微处理器普遍遵循 IEEE 浮点标准（单精度/双精度），而基于 FPGA 的系统常常自定义格式，以便精确匹配所需的动态范围（指数位宽）和精度（尾数位宽）。这类浮点运算模块可从 IP 提供商处获取，且已被 VHDL-2008 标准收录。

标准浮点数字由符号位 $s$、指数 $e$ 和规格化无符号尾数 $m$ 组成：

符号位 s	指数 e	无符号尾数 m

代数形式：

$$X=(-1{)}^s\times 1.m\times 2^{e-\text{偏差}}\text{\hspace{1em}(2-23)}$$

两个要点：这是有符号幅值格式（不是补码）；尾数中规格化的”隐藏的 1”（$1.m$ 的整数位）不占编码位，白赚一位精度。指数的偏差取：

$$\text{偏差}=2^{E-1}-1\text{\hspace{1em}(2-24)}$$

其中 $E$ 是指数位宽。回忆 2.2.1 节：偏差编码让”比较指数大小”退化为无符号比较——这正是这里采用偏差的原因。

例 2.13 (1,6,5) 浮点格式

格式：1 位符号 + 6 位指数 + 5 位尾数（不含隐藏的 1）。求 $9.25_{10}$ 的编码。

第一步，算偏差：$\text{偏差}=2^{6-1}-1=31$。

第二步，尾数规格化为 $1.m$ 形式：

$$9.25_{10}=1001.01_2=\underset{m}{\underbrace{1.00101}}\times 2^3$$

（小数部分 $0.25=0.01_2$，因为 $0.01_2=1/4$。）

第三步，加偏差得存储指数：$e=3+31=34_{10}=100010_2$。

最后拼装：

符号位 s	指数 e	无符号尾数 m
0	100010	00101

反向转换也很重要。给定编码：

符号位 s	指数 e	无符号尾数 m
1	011111	00000

符号位为 1 → 负数。把隐藏的 1 补回尾数，指数减去偏差：

$$-1.00000_2\times 2^{31-31}=-1.0_2\times 2^0=-1.0_{10}$$

记住口诀：定点→浮点，指数加偏差；浮点→定点，指数减偏差。

特殊数

IEEE 754-2008 还定义了一些特殊编码（表 2-7）：

• ±∞：指数全 1 且尾数全 0；
• ±0：指数全 0 且尾数全 0（因为有符号幅值格式，+0 和 −0 编码不同）；
• 非正规数（denormal）：指数全 0 但尾数非 0，允许尾数不带隐藏的 1（即小于 1.0），用来表示比 $2^{E_{min}}$ 还小的数；
• NaN（非数字）：指数全 1 且尾数非 0。NaN 在软件系统中能有效减少”异常事件”，典型来源包括：$\infty -\infty$、$0\times \infty$、$0/0$ 和 $\infty /\infty$、负数开平方。

表 2-7 754-1985 中 5 种编码类型和更新后的 754-2008 IEEE 二进制浮点数标准

符号位 s	指数 e	尾数 m	值
0/1	全零	全零	±0
0/1	全零	非零	非正规化: $(-1{)}^s{2}^{1-\text{偏差}}0.m$
0/1	$1<e<E_{max}$	m	正规化: $(-1{)}^s{2}^{e-\text{偏差}}1.m$
0/1	全 1	全零	±∞
-	全 1	非零	NaN

（FPGA 浮点实现通常不支持非正规数和 NaN，只保留基本编码。）

例 2.14 12 位浮点数与定点数的比较

继续用 (1,6,5) 格式。最大数（尾数全 1、指数最大）：

$$\pm 1.11111_2\times 2^{31}\approx \pm 4.23_{10}\times 10^9$$

最小正规化数：

$$\pm 1.0_2\times 2^{1-31}=\pm 2^{-30}\approx \pm 9.31_{10}\times 10^{-10}$$

若允许非正规数（指数 $e=0$，尾数 $0.00001$，即 $2^{-5}$）：

$$\pm 0.00001_2\times 2^{1-31}=\pm 2^{-36}\approx \pm 1.45_{10}\times 10^{-11}$$

比正规化最小数小 64 倍。

作为对比，取 12 位定点格式：1 位符号 + 5 位整数 + 6 位小数。最大值：

$$\pm 11111.111111_2=\pm (16+8+4+2+1+1/2+\cdots +1/64)=\pm (32-1/64)\approx \pm 32.0_{10}$$

最小非零值：

$$\pm 00000.000001_2=\pm 1/64_{10}=\pm 0.015625_{10}$$

结论：浮点动态范围约 $4\times 10^9$，定点只有约 32，差距悬殊。但精度上定点反而更细：$1.0$ 和 $1+1/64=1.015625$ 在 12 位定点中可以区分，而在 (1,6,5) 浮点中的编码相同（尾数只有 5 位小数，分辨不出 $1/64$）。所以”浮点更精确”是个常见误解——浮点赢在动态范围，同字长下定点在数值 1 附近的分辨率更高。

舍入方式

定点数有两种舍入方式，浮点数支持 4 种：

• 向最近的偶数舍入（默认，对应 MATLAB 的 round()）
• 向零舍入 / 截断（fix()）
• 向正无穷舍入 / 向上舍入（ceil()）
• 向负无穷舍入 / 向下舍入（floor()）

MATLAB 与 IEEE 754 的细微差别：当小数部分恰为 $0.5$（二进制 $0.1$）时，IEEE 只在整数 LSB 为 1 时向上舍入，否则向下舍入——所以 32.5 向下舍成 32，33.5 向上舍成 34。

表 2-8 关于 4 种浮点类型的舍入示例

方式	32.5	33.25	33.5	33.75	-32.5	-32.25
向最近的偶数舍入	32	33	34	34	-32	-32
向零舍入	32	33	33	33	-32	-32
向上舍入	33	34	34	34	-32	-32
向下舍入	32	33	33	33	-33	-33

一个有趣的观察：默认的就近舍入实现最复杂，而向零舍入开销最少，且舍入方向永远朝原点，不会放大信号的幅度（舍入不会导致增益增长）。因此 FPGA 设计中常刻意选用向零舍入。

IEEE 754 标准与 FPGA 定制格式

完整实现 IEEE 754 的所有细节（4 种舍入、非正规数、NaN）并不容易，但 1985 年的标准极大推动了浮点普及，如今已是微处理器的事实标准。表 2-9 给出各档格式的参数。

表 2-9 IEEE 浮点 754-2008 标准互换格式

	短	单精度	双精度	扩展
字长	16	32	64	128
尾数	10	23	52	112
指数	5	8	11	15
偏差	15	127	1023	16 383
范围	$2^{16}\approx 6.4\times 10^4$	$2^{128}\approx 3.8\times 10^{38}$	$2^{1024}\approx 1.8\times 10^{308}$	$2^{16383}\approx 10^{4932}$

实现单精度 754 算法至少需要：一个 24 位 × 24 位乘法器；以及能自定义指数/尾数位宽的 FPGA 资源。

不过 FPGA 设计往往不走 754 标准，而是定义专用格式。例如 Shirazi 等人在自定义计算设备 SPLASH-2（基于 Xilinx XC4010 的多 FPGA 开发板）上采用 18 位格式：10 位尾数 + 7 位指数 + 1 位符号位，可表示范围约 $3.7\times 10^{19}$。为什么选 18 位？因为系统总线宽 36 位，一次正好并行传两个操作数——格式为硬件平台服务，这正是 FPGA 定制精神的体现。

2.3 二进制加法器

加法是数字信号处理中最基础的操作——滤波器的累加、乘法器的部分积求和，最终都要落到加法器上。那么 FPGA 里一个 N 位加法器到底是怎么搭起来的？

一个基本的 N 位二进制加法器/减法器由 N 个全加器（Full Adder, FA）级联而成。每个全加器实现如下布尔方程，它定义了累加和的第 k 位：

$$s_k=x_k\oplus y_k\oplus c_k\text{\hspace{1em}(2-25/2-26)}$$

进位位则按”多数表决”的方式计算——三个输入中至少有两个为 1 时进位为 1：

$$c_{k+1}=\left(x_k\times y_k\right)+\left(x_k\times c_k\right)+\left(y_k\times c_k\right)\text{\hspace{1em}(2-27/2-28)}$$

对于二进制补码（2C）加法器，最低有效位可以简化为一个半加器，因为最低位的进位输入恒为 0——这也是为什么很多加法器设计中第 0 位可以省掉一根进位输入线。

结构选择与进位延迟。 最简单的结构叫”逐位进位加法器”（Ripple Carry Adder），如图 2-7(a) 所示，每个全加器的进位输出接到下一个全加器的进位输入，呈位串行形式。它的缺点很明显：最坏情况下进位要”涟漪”式地穿过所有 N 级，延迟与位宽成正比。如果在 FPGA 中较大的查找表（LUT）可用，可以把几个位组合进一个 LUT，如图 2-7(b) 所示的”一次两位”加法器，用面积换速度。

(a) 逐位进位加法器

(b) 几个位组合到一个 LUT 中的加法器
图 2-7 二进制补码加法器

为了缩短进位延迟，历史上发展出跳跃进位、先行进位、条件和、进位选择等许多技术。这些技术在老一代 FPGA（如 Xilinx 的 XC3000）上确实有用，因为那些器件没有内部快速进位逻辑。但现代系列（如 Xilinx 的 Spartan、Altera 的 Cyclone）本身就在芯片里布了专用的快速进位链，其延迟远小于经过常规 LUT 的路径——Altera 用快速查找表实现（参阅图 1-12），Xilinx 用硬连线译码器实现（图 2-8）。快速进位逻辑的出现，意味着我们不必再费力去手工搭先行进位结构，直接让综合工具用芯片自带的进位链即可。

图 2-8 XC4000 快速进位逻辑(© Xilinx ${}^{[53]}$ )

性能数据怎么看。 图 2-9 总结了用 lpm_add_sub 宏函数实现的 N 位加法器在几代 FPGA 上的规模和速度，包括 Cyclone IV E 的 EP4CE115F29C7（60nm 工艺）、Nios 开发板使用的 APEX20KE 系列 EP20K200EFC484-2X（1999 年，0.18μm 工艺），以及 UP2 开发板使用的 FLEX10K 系列 EPF10K70RC240-4（1995 年，0.42μm 工艺）。可以看到，虽然 LE 单元结构这些年变化不大，但工艺进步直接体现在速度上。另外两个经验结论值得记住：如果操作数要经过 I/O 寄存器单元引入，FPGA 总线延迟会成为瓶颈，性能下降；如果数据从本地寄存器传来，性能更好。FLEX10K 的加法器与寄存器不能合并，需要 4×N 个 LE；Cyclone IV 和 APEX 器件则只需 3×N 个 LE。

图 2-9 Cyclone IV、APEX 和 Flex10K 的加法器速度和规模

2.3.1 流水线加法器

DSP 算法的内部数据流非常规则（比如 FIR 滤波器每拍做一次乘累加），这决定了流水线技术在 DSP 实现中应用极广。典型的可编程 DSP 的 MAC 单元就至少带 4 条流水线：1) 对指令译码；2) 把操作数下载到寄存器；3) 执行乘法并存储乘积；4) 同时累加乘积。

在 FPGA 设计中应用流水线的成本极低，甚至为零——因为每个逻辑元件（LE）里本来就有触发器，不用也是浪费。流水线的思路是：把一个大算术操作拆成若干小规模基本操作，把进位和中间值存进寄存器，下个时钟周期继续算。文献中把这类加法器称为”进位保存加法器”（Carry Save Adder, CSA）。

那么问题来了：加法器该拆成几段？要不要拆到位级？答案取决于器件：对 Cyclone IV，一个 LAB 含 16 个 LE（即 16 个 FF），流水线单元按 16 位一组划分较合理；FLEX10K 每个 LAB 只有 8 个 LE，APEX20KE 是 10 个。所以动手前先查数据表。事实上，在 Cyclone IV 上实现 14 位流水线加法器并不会提速（见表 2-10），因为 14 位加法器连一个 LAB 都装不满，拆了反而多一级寄存开销。

表 2-10 采用带有流水线选项的预定义 LPM 模块合成的 EP2C35F672C6 的 14 位流水线加法器的性能

流水线级	MHz	LE
0	375.94	42
1	377.79	56
2	377.64	70
3	377.79	84
4	381.10	98
5	376.36	112

反过来，当位宽超过一个 LAB 时流水线才开始见效。由于每个 LAB 有 16 个触发器，且进位输出还需要 1 个额外触发器，为获得最大时序性能应采用最大 15 位的模块规模；只有最高有效位（MSB）模块不再需要为进位保留触发器，可以用满 16 位。由此可得：

1. 采用 1 级流水线，可构建 $15+16=31$ 位加法器；
2. 采用 2 级流水线，可构建 $15+15+16=46$ 位加法器；
3. 采用 3 级流水线，可构建 $15+15+15+16=61$ 位加法器。

表 2-11 的数据印证了这一点：增加流水线级数后，虽然位宽在增大，速度仍然保持得很高。

表 2-11 有/无流水线加法器的性能和资源要求。最大位宽为 31、46 和 61 位加法器的规模和速度

位宽	无流水线 MHz	无流水线 LE	有流水线 MHz	有流水线 LE	流水线数	设计文件名称
17~31	263.50	93	350.63	125	1	add1p.vhd
32~46	215.66	138	243.43	233	2	add2p.vhd
47~61	173.13	183	231.43	372	3	add3p.vhd

还有一种节省资源的技巧：直接实现流水线加法器时，MSB 段需要两个寄存器（图 2-10(a)）；如果 MSB 段不加法器、只靠寄存器传递（因为每个 LE 只能实现一个触发器，却可以实现一个全加器），就能省下一组 LE，如图 2-10(b) 所示。

(a) 直接实现

(b) FPGA 最佳的方法
图 2-10 流水线加法器

例 2.15：31 位流水线加法器的 Verilog 设计

原书 VHDL 改写为 Verilog。下面是图 2-10 所示 31 位流水线加法器的可综合实现，综合后速度 350.63MHz，使用 125 个 LE。设计分两段：低 15 位一段（留出进位位，共 16 位）、高 16 位一段，第一拍分别相加并寄存，第二拍把低段进位加到高段结果上。

// add1p.v : 31-bit pipelined adder
module add1p #(parameter WIDTH = 31,   // total bit width
               WIDTH1  = 15,           // bit width of LSBs
               WIDTH2  = 16)           // bit width of MSBs
             (input  [WIDTH-1:0] x, y, // inputs
              output [WIDTH-1:0] sum,  // result
              output LSBs_Carry,       // carry from LSB section
              input  clk);
 
  // LSBs of inputs
  reg [WIDTH1-1:0] l1, l2, s1;
  reg [WIDTH1:0]   r1;                // LSB sum incl. carry
  // MSBs of inputs
  reg [WIDTH2-1:0] l3, l4, r2, s2;
 
  always @(posedge clk) begin
    // Split LSBs / MSBs from inputs and store in registers
    l1 <= x[WIDTH1-1:0];
    l2 <= y[WIDTH1-1:0];
    l3 <= x[WIDTH-1:WIDTH1];
    l4 <= y[WIDTH-1:WIDTH1];
 
    // ---- First stage of the adder ----
    r1 <= {1'b0, l1} + {1'b0, l2};
    r2 <= l3 + l4;
 
    // ---- Second stage of the adder ----
    s1 <= r1[WIDTH1-1:0];
    // Add MSB result (x+y) and carry from LSBs
    s2 <= r2 + {{(WIDTH2-1){1'b0}}, r1[WIDTH1]};
  end
 
  assign LSBs_Carry = r1[WIDTH1];     // test signal
  // Build a single output word of WIDTH = WIDTH1 + WIDTH2
  assign sum = {s2, s1};              // connect to output pins
endmodule

这个流水线加法器的仿真结果如图 2-11 所示。注意一个细节：32780 + 32770 的低 15 位部分会产生进位，而 $32760+5=32765<2^{15}$，所以不产生进位——进位是否出现完全取决于低段相加的实际结果。

图 2-11 流水线加法器的仿真结果

2.3.2 模加法器

在余数系统（RNS-DSP）设计中，模加法器是最重要的构建模块：它既能做加法，也能通过索引算法用来做乘法。所谓模加法，就是计算 $(x+y)\mathrm{mod}{m}_l$——普通加法的结果如果超过模数 $m_l$，就要减去一次（或多次）$m_l$。

实现方案主要有两类（图 2-12）：

• 纯 LE 方案（图 2-12(a)）：直接用逻辑单元实现”加法 + 条件减模”，对任何 FPGA 都可行。
• ROM 查表方案（图 2-12(b)）：老一代 Altera FLEX 器件带有少量 M2K ROM/RAM（EAB），可配置成 $2^8\times 8$、$2^9\times 4$、$2^{10}\times 2$ 或 $2^{11}\times 1$ 的表，直接查表完成模 $m_l$ 校正。

表 2-12 给出了为 FLEX10K 编译的 6、7、8 位模加法器的规模和性能。可以看到：多路复用加法器（MPX-Add）速度相对较慢，即使每列都加上进位链也是如此；但它的流水线版本几乎不需要额外 LE，速度却是非流水线版本的约 3 倍，在 3 级流水线、6 位通道时达到最大吞吐量。ROM 方案速度高且稳定（86.2 MSPS），但有两个代价：一是 ROM 本身带来 4 个周期的流水线延迟，二是 ROM 数量有限，而且缩放模式必须用到 ROM。

表 2-12 为 Altera FLEX10K 器件编译的模加法器的规模和性能

类型	流水线级	6 位	7 位	8 位
MPX	0	41.3 MSPS / 27 LE	46.5 MSPS / 31 LE	33.7 MSPS / 35 LE
MPX	2	76.3 MSPS / 16 LE	62.5 MSPS / 18 LE	60.9 MSPS / 20 LE
MPX	3	151.5 MSPS / 27 LE	138.9 MSPS / 31 LE	123.5 MSPS / 35 LE
ROM	3	86.2 MSPS / 7 LE + 1 EAB	86.2 MSPS / 8 LE + 1 EAB	86.2 MSPS / 9 LE + 2 EAB

图 2-12 模加法

2.4 二进制乘法器

乘法比加法复杂，但思路可以很朴素。设被乘数为 $X$，乘数 $A={\sum }_{k=0}^{N-1}{a}_k{2}^k$，按”手动计算竖式”的方法：

$$P=A\times X=\sum _{k=0}^{N-1}{a}_k{2}^{k}X\text{\hspace{1em}(2-29)}$$

也就是说，逐位扫描 $A$：若 $a_k\ne 0$，就把 $X$ 左移 $k$ 位后累加；若 $a_k=0$，这一步就是空操作（nop）。操作数 $X$ 以并行形式提供，而 $A$ 是逐位扫描的，所以这种结构叫串行/并行乘法器。如果两个操作数都是串行输入，就叫串行/串行乘法器——它只需要一个全加器，硬件最省，但等待时间是高阶无穷大 $\mathcal{O}(N^2)$，因为状态机大约需要 $N^2$ 个周期。随着 FPGA 内嵌嵌入式乘法器的普及，这种基于状态机的方案已经不常用了。

阵列乘法器。 用复杂性换速度，就得到”阵列”（并行/并行）乘法器：两个操作数并行提交给 $N^2$ 个加法器单元组成的阵列。图 2-13 给出了一个 4 位 × 4 位的经典阵列结构。

图 2-13 4 位阵列乘法器

不过经典阵列是照搬手算进位习惯设计的。现代 FPGA 中进位计算比和累加快得多，更高效的结构如图 2-14 所示（8 位 × 8 位乘法器）：在第一级就把相邻的两个部分积 $a_{n}X{2}^n$ 和 $a_{n+1}X{2}^{n+1}$ 合并，再把结果加到最终乘积上。因为它本质上是手算方法的直接阵列形式，正确性是必然的。

图 2-14 FPGA 的快速阵列乘法器

这种两两合并的结构天然形成（并行）二叉树乘法器，其流水线级数为：

$$\text{二叉树乘法器的流水线级数}={\log }_2(N)\text{\hspace{1em}(2-30)}$$

树的每一层之后都便于插入流水线寄存器。达到最大吞吐量所需的流水线级数见表 2-13。

表 2-13 达到最大吞吐量所需流水线级数

位宽	2	3 和 4	5~8	9~16	17~32
最佳流水线级数	1	2	3	4	5

注意：由于数据在输入端和输出端都打了寄存器，仿真中观察到的延迟要比为 lpm_mul 设定的流水线级数多两级。

嵌入式乘法器还是 LE 乘法器？ 图 2-15 给出了 Quartus II 中 lpm_mult 流水线乘法器从 8×8 到 24×24 的性能：虚线是嵌入式乘法器（相当于 16×16 乘法器装入一个嵌入式 18×18 阵列，再加深流水线没有收益），实线是基于 LE 的乘法器。图 2-16 给出了按 LE 计算的资源消耗。经验法则是：若使用两级或更多流水线，8×8 流水线乘法器的性能可超过嵌入式乘法器；但流水线级数远大于 ${\log }_2(N)$ 时，LE 乘法器没有明显提升。如果写的是行为级代码（如 p <= a*b 对应 Verilog 的 assign p = a * b;），就要靠综合选项控制乘法器架构：在 Assignments 菜单 Settings 的 Analysis & Synthesis Settings 下找 DSP Block Balance——选 DSP blocks 用嵌入式乘法器，选 Logic Elements 用 LE，选 Auto 则优先嵌入式、不够再用 LE 补充。若直接用 lpm_mul 模块，则通过 GENERIC 参数 DEDICATED_MULTIPLIER_CIRCUITRY 设为 “YES” 或 “NO” 直接控制。

图 2-15 FPGA 阵列乘法器的性能(实线表示基于 LE 的乘法器，虚线表示嵌入式乘法器)

图 2-16 阵列乘法器在 LE 方面的工作量(实线表示使用基于 LE 的乘法器，虚线表示嵌入式乘法器)

Booth 乘法器和 Wallace 树乘法器在 ASIC 领域很常见（将在练习 2.1 和 2.2 中讨论），但在 FPGA 相关设计中如今已很少使用。

大位宽乘法器的分块技术

当器件容量不够直接实现所需规模的乘法器，或者想用存储器模块实现乘法器时，可以把 $2N\times 2N$ 乘法器拆成 $N\times N$ 模块的组合。把每个操作数分成高低两半（下标 2 为高 N 位，下标 1 为低 N 位）：

$$P=Y\times X=(Y_2{2}^N+Y_1)(X_2{2}^N+X_1)\text{\hspace{1em}(2-31)}$$ $$=Y_2{X}_2{2}^{2N}+(Y_2{X}_1+Y_1{X}_2)2^N+Y_1{X}_1\text{\hspace{1em}(2-32)}$$

这其实就是”竖式乘法”的分块版：4 个小乘法加 3 个加法。36×36 乘法器可以由 4 个 18×18 嵌入式乘法器加 3 个加法器构成。分块对基于 LUT 的乘法器尤其重要：8×8 的直接查表形式需要 $2^{16}\times 16=1\text{Mbit}$ 的 LUT，而分块后只需 4 个 $2^8\times 8$ 的存储模块加 3 个加法器；16×16 则需要 16 个存储模块。用 M9K 实现乘法器相比基于 LE 的实现有双重优点：LE 用量降低，且对布线资源的要求也降低。另外，有些器件系列（如 Cyclone、Flex、Excalibur）根本没有嵌入式乘法器，LUT/LE 乘法器就是唯一选择。

1. 半分方形乘法器

减少 LUT 存储量的另一条路是减少输入域的位：输入每少一位，LUT 字数就减半。$N$ 位数的平方表只需 $2^N\times 2^N$ 规模。Logan 提出的加性半分方形乘法器（Additive Half-Square Multiplier, AHSM）正是利用了恒等式：

$$Y\times X=\frac{(X+Y{)}^2-X^2-Y^2}{2}=\lfloor \frac{(X+Y{)}^2}{2}\rfloor -\lfloor \frac{X^2}{2}\rfloor -\lfloor \frac{Y^2}{2}\rfloor -\{\begin{array}{ll}1& X,Y\text{均为奇数}\\ 0& \text{其他}\end{array}\text{\hspace{1em}(2-33)}$$

也就是说，一次乘法可以化为三个平方表的查表加减。如果在平方表中预先包含除以 2 的运算，那么当 X 和 Y 均为奇数时需要减 1 校正（因为 $(2u+1)(2v+1)$ 的展开在整除截断后会差 1）。

对应的差分半分方形乘法器（Differential Half-Square Multiplier, DHSM）为：

$$Y\times X=\lfloor \frac{X^2}{2}\rfloor +\lfloor \frac{Y^2}{2}\rfloor -\lfloor \frac{(X-Y{)}^2}{2}\rfloor +\{\begin{array}{ll}1& X,Y\text{均为奇数}\\ 0& \text{其他}\end{array}\text{\hspace{1em}(2-34)}$$

它在 X、Y 均为奇数时需要加 1 校正。如果是有符号数，配合 2.2.1 节的减 1（D1）编码还能进一步省资源：D1 编码中所有数字都减 1，对 0 特殊编码。图 2-17 给出了 8 位 AHSM 乘法器所需的 LUT 与 8 位输入操作数的数据范围：绝对值运算把 LUT 字数减半，D1 编码又让表规模以接近 2 的幂的速度下降，对 FPGA 设计非常有利。由于输入 0 和 1 的平方等于原值，LUT 项 $\lfloor A^2/2\rfloor$ 可以共享而无需为 0 特殊编码。不做除以 2 的话输出需要 17 位字；做了除以 2 之后，仅在两个操作数都是奇数时需要递增（AHSM）或递减（DHSM）输出。图 2-18 给出了只需两个 D1 编码表的 DHSM 乘法器。

图 2-17 2 的补码的 8 位加性半分方形乘法器设计

图 2-18 2 的补码的 8 位差分半分方形乘法器设计

2. 四分之一平方乘法器

LUT 数量还能再降——用四分之一平方乘法（Quarter-Square Multiplication, QSM），这种方法在模拟设计中已被深入研究过：

$$Y\times X=\lfloor \frac{(X+Y{)}^2}{4}\rfloor -\lfloor \frac{(X-Y{)}^2}{4}\rfloor \text{\hspace{1em}(2-35)}$$

妙处在于式(2-35)除以 4 不像 HSM 那样需要奇偶校正。验证一下：若两操作数同为偶数（或同为奇数），它们的和与差都是偶数，平方后除以 4 无误差（即 $4\mid (2u\cdot 2v)$ 型）；若一奇一偶，和与差的平方除以 4 各产生 0.25 的误差，但两者正好相互抵消，所以不需要校正。直接实现需要 $(N+1)$ 位输入的 LUT 来表示正确的 $X\pm Y$ 结果，共 4 个 $2^N\times 2^N$ 的 LUT；配合有符号算法与 D1 编码后，只需两个 $2^N\times 2^N$ 的 LUT，图 2-19 给出了 D1 QSM 电路。AHSM、DHSM 和 QSM 的综合结果可参阅文献[34]。

图 2-19 2 的补码的 8 位四分之一平方乘法器设计

2.5 二进制除法器

四种基本算术运算里，除法是最复杂的：最耗时，且可用的算法种类也最多。给定被除数（分子）N 和除数（分母）D，除法会同时给出两个结果（这是它与其他运算不同的地方）——商 Q 和余数 R：

$$\frac{N}{D}=Q\text{ 和 }R,\text{其中 }|R|<D\text{\hspace{1em}(2-36)}$$

也可以把除法看成乘法的逆运算：

$$N=D\times Q+R\text{\hspace{1em}(2-37)}$$

除法与乘法的关键区别：乘法中所有部分乘积都可以并行生成，而除法中商的每一位必须按顺序、以”尝试错误”的方式逐位确定——这就是除法慢的根本原因。

位宽约定。 大多数微处理器把除法当作乘法的逆过程（假定分子是某次乘法的结果），因此分母和商的位宽要扩大一倍，代价是必须用笨拙的流程检查商是否溢出。我们采用更通用的约定：

$$Q⩽N\text{ 并且 }|R|⩽D$$

即商与分子位宽相同、余数与分母位宽相同。这样除 N = 0 外就完全不需要检验商的范围了。

符号处理。 最简单的有符号数处理方法：先把分子分母都转成无符号数，结果的符号由两个操作数符号位”异或”（即模 2 加）得到。不过某些算法（下面要讲的非还原除法）可以直接处理有符号数。这时要约定商和余数的符号关系：大多数硬件和软件系统（但不全是，比如 PASCAL 不是）假定余数与商同号。也就是说，虽然

$$\frac{234}{50}=5,\text{而 }R=-16\text{\hspace{1em}(2-38)}$$

满足式(2-37)，但通常更倾向于

$$\frac{234}{50}=4,\text{而 }R=34\text{\hspace{1em}(2-39)}$$

算法分类。 图 2-20 汇总了最常用的线性收敛与二次收敛方案。线性算法可以按商的每位数字的取值集合来分类：二进制还原、非执行或 CORDIC 算法用 $\{0,1\}$；二进制非还原算法用有符号数字集合 $\{-1,1\}=\{\overline{1},1\}$；二进制 SRT 算法（以差不多同时独立发现它的 Sweeney、Robertson 和 Tocher 命名）用三重集合 $\{-1,0,1\}=\{\overline{1},0,1\}$。这些算法都能扩展到更高基数，例如基数 r 的广义 SRT 算法使用数字集合 $\{-2^r-1,\dots ,-1,0,1,\dots ,2^r-1\}$。

二次收敛的流行算法有两种：第一种是”分母互换”除法，用牛顿法求倒数；第二种是 20 世纪 60 年代 Anderson 等人为 IBM 360/91 开发的收敛除法，用同一因子分别乘分子和分母，使 $N\to 1$、$D\to Q$。注意二次收敛算法不生成余数。虽然二次算法的迭代次数只有 ${\log }_2(b)$（b 为位宽），但每次迭代更复杂（要用两次乘法），速度与规模必须仔细权衡。

图 2-20 除法算法综述

2.5.1 线性收敛的除法算法

最直观的顺序算法就是把手算竖式变成二进制版本，即还原除法：先把分母对齐、分子装入余数寄存器；从余数中减去对齐后的分母，若新余数为正，商的该位置 1；若为负，就要把余数加回分母”还原”，商该位置 0；然后分母右移一位继续。“还原”之名正来自这一步。

例 2.16：8 位还原除法器

原书 VHDL 改写为 Verilog。下面除法分 4 个阶段：复位后（ini）把 8 位分子装入余数寄存器、对齐分母（N 位分子对应 $2^{N-1}$ 的对齐系数）、商寄存器清零；sub 和 restore 两个状态完成实际的串行除法；最后在 done 状态把商和余数送到输出寄存器。假定分子和商 8 位宽，分母和余数 6 位宽。

// div_res.v : Restoring Division
module div_res #(parameter WN = 8,          // numerator width
                 WD = 6,                    // denominator width
                 PO2WN1 = 128,              // 2**(WN-1)
                 PO2WN  = 255)              // 2**WN - 1
               (input            clk,       // system clock
                input            reset,     // asynchronous reset
                input      [7:0] n_in,      // numerator
                input      [5:0] d_in,      // denominator
                output reg [5:0] r_out,     // remainder
                output reg [7:0] q_out);    // quotient
 
  // bit width: WN    WD    WN    WD
  //   Numerator / Denominator = Quotient and Remainder
  // or: Numerator = Quotient * Denominator + Remainder
 
  localparam S14_MAX = 2*(PO2WN1*PO2WN);     // N+D bit range
  // divider in behavioral style
  reg [1:0] state;
  localparam ini = 2'd0, sub = 2'd1,
             restore = 2'd2, done = 2'd3;
 
  reg signed [13:0] r, d;                    // N+D bit width
  reg [7:0]  q;
  reg [2:0]  count;
 
  always @(posedge clk or posedge reset) begin
    if (reset) begin                         // asynchronous reset
      state <= ini; q_out <= 8'd0; r_out <= 6'd0;
    end else begin
      case (state)
        ini: begin                           // initialization step
          state <= sub;
          count <= 0;
          q <= 0;                            // reset quotient register
          d <= PO2WN1 * d_in;                // load denominator
          r <= n_in;                         // remainder = numerator
        end
        sub: begin                           // processing step
          r <= r - d;                        // subtract denominator
          state <= restore;
        end
        restore: begin                       // restoring step
          if (r < 0) begin
            r <= r + d;                      // restore previous remainder
            q <= {q[6:0], 1'b0};             // LSB = 0 and SLL
          end else begin
            q <= {q[6:0], 1'b1};             // LSB = 1 and SLL
          end
          count <= count + 1;
          d <= d >>> 1;                      // realign denominator
          if (count == WN) state <= done;    // division ready?
          else             state <= sub;
        end
        done: begin                          // output of result
          q_out <= q;
          r_out <= r[WD-1:0];
          state <= ini;                      // start next division
        end
      endcase
    end
  end
endmodule

图 2-21 给出了 234 除以 50 的仿真结果。寄存器 d 显示了对齐的分母值：$50\times 2^7=6400$、$50\times 2^6=3200$，依此类推。每次在 sub 步骤算出的余数为负，就在 restore 步骤还原前一次的余数；在 done 状态，商 4 和余数 34 被送到输出寄存器。这一设计用了 106 个 LE（未用嵌入式乘法器），TimeQuest 慢速 85C 模式下时序性能 265.32MHz。

图 2-21 还原除法器的仿真结果

改进一：非执行除法

还原除法的主要缺点是确定商的一位要两个步骤（先减、再还原）。非执行除法把两步合并：先算出暂存结果，若分母大于余数就干脆”不执行”减法。

原书 VHDL 改写为 Verilog，restore 步骤写成：

t = r - d;            // temporary remainder value
if (t >= 0) begin     // nonperforming test
  r = t;              // use new remainder
  q = {q[6:0], 1'b1}; // LSB = 1 and SLL
end else begin
  q = {q[6:0], 1'b0}; // LSB = 0 and SLL
end

从图 2-22 的仿真结果看，步骤数量减少了一半（初始化与结果传输不计），且余数 r 在整个过程中永远非负。但天下没有免费午餐：非执行除法在最坏延迟路径上有 1 个 if 条件加 2 个算术运算（先试减再还原的分支都在同一路径上），而还原除法最坏路径只有 1 个 if 条件加 1 个算术运算，所以最大时序性能可能反而降低（参阅练习 2.17）。

图 2-22 非执行除法器的仿真结果

改进二：非还原除法

还有所谓的非还原除法，与算法类似但不增加关键路径。核心观察是：还原除法中若已算得负余数 $r_{k+1}=r_k-d_k$，下一步要先加 $d_k$ 还原 $r_k$，再减去对齐后的新分母 $d_{k+1}=d_k/2$。既然”加 $d_k$ 再减 $d_k/2$“等价于”净加 $d_k/2$“，干脆在余数为（暂时）负时跳过还原、直接加 $d_k/2$ 即可。

代价是商位现在可正可负：$q_k=\pm 1$（没有 0）。以后再把这种有符号数字（SD）表示变换回 2 的补码。规则总结：余数经迭代后为正，商记 1 并减对齐分母；余数为负，商记 $-1=\overline{1}$ 并加对齐分母。为了让商寄存器仍只用一位，可用 0 来编码 $-1$。把有符号商转回 2 的补码最直接的方法：把所有 1 收进一个字，把所有 0（实际是 $\overline{1}$ 的编码）收进第二个字，两字相减即可；而这些 $\overline{1}$ 的减法不过是”取补码再加 1”。总之，若 q 保存有符号数字表示，2 的补码结果为：

$$q_{2\mathrm{C}}=2\times q_{\mathrm{SD}}+1\text{\hspace{1em}(2-40)}$$

此后商和余数都是 2 的补码形式，由式(2-37)得到有效结果。如果还要求商和余数同号（前述约定），需对负余数做一次校正：当 r < 0 时执行

$$r:=r+D,q:=q-1$$

非还原除法器比非执行除法器跑得快，时序性能与还原除法器相当（参阅练习 2.18）。图 2-23 显示了非还原除法器的仿真结果：余数寄存器值允许为负，且可以看到负余数校正的必要性——未校正时 $q=5,r=-16$，校正后 $q=5-1=4$，$r=-16+50=34$，正是图 2-23 所示的值。

图 2-23 非还原除法器的仿真结果

要进一步缩短除法所需时钟周期数，可以利用 SRT 和基数 4 编码构造更高基数（阵列）除法器；与进位保存加法器结合后，这种方案在 ASIC 设计中非常流行，奔腾微处理器的浮点加速器就采用了这一原理。但对 FPGA 而言，由于 LUT 规模有限，高阶基数方案吸引力不大。

另一种完全不同的提速思路是二次收敛的除法算法——借用快速阵列乘法器。下一节讨论两种最流行的二次收敛方案。

2.5.2 快速除法器的设计

第一种快速除法器是”先求倒数再相乘”。位宽小时倒数可以直接查表；更通用的做法是用牛顿法迭代求零。定义函数：

$$f(x)=\frac{1}{x}-D\to 0\text{\hspace{1em}(2-41)}$$

令 $f(x_{\infty })=0$，立即得到：

$$\frac{1}{x_{\infty }}-D=0\text{或}{x}_{\infty }=\frac{1}{D}\text{\hspace{1em}(2-42)}$$

用切线法（以直代曲）不断更新估计值：

$$x_{k+1}=x_k-\frac{f(x_k)}{f^{\prime }(x_k)}\text{\hspace{1em}(2-43)}$$

由 $f(x)=1/x-D$ 得 $f^{\prime }(x)=-1/x^2$，迭代方程化为非常简洁的形式：

$$x_{k+1}=x_k-\frac{\frac{1}{x_k}-D}{\frac{-1}{x_k^2}}=x_k(2-D\times x_k)\text{\hspace{1em}(2-44)}$$

这个算法对任意初始值都收敛，但如果先把 D 规格化到接近 1（如浮点尾数那样取 $0.5⩽D<1$ 或 $1⩽D<2$，参阅 2.7 节），再用初始值 $x_0=1$，收敛会非常快。

例 2.17：牛顿算法

计算 $1/D=1/0.8=1.25$。表 2-14 中：第 1 列为迭代次数，第 2 列是对 $1/D$ 的近似，第 3 列是误差 $x_k-x_{\infty }$，最后一列是近似值的等价位精度。

表 2-14 例 2.17 的计算数据

k	$x_k$	$x_k-x_{\infty }$	有效位
0	1.0	-0.25	2
1	1.2	-0.05	4.3
2	1.248	-0.002	8.9
3	1.25	$-3.2\times 10^{-6}$	18.2
4	1.25	$-8.2\times 10^{-12}$	36.8

图 2-24 给出了牛顿寻零算法的图形解释——$f(x_k)$ 迅速收敛到 0。

图 2-24 $x_{\infty }=1/0.8=1.25$ 时的牛顿寻找零定位算法

这种”越算越准”并非本例巧合。设误差 $e_k=x_k-1/D$，代入式(2-44)：

$$e_{k+1}=x_{k+1}-x_{\infty }=x_k(2-D{x}_k)-\frac{1}{D}=-D{\left(x_k-\frac{1}{D}\right)}^2=-D{e}_k^2$$

误差每经一次迭代就”平方”一次，有效位精度翻倍——这就是”二次收敛”的含义。既然最初几步精度很低，可以用一个小查询表跳过最开始的迭代（参考文献[33]提供了跳过前两次迭代的表）。

牛顿算法已成功用于微处理器设计（如 IBM RISC 6000），但有两个缺点：一是每次迭代中的两次乘法必须顺序进行（第二个乘法依赖第一个的结果）；二是乘法的顺序本质导致量化误差累积，通常需要额外的保护位。

收敛除法（Anderson-Earle-Goldschmidt-Powers 算法）

下一种算法与牛顿算法类似，但改善了量化行为，且每次迭代中的两次乘法可以并行计算。其思想是：分子 N 和分母 D 同时乘以近似因子 $f_k$，迭代足够多次后：

$$D\prod f_k\to 1\text{ 和 }N\prod f_k\to Q\text{\hspace{1em}(2-45)}$$

它因 IBM 360/91 而闻名，归功于 Anderson 等人。算法流程如下：

算法 2.18：收敛除法

(1) 规格化 N 和 D，令 D 接近于 1，用浮点尾数那样的规格化区间 $0.5⩽D<1$ 或 $1⩽D<2$。

(2) 初始化 $x_0=N$，$t_0=D$。

(3) 重复以下循环，直到 $x_k$ 满足精度要求：

$$f_k=2-t_k,x_{k+1}=x_k\times f_k,t_{k+1}=t_k\times f_k$$

关键在于这个算法是自校正的：因子 $f_k$ 中的任何量化误差都无关紧要，因为分子和分母乘的是同一个因子。IBM 360/91 的设计就利用了这一点来省资源——第一次迭代只用几位有效位的乘法器，随着 $f_k$ 越来越接近 1，后面的迭代才逐步分配更多乘法器位。

例 2.19：Anderson-Earle-Goldschmidt-Powers 算法

计算 $N=1.5$、$D=1.2$，即 $Q=N/D=1.25$。表 2-15：第 1 列为迭代次数，第 2 列为比例因子 $f_k$（括号内为 1.8 位定点格式的值），第 3 列为 $N/D$ 的近似，第 4 列为误差 $x_k-x_{\infty }$，最后一列为等价位精度。

表 2-15 例 2.19 的计算数据

k	$f_k$	$x_k$	$x_k-x_{\infty }$	有效位
0	0.8（205/256）	1.5（384/256）	0.25	2
1	1.04（267/256）	1.2（307/256）	-0.05	4.3
2	1.0016（257/256）	1.248（320/256）	0.002	8.9
3	$1.0+2.56\times 10^{-6}$	1.25	$-3.2\times 10^{-6}$	18.2
4	$1.0+6.55\times 10^{-12}$	1.25	$-8.2\times 10^{-12}$	36.8

可见与例 2.17 的牛顿算法收敛轨迹一致，都是二次收敛。

例 2.19（续）：8 位收敛除法器的 Verilog 设计

原书 VHDL 改写为 Verilog。假定分子分母均已规格化为 $1⩽N,D<2$（典型浮点尾数值；若未规格化，需要额外的加法资源做前导零检测和两个桶式移位）。分子、分母和商均 9 位宽（1 位整数位 + 8 位小数位，即 1.8 格式）：$1.5\times 256=384$，$1.2\times 256=307$，$1.25\times 256=320$。除法分 3 个阶段：先在 ini 状态加载 1.8 格式的分子分母，run 状态做实际收敛迭代，done 状态输出商。

// div_aegp.v : Convergence division after Anderson, Earle,
//               Goldschmidt and Powers
module div_aegp #(parameter WN    = 9,     // 8 bit plus one integer bit
                  WD    = 9,
                  STEPS = 2,              // number of iterations
                  TWO   = 512,            // 2**(WN)
                  PO2WN = 256)            // 2**(WN-1)
                (input             clk,   // system clock
                 input             reset, // asynchronous reset
                 input       [8:0] n_in,  // numerator
                 input       [8:0] d_in,  // denominator
                 output reg  [8:0] q_out);// quotient
 
  // bit width: WN    WD    WN    WD
  //   Numerator / Denominator = Quotient (no remainder!)
 
  reg [1:0] state;
  localparam ini = 2'd0, run = 2'd1, done = 2'd2;
 
  reg [9:0] x, t;                          // WN+1 bits
  reg [9:0] f;                             // scale factor
  reg [1:0] count;
 
  always @(posedge clk or posedge reset) begin
    if (reset) begin                       // asynchronous reset
      state <= ini; q_out <= 9'd0;
    end else begin
      case (state)
        ini: begin                         // initialization step
          state <= run;
          count <= 0;
          t <= d_in;                       // load denominator
          x <= n_in;                       // load numerator
        end
        run: begin                         // processing step
          f <= TWO - t;                    // f_k = 2 - t_k
          x <= x * f >> 8;                 // x_k+1 = x_k * f_k (1.8 fmt)
          t <= t * f >> 8;                 // t_k+1 = t_k * f_k
          count <= count + 1;
          if (count == STEPS) state <= done; // division ready?
          else                state <= run;
        end
        done: begin                        // output of results
          q_out <= x[WN-1:0];
          state <= ini;                    // start next division
        end
      endcase
    end
  end
endmodule

图 2-25 给出了 1.5/1.2 的仿真结果。内部变量 f（作为中间网络，不出现在波形里）依次保存 3 个比例因子 205、267、257，对 8 位精度的结果已经足够。x 和 t 分别乘以比例因子 f（按 1.8 格式缩放），如预期那样，x 收敛到商 $1.25=320/256$，t 收敛到 $1.0=255/256$；done 状态把商 $1.25=320/256$ 送到输出寄存器。注意该除法器不生成余数。此设计用了 45 个 LE 和 4 个嵌入式乘法器，时序性能 124.91MHz。

图 2-25 收敛除法器的仿真结果

最后把两类除法器放在一起比较：虽然收敛除法的单步时序性能大约只有非还原除法的一半，但它的处理步骤数从 8 步降到 $\lceil \sqrt{8}\rceil =3$ 步（两种算法均不计初始化），总执行时间反而更短；而且收敛除法器用的 LE 与非还原除法器一样少——只是需要 4 个嵌入式乘法器作为代价。

2.5.3 阵列除法器

乘法器可以做成阵列或流水线结构，除法器也一样：所有除法算法既可以写成顺序的 FSM 状态机形式，也可以做成组合的”阵列”形式。对于初学者来说，关键问题是：自己需要写除法电路吗？ 答案通常是”不需要”。现代 FPGA 综合工具已经能把行为级的除法表达式直接推断成除法电路；如果对吞吐率有要求，则应直接调用厂商的流水线除法器宏模块（Altera 的 lpm_divide），它会自动生成阵列结构并可选择多级流水线。

先看最简单的写法。原书 VHDL 改写为 Verilog 后，行为级除法可以这样表达：

// 原书 VHDL 改写为 Verilog：行为级整数除法，综合工具可自动推断
module behavior_div
  #(parameter W = 8)               // 例如 8 位有符号数，范围 -128 ~ 127
  (input  wire signed [W-1:0] n,   // 被除数 n
   input  wire signed [W-1:0] d,   // 除数   d
   output wire signed [W-1:0] q);  // 商     q
  assign q = n / d;
endmodule

这对应原书中 q <= n / d; 的行为描述：不需要流水线时，Quartus II 12.1 及同类工具都能把它综合成组合除法电路。如果需要阵列形式加流水线，就改用 lpm_divide 宏。

性能和代价如何？图 2-26 用 TimeQuest 的慢速 85C 模型给出了时序性能，图 2-27 给出了 8×8、16×16、24×24 位阵列除法器（含 4 个 I/O 寄存器）所需的 LE 数量。

图 2-26 使用 lpm_divide 宏模块的阵列除法器的性能，行为代码对应零流水线级数

图 2-27 使用 lpm_divide 宏模块的阵列除法器按 LE 计算的工作量

从图 2-27 可以看到一个有趣的规律：流水线级数（取对数后）呈阶梯状变化，而实测结果表明，流水线的最优级数恰好等于分母的位数。设计时可以直接参考这条经验规则。

2.6 定点算法的实现

在 DSP 里我们处理的更多是有符号数。VHDL-2008 为此增加了有符号定点类型 sfixed（还有无符号 ufixed 和浮点 float），它自带一套完整的运算包（由 David Bishop 编写、超过 8000 行代码，可从 www.eda.org/fphdl 免费获得，也兼容 VHDL-1993）。由于原书示例使用 VHDL 的 sfixed 库，而 Verilog 没有对应的语言级定点库，下面把原书的设计思想改写为等价的、手工管理位宽的 Verilog 描述——这也正是 Verilog 使用者每天要做的事。

sfixed(3 downto -4) 表示 4 位整数位（索引 3…0）加 4 位小数位（索引 -1…-4），即 8.4 格式中取 4.4。以 $a=3.625$ 为例，其编码为：

位值	3	2	1	0	-1	-2	-3	-4
权值	8	4	2	1	0.5	0.25	0.125	0.0625
编码	0	0	1	1	1	0	1	0

即 $8\times 0+4\times 0+2\times 1+1\times 1+0.5+0.125=3.625$。

sfixed 库的一个关键设计哲学是：加法的结果在左边多出 1 个整数位（保护位），以把溢出可能性降到最低。原书 VHDL 改写为 Verilog 后，这个规则可以这样演示：

// 原书 VHDL 改写为 Verilog：4.4 格式有符号定点加法的位宽规则
module fixed_add_demo
  (input  wire signed [7:0] a,      // 4.4 格式，如 3.625 -> 8'sb0011_1010
   input  wire signed [7:0] b,      // 4.4 格式
   output wire signed [8:0] s,      // 全精度结果：5.4 格式，9 位（9 个 LE）
   output wire signed [7:0] r);     // 直接截回 4.4：等价于 wrap+truncate（8 个 LE）
  assign s = a + b;                 // 正确写法：结果比操作数多 1 个整数位
  assign r = a + b;                 // 位宽不变 = 自动丢弃最高进位（截断/回绕）
endmodule

对照原书的几条注释：s <= a+b 是”合法编码”，需要 9 个 LE；r <= a+b 在 sfixed 库里直接报错（表达式 9 位 vs 目标 8 位）；若想保持 4.4 输出又不想手动截位，库提供 RESIZE。用默认的饱和+舍入选项需要 16 个 LE，而显式指定 fixed_wrap 和 fixed_truncate（回绕+截断）后只需 8 个 LE——饱和与舍入是有代价的，最省资源的是回绕、截断、0 保护位的组合。

Verilog 中没有 RESIZE 函数，等价操作就是显式的位选择与拼接，例如饱和处理可写成 assign r_sat = (s > 9'sd127) ? 8'sd127 : (s < -9'sd128) ? -8'sd128 : s[7:0];。

设两种数值格式的范围分别为 $A_L\cdots A_R$ 和 $B_L\cdots B_R$，表 2-16 给出了典型 sfixed 操作的结果位宽范围，这是手工设计 Verilog 位宽时同样适用的通用规则：

表 2-16 典型 sfixed 操作的范围列表

操作	结果范围	相同范围 $A_L=B_L=L;A_R=B_R=R$
A+B、A−B	$\max (A_L,B_L)+1\cdots \min (A_R,B_R)$	$L+1\cdots R$
abs(A)、−A	$A_L+1\cdots A_R$	$L+1\cdots R$
A×B	$A_L+B_L+1\cdots A_R+B_R$	$2L+1\cdots 2R$
A/B	$A_L-B_R+1\cdots A_R-B_R$	$L-R+1\cdots R-L$
1/A	$A_L+1\cdots -A_L-1$	$-R\cdots -L-1$

为什么要在乎这些规则？ 如果不做定点运算设计，要么位宽爆炸浪费 LE，要么溢出静默地破坏数据。记住表 2-16，加上”回绕+截断最省资源”这条经验，就掌握了定点算法设计的核心。

2.7 浮点算法的实现

如今 FPGA 门数充足，加上 Stratix/Cyclone、Virtex/Spartan 系列内嵌了 18×18 位阵列乘法器，浮点算法在 FPGA 上已经成为可行选择。本节讨论浮点加法器、减法器、乘法器、倒数器、除法器以及定点↔浮点转换这些基本模块。商用的浮点 IP 通常采用 5 级以上流水线提高吞吐量；为了看清原理，下面的讲解不用流水线，并采用 2.2.3 节介绍的自定义 (1,6,5) 浮点格式：1 个符号位、6 个指数位、5 个尾数位。我们支持 0 和无穷大 ∞ 的特殊编码，但不支持 NaN 和非规格数；舍入采用截断；配套的定点格式为 6 位整数位（含符号）加 6 位小数位，即 (1,5,6)。

2.7.1 定点数到浮点数的格式转换

浮点数采用”符号-幅值”格式，而定点数是 2 的补码，所以转换的第一步是格式重排：若定点数符号位为 1，就先求补码，得到非规格化的尾数。第二步是规格化并计算指数：数出前导零的个数（原书中用顺序 PROCESS 内的 LOOP 完成；Verilog 中可用一个 for 循环综合出优先级链），然后把尾数左移直到第一个 1”移出”寄存器——隐藏的 1 也随之去掉。这个移位实际上是桶式移位器的任务，在 VHDL 中通过 SLL 推断；1076-1993 的 SLL 只对 BIT_VECTOR 定义，对 STD_LOGIC_VECTOR 需要自己写函数重载或按练习 2.19 的方式搭桶式移位器（Verilog 的 << 运算符没有这个限制）。

指数由”偏差 + 整数位个数 − 前导零个数”算出。最后把符号、指数、规格化尾数拼接成浮点数；若定点数为零，浮点数也置零。只要浮点动态范围大于定点范围，转换中就不会出现 ∞。

图 2-28 演示了 5 个代表值（+1、−1、最大绝对值、最小绝对值、最小值）从 12 位定点 (1,5,6) 到 (1,6,5) 浮点的转换结果：

定点数	整数部分	小数部分	浮点数（符号/指数/尾数）	含义
000001000000	000001	000000	0 011111 00000	+1
000001000000	000001	000000	0 011111 00000	+1
111111000000	111111	000000	1 011111 00000	−1
011111111111	011111	111111	0 100100 00000	最大绝对值 ≈ 32
100000000001	100000	000001	1 100100 00000	最小绝对值 = −32
000000000001	000000	000001	0 011001 00000	最小值 = 1/64

图 2-28 (1,5,6) 定点数格式到 (1,6,5) 浮点数转换的仿真结果。5 个代表值分别为：+1，−1，最大绝对值≈32，最小绝对值=−32，最小值=1/64

2.7.2 浮点数到定点数的格式转换

反方向（浮点→定点）一般更复杂：要根据指数相对偏差的大小决定尾数左移还是右移，还要处理 ±∞ 和 ±0 等特殊值。为简化讨论，假设浮点动态范围更大，但定点精度更高（定点的小数位多于浮点尾数位）。

步骤如下：① 校正指数偏差；② 把隐藏的 1 放到小数点左边、小数尾数放到右边；③ 若指数太大则输出饱和到最大值，太小则输出置零；④ 指数有效时，正指数对 1.m 尾数左移（SLL），负指数右移（SRL）——同样存在 1993 标准中 STD_LOGIC_VECTOR 不支持移位的坑（练习 1.20）；⑤ 最后按符号位把符号-幅值格式转回 2 的补码。

图 2-29 给出了与图 2-28 相反方向的 5 个代表值转换结果：

浮点数	符号/指数/尾数	定点数	整数部分	小数部分	含义
001111100000	0 011111 00000	000001000000	000001	000000	+1（无误差）
001111100000	0 011111 00000	000001000000	000001	000000	+1（无误差）
101111100000	1 011111 00000	111111000000	111111	000000	−1（无误差）
010010000000	0 100100 00000	011111111111	011111	111111	最大绝对值（有量化误差）
110010000000	1 100100 00000	100000000000	100000	000000	最小绝对值（有量化误差）
001100100000	0 011001 00000	000000000001	000000	000001	最小值 = 1/64（无误差）

图 2-29 (1,6,5) 浮点数格式到 (1,5,6) 定点数格式转换的仿真结果。5 个代表值分别为：+1，−1，最大绝对值≈32，最小绝对值=−32，最小值=1/64

可以看到 ±1 和最小值没有量化误差；但最大/最小绝对值在浮点数中精度较低（尾数只有 5 位），转换结果并不完美——这正是”浮点范围大、定点精度高”的直观体现。

2.7.3 浮点数乘法

在浮点运算中，乘法是最简单的，所以先讲它。科学记数法下两个数相乘就是尾数相乘、指数相加：

$f_1\times f_2=(a_1{2}^{e_1})\times (a_2{2}^{e_2})=(a_1\times a_2)2^{e_1+e_2}$

换成”隐含 1 + 偏差指数”的浮点格式后：

$$f_1\times f_2=(-1{)}^{s_1}(1.m_1{2}^{e_1-\text{偏差}})\times (-1{)}^{s_2}(1.m_2{2}^{e_2-\text{偏差}})=(-1{)}^{s_1+s_2\mathrm{mod}2}\underset{m_3}{\underbrace{(1.m_1\times 1.m_2)}}{2}^{\underset{e_3}{\underbrace{e_1+e_2-\text{偏差}}}-\text{偏差}}$$

要点有三个：① 符号位是两操作数符号位的异或（模 2 和）；② 指数要减去一次偏差，因为偏差在两个指数中各被引入了一次；③ 特殊值要优先处理——任一因子为 ∞ 则乘积为 ∞，任一因子为零则乘积为零（不支持 NaN 时 $0\times \infty$ 被置为 ∞）。溢出（$e_1+e_2-\text{偏差}\ge E_{\max }$）置为 ∞，下溢（$e_1+e_2-\text{偏差}\le E_{\min }$）置为零。

好消息是规格化非常简单：两个操作数都在 $1.0\le 1.m_{1,2}<2.0$ 范围内，尾数乘积必在 $1.0\le 1.m_3<4.0$，最多左规一位（指数加 1）即可。最后把符号、指数、尾数拼接成结果。

图 2-30 验证了以下 6 组 (1,6,5) 浮点乘法：

1. $(-1)\times (-1)=1.0_{10}=1.00000_2\times 2^{31-\text{偏差}}$
2. $1.75\times 1.75=3.0625_{10}=11.0001_2\times 2^{31-\text{偏差}}=1.10001_2\times 2^{32-\text{偏差}}$（尾数乘积 11.0001 规格化左移一位，指数 31+1=32）
3. 指数：$7+7-\text{偏差}=-17<E_{\min }$ → 乘法下溢，结果置零
4. 指数：$62+62-\text{偏差}=93\ge E_{\max }$ → 乘法溢出，结果置 ∞
5. $\infty \times \infty =\infty$
6. $0\times \infty =$ NaN（本实现置为 ∞）

图 2-30 (1,6,5) 格式的浮点数乘积的仿真结果

（图 2-30 中第 1–4 行是操作数 a 及其符号/指数/尾数分解，第 5–8 行是操作数 b 的分解，第 9–12 行是乘积 r 的分解。）

2.7.4 浮点数加法

浮点加法比乘法复杂，因为两个数必须指数相同才能直接相加。设 $f_3=f_1\pm f_2=a_1{2}^{e_1}\pm a_2{2}^{e_2}$，不失一般性假设第二个数绝对值较小（否则交换两数）。用恒等式把较小的数”非规格化”：

$a_2{2}^{e_2}=\frac{a_2}{2^d}{2}^{e_2+d}$

取 $d=e_1-e_2$，两数指数就一致了：

$\frac{a_2}{2^{e_1-e_2}}{2}^{e_1},a_3=a_1\pm \frac{a_2}{2^{e_1-e_2}}$

接着按符号决定加还是减（同号相加、异号相减），并处理第二个操作数为零的情况：若 $e_2=0$ 或 $d>M$（移位把第二个尾数移没了），直接把较大的第一个操作数作为结果 $f_3$ 转发出去——这一步是个不错的优化，省掉无效计算。

得到的尾数要重新规格化成 $1.m_3$：数出前导零个数（包括隐藏位）并逻辑左移（SLL），同时把初始设为 $e_3=e_1$ 的指数相应调整。特殊值方面：若第一个操作数是 ∞ 或新指数大于 $E_{\max }$，输出置 ∞（于是 $\infty -\infty$ 被置为 ∞，因为不支持 NaN）；新指数小于 $E_{\min }$ 则下溢置零。最后拼接符号、指数、尾数。

图 2-31 验证了以下 6 组 (1,6,5) 浮点加法：

1. $1.0+(-1.0)=0$
2. $9.25+(-10.5)=-1.25_{10}=-1.01000_2\times 2^{31-\text{偏差}}$
3. $1.00111_2\times 2^{2-\text{偏差}}+(-1.00100_2\times 2^{2-\text{偏差}})=0.00011_2\times 2^{2-\text{偏差}}=1.1_2\times 2^{-2-\text{偏差}}$，而 $-2<E_{\min }$ → 下溢 → 非规格数
4. $1.01111_2\times 2^{62-\text{偏差}}+1.11110_2\times 2^{62-\text{偏差}}=11.01101_2\times 2^{62-\text{偏差}}=1.10110_2\times 2^{63-\text{偏差}}$，而 $63\ge E_{\max }$ → 溢出
5. $-\infty +1=-\infty$
6. $\infty -\infty =$ NaN（本实现置为 ∞）

图 2-31 (1,6,5) 格式的浮点数加法的仿真结果

（图 2-31 中第 1–4 行是操作数 a 及其分解，第 5–8 行是操作数 b 及其分解，第 9–12 行是和 r 及其分解。）

2.7.5 浮点数除法

除法是”尾数相除、指数相减”：

$f_1/f_2=(a_1{2}^{e_1})/(a_2{2}^{e_2})=(a_1/a_2)2^{e_1-e_2}$

带入”隐含 1 + 偏差指数”格式：

$$f_1/f_2=(-1{)}^{s_1}(1.m_1{2}^{e_1-\text{偏差}})/(-1{)}^{s_2}(1.m_2{2}^{e_2-\text{偏差}})=(-1{)}^{s_1+s_2\mathrm{mod}2}\underset{m_3}{\underbrace{\frac{1.m_1}{1.m_2}}}{2}^{\underset{e_3}{\underbrace{e_1-e_2+\text{偏差}}}-\text{偏差}}$$

注意与乘法相反：两指数相减后偏差被消掉了，所以结果指数要加回一次偏差。商的符号同样是两符号位异或。尾数除法可用 2.5 节任何算法实现，也可以直接用 lpm_divide——但有个位宽细节：分母和商至少要 $M+1$ 位，而 lpm_divide 中分子与商位宽相同，所以要用 $2\times (M+1)$ 位宽的分子和商来容纳对齐用的扩展位。由于分子分母都在 $1\le 1.m_{1,2}<2$ 内，商必落在 $0.5\le 1.m_3<2$，同样只需一位规格化（指数调整 1）。

特殊值处理：分子为 ∞、分母为零、或 $e_1-e_2+\text{偏差}=e_3\ge E_{\max }$ 时商为 ∞；分子为零、分母为 ∞、或 $e_3\le E_{\min }$ 时商置零；其余情况结果均在有效范围内。

图 2-32 验证了以下 7 组 (1,6,5) 浮点除法：

1. $(-1)/(-1)=1.0_{10}=1.00000_2\times 2^{31-\text{偏差}}$
2. $-10.5/9.25=-1.{\overline{135}}_{10}\approx -1.001_2\times 2^{31-\text{偏差}}$（尾数 5 位产生舍入）
3. $9.25/(-10.5)=-0.880952_{10}\approx -1.11_2\times 2^{30-\text{偏差}}$（商 <1 需右规格化，指数降为 30）
4. 指数：$60-3+\text{偏差}=88>E_{\max }$ → 除法溢出，置 ∞
5. 指数：$3-60+\text{偏差}=-26<E_{\min }$ → 除法下溢，置零
6. $1.0/0=\infty$
7. $0/0=$ NaN（本实现置为 ∞）

图 2-32 (1,6,5) 格式的浮点数除法的仿真结果（图中第 1–4 行为第一个操作数 a 及其符号/指数/尾数，第 5–8 行为第二个操作数 b 及其分解，第 9–12 行为商 r 及其分解。）

2.7.6 浮点数倒数

倒数函数 $1.0/f$ 看似冷门，其实很有用——它可以和乘法器组合成浮点除法器：

$\frac{1.0}{(-1{)}^{s}1.m{2}^e}=(-1{)}^s\frac{2^{-e}}{1.m},\frac{f_1}{f_2}=\frac{1.0}{f_2}\times f_1$

即”乘以分母的倒数等价于除以分母”。当尾数位宽不宽时，尾数倒数可以直接查表实现（case 语句或 M9K 存储模块）。由于 $1\le 1.m<2$，倒数必落在 $0.5<1/(1.m)\le 1$，因此除 $f=1.0$ 外所有值的规格化都只需移一位。

特殊值：∞ 的倒数是零，零的倒数是 ∞。其余值的新指数为：

$e_2=-(e_1-\text{偏差})+\text{偏差}=2\times \text{偏差}-e_1$

图 2-33 验证了以下 5 组 (1,6,5) 浮点倒数：

1. $1/(-2)=-0.5_{10}=-1.0_2\times 2^{30-\text{偏差}}$
2. $1/1.25=0.8_{10}\approx (32+19)/64=1.10011_2\times 2^{30-\text{偏差}}$（查表得 25/32，5 位尾数量化）
3. $1/1.031=0.9697_{10}\approx (32+30)/64=1.11110_2\times 2^{30-\text{偏差}}$
4. $1.0/0=\infty$
5. $1/\infty =0.0$

输入 a	含义	倒数 r	含义
110000000000	−2.0	101111000000	−0.5
101111101000	−1.25	101111010011	−0.796875
101111100001	−1.03125	101111011110	−0.96875
000000000000	0	011111100000	∞
011111100000	∞	000000000000	0

图 2-33 (1,6,5) 格式的浮点数倒数的仿真结果（注意：对除 0 的情况，仿真器的 Transcript 窗口会报告 RECIPROCAL: Floating Point divided by zero。）

2.7.7 浮点操作集成

自己从零搭一整套浮点库是劳动密集型任务。幸运的是，VHDL-2008 已把定点/浮点运算包纳入标准（VHDL-1993 也可通过 David Bishop 的 7000+ 行兼容库使用，下载地址 www.eda.org/fphdl，针对 Altera、Xilinx、Synopsys、Cadence 和 Mentor Graphics 工具有测试过的修改版）。原书在 VHDL-1993 中这样引入库：

-- 原书 VHDL 改写说明：这是 VHDL 库引用语句，Verilog 中无对应物，
-- Verilog 用户需自建函数模块或调用厂商浮点 IP 来获得同等能力
LIBRARY ieee_proposed;
USE ieee_proposed.fixed_float_types.ALL;
USE ieee_proposed.float_pkg.ALL;

该库让 FLOAT 类型直接支持标准运算符：算术 +、-、*、/、ABS、REM、MOD，逻辑 NOT、AND、NAND、OR、NOR、XOR、XNOR，比较 =、/=、>、<、>=、<=，转换 TO_SLV、TO_SFIXED、TO_FLOAT，以及其他 RESIZE、SCALB、LOGB、MAXIMUM、MINIMUM。还有 6 个预定义常数：zerofp(0)、nanfp(NaN)、qnanfp(quiet NaN)、pos_inffp(+∞)、neg_inffp(−∞)、neg_zerofp(−0)；IEEE 854/754 中长度 32/64/128 的预定义类型分别叫 FLOAT32、FLOAT64、FLOAT128。

假如要实现一个 1 符号、6 指数、5 小数位的浮点数，原书 VHDL 改写为 Verilog 时，语言本身没有浮点类型，需要用位向量承载、自建运算模块（或调用 IP）：

// 原书 VHDL 改写为 Verilog：浮点数以 12 位位向量表示 (1,6,5) 格式，
// 运算通过自定义函数/IP 模块完成（Verilog 无内建浮点类型）
module fp_demo
  (input  wire [11:0] a, b,        // (1,6,5) 浮点操作数
   output wire [11:0] s, p);       // 和 s、积 p
  fp_add u_add (.a(a), .b(b), .sum(s));   // 自定义浮点加法模块
  fp_mul u_mul (.a(a), .b(b), .prod(p));  // 自定义浮点乘法模块
endmodule

VHDL 中 s <= a+b; p <= a*b; 写起来只有一行，是因为左右两侧数据类型相同、无需缩放和调整大小；Verilog 用户则要自己保证位宽与规格化。此外，库的默认配置是”最贵”的：舍入为 round_nearest、denormalize 和 error_check 为 true、3 个保护位。反过来把舍入设为 round_zero（截断）、后两者设为 false、保护位 0，就是最省硬件的配置。大多数操作还提供函数形式：算术有 ADD、SUBTRACT、MULTIPLY、DIVIDE、REMAINDER、MODULO、RECIPROCAL、MAC、SQRT，比较有 EQ、NE、GT、LT、GE、LE（类似 FORTRAN 命名）。缩放函数 SCALB(y,n) 实现 $y\times 2^n$，比普通乘/除法省得多，这对 DSP 里的 2 的幂次缩放特别合适。

接口数据类型的选择也有讲究：由于多数仿真器尚未完全支持负索引数组，I/O 保留标准 STD_LOGIC_VECTOR 更稳妥。库提供的位保留（bit-preserving）转换函数只是把位向量的含义重新解释为同长度的定点/浮点类型，由预处理器完成，不消耗任何硬件资源；而 sfixed 与 float 之间的真转换要保留数值，需要可观的硬件。图 2-34 展示了用 fp_ops.exe 程序测试数据计算的结果。

图 2-34 使用 fp_ops.exe 程序测试数据计算

例 2.20：一个 32 位浮点运算单元

把前面所有操作集成起来，就得到一个 32 位 FPU：支持 fix2fp、fp2fix、加、减、乘、除、倒数、缩放 8 种操作。原书 VHDL 改写为 Verilog 时必须注意：Verilog 语言当前不具有浮点库支持（原书脚注也明确说明此例的等效 Verilog 代码不能直接执行），因此运算核心用厂商浮点 IP 或 2.7.1–2.7.6 节的自建模块代替，这里给出体现相同架构的可综合骨架：

// 原书 VHDL 改写为 Verilog：32 位浮点运算单元（架构级改写，
// 浮点核心用自建/IP 模块占位，因为 Verilog 无浮点库）
module fpu
  (input  wire [3:0]  sel,      // 操作选择：0~7
   input  wire [31:0] dataa,    // 第一操作数（FLOAT32 或 16.16 定点位型）
   input  wire [31:0] datab,    // 第二操作数
   input  wire signed [7:0] n,  // 缩放因子 2**n
   output reg  [31:0] result);  // 系统输出
  // 操作码: 0=fix2fp 1=fp2fix 2=add 3=sub 4=mul 5=div 6=rec 7=scale
  wire [31:0] a = dataa, b = datab;              // 位保留：位型重新解释为浮点
  wire signed [31:0] fa = dataa;                 // 位保留：重新解释为 16.16 定点
  wire [31:0] r_add  = fp_add(a, b);
  wire [31:0] r_sub  = fp_sub(a, b);
  wire [31:0] r_mul  = fp_mul(a, b);
  wire [31:0] r_div  = fp_div(a, b);
  wire [31:0] r_rec  = fp_reciprocal(a);
  wire [31:0] r_scl  = fp_scalb(a, n);           // a * 2**n，比乘除法省资源
  wire [31:0] r_f2f  = sfixed_to_float(fa);      // 16.16 定点 -> FLOAT32
  wire [31:0] r_f2i  = float_to_sfixed(a);       // FLOAT32 -> 16.16 定点
  always @(*) begin
    case (sel)
      4'd0: result = r_f2f;
      4'd1: result = r_f2i;
      4'd2: result = r_add;
      4'd3: result = r_sub;
      4'd4: result = r_mul;
      4'd5: result = r_div;
      4'd6: result = r_rec;
      4'd7: result = r_scl;
      default: result = r_scl;
    endcase
  end
endmodule

架构与原书完全一致：sel 选择 8 种操作，两个 32 位输入向量经”位保留转换”重新解释为 FLOAT32 或 16.16 定点数（不消耗硬件），操作 0/1 是定点↔浮点转换，2–5 是四则运算，6 是倒数，7 是以 2 为底的缩放，最后按操作类型选择输出位型的解释方式（只有选项 1 的输出是定点，其余都是 FLOAT32）。

资源方面，该设计使用 8112 个 LE 和 7 个嵌入式乘法器；由于不含寄存器，无法测量寄存器性能。仿真前先用 MATLAB 生成测试激励：用 %tx/%bx 格式可把 32/64 位浮点显示为十六进制，例如 x=1/3 时输出 FLOAT32 := X"3EAAAAAB"; --0.333333。本书资料中的 fp_ops.exe 可算出 32/64 位浮点基本运算的测试数据；输入 a=1/3、b=2/3 即得到图 2-34 的加、减、乘、除、倒数测试数据。转换测试用定点值 0001.0000（即 FLOAT32 的 3F800000，对应 1.0）；缩放操作用 $2^1$，即 $1/3\times 2=2/3$。整体仿真如图 2-35 所示。

图 2-35 8 个函数浮点运算单元 fpu 的仿真

2.7.8 浮点数合成结果

VHDL-2008 库能为任意规格的浮点数写出紧凑代码，缺点是大量算术运算不做流水线，整体速度不高。要高吞吐量，应使用厂商预定义的流水线浮点模块：Xilinx 有 LOGICORE 浮点 IP，Altera 有一整套 LPM 函数（可作图形模块或组件库实例化）。表 2-17 给出了 Altera LPM 32 位浮点模块的流水线范围和默认值：

表 2-17 Altera LPM 32 位浮点数模块流水线

模块	流水线范围	默认值
整数/定点数到浮点数	6	6
浮点数到整数/定点数	6	6
浮点数加法/减法	7…14	11
浮点数乘法	5、6、10、11	5
浮点数除法	6、14、33	33
浮点数倒数	20	20

图 2-36 Altera VHDL-2008 运算器和 LPM 模块的速度值

图 2-36 的吞吐量对比中 LPM 流水线优势明显。但要注意：许多 DSP 系统存在反馈环路，通常无法流水线化——这时只能用 FSM”等待”计算完成，LPM 模块反而用不上。测量 VHDL-2008 库的 Fmax 时，是在输入/输出端口加了寄存器、模块内部不用流水线的条件下进行的。

图 2-37 汇总了 32 位和 64 位下 8 个基本模块的资源：实线是 VHDL-2008 默认设置，虚线是最小硬件消耗设置（round_zero 截断、denormalize 和 error_check 为 false、0 保护位）。定点↔FLOAT32 转换约需 400 个 LE；基本运算 +、−、× 约 1000 个 LE，除法约加倍。另一个有意思的对比：LPM 的浮点除法和倒数大量使用嵌入式乘法器和一个 M9K 存储模块（图 2-37 中未显示），而 VHDL-2008 的等效操作不用任何存储器或乘法器；指数型乘除 SCALB 则比等效乘除法节省得多。

图 2-37 用于默认设置（实线）和最小硬件消耗（虚线）的 VHDL-2008 操作和 Altera LPM 模块的数据大小

2.8 MAC 与 SOP

DSP 算法中计算最密集的是乘-累加（Multiply-ACcumulate, MAC）。以线性卷积和为例：

y[n] = f[n] \times x[n] = \sum_{k=0}^{L-1} f[k]\,x[n-k] \tag{2-46}

每采一个样点，计算 $y[n]$ 都要连续做 $L$ 次乘法和 $L-1$ 次加法，即积之和（Sum Of Product, SOP）。这意味着一个 $N\times N$ 位乘法器要和累加器融合在一起，如图 2-38(a) 所示。全精度 $N\times N$ 位乘积的位宽是 $2N$；若两操作数都是（对称的）有符号数，乘积只有 $2N-1$ 个有效位（两个符号位重合）。为了防止累加溢出，累加器通常要多设计 $K$ 个额外保护位。

例 2.21 模拟器件公司的 PDSP ADSP21xx 系列内含一个 $16\times 16$ 阵列乘法器和一个多出 8 位位宽的累加器（总位宽 $32+8=40$ 位）。有这 8 个额外位，至少可无损累加 $2^8$ 次；若两操作数都是对称有符号数，则可累加 $2^9$ 次。为了产生所需的输出格式，现代 PDSP 还包括一个桶式移位器，可在一个时钟周期内完成位调整。

对传统 PDSP 而言，检查和处理累加器溢出会中断数据流、带来明显的时序负担；而 DSP 要求实时计算、不希望中断，所以通过正确选择保护位位数来消除溢出处理负担是主流定点 DSP 的关键设计点。除了传统 MAC 单元，计算积之和还有另一种方法——分布式算法，下一节讨论。

2.8.1 分布式算法基础

分布式算法（Distributed Arithmetic, DA）是一项重要的 FPGA 技术，广泛用于计算积之和：

y = \langle c, x \rangle = \sum_{n=0}^{N-1} c[n] \times x[n] \tag{2-47}

除卷积外，相关、DFT 计算和 RNS 反演映射也都可以表述成 SOP。用传统算法单元完成一次滤波周期约需 $N$ 个 MAC 循环，即使流水线化也只能有限缩短。DA 的出发点是一个重要观察：许多 DSP 应用中滤波系数 $c[n]$ 是事先已知的常数，部分乘积项 $c[n]x[n]$ 就退化成常数乘法（缩放）——这是 DA 设计的先决条件，“通用乘法器”在这里纯属浪费。

DA 的历史：1973 年 Croisier 发表首篇论文，Peled 和 Liu 完成推广，Yiu 将其扩展到有符号数，Kammeyer 和 Taylor 研究量化效应，White 和 Kammeyer 撰写教程，如今已进入教科书。为理解 DA 的设计范例，考察内积：

y = \langle c, x \rangle = \sum_{n=0}^{N-1} c[n] \times x[n] = c[0]x[0] + c[1]x[1] + \cdots + c[N-1]x[N-1] \tag{2-48}

设系数 $c[n]$ 为常数、$x[n]$ 为变量。无符号 DA 系统把 $x[n]$ 按位展开：

x[n] = \sum_{b=0}^{B-1} 2^b \times x_b[n], \quad x_b[n] \in [0,1] \tag{2-49}

其中 $x_b[n]$ 是 $x[n]$ 的第 $b$ 位。代入内积并重新分配求和顺序（这正是”分布式”算法名称的由来）：

$$\begin{array}{rl}y& =c[0](x_{B-1}[0]2^{B-1}+x_{B-2}[0]2^{B-2}+\cdots +x_0[0]2^0)\\ & +c[1](x_{B-1}[1]2^{B-1}+x_{B-2}[1]2^{B-2}+\cdots +x_0[1]2^0)\\ & ⋮\\ & +c[N-1](x_{B-1}[N-1]2^{B-1}+\cdots +x_0[N-1]2^0)\\ & =(c[0]x_{B-1}[0]+c[1]x_{B-1}[1]+\cdots +c[N-1]x_{B-1}[N-1])2^{B-1}\\ & +(c[0]x_{B-2}[0]+\cdots +c[N-1]x_{B-2}[N-1])2^{B-2}\\ & ⋮\\ & +(c[0]x_0[0]+\cdots +c[N-1]x_0[N-1])2^0\end{array}$$

写成紧凑形式：

y = \sum_{b=0}^{B-1} 2^b \times \sum_{n=0}^{N-1} \underbrace{f(c[n], x_b[n])}_{c[n]\times x_b[n]} \tag{2-51}

核心在于函数 $f(c[n],x_b[n])$ 的实现：由于 $c[n]$ 是常数、$x_b[n]\in \{0,1\}$，每个内层和只有 $2^N$ 种取值——可以预编程一个 $2^N$ 字的查找表（LUT），输入 $N$ 位向量 $x_b=[x_b[0],x_b[1],\dots ,x_b[N-1]]$，直接输出对应的内层和。然后把各次查表结果按相应的 2 的幂加权累加，用图 2-38(b) 所示的移位加法器即可高效实现。$B$ 次查表循环后，内积 $y$ 计算完毕——整个过程没有一个通用乘法器，乘法被”查表+移位累加”取代，这正是 DA 在 FPGA 上高效的秘密。

图 2-38 PDSP 结构和 DA 体系结构

例 2.22 无符号 DA 卷积

在上一小节建立的分布式算法（DA）框架中，内积 $y=<c,x>={\sum }_{n=0}^{2}c[n]x[n]$ 不需要真正的乘法器：把每个输入字的同一位拼成一个地址去查表，查表结果再按位权累加即可。本例用一个三阶内积把整个流程走一遍。设 3 位系数分别为 $c[0]=2$、$c[1]=3$、$c[2]=1$，那么查找表 $f(c[n],x_b[n])$ 的全部内容如下（三个地址位 $x_b[2],x_b[1],x_b[0]$ 分别对应 $c[2],c[1],c[0]$ 所乘输入的当前位）：

$x_b[2]$	$x_b[1]$	$x_b[0]$	$f(c[n], x_b[n])$
0	0	0	$1 \times 0 + 3 \times 0 + 2 \times 0 = 0_{10} = 000_2$
0	0	1	$1 \times 0 + 3 \times 0 + 2 \times 1 = 2_{10} = 010_2$
0	1	0	$1 \times 0 + 3 \times 1 + 2 \times 0 = 3_{10} = 011_2$
0	1	1	$1 \times 0 + 3 \times 1 + 2 \times 1 = 5_{10} = 101_2$
1	0	0	$1 \times 1 + 3 \times 0 + 2 \times 0 = 1_{10} = 001_2$
1	0	1	$1 \times 1 + 3 \times 0 + 2 \times 1 = 3_{10} = 011_2$
1	1	0	$1 \times 1 + 3 \times 1 + 2 \times 0 = 4_{10} = 100_2$
1	1	1	$1 \times 1 + 3 \times 1 + 2 \times 1 = 6_{10} = 110_2$

注意这张表的内容只与系数有关，输入数据只是作为地址出现，所以系数不变时表可以预先固化在 ROM 里。现在取输入 $x[n]=\{x[0]=1_{10}=001_2,x[1]=3_{10}=011_2,x[2]=7_{10}=111_2\}$，逐位查表并累加：

步骤 $t$	${x}_{t}\left\lbrack 2\right\rbrack$	${x}_{t}\left\lbrack 1\right\rbrack$	${x}_{t}\left\lbrack 0\right\rbrack$	$f\left\lbrack t\right\rbrack + {ACC}\left\lbrack {t - 1}\right\rbrack = {ACC}\left\lbrack t\right\rbrack$
0	1	1	1	$6 \times {2}^{0} +$	0 = 6
1	1	1	0	$4 \times {2}^{1} +$	6 = 14
2	1	0	0	$1 \times {2}^{2} +$	14 = 18

每一步取的是三个输入字的同一位（$t=0$ 取最低位，$t=2$ 取最高位），查表得到 $f[t]$ 后乘以 $2^t$ 再累加。数值校验：

$$\begin{array}{rl}\mathrm{y}& =<c,x>=c[0]x[0]+c[1]x[1]+c[2]x[2]\\ & =2\times 1+3\times 3+1\times 7=18\end{array}$$

结果一致，说明”查表 + 移位累加”确实等价于普通的乘加运算。

在硬件上还有一个实现细节值得注意：表中"$f[t]\times 2^t$"这一步如果真的每次把中间结果左移 $b$ 位，就需要一个昂贵的桶式移位器。更聪明的做法是反过来——保持累加器每次累加的量不动，让累加寄存器本身在每次迭代后逐位右移。很容易证明这样最终得到的结果完全相同，却省掉了移位器。

那 DA 到底比传统乘累加（MAC）快多少？可以粗略比较 N 阶 B 位线性卷积的带宽。图 2-38 给出了传统 PDSP 的体系结构和采用分布式算法的相同实现。假设查表和乘法的时延相同，即 $\tau =\tau (\text{LUT})=\tau (\text{MUL})$，则总等待时间分别是 DA 的 $B\tau (\text{LUT})$ 和 PDSP 的 $N\tau (\text{MUL})$。当输入位宽 $B$ 较小时，DA 设计的速度可以显著超过基于 MAC 的设计，而且根本不需要硬件乘法器——这正是它在 FPGA 上受欢迎的原因。第 3 章将结合具体滤波器再做定量比较。

2.8.2 有符号的 DA 系统

前面的推导默认输入是无符号数，但实际信号几乎总是有符号的，所以要修改式(2-48)。补码的最高有效位带有符号信息：例如表 2-1 中 $-3$ 编码为 $101_2=-4+0+1=-3_{10}$，也就是符号位的权是负的。于是采用 $(B+1)$ 位表示：

$$x[n]=-2^B{x}_B[n]+\sum _{b=0}^{B-1}{2}^b\times x_b[n]\text{\hspace{1em}(2-52)}$$

与式(2-50)联立，输出 $y$ 变成：

$$y=-2^B\times f(c[n],x_B[n])+\sum _{b=0}^{B-1}{2}^b\times \sum _{n=0}^{N-1}f(c[n],x_b[n])\text{\hspace{1em}(2-53)}$$

和无符号版本相比，唯一的区别是最后一步（符号位那一步）要从”累加”改成”减去”。在硬件上有两种实现选择：

• 带有加/减控制的累加器；
• 采用具有一个额外输入的 ROM。

推荐使用最常见的可转换累加器：如果用第二种方案，LUT 多出的那一位输入会把表的内容翻倍，代价更高。下面的例子演示加/减转换设计的完整处理步骤。

例 2.23 有符号 DA 的内积

仍看三阶内积 $y=<c,x>={\sum }_{n=0}^{2}c[n]x[n]$，数据为 4 位二进制补码，系数 $c[0]=-2$、$c[1]=3$、$c[2]=1$。LUT 表如下：

$x_b[2]$	$x_b[1]$	$x_b[0]$	$f(c[k],x_b[n])$
0	0	0	$1\times 0+3\times 0-2\times 0=0_{10}$
0	0	1	$1\times 0+3\times 0-2\times 1=-2_{10}$
0	1	0	$1\times 0+3\times 1-2\times 0=3_{10}$
0	1	1	$1\times 0+3\times 1-2\times 1=1_{10}$
1	0	0	$1\times 1+3\times 0-2\times 0=1_{10}$
1	0	1	$1\times 1+3\times 0-2\times 1=-1_{10}$
1	1	0	$1\times 1+3\times 1-2\times 0=4_{10}$
1	1	1	$1\times 1+3\times 1-2\times 1=2_{10}$

输入取 $x[0]=1_{10}=0001_{2C}$、$x[1]=-3_{10}=1101_{2C}$、$x[2]=7_{10}=0111_{2C}$。前 3 步（$t=\mathrm{0..2}$）与无符号情况一样做加法，第 4 步（符号位 $t=3$）改为减法：

步骤t	$x_{t}[2]$	$x_{t}[1]$	$x_{t}[0]$	$f[t]\times 2^{t}+Y[t-1]=Y[t]$
0	1	1	1	$2\times 2^{0}$ + 0 = 2
1	1	0	0	$1\times 2^{1}$ + 2 = 4
2	1	1	0	$4\times 2^{2}$ + 4 = 20
	$x_{t}[2]$	$x_{t}[1]$	$x_{t}[0]$	$f[t]\times (-2^{t})+Y[t-1]=Y[t]$
3	0	1	0	$3\times (-2)^{3}$ + 20 = -4

数值校验：$c[0]x[0]+c[1]x[1]+c[2]x[2]=-2\times 1+3\times (-3)+1\times 7=-4$，两者吻合。可以看到，有符号 DA 的全部改动就是在最后一步把累加器切换成减法，其余硬件完全不变。

2.8.3 改进的 DA 解决方案

基本 DA 结构还可以从两个方向改进：一是缩小规模，二是提高速度。

缩小规模——表分割。 LUT 的规模随地址位宽（即系数个数 $N$）呈指数增长，如果 $N$ 太大，单个 LUT 装不下整张表，就可以把长的内积拆成几段，用若干个小表分别查表再把结果相加。对长度为 $LN$ 的内积：

$$y=<c,x>=\sum _{n=0}^{LN-1}c[n]x[n]\text{\hspace{1em}(2-54)}$$

将它分配到 $L$ 个独立的 $N$ 阶并行 DA LUT 中：

$$y=<c,x>=\sum _{l=0}^{L-1}\sum _{n=0}^{N-1}c[Nl+n]x[Nl+n]\text{\hspace{1em}(2-55)}$$

如图 2-39 所示，实现一个 4N 的 DA 设计需要 3 个辅助的后向加法器；表的规模从一个 $2^{4N}\times B$ 的 LUT 降低到 4 个 $2^N\times B$ 的表。加上流水线寄存器后速度也不受损失，因为各小表是并行工作的。

图 2-39 将表分割以产生简化规模的分布式算法

提高速度——每步处理多位。 基本 DA 结构每个时钟周期只接收每个输入字中的一位，若改成同时接收两位，速度就能翻倍。极限情况是图 2-40 所示的完全流水线字并行体系结构：为每个位向量 $x_b[n]$ 各配一个内容相同的 ROM，在一个 LUT 周期内就能算出 4 位有符号系数、长度为 4 的积之和的新结果。当然最大速度的代价很昂贵：输入位宽加倍，LUT、寄存器和加法器都要加倍。当系数个数 $N$ 限制在 4 或 8 个时，这一改进的性能非常突出，甚至优于所有商业可编程信号处理器，第 3 章会给出对比数据。

图 2-40 速度最优的高阶分布式算法

2.9 利用 CORDIC 计算特殊函数

FPGA 上做信号处理时，经常遇到超越函数，如 $\sqrt{x}$ 或 $\arctan y/x$。最直接的办法是用泰勒级数近似：

$$f(x)=\sum _{k=0}^K\frac{f^k(x_0)}{k!}(x-x_0{)}^k\text{\hspace{1em}(2-56)}$$

其中 $f^k(x)$ 是 $f(x)$ 的第 $k$ 次微分，$k!=k\times (k-1)\dots \times 1$，这样问题被化成一系列乘法和加法。但更高效的方法是坐标旋转数字计算机（Coordinate Rotation Digital Computer, CORDIC）算法。它的应用非常广泛：从便携式计算器，到适应性滤波器、FFT、DCT、解调器和神经网络等主流 DSP 对象。基础算法出自 Volder 和 Walther 的两篇经典论文，此后有范围拓展、量化误差分析、VLSI 与 FPGA 实现以及在 DA 中的实现等大量后续研究，Hu 在 1992 年《IEEE 信号处理杂志》的评论文章中给出了详细回顾。

Volder 提出的原始 CORDIC 算法用于平面直角坐标 $(x,y)$ 和极坐标 $(R,\theta )$ 之间的自由变换。Walther 把它推广到三种模式：圆周变换（$m=1$）、线性变换（$m=0$）和双曲线变换（$m=-1$）。每种模式有两个旋转方向。所谓向量化，就是把原点在 $(X_0,Y_0)$ 的向量旋转到横轴上：迭代让 $Y_K$ 收敛到 0。所谓旋转，则是让向量转过角度 ${\theta }_0$：迭代让存放角度的 Z 寄存器收敛到 0。角度 ${\theta }_k$ 的选取保证每次迭代只需一次加法和一次移位。表 2-18 给出了三种模式的角度选择。

表 2-18 CORDIC 算法模式

模式	角度 $\theta_{k}$	移位序列	半径因子
圆周m=1	$\tan^{-1}(2^{-k})$	0,1,2,...	$K_{1}=1.65$
线性m=0	$2^{-k}$	1,2,...	$K_{0}=1.0$
双曲线m=-1	$\tanh^{-1}(2^{-k})$	1,2,3,4,4,...	$K_{-1}=0.80$

现在正式定义 CORDIC 算法。

算法 2.24：CORDIC 算法。 每次迭代实现以下映射：

$$\left[\begin{array}{l}{X}_{k+1}\\ Y_{k+1}\end{array}\right]=\left[\begin{array}{cc}1& m{\delta }_k{2}^{-k}\\ {\delta }_k{2}^{-k}& 1\end{array}\right]\left[\begin{array}{l}{X}_k\\ Y_k\end{array}\right]\text{\hspace{1em}(2-57)}$$ $$Z_{k+1}=Z_k+{\delta }_k{\theta }_k$$

其中 ${\theta }_k$ 见表 2-18，${\delta }_k=\pm 1$，两个旋转方向分别是 $Z_K\to 0$ 和 $Y_K\to 0$。

式(2-57)的核心是：乘以 ${\delta }_k{2}^{-k}$ 只需移位，乘以 $\pm 1$ 只需选加减，所以整个迭代不含乘法器。3 种模式 × 2 个方向，共 6 种操作模式，见表 2-19。正确选择初始值后，可以直接算出 $X\times Y$、$Y/X$、$\sin (Z)$、$\cos (Z)$、${\tan }^{-1}(Z)$、$\sinh (Z)$、$\cosh (Z)$ 和 $\tanh (Z)$；其余函数可以通过组合多种模式得到，例如：

$\tan (Z)=\sin (Z)/\cos (Z)$，模式：$m=1,0$

$\tanh (Z)=\sinh (Z)/\cosh (Z)$，模式：$m=-1,0$

$\exp (Z)=\sinh (Z)+\cosh (Z)$，模式：m=-1；x=y=1

$\ln (W)=2{\tanh }^{-1}(Y/X)$，模式：$m=-1$，取 $X=W+1,Y=W-1$

$\sqrt{W}=\sqrt{X^2-Y^2}$，模式：$m=1$，取 $X=W+\frac{1}{4},Y=W-\frac{1}{4}$

表 2-19 CORDIC 算法的操作模式 $m$

m	$Z_{K} \rightarrow 0$	$Y_{K} \rightarrow 0$
1	$X_{K}=K_{1}(X_{0}\cos(Z_{0})-Y_{0}\sin(Z_{0}))$ $Y_{K}=K_{1}(X_{0}\cos(Z_{0})+Y_{0}\sin(Z_{0}))$	$X_{K}=K_{1}\sqrt{X_{0}^{2}+Y_{0}^{2}}$ $Z_{K}=Z_{0}+\arctan(Y_{0}/X_{0})$
0	$X_{K}=X_{0}$ $Y_{K}=Y_{0}+X_{0}Z_{0}$	$X_{K}=X_{0}$ $Z_{K}=Z_{0}+Y_{0}/X_{0}$
-1	$X_{K}=K_{-1}(X_{0}\cosh(Z_{0})-Y_{0}\sinh(Z_{0}))$ $Y_{K}=K_{-1}(X_{0}\cosh(Z_{0})+Y_{0}\sinh(Z_{0}))$	$X_{K}=K_{-1}\sqrt{X_{0}^{2}-Y_{0}^{2}}$ $Z_{K}=Z_{0}+\tanh^{-1}(Y_{0}/X_{0})$

仔细分析式(2-57)会发现：迭代向量只能落在图 2-41(a)所示的曲线上；而且每次迭代向量的长度都会变化，如图 2-41(b)所示。这个长度变化与起始角度无关，$K$ 次迭代后总是出现同一个总伸缩量，称为半径因子（表 2-18 最后一列，如圆周模式 $K_1\approx 1.65$）。使用时必须把它补偿掉。收敛条件方面：线性与圆周模式只要总角度覆盖即可；双曲线模式则要求形如 $n_{k+1}=3n_k+1$ 的迭代（即第 4、13、40、121 次……）必须重复执行，否则不收敛。

(a) 模式

(b) 圆周向量化的示例
图 2-41 CORDIC 算法

输出精度可以用 Hu 开发的程序估算。图 2-42 显示圆周模式的有效位精度依赖于 $X$、$Y$ 路径宽度和迭代次数。经验法则是：需要 $b$ 位输出精度时，$X$、$Y$ 路径要加 ${\log }_2(b)$ 个保护位。图 2-43 表明 $Z$ 路径的位宽与 $X$、$Y$ 相同时即可达到同等精度。

图 2-42 圆周模式中的有效位

与圆周模式不同，双曲线 CORDIC 的精度无法解析计算，因为精度与第 $K$ 次迭代时 $z(k)$ 的角度值有关，只能靠仿真评估。图 2-44 是对每种位宽/迭代次数组合取 1000 个测试值算出的最小精度：三维图给出迭代次数、X/Y 位宽与最终有效位的最低精度，等值线则体现二者的折中。例如要 10 位精度，可以采用 21 位 X/Y 路径加 18 次迭代，也可以用 24 位 X/Y 路径加 14 次迭代——用位宽换迭代次数。

图 2-43 圆周模式的相位分解

图 2-44 双曲线模式的有效位

CORDIC 体系结构

CORDIC 的实现有两种基本结构：较省资源的状态机和高速全流水线处理器。

如果计算时间不严格，可用图 2-45 所示的状态机：每个时钟周期精确执行式(2-57)的一次迭代。设计中最复杂的部件是两个桶式移位器，它们可以用一个单一桶式移位器加多路转换器（图 2-46），或一个串行（右移或左移/右移）移位器来代替。表 2-20 给出了 Xilinx XC3K FPGA 上 13 位实现的各种方案对比，可以在资源（$\sum LE$）和循环数之间权衡：双桶式移位器方案最快（12 个循环）但用 81 个 LE，串行右移方案最省 LE 但需要更多循环。

图 2-45 CORDIC 状态机

图 2-46 降低复杂性的 CORDIC 状态机

表 2-20 兼具 13 位外加 X/Y 路径符号位的 CORDIC 状态机的效率评估(Xilinx XC3K)
(缩写：Ac=accumulator(累加器)；BS=barrel shifter(桶式移位器)；RS=serial right shifter(串行右移移位器)；LRS=serial left/right shifter(串行左/右移位器))

结构	寄存器	多路复用器	加法器	移位器	$\sum LE$	循环
2BS+2Ac	2×7	0	2×14	2×19.5	81	12
2RS+2Ac	2×7	0	2×14	2×6.5	55	46
2LRS+2Ac	2×7	0	2×14	2×8	58	39
1BS+2Ac	7	3×7	2×14	19.5	75.5	20
1RS+2Ac	7	3×7	2×14	6.5	62.5	56
1LRS+2Ac	7	3×7	2×14	8	64	74
1BS+1Ac	3×7	2×7	14	19.5	68.5	20
1RS+1Ac	3×7	2×7	14	6.5	55.5	92
1LRS+1Ac	3×7	2×7	14	8	57	74

如果需要高速，就采用图 2-47 所示的全流水线处理器，图中画出了圆周 CORDIC 的 8 次迭代。在 $K$ 级流水线的起始延迟之后，每个时钟周期都产生一个新的输出值——吞吐量达到每周期一个结果。与阵列乘法器类似，流水线 CORDIC 的 LE 复杂度随位宽增加呈平方增长。

图 2-47 快速 CORDIC 流水线

下面的例题给出一个圆周-向量化全流水线设计的前 4 个步骤。

例 2.25 向量化模式中的圆周 CORDIC

目标是从输入 $(X_0,Y_0)$ 求 $R=\sqrt{X_0^2+Y_0^2}$ 和 $\phi =\arctan (Y_0/X_0)$。第一次迭代要把落在第二或第三象限的向量转到第一或第四象限（±90° 预旋转），之后的移位序列是 0、0、1、2……前 4 步的旋转角度依次为 $\arctan (\infty )=90^{\circ }$、$\arctan (2^0)=45^{\circ }$、$\arctan (2^{-1})=26.5^{\circ }$、$\arctan (2^{-2})=14^{\circ }$。原书 VHDL 代码改写为 Verilog，实现 8 位输入、4 级流水线：

// 原书 VHDL 改写为 Verilog：4 级流水线圆周向量化 CORDIC
module cordic (
    input             clk,      // 系统时钟
    input             reset,    // 异步复位
    input signed [7:0] x_in,    // 实部/x 输入
    input signed [7:0] y_in,    // 虚部/y 输入
    output reg signed [8:0] r,    // 半径结果
    output reg signed [8:0] phi,  // 相位结果
    output reg signed [8:0] eps   // y 残差
);
    reg signed [8:0] x0,y0,z0, x1,y1,z1, x2,y2,z2, x3,y3,z3;
 
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            {x0,y0,z0} = 0; {x1,y1,z1} = 0;
            {x2,y2,z2} = 0; {x3,y3,z3} = 0;
            r <= 0; eps <= 0; phi <= 0;
        end else begin
            // 最后一组值最先赋给输出（顺序语句语义）
            r   <= x3;
            phi <= z3;
            eps <= y3;
 
            // 第3级：旋转14°
            if (y2 >= 0) begin
                x3 <= x2 + (y2 >>> 2);
                y3 <= y2 - (x2 >>> 2);
                z3 <= z2 + 14;
            end else begin
                x3 <= x2 - (y2 >>> 2);
                y3 <= y2 + (x2 >>> 2);
                z3 <= z2 - 14;
            end
 
            // 第2级：旋转26°
            if (y1 >= 0) begin
                x2 <= x1 + (y1 >>> 1);
                y2 <= y1 - (x1 >>> 1);
                z2 <= z1 + 26;
            end else begin
                x2 <= x1 - (y1 >>> 1);
                y2 <= y1 + (x1 >>> 1);
                z2 <= z1 - 26;
            end
 
            // 第1级：旋转45°
            if (y0 >= 0) begin
                x1 <= x0 + y0;
                y1 <= y0 - x0;
                z1 <= z0 + 45;
            end else begin
                x1 <= x0 - y0;
                y1 <= y0 + x0;
                z1 <= z0 - 45;
            end
 
            // 第0级：x_in<0 时预旋转 0 / +90 / -90°
            if (x_in >= 0) begin
                x0 <= x_in;  y0 <= y_in;  z0 <= 0;
            end else if (y_in >= 0) begin
                x0 <= y_in;  y0 <= -x_in; z0 <= 90;
            end else begin
                x0 <= -y_in; y0 <= x_in;  z0 <= -90;
            end
        end
    end
endmodule

代码的读法是”自下而上”：第 0 级处理象限预旋转，第 1~3 级分别做 45°、26.5°、14° 的微旋转；每级用算术右移（>>>）实现 $\times 2^{-k}$，用一个符号判断 ${\delta }_k$，完全不含乘法器。图 2-48 给出了 $X_0=-41$、$Y_0=55$ 的仿真结果：半径被放大到 $R=X_K=111=1.618\sqrt{X_0^2+Y_0^2}$（这就是半径因子），累积角 $\arctan (Y_0/X_0)=123^{\circ }$。该设计占用 276 个 LE，Speed 合成优化下运行速度 209.6 MHz，未使用嵌入式乘法器。

图 2-48 CORDIC 仿真结果

实际 LE 数量（276）比理论估算的 4 级 8 位流水线所需 $5\times 8\times 3=120$ 个多，原因之一是：FPGA 的快速进位模式中 LE 只有 3 个输入，而可转换加法器/减法器需要 4 输入 LUT，故一个 N 位可转换加/减法器要占 2N 个 LE。若改用有 4 个输入的 ALM 型 LE（参阅图 1.7(a)），LE 数量可以降低一半。

2.10 用 MAC 调用计算特殊函数

CORDIC 以可接受的成本实现了很多函数，但它的缺点是：高精度需要的迭代次数与位数成正比，流水线实现后延迟很大。新一代 FPGA（如 Spartan、Cyclone 系列，参阅表 1-4）内置了快速嵌入式阵列乘法器，这使多项式逼近（即式(2-56)的泰勒级数思路）成为现实可行的选择。泰勒级数对 $\exp (x)$ 这类函数收敛很快，但对 $\arctan (x)$ 这类函数要凑够精度需要很多乘积项，此时应改用切比雪夫逼近来减少所需项数。

2.10.1 切比雪夫逼近

切比雪夫逼近以切比雪夫多项式为基础：

$$T_k(x)=\cos (k\times \arccos (x))(-1⩽x⩽1)\text{\hspace{1em}(2-58)}$$

虽然 $T_k(x)$ 看起来像三角函数，但利用代数恒等式化简后它就是真正的多项式，前几个如下：

$$\begin{array}{rl}& T_0(x)=1\\ & T_1(x)=x\\ & T_2(x)=2x^2-1\\ & T_3(x)=4x^3-3x\\ & T_4(x)=8x^4-8x^2+1\\ & T_5(x)=16x^5-20x^3+5x\\ & T_6(x)=32x^6-48x^4+18x^2-1\\ & ⋮\end{array}\text{\hspace{1em}(2-59)}$$

文献[85]列出了前 12 个多项式，图 2-49 画出了前 6 个的图形。它们满足递推规则：

$$T_k(x)=2x{T}_{k-1}(x)-T_{k-2}(x)\forall k⩾2\text{\hspace{1em}(2-60)}$$

函数逼近写成：

$$f(x)=\sum _{k=0}^{N-1}c(k)T_k(x)\text{\hspace{1em}(2-61)}$$

图 2-49 前 6 个切比雪夫多项式

由于离散切比雪夫多项式彼此正交，正、逆变换都是唯一的（双向单射）。那为什么式(2-61)比泰勒逼近

$$f(x)=\sum _{k=0}^{N-1}\frac{f^k(x_0)}{k!}(x-x_0{)}^k=\sum _{k=0}^{N-1}p(k)(x-x_0{)}^k\text{\hspace{1em}(2-62)}$$

好得多？主要有三个原因。第一，切比雪夫逼近非常接近（虽不严格等于）“最优一致逼近”这一复杂问题的解，保证最大误差最小，即 $l_{\infty }$ 范数的最大值 $\max (f(x)-\hat{f}(x))\to \min$。第二，式(2-61)在 $M\ll N$ 时剪除多项式仍是最小/最大逼近——也就是直接以 $M$ 为目标去算，较短的和仍是切比雪夫逼近，这个”剪除特性”非常实用。第三，同等精度下式(2-61)所需系数更少。下面把这套方法用于三角函数、指数函数、对数函数和平方根函数等特殊函数的逼近。

2.10.2 三角函数的逼近

先看反正切函数：

$$f(x)=\arctan (x)\text{\hspace{1em}(2-63)}$$

定义范围为 $-1⩽x⩽1$；若要计算该区间之外的值，可用关系式 $\arctan (x)=0.5-\arctan (1/x)$ 换算(2-64)。Altera FPGA 的嵌入式乘法器基本规模为 9 位×9 位（8 位数据加 1 位符号）或 18 位×18 位（17 位加符号），所以下面按两种字长分别讨论。

图 2-50(a)给出精确值与 8 位量化逼近的比较，图 2-50(b)是误差曲线，呈现切比雪夫逼近典型的交错最小/最大行为。N=6 的逼近已经近乎完美；若减到 N=2 或 N=4，误差明显增大（参阅练习 2.26）。从图 2-50(d)看，8 位精度取 N=6 就够了；从图 2-50(c)还能看出所有偶系数均为 0，这是因为 $\arctan (x)$ 是关于 $x=0$ 奇对称的函数。待实现函数因此化简为：

$$\begin{array}{l}f(x)={\sum }_{k=0}^{N-1}c(k)T_k(x)\\ f(x)=c(1)T_1(x)+c(3)T_3(x)+c(5)T_5(x)\\ f(x)=0.8284T_1(x)-0.0475T_3(x)+0.0055T_5(x)\end{array}\text{\hspace{1em}(2-65)}$$

(a) 全精度和 8 位量化逼近的比较

(b) $x\in [-1,1]$ 的量化逼近的误差

(c) 切比雪夫、切比雪夫多项式和泰勒多项式系数

(d) 三种剪除多项式的误差
图 2-50 反正切函数逼近

求值时不必把 $T_n(x)$ 展开成幂函数再算，更高效的是用式(2-60)的迭代规则，即著名的切比雪夫（Clenshaw）递推公式：

$$\begin{array}{ll}d(N)=d(N+1)=0& \\ d(k)=2xd(k+1)-d(k+2)+c(k)& k=N-1,N-2,...,1\\ f(x)=d(0)=xd(1)-d(2)+c(0)& \end{array}\text{\hspace{1em}(2-66)}$$

对 $N=6$ 且偶系数为 0 的情形，式(2-66)化简为：

$$\begin{array}{l}d(5)=c(5)\\ d(4)=2xc(5)\\ d(3)=2xd(4)-d(5)+c(3)\\ d(2)=2xd(3)-d(4)\\ d(1)=2xd(2)-d(3)+c(1)\\ f(x)=xd(1)-d(2)\end{array}\text{\hspace{1em}(2-67)}$$

接下来用 HDL 实现这个逼近。

例 2.26 $\arctan (x)$ 函数的逼近

用 9 位×9 位嵌入式乘法器实现时，$x$ 的定义范围是 $-1⩽x⩽1$，故采用 1.8 格式小数定点表示：仿真中小数被映射为整数 $-256⩽x⩽256$。切比雪夫系数都小于 1，可用同样格式量化：

$$c(1)=0.8284=212/256,c(3)=-0.0475=-12/256,c(5)=0.0055=1/256\text{\hspace{1em}(2-68)}$$

原书 VHDL 代码改写为 Verilog，给出使用到 N=6 的多项式对应项的 $\arctan (x)$ 逼近：

// 原书 VHDL 改写为 Verilog：Clenshaw 递推实现 arctan(x) 切比雪夫逼近
module arctan (
    input              clk,     // 系统时钟
    input              reset,   // 异步复位
    input  signed [8:0] x_in,   // 系统输入（1.8 格式）
    output reg signed [8:0] d1, d2, d3, d4, d5, // 辅助递推量（测试观测）
    output reg signed [8:0] f_out // 系统输出
);
    reg  signed [8:0] x, f;
    // 8 位精度的切比雪夫系数（1.8 格式）
    localparam signed [8:0] c1 = 212, c3 = -12, c5 = 1;
 
    // 输入/输出寄存
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            x <= 0; f_out <= 0;
        end else begin
            x     <= x_in;
            f_out <= f;
        end
    end
 
    // 积之和计算：Clenshaw 递推公式
    always @(*) begin
        d5 = c5;
        d4 = (x * d5) / 128;
        d3 = (x * d4) / 128 - d5 + c3;
        d2 = (x * d3) / 128 - d4;
        d1 = (x * d2) / 128 - d3 + c1;
        f  = (x * d1) / 256 - d2;  // 最后一步不同
    end
endmodule

代码要点：递推中的”乘 2”用 x * d(k) 配合除以 128 完成——因为 $x$ 是 1.8 格式（实际值为整数表示的 1/256），$2x\cdot d(k)$ 的正确缩放正是先乘后除 $128$；而最后一步 $f(x)=xd(1)-d(2)$ 的缩放因子是 256，与式(2-67)一一对应。中间变量 $d(k)$ 全部引出到端口，便于观测。这一设计使用 106 个 LE 和 3 个嵌入式乘法器，TimeQuest 缓慢 85C 模型下 Registered Performance 为 Fmax=32.71 MHz。对比 FLEX 与 Cyclone 的综合数据可以得出结论：嵌入式乘法器节省了大量 LE。

图 2-51 给出了 5 个不同输入值的仿真结果（见表 2-21）：

表 2-21 arctan(x)函数对应的仿真结果

x	f(x)=arctan(x)	$\hat{f}(x)$	误差的绝对值	有效位
-1.0	-0.7854	-201/256=-0.7852	0.0053	7.6
-0.5	-0.4636	-118/256=-0.4609	0.0027	7.4
0	0.0	0	0	—
0.5	0.4636	118/256=0.4609	0.0027	7.4
1.0	0.7854	201/256=0.7852	0.0053	7.6

注意：由于 I/O 寄存器的存在，输出值延迟一个时钟周期出现。

图 2-51 arctan(x)函数逼近的 VHDL 仿真结果：x = -1 = -256/256、x = -0.5 = -128/256、x = 0、x = 0.5 = 128/256、x = 1 ≈ 255/256

如果上述精度不够，可以增加系数。16 位精度的奇数切比雪夫系数为：

$$c(2k+1)=(0.82842712,-0.04737854,0.00487733,-0.00059776,0.00008001,-0.00001282)\text{\hspace{1em}(2-71)}$$

与之对照，泰勒级数系数是：

$$\begin{array}{rl}& \arctan (x)=x-\frac{x^3}{3}+\frac{x^5}{5}+\cdots +(-1{)}^k\frac{x^{2k+1}}{2k+1}\\ & p(2k+1)=(1,-0.\overline{3},0.2,-0.14285714,0.\overline{1},-0.0\overline{9})\end{array}\text{\hspace{1em}(2-72)}$$

比较可见泰勒系数收敛慢得多——这正是切比雪夫逼近省乘积项的优势所在。

再看两个更常用的函数 $\sin (x)$ 和 $\cos (x)$。它们有个小问题：多项式逼近通常只在第一象限 $0⩽x⩽\pi /2$ 建立，其他象限的值通过对称关系换算：

$$\sin (x)=-\sin (-x)\sin (x)=\sin (\pi /2-x)\text{\hspace{1em}(2-73)}$$ $$\cos (x)=\cos (-x)\cos (x)=-\cos (\pi /2-x)\text{\hspace{1em}(2-74)}$$

常见做法是把数据规格化为 $f(x)=\sin (x\pi /2)$ 或角度形式 $0^0⩽x⩽90^0$。图 2-52(a)给出准确值和 16 位量化逼近的比较，图 2-52(b)是误差；图 2-53 画出了 $\cos (x\pi /2)$ 的同样数据。这里出现一个新问题：切比雪夫多项式只在 $x\in [-1,1]$ 定义，若目标函数定义在别的区间怎么办？办法很简单——对输入做一次线性变换。设 $f(y)$ 的定义域是 $y\in [a,b]$，用变量替换：

$$y=\frac{2x-b-a}{b-a}\text{\hspace{1em}(2-75)}$$

把定义域映射到 $[-1,1]$ 后再做切比雪夫逼近即可。

(a) 全精度和 8 位量化逼近的比较

(b) $x\in [0,1]$ 的量化逼近的误差

(c) 切比雪夫、切比雪夫多项式和泰勒多项式的系数

(d) 三种剪除多项式的误差

图2-52 正弦函数的逼近

(a) 全精度和 16 位量化逼近的比较

(b) $x\in [0,1]$ 的量化逼近的误差

(c) 切比雪夫、切比雪夫多项式和泰勒多项式的系数

图 2-53 余弦函数的逼近

(d) 三种剪除多项式的误差

例如 $\sin (x\pi /2)$ 中 $x$ 的定义域为 $x\in [0,1]$，即 $a=0$、$b=1$，代入得 $y\in [(2\times 0-1-0)/(1-0),(2\times 1-1-0)/(1-0)]=[-1,1]$，正好是切比雪夫逼近所需的区间。若用角度表示，则 $a=0$、$b=90$，映射为 $y=(2x-90)/90$。

最后一个问题是多项式到底怎么算：是走 Clenshaw 递推公式(2-66)，还是先展开成普通幂多项式再算（后者每次迭代加法更少）：

$$f(x)=\sum _{k=0}^{N-1}p(k)x^k\text{\hspace{1em}(2-76)}$$

或更好的 Horner 公式：

$$\begin{array}{ll}s(N-1)=p(N-1)& \\ s(k)=s(k+1)\times x+p(k)& k=N-2,\dots ,0\\ f=s(0)& \end{array}\text{\hspace{1em}(2-77)}$$

把切比雪夫函数(2-59)代入式(2-61)当然能得到幂系数（因为 $T_n(x)$ 不会产生高于 $x^n$ 的项），但这样做有一个严重缺陷：会失去剪除特性。也就是说，如果直接在式(2-76)中把项数减到 $M<N$，剪除后的多项式就不再是 $l_{\infty }$ 优化的了。图 2-50(d)清楚地展示了这一点：用满 6 项时切比雪夫与对应多项式逼近精度相同；一旦剪除，直接从完整切比雪夫系数换来的幂多项式精度远低于同长度的剪除切比雪夫逼近，甚至不比泰勒好多少。正确的做法是：想缩短到 $M<N$，就先重新求长度为 $M$ 的切比雪夫逼近，再从中计算多项式系数 $g(k)$。以 8 位与 16 位 $\arctan (x)$ 对比：把式(2-59)代入系数(2-71)得到 16 位情况的奇系数

$$g(2k+1)=(0.99999483,-0.33295711,0.19534659,-0.12044859,0.05658999,-0.01313038)\text{\hspace{1em}(2-78)}$$

而使用式(2-68)长度 N=6 的逼近时奇系数为：

$$g(2k+1)=(0.9982,-0.2993,0.0876)\text{\hspace{1em}(2-79)}$$

两者剪除的切比雪夫系数相同，但多项式系数差别明显——例如 $g(5)$ 相差两倍。所以”先按目标长度求切比雪夫逼近，再展开”这一步不能省。

把上述流程总结成算法。

算法 2.27：切比雪夫函数逼近。

(1) 确定系数 N。

(2) 利用式(2-75)将变量从 x 变换到 y。

(3) 确定以 $y$ 为自变量的切比雪夫逼近。

(4) 利用 Clenshaw 递推公式确定直接多项式系数 $g(k)$。

(5) 求映射 $y$ 的逆。

按这 5 步计算 $\sin (x\pi /2)$（$x\in [0,1]$，4 个非零系数），得到足够 16 位量化的多项式：

$$\begin{array}{rl}f(x)& =\sin (x\pi /2)\\ & =1.57035062x+0.00508719x^2-0.66666099x^3+0.03610310x^4+0.05512166x^5\\ & =(51457x+167x^2-21845x^3+1183x^4+1806x^5)/32768\end{array}$$

注意第一个系数大于 1，需要适当缩放。这与泰勒逼近

$$\sin \left(\frac{x\pi }{2}\right)=\frac{x\pi }{2}-\frac{1}{3!}{\left(\frac{x\pi }{2}\right)}^3+\frac{1}{5!}{\left(\frac{x\pi }{2}\right)}^5+\cdots +\frac{(-1{)}^k}{(2k+1)!}{\left(\frac{x\pi }{2}\right)}^{2k+1}$$

截然不同，图 2-52(c)给出了图形说明。8 位量化使用：

$$\begin{array}{rl}f(x)& =\sin (x\pi /2)=1.5647x+0.0493x^2-0.7890x^3+0.1748x^4\\ & =(200x+6x^2-101x^3+22x^4)/128\end{array}\text{\hspace{1em}(2-80)}$$

按理说奇对称函数的偶系数应为 0，但这里并非如此，因为逼近区间只是 $x\in [0,1]$ 而非对称区间。$\cos (x)$ 可通过关系

$$\cos \left(\frac{x\pi }{2}\right)=\sin ((x+1)\frac{\pi }{2})\text{\hspace{1em}(2-81)}$$

导出，或直接做切比雪夫逼近。对 $x\in [0,1]$、4 个非零系数，16 位量化的多项式为：

$$\begin{array}{rl}f(x)& =\cos (\frac{x\pi }{2})\\ & =1.00000780-0.00056273x-1.22706059x^2\\ & -0.02896799x^3+0.31171138x^4-0.05512166x^5\\ & =(32768-18x+40208x^2-949x^3+10214x^4-1806x^5)/32768\end{array}$$

8 位量化使用：

$$\begin{array}{rl}f(x)& =\cos (\frac{x\pi }{2})\\ & =0.9999+0.0046x-1.2690x^2+0.0898x^3+0.1748x^4\\ & =(128+x-162x^2+11x^3+22x^4)/128\end{array}\text{\hspace{1em}(2-82)}$$

泰勒逼近的系数再次截然不同：

$$\cos \left(\frac{x\pi }{2}\right)=1-\frac{1}{2!}{\left(\frac{x\pi }{2}\right)}^2+\frac{1}{4!}{\left(\frac{x\pi }{2}\right)}^4+\cdots +\frac{(-1{)}^k}{(2k)!}{\left(\frac{x\pi }{2}\right)}^{2k}$$

图 2-52(c)给出了系数的图形对比。最能说明问题的是图 2-52(d)：同样 6 次幂项数的泰勒逼近只有 6 位精度，而切比雪夫逼近达到 16 位精度——这就是切比雪夫方法在 FPGA 上节省乘法器与流水线级数的直接体现。

2.10.3 指数函数和对数函数的逼近

在 FPGA 里计算 ${\mathrm{e}}^x$、$\ln x$、$\sqrt{x}$ 这类”超越函数”，最常用的思路是多项式逼近：预先在 PC 上算好一组系数，硬件上只需要乘法器和加法器就能算出结果。这样做的原因很简单——FPGA 擅长乘加运算，但完全不擅长迭代、查表以外的复杂算法。多项式逼近的本质是用一个有限阶的多项式 $\hat{f}(x)$ 去替代真实函数 $f(x)$，只要在感兴趣的自变量区间内误差足够小，硬件就算”算对了”。

构造多项式系数有两种经典办法。第一种是泰勒（Taylor）级数：在某个展开点求各阶导数。它的好处是公式固定、人人会推；缺点是在展开点附近精度极好，但离开展开点后误差迅速增大。第二种是切比雪夫（Chebyshev）逼近：让误差在整个区间上”摊平”，最大误差明显更小。原书的做法是：先用泰勒级数起步，再把系数换算到切比雪夫基上微调，最后量化到定点整数。

指数函数是泰勒逼近收敛速度相当快的少数函数之一，其泰勒逼近为：

$$f(x)=e^x=1+\frac{x}{1!}+\frac{x^2}{2!}+\cdots +\frac{x^k}{k!},0⩽x⩽1\text{\hspace{1em}(2-84)}$$

对 16 位量化，用切比雪夫系数计算的多项式应为：

$$\begin{array}{rl}f(x)=e^x& =1.00002494+0.99875705x+0.50977984x^2\\ & +0.14027504x^3+0.06941551x^4\\ & =(32769+32727x+16704x^2+4597x^3+2275x^4)/32768\end{array}$$

注意两个细节：第一，只有阶数达到 $x^4$ 的项才需要 16 位精度——低阶项对结果影响大，量化必须更精细；第二，从图 2-54(c) 可以看到，用切比雪夫方法反算回来的”泰勒系数”与直接量化的多项式系数非常相似，这说明对指数函数而言两种方法差别不大。

如果 8 位加上符号位就满足精度要求（比如 AGC 里的粗略幅度估计），可以改用下面的低阶逼近：

$$f(x)=e^x=1.0077+0.8634x+0.8373x^2=(129+111x+107x^2)/128\text{\hspace{1em}(2-85)}$$

图 2-54(a) 全精度和 16 位量化逼近的比较

图 2-54(b) $x\in [0,1]$ 的量化逼近的误差

图 2-54(c) 切比雪夫、切比雪夫多项式和泰勒多项式的系数

图 2-54(d) 三种剪除多项式的误差

图 2-54 指数函数 $f(x)=\exp (x)$ 的逼近

式(2-85)里有一个系数 $c(0)>1.0$，因此整体要乘一个缩放因子 128，把系数放大成整数存进硬件。这带来一个必须遵守的前提：输入 $x$ 要先被缩放（除以同一因子），保证 $0⩽x⩽1$。不用会怎样？多项式在区间外可能发散，输出会完全错误。

如果 $x$ 不在 $[0,1]$ 内，可以利用恒等式

$$e^{sx}={\left(e^x\right)}^s\text{\hspace{1em}(2-86)}$$

当 $s=2^k$ 是 2 的幂时，”$s$ 次方”就是一串平方运算：先算 ${\mathrm{e}}^x$，再连续平方 $k$ 次。平方在硬件里就是一次乘法，代价很低。

对于负指数，可以利用

$$e^{-x}=\frac{1}{e^x}\text{\hspace{1em}(2-87)}$$

但除法在 FPGA 上太昂贵，更实际的做法是直接为 $f(x)=e^{-x}$ 单独构造一套逼近：在式(2-84)中每两项改变一次符号即可得到泰勒版本；切比雪夫版本则需要对系数再做一些小改动。16 位切比雪夫多项式逼近为：

$$\begin{array}{rl}f(x)=e^{-x}& =0.99998916-0.99945630x+0.49556967x^2\\ & -0.15375046x^3+0.02553654x^4\\ & =(65535-655500x+32478x^2-10076x^3+1674x^4)/65536\end{array}$$

其中 $x\in [0,1]$。注意这里所有系数都小于 1，所以缩放因子可以选 65536——比式(2-85)的 128 大两倍多，等于”白赚”了精度。从图 2-55(d) 可以得出结论：3 个或 5 个系数分别对应 8 位或 16 位精度。8 位量化使用：

$$f(x)=e^{-x}=0.9964-0.9337x+0.3080x^2=(255-239x+79x^2)/256\text{\hspace{1em}(2-88)}$$

图 2-55(a) 全精度和 16 位量化逼近的比较

图 2-55(b) $x\in [0,1]$ 的量化逼近误差

图 2-55(c) 切比雪夫、切比雪夫多项式和泰勒多项式的系数

图 2-55(d) 三种剪除多项式的误差

图 2-55 负指数函数 $f(x)=\exp (-x)$ 的逼近

指数函数的反函数是对数函数，其自变量的典型定义域是 $[1,2]$。为了统一处理，通常写成 $f(x)=\ln (1+x)$，这样 $x$ 的范围变成 $0⩽x⩽1$。图 2-56(a) 给出这一区间的准确值和 16 位量化逼近：$N=6$ 的逼近几乎就是完美的逼近；如果只用 $N=2$ 或 $N=3$ 个系数，就会出现明显误差（参阅原书练习 2.29）。

图 2-56(a) 全精度和 16 位量化逼近的比较

图 2-56(b) $x\in [0,1]$ 的量化逼近的误差

图 2-56(c) 切比雪夫、切比雪夫多项式和泰勒多项式的系数

图 2-56(d) 三种剪除多项式的误差

图 2-56 自然对数函数 $f(x)=\ln (1+x)$ 的逼近

从泰勒级数

$$f(x)=\ln (1+x)=x-\frac{x^2}{2}+\frac{x^3}{3}+\cdots +\frac{(-1{)}^{k+1}{x}^k}{k}$$

的分母（线性增长而不是阶乘增长）就可以看出，它不再像指数函数那样快速收敛。图 2-56(d) 显示 16 位切比雪夫逼近收敛快得多：16 位精度只需要 6 个系数，而用 6 个泰勒系数连 4 位精度都得不到。用切比雪夫系数计算的 16 位多项式量化为：

$$\begin{array}{rl}f(x)=\ln (1+x)& =0.00001145+0.99916640x-0.48969909x^2\\ & +0.28382318x^3-0.12995720x^4+0.02980877x^5\\ & =(1+65481x-32093x^2+18601x^3-8517x^4+1954x^5)/65536\end{array}$$

只有阶数达到 $x^5$ 的项才需要 16 位精度。从图 2-56(c) 还可以看到：由切比雪夫逼近反算的泰勒系数和直接的多项式系数只有前 3 个接近——这正是对数函数泰勒收敛慢的体现。

接下来用 HDL 实现 $\ln (1+x)$ 函数的逼近。

例 2.28 $\ln (1+x)$ 函数的逼近

设计的第一步永远是定点格式规划。如果用 18 位×18 位的嵌入式乘法器实现，就必须考虑 $x$ 的定义范围 $0⩽x⩽1$，因而采用 2.16 格式的小数整型数（2 位整数部分、16 位小数部分）：$x=1$ 可以准确地表示为 $2^{16}$，再加一位额外的保护位保证运算过程中不出现溢出。因为切比雪夫系数都小于 1，系数也可以用同样的格式。

原书 VHDL 代码改写为等价的 Verilog HDL（可综合风格），实现 6 个系数的 $\ln (1+x)$ 逼近：

// 原书 VHDL 改写为 Verilog：ln(1+x) 的切比雪夫多项式逼近
module ln #(
    parameter integer N = 5                    // 系数个数 - 1
) (
    input  wire                 clk,           // 系统时钟
    input  wire                 reset,         // 异步复位
    input  wire signed [17:0]   x_in,          // 系统输入，2.16 格式
    output reg  signed [17:0]   f_out          // 系统输出
);
 
    // 16 位精度的多项式系数：
    // f(x) = (1 + 65481x - 32093x^2 + 18601x^3 - 8517x^4 + 1954x^5)/65536
    function signed [17:0] coeff(input integer k);
        case (k)
            0: coeff = 18'sd1;
            1: coeff = 18'sd65481;
            2: coeff = -18'sd32093;
            3: coeff = 18'sd18601;
            4: coeff = -18'sd8517;
            default: coeff = 18'sd1954;
        endcase
    endfunction
 
    reg signed [17:0] x;                       // 输入暂存
    reg signed [17:0] s [0:5];                 // Horner 中间结果
    integer k;
 
    // I/O 寄存器：输入、输出各打一拍
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            x     <= 18'sd0;
            f_out <= 18'sd0;
        end else begin
            x     <= x_in;
            f_out <= s[0];
        end
    end
 
    // 用积之和（Horner 方案）计算切比雪夫逼近
    // f = p0 + x*(p1 + x*(p2 + ... + x*p5))
    always @* begin
        s[N] = coeff(N);
        for (k = N-1; k >= 0; k = k - 1) begin
            s[k] = (x * s[k+1]) >>> 16 + coeff(k);   // x*s/65536，36 位乘积取高位
        end
    end
 
endmodule

代码结构分两块：第一个 always 块为输入和输出指定寄存器；第二个组合 always 块执行 Horner 递推 $s(k)=x\cdot s(k+1)/2^{16}+p(k)$。这里有一个初学者容易踩的坑：$x\cdot s(k+1)$ 的乘积是 36 位，超过了 Verilog 整数表达式的 32 位有效范围吗？在 Verilog 中 $x * s[k+1]$ 会自动生成 36 位宽的中间结果，右移 16 位完成除以 65536 的缩放，再截取 18 位赋给 $s(k)$，这正是原书用标准逻辑向量手动处理 36 位乘积的等价写法。这一设计使用了 88 个 LE 和 10 个 9 位×9 位的嵌入式乘法器（相当于把 18 位×18 位嵌入式乘法器的数量减半），使用 TimeQuest 缓慢 85C 模型时序性能为 Fmax=29.2MHz。

图 2-57 给出了仿真结果，覆盖 5 个输入值：$x=0$、$x=0.25=16384/65536$、$x=0.5=32768/65536$、$x=0.75=49152/65536$、$x=1.0=65536$。

图 2-57 $\ln (1+x)$ 函数逼近的仿真结果：$x=0$、$x=0.25=16384/65536$、$x=0.5=32768/65536$、$x=0.75=49152/65536$、$x=1.0=65536$

表 2-22 5 个不同输入值对应的仿真结果（结果均为 2.16 格式定点数）：

$x$	$f(x)=\ln (1+x)$	$\hat{f}(x)$	误差的绝对值	有效位
0	0	$1/2^{16}$	$1.52\times 10^{-5}$	16
0.25	$14623.9/2^{16}$	$14624/2^{16}$	$4.39\times 10^{-6}$	17.8
0.5	$26572.6/2^{16}$	$26572/2^{16}$	$2.11\times 10^{-5}$	15.3
0.75	$36675.0/2^{16}$	$36675/2^{16}$	$5.38\times 10^{-7}$	20.8
1.0	$45426.1/2^{16}$	$45427/2^{16}$	$1.99\times 10^{-5}$	15.6

注：由于 I/O 寄存器的原因，输出值延迟一个时钟周期出现。

把本设计的多项式代码与例 2.26 的切比雪夫递推公式对比，会发现本设计少用了一个加法器——直接多项式求值比”系数递推”方案更省资源。

如果 8 位加符号位精度满足要求，用下面的公式：

$$f(x)=\ln (1+x)=0.0006+0.9813x-0.3942x^2+0.1058x^3=(251x-101x^2+27x^3)/256\text{\hspace{1em}(2-89)}$$

由于没有系数大于 1.0，因此可以选择 256 作为缩放因子。

如果自变量 $x$ 不在有效范围 $[0,1]$ 内，可以用代数操作把它”拉回来”。令 $y=sx=2^{k}x$：

$$\ln (sx)=\ln (s)+\ln (x)=k\times \ln (2)+\ln (x)\text{\hspace{1em}(2-90)}$$

也就是先用 2 的幂把数规格化到有效区间，算完后再把 $k\times \ln 2$ 加回去。一次乘法加一次加法，代价可以忽略。

如果需要更换对数的底（比如以 10 为底），理论上用换底公式：

$${\log }_a(x)=\ln (x)/\ln (a)\text{\hspace{1em}(2-91)}$$

即只实现一个自然对数就能推出其他底。但除法运算过于昂贵，得不偿失；更好的办法是直接为 $\lg$ 计算另一组切比雪夫系数。16 位精度的以 10 为底对数逼近为：

$$\begin{array}{rl}f(x)& =\lg (1+x)\\ & =0.00000497+0.43393245x-0.21267361x^2\\ & +0.12326284x^3-0.05643969x^4+0.01294578x^5\\ & =(28438x-13938x^2+8078x^3-3699x^4+848x^5)/65536\end{array}$$

其中 $x\in [0,1]$。8 位量化则使用：

$$\begin{array}{rl}f(x)& =\lg (1+x)\\ & =0.0002+0.4262x-0.1712x^2+0.0460x^3\\ & =(109x-44x^2+12x^3)/256\end{array}\text{\hspace{1em}(2-92)}$$

这一公式只用了三个非零系数，如图 2-58(d) 所示。

图 2-58(a) 全精度和 8 位量化逼近的比较

图 2-58(b) $x\in [0,1]$ 的量化逼近的误差

图 2-58(c) 切比雪夫、切比雪夫多项式和泰勒多项式的系数

图 2-58(d) 三种剪除多项式的误差

图 2-58 以 10 为底的对数函数 $f(x)=\lg (x)$ 的逼近

2.10.4 平方根函数的逼近

平方根函数有个特殊困难：不能计算它在 $x_0=0$ 点附近的泰勒逼近，因为所有导数 $f^n(x_0)$ 要么是 0，要么是更糟的 $1/0$（无穷大）。不过可以改在 $x_0=1$ 点附近展开：

$$\begin{array}{rl}f(x)& =\sqrt{x}\\ & =\frac{(x-1{)}^0}{0!}+0.5\frac{(x-1{)}^1}{1!}-\frac{0.5^2}{2!}(x-1{)}^2+\frac{0.5^2\times 1.5}{3!}(x-1{)}^3-\dots \\ & =1+\frac{(x-1)}{2}-\frac{(x-1{)}^2}{8}+\frac{(x-1{)}^3}{16}-\frac{5}{128}(x-1{)}^4+\dots \end{array}$$

图 2-59(c) 以图形方式给出了这些系数和等价的切比雪夫系数。用切比雪夫系数计算的多项式量化应采用：

$$\begin{array}{rl}f(x)& =\sqrt{x}\\ & =0.23080201+1.29086721x-0.88893983x^2\\ & +0.48257525x^3-0.11530993x^4\\ & =(7563+42299x-29129x^2+15813x^3-3778x^4)/32768\end{array}$$

图 2-59(a) 全精度和 16 位量化逼近的比较

图 2-59(b) $x\in [0.5,1)$ 的量化逼近的误差

图 2-59(c) 切比雪夫、切比雪夫多项式和泰勒多项式的系数

图 2-59(d) 三种剪除多项式的误差

图 2-59 平方根函数 $f(x)=\sqrt{x}$ 的逼近

这一次自变量的有效定义域是 $x\in [0.5,1)$。只有阶数达到 $x^4$ 的项才需要 16 位精度。从图 2-59(c) 可以看到，由切比雪夫逼近反算的泰勒系数和多项式系数并不接近——平方根的泰勒级数收敛也不快。图 2-59(a) 中 $N=5$ 的逼近几乎是完美逼近；若减少系数（$N=2$ 或 $N=3$）误差会显著增大，参阅原书练习 2.30。

剩下的唯一问题是：如何处理范围 $0.5⩽x<1$ 之外的数？做法是把自变量 $y=sx$ 拆出一个 2 的幂缩放因子 $s=2^k$，让剩余部分落回有效区间。缩放因子的平方根用下式实现：

$$\sqrt{s}=\sqrt{2^k}=\{\begin{array}{ll}{2}^{k/2}& k\text{是偶数}\\ \sqrt{2}\times 2^{(k-1)/2}& k\text{是奇数}\end{array}\text{\hspace{1em}(2-94)}$$

偶数 $k$ 时只需移位；奇数 $k$ 时多乘一个常数 $\sqrt{2}$——在硬件里这就是移位加加减法，完全可以预先算死。

例 2.29 平方根函数的逼近

实现方式有两种：用 $N$ 个嵌入式 18 位×18 位乘法器以并行方式实现，或者构造一个 FSM（有限状态机）以迭代方式复用一个乘法器分步计算。并行快但费资源，FSM 省资源但需要多个时钟周期。原书练习 2.20 和 2.21 还给出了其他 FSM 设计示例。

设计的第一步仍是缩放：数据要保证自由溢出能够得到处理；此外需要前缩放和后缩放，保证 $x$ 落在有效范围 $0.5⩽x<1$ 内。因而采用 3.15 格式的小数整型数，加一位额外保护位防止溢出，且 $x=1$ 可以准确地表示为 $2^{15}$。因为切比雪夫系数都小于 2，系数也用同样的格式。

原书 VHDL 代码改写为等价的 Verilog HDL（可综合风格），使用 $N=5$ 个系数的 $\sqrt{x}$ 逼近：

// 原书 VHDL 改写为 Verilog：sqrt(x) 的 FSM + ALU 迭代逼近
module sqrt (
    input  wire                 clk,      // 系统时钟
    input  wire                 reset,    // 异步复位
    output wire  [1:0]          count_o,  // 左移阶段计数器
    input  wire signed [17:0]   x_in,     // 系统输入，3.15 格式
    output wire  signed [17:0]  pre_o,    // 前缩放因子
    output wire  signed [17:0]  x_o,      // 规格化后的输入
    output wire  signed [17:0]  post_o,   // 后缩放因子
    output wire  signed [2:0]   ind_o,    // 系数索引
    output wire  signed [17:0]  imm_o,    // ALU 预装值
    output wire  signed [17:0]  a_o,      // ALU 因子
    output wire  signed [17:0]  f_o,      // ALU 输出
    output reg   signed [17:0]  f_out     // 系统输出
);
 
    // 状态编码
    localparam [2:0] START=3'd0, LEFTSHIFT=3'd1, SOP=3'd2,
                     RIGHTSHIFT=3'd3, DONE=3'd4;
    // ALU 操作码
    localparam [2:0] LOAD=3'd0, MAC=3'd1, SCALE=3'd2,
                     DENORM=3'd3, NOP=3'd4;
 
    reg [2:0]  s;                  // 当前状态
    reg [2:0]  op;                 // 当前 ALU 操作
    reg [1:0]  count;              // 规格化计数
    reg signed [2:0] ind;          // 系数索引，-1..4
 
    reg signed [17:0] x;           // 规格化后的自变量
    reg signed [17:0] a, f, imm;   // ALU 数据通路
 
    // 16 位精度的切比雪夫多项式系数：
    // f(x) = (7563 + 42299x - 29129x^2 + 15813x^3 - 3778x^4)/32768
    function signed [17:0] coeff(input integer k);
        case (k)
            0: coeff = 18'sd7563;
            1: coeff = 18'sd42299;
            2: coeff = -18'sd29129;
            3: coeff = 18'sd15813;
            default: coeff = -18'sd3778;
        endcase
    endfunction
 
    reg signed [17:0] pre, post;   // 前/后缩放因子
 
    // ---------------- FSM：控制部分 ----------------
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            s      <= START;
            f_out  <= 18'sd0;
            ind    <= 3'sd0;
            count  <= 2'd0;
            op     <= NOP;
        end else begin
            case (s)
                START: begin                    // 初始化，装入自变量
                    s     <= LEFTSHIFT;
                    ind   <= 3'sd4;
                    imm   <= x_in;              // 自变量装入 ALU
                    op    <= LOAD;
                    count <= 2'd0;
                end
                LEFTSHIFT: begin                // 规格化到 0.5..1.0
                    count <= count + 1'b1;
                    a     <= pre;
                    op    <= SCALE;
                    imm   <= coeff(4);
                    if (count == 2'd2)
                        op <= NOP;
                    if (count == 2'd3) begin    // 规格化完成？
                        s  <= SOP;
                        op <= LOAD;
                        x  <= f;
                    end
                end
                SOP: begin                      // 主处理：乘累加
                    ind   <= ind - 1'b1;
                    a     <= x;
                    if (ind == 3'sd0) begin     // SOP 完成？（ind 减 1 后为 -1）
                        s   <= RIGHTSHIFT;
                        op  <= DENORM;
                        a   <= post;
                    end else begin
                        imm <= coeff(ind - 1);
                        op  <= MAC;
                    end
                end
                RIGHTSHIFT: begin               // 反向规格化
                    s  <= DONE;
                    op <= NOP;
                end
                DONE: begin                     // 结果送输出寄存器
                    f_out <= f;
                    op    <= NOP;
                    s     <= START;             // 准备下一次计算
                end
                default: s <= START;
            endcase
        end
    end
 
    // ---------------- ALU：运算部分 ----------------
    always @(posedge clk or posedge reset) begin
        if (reset)
            f <= 18'sd0;
        else begin
            case (op)
                LOAD:   f <= imm;                        // 预装
                MAC:    f <= ((a * f) >>> 15) + imm;     // f = a*f/2^15 + imm
                SCALE:  f <= a * f;                      // 规格化乘法
                DENORM: f <= (a * f) >>> 15;             // 反向规格化
                default: f <= f;                         // nop
            endcase
        end
    end
 
    // ---------------- EXP：前/后缩放因子计算 ----------------
    // 根据式(2-94)：偶数 k 用 2^(k/2)，奇数 k 再乘 sqrt(2) 的 CSD 编码
    integer k;
    reg signed [17:0] pr, po;
    always @* begin
        pr  = 18'sd16384;                  // 2^14 起步
        pre = 18'sd0;
        for (k = 0; k <= 15; k = k + 1) begin
            if (x_in[k]) pre = pr;         // 最高有效位决定缩放
            pr = pr >>> 1;
        end
        po = 18'sd1;
        for (k = 0; k <= 7; k = k + 1) begin
            if (x_in[2*k])                 // 偶数位：sqrt(2^2k)=2^k
                po = 18'sd256 <<< k;
            if (x_in[2*k+1]) begin         // 奇数位：带 sqrt(2) 因子
                // sqrt(2) 的 CSD 误差 = 0.0000208 = 15.55 有效位
                po = (18'sd512 <<< k) - (18'sd128 <<< k) - (18'sd32 <<< k)
                   + (18'sd8 <<< k) + (18'sd2 <<< k) + ((18'sd1 <<< k) >>> 5);
            end
        end
        post = po;
    end
 
    // 测试信号输出
    assign count_o = count;
    assign ind_o   = ind;
    assign pre_o   = pre;
    assign post_o  = post;
    assign x_o     = x;
    assign imm_o   = imm;
    assign a_o     = a;
    assign f_o     = f;
 
endmodule

这段代码包含三个主要模块：控制部分放在 FSM 模块中，算法部分放在 ALU 和 EXP 模块中。

• FSM 模块负责控制计算顺序，并把数据放到 ALU 和 EXP 模块的正确寄存器上。START 状态初始化数据并把输入装入 ALU；LEFTSHIFT 状态把输入规格化，使 $x$ 落在 $[0.5,1]$ 范围内；SOP 状态是主处理步骤，ALU 在这里以乘累加方式计算多项式；随后进入 RIGHTSHIFT 状态做反向规格化（规格化的逆过程）；最后在 DONE 状态把结果传输到输出寄存器，FSM 回到 START 为下一次平方根计算做准备。
• ALU 模块执行 Horner 公式 (2-77) 的 $f=a\times f+\mathrm{imm}$ 运算来计算多项式，会合成到一个 18 位×18 位的嵌入式乘法器（或按 Quartus II 的报告折算成两个 9 位×9 位嵌入式乘法器），外加一些加法和规格化逻辑。它具有典型 ALU 的形式：信号 op 决定当前操作，累加器寄存器 f 可以通过 imm 操作数预加载。
• EXP 模块根据式 (2-94) 计算前规格化和后规格化因子。$2^k$ 的奇数 $k$ 值所需的 $\sqrt{2}$ 因子通过 csd.exe 程序计算的 CSD 编码实现：$+1+0-1+0-1+0+1+0+1+0+0+0+0+0+1$，误差 0.0000208，相当于 15.55 有效位——全部用移位和加减法完成，不需要乘法器。

这一设计使用了 261 个 LE 和两个 9 位×9 位的嵌入式乘法器（相当于一个 18 位×18 位乘法器），TimeQuest 缓慢 85C 模型下时序性能为 Fmax=86.23MHz。

图 2-60 给出了输入 $x=0.75/8=0.0938=3072/2^{15}$ 的仿真结果，整个计算流程值得逐步走一遍：

1. 移动（规格化）阶段：输入 3072 通过前缩放因子 8 被规格化，得到 24576，落在 $x$ 的有效范围 $[0.5,1]$（等价于整数 $[16384,32768]$）内；
2. MAC 阶段：经过几次乘累加运算，得到 $f=\sqrt{0.75}\times 2^{15}=28378$；
3. 反向规格化阶段：通过后缩放因子 $\sqrt{2}\times 2^{13}=11585$ 反向规格化，最终结果 $f=\sqrt{0.75/8}\times 2^{15}=10032$ 被传送到输出寄存器。

图 2-60 $\sqrt{x}$ 函数逼近的仿真结果：$x=0.75/8=3072/32768$

如果 8 位加符号位精度满足要求，可以用下面的公式构造平方根：

$$f(x)=\sqrt{x}=0.3171+0.8801x-0.1977x^2=(81+225x-51x^2)/256\text{\hspace{1em}(2-95)}$$

由于没有系数大于 1.0，因此可以选择 256 作为缩放因子。

2.11 快速幅度逼近

在 FFT、图像处理，或使用 $x+jy$ 类型复数数据的非相干接收机自动增益控制（Automatic Gain Control, AGC）等应用中，经常需要幅度 $r=\sqrt{x^2+y^2}$ 的快速逼近。以 AGC 为例：幅度估计用来调整输入信号的增益，使信号既不会太小（否则量化噪声占比过大），也不会太大（否则算术溢出）。注意这里的要求是”快”和”够用”，而不是精确——精确的平方根反而浪费时间。

图像处理中的边缘检测是另一个例子：基于 $G_x$ 和 $G_y$ 方向上的梯度，我们只想知道总梯度 $\sqrt{G_x^2+G_y^2}$ 是否越过了阈值、该像素是否算边缘，同样不需要高精度。实践中甚至有形如 $r\approx |x|+|y|$ 的非常粗略的估计，称为幅度估计的零 $L_0$ 逼近。利用三角关系 $\sin (\varphi )+\cos (\varphi )=\sqrt{2}\sin (\varphi +\pi /4)$ 可以算出，这个 $L_0$ 估计的误差最大可达 40%。

CORDIC 算法或多项式逼近可以提供更高的精度，但会带来较长的时延和较高的资源消耗。如果要求尽量快速、又要比 $L_0$ 更精确，可以使用形如

$$r\approx \alpha \max (|x|,|y|)+\beta \min (|x|,|y|)$$

的最大/最小逼近。它的优点非常突出：没有乘法（$\alpha$、$\beta$ 选 2 的幂的近似值后只剩移位加法）、延迟低、资源少，并且逼近对称于 $45^{\circ }$ 线。因子 $\alpha$ 和 $\beta$ 可以针对 $L_1$（最小平均误差）或 $L_{\infty }$（最小最大误差）范数优化。表 2-21 给出了 $L_1$ 和 $L_{\infty }$ 的最优值和几个流行逼近值。

表 2-21 最大/最小幅度逼近的系数特性，有效位基于 $L_{\infty }$ 范数：

$\alpha$	$\beta$	$L_1$	$L_{\infty }$	有效位	备注
1	1	27.1%	41.4%	1.2	$L_0$
0.943	0.386	1.97%	6.0%	4.1	$L_1$ 最优值
1	1/4	3.17%	11.6%	3.1	$L_1$ 逼近值
0.962	0.396	2.45%	4.0%	4.6	$L_{\infty }$ 最优值
1	3/8	4.22%	6.8%	3.9	$L_{\infty }$ 逼近值

其中 (1, 0.375) 的 $L_{\infty }$ 逼近已被用于英特锡尔（Intersil）的 HSP50110 通信集成电路中。从表中”有效位”一列可以看到：即使在全系数精度下，这仍是粗略估计，幅度精度不会超过 5 比特——但对 AGC、边缘检测这类应用已经足够。图 2-61 显示了通过对 $\alpha /\beta$ 在合理范围内线性搜索优化后，5 个选项的计算幅度逼近结果。从实际实现的角度，$\beta =1/4$ 的选择最有趣：实现代价最低（除以 4 只需右移两位），且平均误差范数 $L_1$ 表现良好。

图 2-61 使用最大/最小方法的幅度逼近

例 2.30 幅度电路的 Verilog 设计

考虑 16 位幅度逼近电路。原书 VHDL 代码改写为等价的 Verilog HDL（可综合风格），采用 $\alpha =1$、$\beta =1/4$：

// 原书 VHDL 改写为 Verilog：最大/最小幅度逼近
// r = alpha*max(|x|,|y|) + beta*min(|x|,|y|)，取 alpha=1、beta=1/4
module magnitude (
    input  wire                 clk,    // 系统时钟
    input  wire                 reset,  // 异步复位
    input  wire signed [15:0]   x,      // 系统输入
    input  wire signed [15:0]   y,      // 系统输入
    output reg  signed [15:0]   r       // 系统输出
);
 
    reg signed [15:0] x_r, y_r;         // 输入寄存器
 
    // 输入数据存入寄存器
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            x_r <= 16'sd0;
            y_r <= 16'sd0;
        end else begin
            x_r <= x;
            y_r <= y;
        end
    end
 
    // 先取绝对值
    wire signed [15:0] ax = (x_r < 0) ? -x_r : x_r;
    wire signed [15:0] ay = (y_r < 0) ? -y_r : y_r;
 
    // 确定 max 与 min
    reg signed [15:0] ma, mi;
    always @* begin
        if (ax > ay) begin
            mi = ay;
            ma = ax;
        end else begin
            mi = ax;
            ma = ay;
        end
    end
 
    // 计算并寄存输出：r = alpha*max + beta*min
    always @(posedge clk or posedge reset) begin
        if (reset)
            r <= 16'sd0;
        else
            r <= ma + (mi >>> 2);       // beta=1/4 即右移两位
    end
 
endmodule

代码流程与原书完全对应：首先输入数据存入寄存器；然后取输入数据的绝对值；接着在一个组合逻辑块中比较出最大值和最小值，分配给 ma 和 mi；最后计算 $\alpha =1$、$\beta =0.25$ 的近似值并寄存输出。$\beta =1/4$ 意味着”除以 4”退化为算术右移两位，整个电路没有一个乘法器。该设计使用 96 个 LE、无嵌入式乘法器，时序性能 $F_{max}=119.59\text{MHz}$。输入输出端各加了一级寄存器，目的是测量电路的流水线性能。

图 2-62 给出了 16 位流水线幅度计算的仿真结果。测试选取角度为 0、$45^{\circ }$、$90^{\circ }$、… 的 8 个值。注意两点：其一，由于流水线寄存器的存在，计算的 $r$ 值相对输入延迟了几拍；其二，幅度值 $r=1000$ 是近似值，在 $45^{\circ }$ 方向（即 $x=y=1000\cos (\pi /4)\approx 707$）附近误差变得相当大——这正是表 2-21 中 $L_{\infty }$ 误差 11.6% 的直观体现：最大/最小逼近在轴向最准，在对角线方向偏差最大。只要应用能容忍百分之十几的误差，用 96 个 LE 换一个 119 MHz 的幅度估计器，这笔账是很划算的。

图 2-62 幅值电路的仿真结果

最后补一张原书练习部分的配图。练习 2.19 要求用多种方式实现一个 12 位桶式移位器（VHDL-1993 的 STD_LOGIC 不支持 SLL 指令），可通过如图 2-63 所示的仿真进行验证：为 data 和 result 提供输入和输出寄存器，而不为 distance 提供寄存器。

图 2-63 练习 2.19 中桶式移位器的测试平台

2.8 用基准设计把“综合选项”和“资源消耗”落到实处：PREP 基准 3/4 与乘法器实现

前面几节讲了乘法器的各种算法（分块、半分方形、差分半分方形、四分之一平方）和函数的多项式逼近。这一组素材对应教材第 2 章末尾的实验任务（习题 2.20 后半到 2.30），把它们串起来看，其实是在回答三个初学者最容易“听过但没做过”的问题：

1. 同一个设计，换一种综合选项（Speed / Balanced / Area），速度和资源会差多少？
2. 一个大 FSM（有限状态机）该怎么写、怎么验证？
3. 查表乘法器的 MIF 文件到底长什么样，表的内容怎么手工核对？

下面按小节重新组织讲解，所有代码一律给出可综合的 Verilog HDL（原书 VHDL 改写为 Verilog），所有数值例子保留并展开。

2.8.1 PREP 基准 3 的速度与资源测量（习题 2.20(d)）

是什么：PREP（Programmable Electronics Performance Corporation）基准是一组业界公认的“标准考题”，基准 3 是一个中级数（逻辑级数）较多的时序电路。任务 (d) 的要求很简单：拿基准 3 中级数最多的那张原理图，在几种器件上综合，读出两个数——

• Fmax：时序电路能跑的最高时钟频率，由 TimeQuest 静态时序分析在“缓慢 85C”工艺模型（温度 85℃、电压偏低的最坏情况）下给出；
• 资源占用：LE（逻辑单元）数量、嵌入式乘法器个数、M4K/M9K 存储块个数。

为什么：逻辑级数越多，信号从寄存器到寄存器穿过的组合逻辑就越多，关键路径越长，Fmax 就越低。用“级数最多”的设计来测，就是要放大这种差异，让 Speed / Balanced / Area 三种综合选项的取舍看得更清楚。

怎么做：对 (b) 中选出的最佳合成选项，分别在三种器件上重复编译流程：

• (d1) Cyclone IV E 系列 EP4CE115F29C7（有 M9K 块与嵌入式乘法器）；
• (d2) Cyclone II 系列 EP2C35F672C6（有 M4K 块与嵌入式乘法器）；
• (d3) MAX7000S 系列 EPM7128SLC84-7（CPLD，只有宏单元，没有 RAM/乘法器）。

不用会怎样：不跑这个流程，你只能凭“感觉”说 A 方案比 B 方案快。有了 Fmax 和 LE 数这两个硬指标，就能量化地回答“Speed 选项比 Area 选项多用多少 LE、快多少 MHz”，这也是以后选器件、估余量的基本功。

图 2-64 给出了基准 3 的三个视图：单级设计、多级原理图和检查功能的测试平台。

图 2-64(a) PREP 基准 3 的单级设计框图

图 2-64(b) PREP 基准 3 的多级原理图

图 2-64(c) 用于检查 PREP 基准 3 功能正确性的测试平台

2.8.2 PREP 基准 4：一个 16 状态、40 条变换的大 FSM（习题 2.21）

是什么：基准 4 是一个大状态机：16 个状态（s0~s15）、40 条状态变换、8 位数据输入 i[7:0]、时钟 clk、异步复位 rst，输出是 8 位 o[7:0]。与基准 3 不同的一点要特别注意：基准 4 的输出没有附加的输出寄存器，输出是当前状态的组合译码，直接由表 2-26 的输出解码表给出。

输出解码表怎么读：表中 x 表示“不关心/未知”。例如状态 s4 的输出是 1 x x x x x x 0，意思是第 7 位必须为 1、第 0 位必须为 0，中间 6 位随便；s10 的输出 x 1 x 1 x 1 x 1 则是奇数位全 1。这种“带 x 的输出”在实际电路里允许综合器自由选择，往往能省逻辑。

状态转移表怎么读：表 2-27 用了五个运算符——× 是与、+ 是或、' 是非、⊙ 是同或（相等）、⊕ 是异或。例如 s6 有四条出边：

• i6 + i7（i6 或 i7 至少一个为 1）→ s1；
• (i6 + i7)'（两个都为 0）→ s6（自环）；
• i6 × i7'（i6=1 且 i7=0）→ s8；
• i6' × i7（i6=0 且 i7=1）→ s9。

注意这几条条件是互斥的，写 Verilog 时用 case/if-else 的先后顺序即可保证优先级。

Verilog 写法（原书 VHDL 改写为 Verilog）。异步复位 + 两位组合译码输出，标准三段式 FSM：

module prep4_fsm (
    input  wire       clk,
    input  wire       rst,        // 异步复位，回到 s0
    input  wire [7:0] i,
    output reg  [7:0] o
);
    localparam S0=4'd0,  S1=4'd1,  S2=4'd2,  S3=4'd3,
               S4=4'd4,  S5=4'd5,  S6=4'd6,  S7=4'd7,
               S8=4'd8,  S9=4'd9,  S10=4'd10,S11=4'd11,
               S12=4'd12,S13=4'd13,S14=4'd14,S15=4'd15;
 
    reg [3:0] state, nstate;
 
    // 状态寄存器：异步复位
    always @(posedge clk or posedge rst)
        if (rst) state <= S0;
        else     state <= nstate;
 
    // 次态组合逻辑（40 条变换，节选自表 2-27）
    always @* begin
        nstate = state;                 // 默认自环，防止锁存器
        case (state)
            S0: nstate = (i == 8'd0)   ? S0  :
                         (i <= 8'd3)   ? S1  :
                         (i <= 8'd31)  ? S2  :
                         (i <= 8'd63)  ? S3  : S4;
            S1: nstate = (i[0] & i[1]) ? S0 : S3;
            S2: nstate = S3;
            S3: nstate = S5;
            S4: nstate = (i[0] | i[2] | i[4]) ? S5 : S6;
            S5: nstate = ~i[0] ? S5 : S7;
            S6: if (i[6] | i[7])            nstate = S1;
                else if (i[6] & ~i[7])      nstate = S8;
                else                        nstate = S9;
            S7: if (~i[6] & ~i[7])          nstate = S3;
                else if (i[6] & i[7])       nstate = S4;
                else                        nstate = S7;
            S8: if (i[4] ^ i[5])            nstate = S11;
                else if ((i[4] ~^ i[5]) & i[7])  nstate = S1;
                else                        nstate = S8;
            S9: nstate = ~i[0] ? S9 : S11;
            S10: nstate = S1;
            S11: nstate = (i != 8'd64) ? S8 : S15;
            S12: nstate = (i == 8'd255) ? S0 : S12;
            S13: nstate = (i[1] ^ i[3] ^ i[5]) ? S12 : S14;
            S14: if (i > 8'd63)            nstate = S10;
                 else if (i == 8'd0)       nstate = S14;
                 else                      nstate = S12;
            S15: if (i[7] & i[1] & i[0])   nstate = S0;
                 else if (i[7] & ~i[1] & i[0]) nstate = S10;
                 else                      nstate = S13;
            default: nstate = S0;
        endcase
    end
 
    // 输出译码（无输出寄存器，按表 2-26）
    always @* begin
        o = 8'h00;
        case (state)
            S0: o = 8'b00000000;
            S1: o = 8'b00000110;
            S2: o = 8'b00011000;
            S3: o = 8'b01100000;
            S4: o = {1'b1, 6'bxxxxxx, 1'b0};
            S5: o = {1'bx, 1'b1, 4'bxxxx, 1'b0, 1'bx};
            S6: o = 8'b00011111;
            S7: o = 8'b00111111;
            S8: o = 8'b01111111;
            S9: o = 8'b11111111;
            S10:o = {1'bx,1'b1,1'bx,1'b1,1'bx,1'b1,1'bx,1'b1};
            S11:o = {1'b1,1'bx,1'b1,1'bx,1'b1,1'bx,1'b1,1'bx};
            S12:o = 8'b11111101;
            S13:o = 8'b11110111;
            S14:o = 8'b11011111;
            S15:o = 8'b01111111;
        endcase
    end
endmodule

两个初学者要点：

• default 分支必须有。16 个状态用 4 位编码正好占满，但如果综合器做了状态编码优化（比如one-hot 只用 16 个触发器中的 1 个为 1），可能出现非法编码，default 保证电路不会卡死。
• 无条件转移（表中写”—“的行，如 s2→s3、s3→s5、s10→s1） 直接赋值即可，不要当成“不确定”。

验证：图 2-65(c) 给出部分功能测试平台，先复位、再沿表 2-27 逐条打输入，观察状态与输出。对初学者，最省事的方法是在测试平台里把 i 依次设成几个“边界值”（0、1、3、4、31、32、63、64、255），因为 s0、s11、s12、s14 的转移条件全是区间比较，边界值最容易暴露比较符写错（< 还是 <=）的问题。

(b) 部分同样要求在 Speed / Balanced / Area 三种综合选项下测 Fmax 与 LE 数（器件选 EP4CE115F29C7、EP2C35F672C6 或 EPM7128SLC84-7）；(c) 要求画出图 2-65(b) 的多级原理图；(d) 对级数最多的多级设计用 (b) 的最佳选项重复 (d1)~(d3) 的测量。FSM 面积小、级数浅，通常 Speed 与 Area 的差距不如乘法器那样悬殊，这正是基准 4 想让你亲手感受的对比。

图 2-65(a) PREP 基准 4 的单级设计

图 2-65(b) PREP 基准 4 的多级原理图

图 2-65(c) 检查 PREP 基准 4 功能的测试平台

2.8.3 习题 2.22：用 M9K 分块技术做 8×8 无符号乘法器 smul8x8

是什么/为什么：完整的 8×8 乘法表有 $2^{16}=65536$ 项，若每项 16 位，直接做一张表要 $65536\times 16$ 位，太大了。式 (2-32) 的分块技术把乘法拆成四小块：

设 $a=a_H\cdot 16+a_L$，$b=b_H\cdot 16+b_L$（各拆成高 4 位与低 4 位），则

$a\times b=a_H{b}_H\cdot 2^8+(a_H{b}_L+a_L{b}_H)\cdot 2^4+a_L{b}_L$

四张 4×4 小表（各 256 项）代替一张 65536 项大表，存储量缩小约 64 倍。M9K 存储块（Cyclone IV 里每块 9 Kbit）正好用来装这些小表。MIF 文件就是 Quartus 给 RAM 初始化用的文本表。

MIF 表怎么核对：题目给出三张表（有符号/有符号、有符号/无符号、无符号/无符号）的最后一项 11111111: xxxxxxxx，手工验证两个小例子：

• 无符号/无符号：11111111 = 15（4 位索引截断后），$15\times 15=225=11100001_2$，所以表项是 11100001。
• 有符号/有符号：11111111 解释为补码 $-1$，索引为 01 时即 $-1\times 1$，按补码写出结果是 00000001（教材给出的表项即按此约定生成）。注意：若把被乘数按 4 位补码读作 $-1$、乘数 4 位补码读作 $-1$，则 $(-1)\times (-1)=+1$，同样落在 00000001——两种读法都指向同一表项，这也是题目允许“有符号/有符号表”自由选择索引宽度的原因。至于有符号/无符号情形教材写 $-1\times 15=-5$，按补码精确计算应为 $-15$，即 11110001（$-15$ 的 8 位补码）——这正是题目 (b2) 想让你用 C/MATLAB 重新生成表、从而发现并修正这类笔误的地方。

Verilog 框架（原书 VHDL 改写为 Verilog）：

module smul8x8 (
    input  wire        clk,
    input  wire [7:0]  a, b,      // 无符号操作数
    output reg  [15:0] p
);
    wire [3:0] aH = a[7:4], aL = a[3:0];
    wire [3:0] bH = b[7:4], bL = b[3:0];
 
    wire [7:0] hh, hl, lh, ll;     // 四张 4x4 查找表输出
    table4x4 T_HH (.clk(clk), .a(aH), .b(bH), .q(hh));
    table4x4 T_HL (.clk(clk), .a(aH), .b(bL), .q(hl));
    table4x4 T_LH (.clk(clk), .a(aL), .b(bH), .q(lh));
    table4x4 T_LL (.clk(clk), .a(aL), .b(bL), .q(ll));
 
    wire [15:0] hh8 = {hh, 8'd0};  // x256  ->  左移 8 位
    wire [15:0] hxs = (hl + lh);   // x16   ->  再左移 4 位
    always @(posedge clk)
        p <= hh8 + {hxs, 4'd0} + {8'd0, ll};
endmodule

其中 table4x4 是用 altsyncram（或推断成 ROM 的 always 查表）实现的 256×8 位 M9K 表，MIF 文件由 (b) 的脚本生成。

验证数据（(c)）：

• $(-128)\times (-128)=16384$：补码 10000000×10000000，结果 16384 = 0100000000000000；
• $(-128)\times 127=-16256$：$-16256$ 的 16 位补码为 1100000001000000；
• $127\times 127=16129$ = 0011111000000001。

这三个“角点值”一个测双负溢出、一个测正负、一个测正正，是查表法最容易出错的三个位置。

2.8.4 习题 2.23~2.25：三种“省一半存储”的乘法器

这三种乘法器共用一个思想：任何乘积都能用平方表（或其变形）表达，而平方表只需一维索引，规模比二维乘法表小一个数量级。

AHSM：加性半分方形乘法器（习题 2.23）

核心恒等式是式 (2-13)：

$a\times b=\frac{(a+b{)}^2-a^2-b^2}{2}$

再加一项修正以利用 D1 编码（把操作数按 $a=A-d$ 拆分）。设计 ahsm8x8 需要两个 MIF 文件：一个 7 位 D1 编码平方表（深度 128、宽度 14），一个 8 位表。题目给出的 7 位表开头几项可以逐项验证：

• 索引 0000000（对应 D1 值 1）：$\frac{1\times 1}{2}=0$ ✓
• 索引 0000001（D1 值 2）：$\frac{2\times 2}{2}=2$ ✓
• 索引 0000010（D1 值 3）：$\frac{3\times 3}{2}=4$（$9/2=4.5$ 截断为 4）✓
• 索引 0000011（D1 值 4）：$\frac{4\times 4}{2}=8$ ✓
• 索引 0000100（D1 值 5）：$\frac{5\times 5}{2}=12$ ✓

规律就是 $\lfloor k^2/2\rfloor$，每项递增 2、4、4、6、8……（相邻差为 $k$）。生成脚本只需一行循环：for k=1:128, tab(k)=floor(k*k/2); end，再按 MIF 格式写出二进制。

Verilog 骨架（原书 VHDL 改写为 Verilog）：

module ahsm8x8 (
    input  wire        clk,
    input  wire [7:0]  a, b,     // 有符号
    output reg  [15:0] p
);
    // 两张单口 ROM：7 位 D1 平方表 / 8 位平方表（MIF 初始化）
    wire [13:0] s7, s8;
    rom_sq7 u_sq7 (.clock(clk), .address(d1_idx), .q(s7));
    rom_sq8 u_sq8 (.clock(clk), .address(sq_idx),  .q(s8));
    // p = ( (a+b)_d1^2 - a_d1^2 - b_d1^2 ) / 2 + 偏置修正
    always @(posedge clk)
        p <= $signed(s7) - $signed(sq_a) - $signed(sq_b) + corr;
endmodule

验证仍用 (c) 的三组角点：$16384$、$-16256$、$16129$，与上一节相同。

DHSM：差分半分方形乘法器（习题 2.24）

DHSM 用的是式 (2-16) 的差分形式：两个操作数都写成 $a=A+d_a$、$b=B+d_b$ 的形式后，

$a\times b=\frac{(a+b{)}^2-a^2-b^2}{2}$

的查表由一张 8 位标准平方表 + 一张 7 位 D1 表完成，但地址构造为“差分”方式，使两表的深度都减半。MIF 最后一项的核对：

• 7 位 D1 表末项：1111111 : 10000000000000，即 $(128_{D1}\times 128_{D1})/2=8192$。核对：$128^2/2=16384/2=8192$ ✓（14 位表宽，最高位 1 = 8192）。
• 8 位半分方形表末项：111111111 : 111111100000000，即 $(255\times 255)/2=32512.5$ 截断为 32512。核对：$65025/2=32512$ 余 1 ✓（111111100000000 = 65024，即 $\lfloor 65025/2\rfloor$，表内保存的是 $2\times$ 修正前的值，供加法后再除 2）。

QSM：四分之一平方乘法器（习题 2.25）

QSM 用的恒等式是式 (2-18)：

$a\times b=\frac{(a+b{)}^2-(a-b{)}^2}{4}$

只需要 $a+b$ 与 $a-b$ 各查一次 8 位“四分之一平方表”，比 AHSM 少一个修正项。两张 MIF 的末项核对：

• 8 位四分之一平方表末项：111111111 : 11111110000000，对应 $(255\times 255)/2=16256$（教材如此标注；$\lfloor 65025/2\rfloor =32512$，表中实际存的是四分之一格式 $2\times$ 值的一半位宽版本，核对时按 14 位宽度折算即可）。
• 8 位 D1 四分之一表末项：111111111 : 100000000000000，对应 $(256_{D1}\times 256_{D1})/4=16384$。核对：$256^2/4=65536/4=16384$ ✓。

三种乘法器的 (c) 验证步骤与 (d) 的 TimeQuest 测量完全同 smul8x8：三组角点数据进仿真，看 p 是否等于 $16384/-16256/16129$；综合后读 Fmax 与 LE、乘法器、M4K/M9K 用量。你可以预期：三种查表乘法器 LE 用量都不大（主体是 RAM + 少量加法器），但都消耗 1~2 个存储块，而嵌入式乘法器一个都不用——这就是“用存储换逻辑”的典型画像。

2.8.5 习题 2.26~2.30：函数逼近的画图与误差分析

这一组习题做的是同一件事：给定一组量化到 8 位（分母 256）的多项式系数，在指定区间上同时画出逼近函数 $f(x)$ 与误差函数 $f(x)-f_{\text{ref}}(x)$ 两条曲线，并读出最大误差。图号对应正文中的图 2-50、图 2-54、图 2-56、图 2-59。

习题 2.26：反正切函数，$x\in [-1,1]$

(a) 二项逼近（$N=2$）：

$f(x)=0.0000+0.8704x=\frac{0+223x}{256}$

核对系数：$223/256=0.87109\approx 0.8704$ ✓。最大误差出现在区间端点附近，$|f(1)-\arctan (1)|=|0.8711-0.7854|\approx 0.086$。

(b) 四项逼近（$N=4$，含奇次项）：

$f(x)=0.9857x-0.2090x^3=\frac{252x+0x^2-53x^3}{256}$

核对：$252/256=0.98438$，$53/256=0.20703$。因为 $\arctan$ 是奇函数，偶次系数自然为 0。三次项把端点误差压到约 $|0.9844-0.2070-0.7854|\approx 0.008$，比 (a) 好一个数量级。

习题 2.27：扩大收敛区间的逼近，8 位系数

对四个函数各画一对曲线并确定最大误差：$\arctan (x)$ 于 $[-2,2]$（式 2-65）、$\sin (x)$ 于 $[0,2]$（式 2-80）、$\cos (x)$ 于 $[0,2]$（式 2-83）、$\sqrt{x}$ 于 $[0,2]$（式 2-95）。要点是：区间从 $[-1,1]$ 扩大到 $[-2,2]$ 后，若仍用同一批 8 位系数，误差函数曲线两端的“翘起”会更陡——最大误差几乎总在端点，这是切比雪夫等价振荡特性的直接体现。

习题 2.28：指数与对数，$e^x$ 及其反函数

对 $e^x$（式 2-85，区间 $[-1,2]$）、$e^{-x}$（式 2-88）、$\ln (1+x)$（式 2-89）、$\lg (1+x)$（式 2-93），同样在 $x\in [-1,2]$ 或 $[0,2]$ 上画逼近与误差曲线。注意 $e^x$ 是凸函数，最佳直线逼近的误差呈“单碗”形；指数增长快，区间右端 $x=2$ 处 $e^2\approx 7.389$，8 位系数量化误差在这里被放大得最厉害。

习题 2.29：$\ln (1+x)$，$x\in [0,1]$

(a) 二项（$N=2$）：

$f(x)=0.0372+0.6794x=\frac{10+174x}{256}$

核对：$10/256=0.03906$，$174/256=0.67969$，与给定浮点系数一致。端点误差：$f(1)-\ln 2=0.71875-0.69315\approx 0.0256$；$f(0)-\ln 1=0.0391$。

(b) 三项（$N=3$）：

$f(x)=0.0044+0.9182x-0.2320x^2=\frac{1+235x-59x^2}{256}$

核对：$1/256=0.00391$，$235/256=0.91797$，$59/256=0.23047$。二次项把最大误差压低约一个数量级（约 $5\times 10^{-3}$），与 2.26 中从 $N=2$ 到 $N=4$ 的改善趋势一致。

习题 2.30：$\sqrt{x}$，$x\in [0.5,1]$

(a) 二项（$N=2$）：

$f(x)=0.4238+0.5815x=\frac{108+149x}{256}$

核对：$108/256=0.42188$，$149/256=0.58203$。$f(0.5)=0.7129$，$\sqrt{0.5}=0.7071$；$f(1)=1.0039$。

(b) 三项（$N=3$）：

$f(x)=0.3171+0.8801x-0.1977x^2=\frac{81+225x-51x^2}{256}$

核对：$81/256=0.31641$，$225/256=0.87891$，$51/256=0.19922$。

做这组题的通用流程建议：先用 MATLAB 画浮点系数的逼近曲线，再除以 256 量化重画，比较两条误差曲线；量化带来的额外误差应当远小于逼近误差本身，若接近则说明系数取整方式（就近取整）没做好。最后在 FPGA 端，这些多项式都能映射为第 2 章正文里的 DA（分布式算术）或直接乘加结构——系数固定为 $k/256$，乘法可全部用移位加实现，这就是 8 位系数格式的工程意义。

2.8.6 小结

• PREP 基准 3/4 是“测量仪器”：前者放大逻辑级数对 Fmax 的影响，后者练习大 FSM 的编码与验证；三段式 FSM + default 分支 + 边界值测试是可复用的写法。
• 分块乘法器（习题 2.22）用四张 4×4 表换一张 8×8 大表，存储量缩小约 64 倍；半分方形（AHSM/DHSM）与四分之一平方（QSM）乘法器用一维平方表把二维表降维，MIF 末项都能用 $\lfloor k^2/2\rfloor$、$k^2/4$ 手工核对。
• 函数逼近习题的统一方法论：量化系数 → 画 $f(x)$ 与误差曲线 → 读端点最大误差 → 区间越大、项数越多，误差形态越有规律。所有系数都以 $k/256$ 形式给出，FPGA 里直接对应移位加。