基于FPGA的数字信号处理3.7开平方运算分析

开平方运算（Non-Restoring）

在书《基于FPGA的数字信号处理》中的开平方运算一节中，作者给出了一种基于Non-Restoring的算法，但是没有给出具体推导过程，本文试图给出算法的简单证明。熟悉十位数的开方运算的读者可以跳转到第二部分

十位数的开平方运算

对于十位数的开方运算，我们有许多算法，比如泰勒展开、牛顿迭代法、利用常数0x5f375a86加速收敛的牛顿迭代法等。在书中使用的是Non-Restoring算法，本文就专门讨论这个算法。
对于一个6位的十位数D，它的每一位用dn表示，即
$$D=d_5{d}_4{d}_3{d}_2{d}_1{d}_0$$
它的n次根号的位数为 6 / n位，不严谨的简单证明如下：
一个数字a的整数部分长度为n，小数部分长度为无限长，它的平方的取值范围即为$[10^n,10^{n+2}]$。

类比除法的运算规则，我们可以利用迭代的减法操作，从高位到低位逐级算出解。
回顾除法运算，例如除法算式 $1354/3$，我们有如下操作
$$\begin{array}{lrc}& 151& \\ 3& \overline{)1354}& \\ & \underline{12}& \\ & 15& \\ & \underline{15}& \\ & 4& \\ & \underline{3}& \\ & 1& \end{array}$$
为了取得每一位的解，在每一位上，我们都需要找出数字n，n满足两个条件：

1. n是正整数且只有一位，对于十进制取值范围为{0,1,2,3,4,5,6,7,8,9}。
2. {上一次运算的余项，本次运算的被除数对应位} - n * 除数 的值 为最小值。

我们来看看根号的手解法。
$$D=d_5{d}_4{d}_3{d}_2{d}_1{d}_0$$
$$Q=\sqrt{D}=q_2{q}_1{q}_0=q_2\ast 10^2+q_1\ast 10+q_0$$
当我们要计算Q时，$q_2$的值可以由$d_5{d}_4$直接平方后试减得到，如$\sqrt{598213}$的百位为7。
可以把算出来的数字设置为a（上句话运算结果为a=700），下一位待计算的值为b。我们可以利用$(a+b{)}^2-a^2=2ab+b^2$试出数字b，b满足两个条件：

1. b是正整数，对于十进制，最高位取值范围为{0,1,2,3,4,5,6,7,8,9}，其余位为0。
2. $2ab+b^2$ - {上一次运算的余项，本次运算的被除数对应位} 的值 为最小值。

如果a、b的值不包含未计算部分,则上述条件将会有些许变化：

1. b是正整数且只有一位，对于十进制，取值范围为{0,1,2,3,4,5,6,7,8,9}
2. {上一次运算的余项，本次运算的被除数对应位} $20ab+b^2$ 的值 为最小值。

一个例子是
$$\begin{array}{lrc}& 773& \\ & \sqrt{598213}& \\ & \underline{49}& \\ & 1082& \\ & \underline{1029}& \\ & 5313& \\ & \underline{4629}& \\ & 684& \end{array}$$
$$\sqrt{795213}=773^2+684$$

a	前次余项+本次运算块	$20ab+b^2$	b
0	59	49	7
7	1082	1029	7
77	5313	4629	3
773	684	-	-

对于n次根号，我们需要将被n次开方处理的数字每n位分割。将上述迭代减法的计算式换成$(a+b{)}^n-a^n$即可。

fpga上的开平方运算

针对上述算法我们可以进行电路实现。首先面临的就是进制转换问题。二进制数字的平方与十进制数有很多相似之处，比如：

1. 一个定点数二进制数字a的整数部分（无符号位）长度为n，小数部分长度为无限长，它的平方的取值范围即为$[2^n,2^{n+2}]$
2. 它的n次根号的位数为 6 / n位
3. 迭代计算中用到的式子同样是$(a+b{)}^2-a^2$（a、b包括未计算部分），第二种形式为$4ab+b^2$（a、b不包括未计算部分）

不过，二进制数有一些有利的性质，我们可以利用他们来简化大量计算，如迭代用的式子$4ab+b^2$，用二进制表示可以优化为$(100{)}_{b}ab+(b{)}^2=(ab0b{)}_b$ 其中a为一串二进制数，b为一位二进制数。

当迭代的减法式子写作$(ab0b{)}_b$时：
$$\begin{array}{lrc}& 1101& \\ & \sqrt{10101001}& \\ & \underline{01}& \\ & 0110& \\ & \underline{101}& \\ & 0110& \\ & \underline{0000}& \\ & 11001& \\ & \underline{11001}& \\ & 0& \end{array}$$
$$\sqrt{(10101001{)}_b}=(1101{)}_b^2$$

a	前次余项+本次运算块	$(ab0b{)}_b$	b
0	010	001	1
1	0110	101	1
11	0110	0000	0
110	11001	11001	1
1101	0	-	-

此计算过程可以用恢复余数算法表示。

当相减后b<0时，需要恢复相减之前的值。

Restoring 算法很直观，但存在的问题是：若当前部分余数为负时，需要额外的加法操作使其恢复到前一次迭代时产生的部分余数。这不仅需要额外的资源，还增加了Latency。是否可以不恢复余数呢？这需要对算法进行一定的修改。
由式（3.62）可知，式（3.63）可重写为$$r_0=r_{1}D(1)D(0)+q_{2}0100-q_{1}01$$

书中的开平方运算

在《基于FPGA的数字信号处理》一书中，用到的算法（Non Restoring）如下：

简单分析可知，书中算法比我们的算法大致相同，不过为了保持运算的一致性，将减法操作统一为加法。图片中有一个错误，即$0S_3{S}_{2}11$应该为$0S_3{S}_2{S}_{1}11$
看懂这个算法需要两个小知识点：

1. A - B = A + ( -B) ,等于A的补码加(-B)的补码
2. if (A >= B) A+补(B)将会产生进位
   else if(A < B) A+补(B)将不会产生进位

我们可以根据 $补(余数:当前项)+迭代式$ 是否产生进位来决定下一步迭代的迭代式选取。Q(k+1)代表前一次运算的进位情况，补码情况下
由前文可知，
$$r_k=\{\begin{array}{ll}{r}_{k+1}D(2k+1)D(2k)-q_{k+1}01& \text{if Q(k+1)==1 }\\ r_{k+1}D(2k+1)D(2k)+q_{k+1}11(& \text{if Q(k+1)==0}\end{array}$$
如果使用补码运算，式子将会修改如下+
$$r_k=\{\begin{array}{ll}{r}_{k+1}D(2k+1)D(2k)+0q_{k+1}11& \text{if Q(k+1)==0}\\ r_{k+1}D(2k+1)D(2k)+1{\tilde{q}}_{k+1}11(& \text{if Q(k+1)==1}\end{array}$$

$q_k$	$r_{k+1}D(2k+1)D(2k)$	$Q(k+1)?0q_{k+1}11:0q_{k+1}11$	进位情况
-	010	111	1
1	0110	$1{\tilde{S}}_{3}11=1011$	1
11	00110	$1{\tilde{S}}_3{\tilde{S}}_{2}11=10011$	0
110	100101	$0S_3{S}_2{S}_{1}11=011011$	1
1101	0	-	-