定点数平方根的VHDL实现

定点数平方根的VHDL实现

在嵌入式系统中处理数字时，最常用的数学算法是加/减法、乘法、除法。当我们从标量扩展到矢量和复数时，在处理矢量长度时不可避免地会遇到平方根运算。平方根还出现在数字的几何平均值以及测量的各种能量和RMS计算中。因此，就基本数学运算而言，它是下一个最重要的运算。

在本文中，我们将使用Newton-Raphson逆平方根算法开发一种参数化字长和精度的定点数平方根基本算法。然后，该算法将针对real_math VHDL包实现进行测试，在GHDL中仿真，最后综合到Efinix Titanium FPGA上。Efinix Titanium的源码可在示例项目的GitHub页面上获取，此处展示的乘法和定点数平方根算法可在hVHDL定点数库中获取。项目中所有的测试平台由VUnit管理，仿真可以通过以下命令启动

python vuint_run.py -p 32 --gtkwave-fmt ghw 然后可以在vunit_out/test_output/test_name /ghdl/wave.ghw中找到波形文件

VHDL模块的高层设计

尽管VHDL用于描述数字电路的功能，但它具有若干高级语言特性，如自定义类型和复合类型以及子程序重载。借助这些特性，我们可以为代码创建接口，从而使模块更易于使用。由于数字电路是大规模并行的，且我们实现的是定点算法，因此抽象的要点在于时序和基数处理。将这些处理封装为子程序，使我们能够专注于使用算法本身，而不必担心底层细节。

以下代码片段展示了如何使用记录和子程序为VHDL模块创建接口。create_multiplier和create_square_root过程围绕记录中定义的寄存器建立逻辑，随后的request_square_root过程则修改状态以触发计算。请注意，平方根外部没有计数器，取而代之的是记录-过程对，它将所有状态机和时序包含在平方根内部，我们通过轮询square_root_is_ready函数来读取计算就绪的时刻。

-- 一个模块的架构 signal multiplier : multiplier_record := init_multiplier signal square_root : square_root_record := init_square_root; begin process_with_fixed_point_square_root : process(clock) begin if rising_edge(clock) then create_multiplier(multiplier); create_square_root(square_root, multiplier); if square_root_calculation_is_needed then request_square_root(square_root, input_number); end if; if square_root_is_ready(square_root) then result <= get_square_root_result(square_root, radix); end if; end process; 平方根模块的接口具有握手接口。create_过程创建模块或初始化模块，然后通过request_square_root过程启动计算，最后我们等待模块就绪，然后以期望的基数调用get_square_root_result获取结果。

通过使用模块接口以期望的基数获取结果，我们大幅简化了定点数计算。借助该接口，我们无需在代码外部跟踪基数，而是选择所使用的基数，模块便以期望的基数返回其输出。通过square_root_is_ready函数，我们将时序信息对应用代码隐藏，因为模块自身会跟踪其时序。这有助于将变更封装在接口内部——例如，如果我们使用该接口并更改延迟（比如添加寄存器以满足时序），应用代码将自动适应。

我们使用记录和子程序创建接口，平方根功能无需单独的实体即可实现；相反，我们直接使用记录和子程序。如有需要，我们可以将其封装在一个实体内，从而通过独立的记录结构来管理输入和输出数据。下面的代码片段展示了如何在VHDL模块中使用这种高层次抽象。我们用记录和子程序编写接口，此处的平方根设计无需实体，直接使用记录和子程序实现。如果将其放入一个实体内，那么我们将分别拥有用于输入和输出数据的记录。--> <!--

如代码片段所示，平方根功能将直接通过记录和子程序实现，无需绑定实体。这使得平方根计算只需几行代码即可添加到任何需要它的进程中，并且可以以最小的依赖进行测试和设计。如有需要，我们可以将平方根模块封装在一个实体内，在这种情况下，我们将修改接口，通过独立的记录来使用输入和输出数据。

平方根算法

对于快速且精确的平方根实现，有两种相似的算法：Goldschmidt算法和基于Newton-Raphson迭代的倒数根求解算法。Michael Morris撰写的一篇关于Goldschmidt算法的文章此前已发表在embeddedrelated上，而使用NR进行倒数平方根的快速计算更是颇具传奇色彩，以至于拥有自己的维基百科条目。这两种算法都是迭代式的，这意味着我们可以直接通过循环运行它们，利用前一次运行的结果在下一次获得更精确的结果。我们将快速审视这两种算法的相似之处。

Goldschmidt平方根

Goldschmidt算法同时估算$\sqrt{x}$和$\frac{1}{\sqrt{x}}$，该算法允许两次乘法并行或流水线计算，并且需要三个初始值来启动。

$$x_{0} = S\left(\frac{1}{\sqrt{S}}\right)_{est}$$

$$h_{0} = \frac{1}{2}\left(\frac{1}{\sqrt{S}}\right)_{est}$$

$$r_{0} = \frac{1-2x_0 h_0}{2}$$

收敛到$\sqrt{S} = x_i$和$\frac{1}{\sqrt{S}} = 2h_i$的实际算法为

$$x_{i} = x_{i-1}\cdot (1+r_{i-1})$$

$$h_{i} = h_{i-1}\cdot (1+r_{i-1})$$

$$r_{i} = \frac{1-2x_i \cdot h_i}{2}$$

$x_i$和$h_i$的计算可以并行进行，它们的结果用于计算$r_i$。算法中明确写出了必需的乘法运算，以表明何处必须使用实际的乘法器。

Newton-Raphson倒数平方根

使用 NR 求根算法 $x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$ 作用于 $\frac{1}{x^2}-S = 0$（其根位于 $\frac{1}{\sqrt{x}}$），可得到经典的递归 NR 平方根算法

$$x_{n+1} = 1.5x_n - 0.5Sx_n^3$$

经过一些改写，我们可以看出 Newton-Raphson 与 Goldschmidt 算法之间的相似之处

$$S_{xi} = S\cdot x_{i-1}$$

$$x_{xi} = x_{i-1} \cdot x_{i-1}$$

$$x_{i} = 1.5x_{i-1} - 0.5x_{xi}\cdot S_{xi}$$

要启动 NR 迭代，我们只需要 $S$ 的倒数平方根的一个估计值

$$x_{0} = \left(\frac{1}{\sqrt{S}}\right)_{est}$$

NR 算法给出倒数平方根，由此我们可以很容易地得到 $\sqrt{x}$，即 $x \cdot \frac{1}{\sqrt{x}} = \sqrt{x}$。从改写后的方程可以看出，两种算法具有非常相似的结构。两者都有两个可流水化的乘法，并需要第三个乘法将前两个乘法的结果相乘。由于该算法的初始化更为简便，我们将使用 NR 来实现平方根。hVHDL 库还有一个除法算法，其实现实际上完全相同。

定点实现

这里给出的算法同样适用于定点或浮点实现。定点与浮点之间的区别在于，在定点中我们将指数（或等价地，整数位数）编码到实现中，因为它不像浮点数那样随数字表示动态携带。

选择NR算法的基数

为逆平方根算法的实现选择整数位数相当直接。只要初始猜测足够接近，收敛就会很快，因此我们将输入范围限制在$[1.0, 2.0)$内，而该范围的逆平方根映射到$[0.707, 1.0]$。我们将为算法使用2个整数位，这是覆盖两个范围所需的最少整数位数，因此我们的基数将由字长减去2来确定。我们将在整个计算过程中使用该基数，因此在迭代实现期间无需调整任何基数或字长。

采用此选择后，表示最小值和最大值的位向量分别为$1.00000...$和$1.111111...$。因此，我们可以通过左移输入的任何位向量直到其没有前导零来进行映射，然后再次左移输入中前导零数量的一半进行后缩放。

NR迭代模块

为了开始实现，我们可以将NR迭代分解为加法和乘法的单一操作，以获得

步骤 0：x <= initial_value；循环：步骤 1：ax <= a * x；three_halfs <= x + x/2；步骤 2：x2 <= x * x；步骤 3：ax3 <= ax * x2；步骤 4：x <= three_halfs - ax3；跳转至循环；步骤 5：result <= x；步骤 1 至 4 可在循环中运行以获得更高精度。我们将使用定点乘法器模块来执行乘法运算。该模块同样以记录和过程的形式编写，具有可参数化的输入和输出流水线深度。该模块通过配置包对位宽和流水线深度进行参数化。我们本可以使用包泛型来实现参数化，但我们仍在用于 FPGA 的部分工具——主要是 Intel Quartus 和旧版本的 Vivado——不允许使用 VHDL2008 特性。

构建逆平方根所需的寄存器定义在 isqrt_record 中。两个状态计数器寄存器 state_counter 和 state_counter2 用于对计算过程进行时序控制，loop_value 寄存器使我们能够跟踪希望运行 NR 迭代循环的次数。完整的逆平方根记录如下所示

type isqrt_record is record x_squared : signed(int_word_length-1 downto 0); threehalfs : signed(int_word_length-1 downto 0); x : signed(int_word_length-1 downto 0); result : signed(int_word_length-1 downto 0); sign_input_value : signed(int_word_length-1 downto 0); state_counter : natural range 0 to 7; state_counter2 : natural range 0 to 7; isqrt_is_ready : boolean; loop_value : natural range 0 to 7; end record; 我们通过 create_isqrt 过程构建实现平方根倒数算法的实际逻辑。该过程通过参数列表获取乘法器，这样便可在 create_isqrt 过程之外，将同一硬件乘法器复用于平方根算法的其余部分。

procedure create_isqrt ( signal self : inout isqrt_record; signal multiplier : inout multiplier_record ) is variable mult_result : signed(int_word_length-1 downto 0); variable inverse_square_root_result : signed(int_word_length-1 downto 0); begin CASE self.state_counter is WHEN 0 => multiply_and_increment_counter(multiplier,self.state_counter, self.x, self.x); self.threehalfs <= self.x + self.x/2; WHEN 1 => multiply_and_increment_counter(multiplier,self.state_counter, self.x, self.sign_input_value); WHEN others => --do nothign end CASE; self.isqrt_is_ready <= false; CASE self.state_counter2 is WHEN 0 => if multiplier_is_ready(multiplier) then self.x_squared <= get_multiplier_result(multiplier, isqrt_radix); self.state_counter2 <= self.state_counter2 + 1; end if; WHEN 1 => if multiplier_is_ready(multiplier) then multiply(multiplier, self.x_squared, get_multiplier_result(multiplier,isqrt_radix)); self.state_counter2 <= self.state_counter2 + 1; end if; WHEN 2 => if multiplier_is_ready(multiplier) then mult_result := get_multiplier_result(multiplier,isqrt_radix+1); inverse_square_root_result := self.threehalfs - mult_result; self.x <= inverse_square_root_result; self.state_counter2 <= self.state_counter2 + 1; if self.loop_value > 0 then request_isqrt(self , self.sign_input_value , inverse_square_root_result , self.loop_value); else self.isqrt_is_ready <= true; self.result <= inverse_square_root_result; end if; end if; WHEN others => --do nothign end CASE; end create_isqrt; 如果我们的乘法器只有一个流水线周期，那么第一次乘法的结果就绪时间与第二次乘法被流水化送入乘法器的时间相同。

通过使用两个状态计数器，我们可以处理这种情况。我们还在第一个 state_counter 中安排 $1.5\cdot x$ 的计算。

isqrt 计算通过设置输入值并使用 request_isqrt 过程复位状态计数器来启动。两个状态机都按线性顺序逐步运行，并在运行完成后停止。该请求将两个状态机都置为零状态，这同时将状态复位到已知且已定义的位置。

procedure request_isqrt ( signal self : inout isqrt_record; input_number : signed; guess : signed; number_of_loops : natural range 1 to 7 ) is begin self.sign_input_value <= input_number; self.state_counter <= 0; self.state_counter2 <= 0; self.x <= guess; self.loop_value <= number_of_loops-1; end request_isqrt; 逆平方根的完整实现及其接口函数可在 fixed_isqrt_pkg.vhd 中找到。

该算法在一个测试平台中开发，测试平台可在此处找到。若初始值接近真实根，算法便会收敛。如前所述，该算法具有二次收敛性，这大致意味着每次迭代我们都会将正确位数翻倍。这种收敛速度非常快，以至于对于1.0到2.0之间的输入数值，一个常数初始值就足以在仅4次迭代中给出大约32位的精度。4次迭代的误差行为如图1所示。误差是在定点平方根与参考值之间计算的，参考值通过使用real_math库计算平方根获得。图中显示的是最大误差，在4次迭代时约为$1.8\cdot 10^{-10}$。

多个初始值

由于收敛非常快，且该算法只需要3次乘法，其中两次可以流水线化，我们在初始值方面面临一个有趣的情况。获得更好的初始猜测并不会带来好处，除非获得它的速度比再运行一次迭代更快。即使是非常粗糙的猜测也能让我们快速收敛，而我们期望的精度越高，就需要越高精度的初始猜测来避免通过算法进行一次迭代。图2显示了初始值为16种可能性之一时的误差行为。这使我们能够仅使用2次迭代就将误差降至低于$2.0\cdot 10^{-7}$，这接近32位浮点数的精度。