无乘法器的多常数乘法

Yevgen Voronenko and Markus Püschel. 2007. Multiplierless multiple constant multiplication. ACM Trans. Algorithms 3, 2 (May 2007), 11–es.

第一章 引言与问题定义

在数字信号处理（DSP）和计算机算术领域，一个核心问题是如何高效地计算变量x与多个已知定点常数t₁, t₂, ..., tₙ的乘积。传统的硬件乘法器占用大量芯片面积且功耗较高，因此研究者们探索了仅使用加法、减法和移位操作来实现乘法的方法。这种方法在硬件实现中特别重要，因为移位操作可以通过简单的布线实现，几乎不占用芯片面积。

本文研究的核心问题——多常数乘法（Multiple Constant Multiplication, MCM）问题——可以形式化定义为：给定一组正整数目标常数T = {t₁, ..., tₙ}，找到最小的基本数集合R = {r₀, r₁, ..., rₘ}，使得T ⊂ R，其中r₀ = 1，且每个rₖ（1 ≤ k ≤ m）都可以通过之前的基本数通过一次A操作得到。

第二章 单常数乘法的演化

2.1 基本方法与CSD表示

考虑最简单的例子，将变量x乘以常数71。二进制表示71 = 1000111₂，直接的分解方法是：

$$71x = 2^6x + 2^2x + 2^1x + x$$

这需要3次加法操作。使用减法的另一种分解是：

$$71x = 2^7x - x - 2^5x - 2^4x - 2^3x = (2^7 - 1)x - (2^5 + 2^4 + 2^3)x$$

规范符号数字（CSD）表示法引入了负数位的概念，允许使用1̄来表示-1。在CSD表示中，71 = 100100̄1，因此：

$$71x = 2^6x + 2^3x - x$$

这将操作数减少到2次。CSD表示的关键性质是任意两个非零位之间至少有一个零位，这保证了表示的唯一性。

2.2 图形表示与拓扑结构

图2描述：此图展示了常数45的两种不同分解方式。上部显示CSD分解，需要3次加减操作；下部显示最优分解，仅需2次操作。图中每个圆形节点代表一次加减操作，节点内的数字表示该操作的输出值。连接节点的有向边标注了移位量（2的幂次），负值表示该边对应的是减法操作。值得注意的是，最优分解使用了不同于CSD的图拓扑结构——它先计算了中间值3 = x + 2x，然后利用这个中间值计算45 = 3 + 16×3 - 4x。

第三章 A操作的数学框架

3.1 A操作的完整定义

A操作是本文算法的核心构建块。其通用形式定义为：

$$A_p(u, v) = \left|2^{l_1}u + (-1)^s 2^{l_2}v\right| \cdot 2^{-r}$$

其中参数集p = (l₁, l₂, r, s)必须满足以下约束：

• l₁, l₂ ≥ 0（左移量）
• r ≥ 0（右移量）
• s ∈ {0, 1}（符号位，决定加法或减法）
• 2^r 必须整除 2^{l₁}u + (-1)^s2^{l₂}v（保证不丢失有效位）

图6描述：A操作的图形表示显示两个输入基本数u和v通过各自的移位后进行加法或减法运算，产生输出基本数w。图中标注了具体的移位参数和符号。直接连接到乘法器块输入的A操作具有u = v = 1的特殊情况。

3.2 奇数基本数图的约简

任何A图都可以转换为等价的仅包含奇数基本数的A图，这种转换不会增加操作数。转换的关键在于限制A配置：

• 最多允许一个非零左移（l₁或l₂）
• 当l₁ = l₂ = 0时，r必须是产生奇数值的唯一右移

这种约简显著减少了搜索空间，因为对于给定的u和v，A_{odd}的自由参数仅为s和非零左移，而通用A操作中l₁、l₂、r和s都可以变化。

第四章 多常数乘法问题

4.1 问题复杂度与共享优化

图3描述：此图展示了一个典型的乘法器块结构，输入x通过一系列A操作产生多个输出t₁, t₂, ..., tₙ。内部节点表示中间计算结果，这些中间结果的共享是MCM优化的关键。

MCM问题的NP完全性可以通过归约SCM问题来证明。关键观察是，如果我们能够高效地计算A距离，那么算法2对于单个常数就是最优的。由于SCM问题是NP完全的，A距离计算也必然是NP完全的。

图4描述：此图显示了随着目标常数集合大小n的增加，平均每个12位常数所需的加减操作数。当使用独立的最优SCM分解时（上方曲线），操作数线性增长；而使用RAG-n算法（下方曲线）时，由于中间结果的共享，增长率显著降低。

4.2 实例分析：23和81的联合优化

图5描述：这是MCM优化效果的一个具体例子。图中展示了同时计算23x和81x的乘法器块。通过巧妙地选择中间节点9（= x + 8x），算法首先计算23 = 32×9 - 9x，然后利用已有的9和23计算81 = 8×(9x) + 9x。整个过程仅需3次加减操作，而独立的最优分解需要4次（每个常数各需2次）。

第五章 算法框架与现有方法

5.1 通用图基算法框架

所有图基MCM算法共享以下高层结构：

算法1：图基MCM算法通用框架
输入：目标常数集T
输出：基本数集R，满足T ⊂ R

1. R ← {1}  // 初始化就绪集
2. while T ≠ ∅ do
3.     计算R的后继集S = {s | dist(R,s) = 1}
4.     根据启发式从S中选择s
5.     合成s：R ← R ∪ {s}, T ← T \ {s}

5.2 现有算法的详细分析

Bull-Horrocks算法（BHA）

BHA使用误差最小化策略：
$$\epsilon = \min(T) - \max(R)$$

算法按升序合成目标，每次迭代选择两个后继s₁和s₂：
$$s_1 = \arg\min_{s \in S, s \leq \epsilon}(\epsilon - s)$$
$$s_2 = \max(R) + s_1$$

RAG-n算法的三种情况

图7描述：RAG-n算法考虑三种不同的情况。左侧的"最优情况"显示目标t直接在后继集S中（距离1）。中间的"启发式情况A"显示目标在S₂中（距离2），虚线圆表示S₂不是显式计算的，而是使用启发式测试。右侧的"启发式情况B"显示目标距离超过2，需要使用预计算的查找表。

第六章 新算法设计

6.1 算法结构与优化

新算法的核心创新在于其启发式函数的设计。算法维护三个关键集合：

• R：就绪集（已合成的基本数）
• S：后继集（距离R为1的所有常数）
• W：工作列表（临时存储新合成的常数）

后继集的增量更新是算法效率的关键。设新合成的常数集为W，则更新公式为：

$$S_{new} = A^*(R_{new}, R_{new}) = A^*(R \cup W, R \cup W)$$

通过集合运算的性质，可以推导出：

$$S_{new} = (S \cup A^*(R_{new}, W)) - W$$

这个公式避免了每次迭代重新计算整个后继集。

6.2 效益函数与启发式

定义效益函数量化添加后继s对到达目标t的改进：

$$B(R, s, t) = \text{dist}(R, t) - \text{dist}(R + s, t)$$

考虑到距离估计的准确性随距离增加而降低，引入指数衰减的权重：

$$\bar{B}(R, s, t) = 10^{-\text{dist}(R+s, t)} \cdot B(R, s, t)$$

累积效益启发式通过最大化所有目标的总效益来实现联合优化：

$$H_{cub}(R, S, T) = \arg\max_{s \in S} \sum_{t \in T} \bar{B}(R, s, t)$$

第七章 A距离的精确计算与估计

7.1 距离1-3的精确测试

图8描述：此图展示了算法处理的四种特殊距离情况。(a)显示距离1的情况，目标直接在后继集S中。(b)显示距离2，目标在S²中但不在S中。(c)显示距离3，目标在S³中。(d)显示距离≥4，此时S³集合为空（用虚线表示未计算）。

图9描述：这组图详细展示了用于精确距离测试的所有图拓扑。(a)显示唯一的距离1拓扑。(b)显示两种距离2拓扑，分别对应不同的中间节点连接方式。(c)显示五种距离3拓扑，包括串联、并联和混合结构。每个拓扑上标注了相应的A方程，这些方程用于确定是否存在满足条件的中间值。

对于距离2的情况1，A方程为t = c₁s，其中c₁ ∈ C₁。解的存在性测试为：

$$\frac{t}{C_1} \cap S \neq \emptyset$$

对于距离2的情况2，使用引理5.2，从t = Aₚ(s, r₂)推导出：

$$s \in A^*(t, r_2) \Rightarrow A^*(t, R) \cap S \neq \emptyset$$

7.2 高距离的估计方法

图10描述：此图展示了用于距离估计的部分图拓扑。每个拓扑包含一个未合成的输入z（用虚线圆表示）。估计过程首先确定所有可能的z值，然后选择CSD代价最小的值，最后将该代价加到拓扑中的操作数（减1，因为假设s已可用）。

距离估计基于以下不等式：

$$\text{dist}(R + s, t) \leq \text{dist}(R + s, z) + \text{dist}(R + s + z, t)$$

使用CSD代价作为辅助度量：

$$\text{Est}(z) = \text{CSD-Cost}(z) = \text{CSD表示中非零位的个数}$$

对于可能值集合Z：

$$\text{Est}(Z) = \min_{z \in Z} \text{Est}(z)$$

最终的距离估计为：

$$\text{dist}(R + s, t) \lesssim \min(\text{dist}(R, t), E_1, E_2, E_3, E_4)$$

其中各估计值为：

• $E_1 = 1 + \text{Est}(A^*(s, t))$
• $E_2 = 2 + \text{Est}(A^*(s, t/C_1))$
• $E_3 = 2 + \text{Est}(A^*(C_1 \cdot s, t))$
• $E_4 = 2 + \text{Est}(A^(R, A^(s, t)))$

第八章 实验评估

8.1 固定常数数量的实验

图11描述：这组图展示了四种不同算法（BHM、Lefèvre、Hcub、RAG-n）在固定常数数量（n = 1, 2, 10, 20）下的性能比较。横轴是常数位宽b，纵轴是平均加减操作数。对于单个常数（n=1），RAG-n使用最优SCM分解，因此性能最好。随着常数数量增加，Hcub的优势逐渐显现，特别是在n=20时，Hcub比RAG-n减少了17%的操作，比BHM减少了25%的操作。

8.2 固定位宽的实验

图12描述：这组六个子图展示了在不同固定位宽（b = 12, 16, 19, 22, 28, 32）下，操作数随常数数量n的变化。对于12位常数，所有算法都快速收敛到理论下界n。随着位宽增加，收敛速度减慢。在28位和32位的情况下，Hcub相比BHM的优势达到26%。

8.3 相对性能分析

图13和图14描述：这两组图展示了Hcub和BHM相对于RAG-n的性能比率。图13固定位宽，变化常数数量；图14固定常数数量，变化位宽。比率小于1表示优于RAG-n。Hcub在大多数情况下都显著优于RAG-n，特别是在b=19、n=80时，改进达到20%。

8.4 距离测试的影响

图15描述：此图比较了使用距离2测试和使用距离2+3测试的Hcub算法性能。对于单个常数，仅使用距离2测试会导致6-15%的性能下降。随着常数数量增加，两种版本的差异逐渐减小，这是因为联合优化的效果超过了精确距离测试的重要性。

第九章 复杂度分析

9.1 最坏情况界限

集合大小的最坏情况界限：

• $|A^*(u, v)| = O(b)$：每对输入最多产生O(b)个输出
• $|C_1| = O(b)$：复杂度1的常数数量
• $|C_2| = O(b^2)$：复杂度2的常数数量
• $|R| = O(nb)$：解的大小
• $|S| = O(n^2b^3)$：后继集大小

9.2 运行时分析

算法各部分的运行时间：

如果仅使用距离2测试，运行时间降至O(n³b⁵)。

附录A：关键引理与定理

A.1 引理5.2的证明（A操作的可逆性）

引理5.2：如果w = Aₚ(u, v)，则存在A配置p'使得u = Aₚ'(w, v)。

证明：从A操作的定义出发：
$$w = |2^{l_1}u + (-1)^s 2^{l_2}v| \cdot 2^{-r}$$

两边乘以2^r：
$$2^r w = |2^{l_1}u + (-1)^s 2^{l_2}v|$$

考虑符号，我们有：
$$2^r w = 2^{l_1}u + (-1)^s 2^{l_2}v \quad \text{或} \quad 2^r w = -(2^{l_1}u + (-1)^s 2^{l_2}v)$$

对于第一种情况，解出u：
$$2^{l_1}u = 2^r w - (-1)^s 2^{l_2}v = 2^r w + (-1)^{s+1} 2^{l_2}v$$
$$u = 2^{-l_1}(2^r w + (-1)^{s'} 2^{l_2}v) = |2^r w + (-1)^{s'} 2^{l_2}v| \cdot 2^{-l_1}$$

其中s' = (s+1) mod 2。这正是Aₚ'(w, v)的形式，其中p' = (r, l₂, l₁, s')。

A.2 定理4.2的证明（算法终止性）

定理4.2：如果距离函数dist是可容许的，则算法2以Hmaxb或Hcub作为启发式函数时必然终止。

证明：定义未合成目标的距离和：
$$D = \sum_{t \in T} \text{dist}(R, t)$$

根据可容许性条件：

1. D是有限的非负整数（条件1）
2. T = ∅当且仅当D = 0（条件2-3）
3. 添加元素到R不会增加距离（条件4）
4. 总存在正效益的后继（条件5）

在每次迭代中，启发式选择具有正效益的后继s，因此：
$$\text{dist}(R + s, t) < \text{dist}(R, t)$$

这意味着D至少减少1。由于D是有限的且每次迭代都减少，算法必然在有限步内使D = 0，此时T = ∅，算法终止。□

A.3 定理5.1的证明（A距离计算的NP完全性）

定理5.1：在约束Aₚ(u, v) ≤ 2^{b+1}下计算A距离是NP完全问题。

证明：通过多项式时间归约SCM问题来证明。

首先，A距离问题在NP中，因为给定一个声称的距离d和相应的A操作序列，我们可以在多项式时间内验证该序列确实产生目标常数。

其次，我们将NP完全的SCM问题归约到A距离计算：

1. 如果dist函数是精确的，根据定理4.3，算法2对单个常数是最优的
2. 算法中，启发式每次迭代调用一次，共O(b)次迭代
3. 每次迭代计算O(|S|) = O(b³)个效益值
4. 因此A距离计算O(b⁴)次

这表明最优SCM分解可以在多项式时间内归约到A距离计算，因此A距离计算是NP完全的。□

附录B：A距离估计的可容许性证明

B.1 定理5.6的证明

定理5.6：设dist(R+s, t) ≃ d是从方程(23)获得的距离估计，且d < dist(R, t)。则在下一次迭代中，存在s' ∈ S_new使得dist(R+s+s', t) = d-1。

证明：我们需要对四种估计情况分别证明。以E₁情况为例进行详细推导。

情况1（d = E₁）：
根据估计器，t ∈ A*(s, z)对某个z成立，且：
$$d = 1 + \text{Est}(z)$$

由于d ≥ 3，有Est(z) ≥ 2，即z的CSD代价至少为2。

从z的CSD表示中移除一个非零位，得到z'，满足：
$$\text{Est}(z') = \text{Est}(z) - 1$$

存在A配置p使得：
$$z = A_p(1, z')$$

将此代入原A方程：
$$t \in A^*(s, A_p(1, z'))$$

通过A操作的结合性（这需要仔细的代数操作），存在配置p'和p''使得：
$$t \in A^*(A_{p'}(s, 1), z')$$

注意到$A{p'}(s, 1) \in S{new}$（因为s ∈ S，1 ∈ R），设s' = $A_{p'}(s, 1)$，则：
$$t \in A^*(s', z')$$

应用E₁估计器到新配置：
$$\text{dist}((R+s)+s', t) \leq 1 + \text{Est}(z') = 1 + (\text{Est}(z) - 1) = d - 1$$

因此B(R+s, s', t) > 0，满足定理要求。

情况2-4的证明采用类似方法，关键是利用A操作的代数性质和CSD表示的可分解性。完整证明需要对每种拓扑结构进行详细的代数推导。□

附录C：实验数据补充

C.1 标准偏差分析

除了平均操作数，我们还测量了不同算法在随机常数集上的标准偏差：

Hcub显示出最小的标准偏差，表明其性能更加稳定。

C.2 特殊常数集的性能

对于某些具有特殊结构的常数集（如FIR滤波器系数），算法表现出不同的特性：

对称FIR滤波器系数：由于系数的对称性，中间结果共享的机会增加，Hcub的优势更加明显，相比RAG-n可减少25-30%的操作。

2的幂次附近的常数：如{31, 33, 63, 65, 127, 129}，这类常数的CSD表示特别简单，所有算法的性能都接近最优。

素数常数集：素数通常具有较复杂的二进制模式，是MCM算法的挑战性输入。Hcub在这类输入上的改进最为显著。