CN103092561B - 一种基于除数映射的Goldschmidt除法实现方法 - Google Patents

Abstract

本发明公开了一种基于除数映射的Goldschmidt除法实现方法，首先，把浮点形式的被除数N_f和除数D_f规格化为f×2^e的形式，规格化后的被除数和除数记为N和D；根据给定的最小相对误差E以及迭代次数M求出分界值p；若规格化后的除数D落在[1,p]区间内，则直接进行M次迭代；如果除数D落在[p,2)区间内，则将D映射到[1,p]区间内，然后，再进行M次迭代。迭代时，初始值F₀=2-D₀。M次迭代得到f部分的相除结果，最后将f部分的相除结果与2^e部分的相减组合起来，得到最终的除法运算结果。该方法不需要初始估计值，从而能够节省大量的存储资源。

Description

一种基于除数映射的Goldschmidt除法实现方法

技术领域

本发明属于信号处理领域，涉及一种基于除数映射的Goldschmidt除法实现方法。

背景技术

浮点除法被越来越多地应用到处理器设计当中，成为信号处理的一个重要部分。虽然除法器没有加减法器的使用频率高，但是除法器的性能直接影响了信号处理的效果，所以加快除法器的运算速度和减少硬件资源成为了一个关键问题。

现代商业处理器大都采用并行乘法器来提高除法的速度，利用乘法的迭代来实现除法运算，且所需精度越高，迭代的次数越多。现在主要的基于乘法的迭代式除法有M.J.Flynn在Ondivisionbyfunctionaliteration一文中提出的Newton-Raphson算法和R.E.Goldschmidt在Applicationsofdivisionbyconvergence一文中提出的Goldschmidt算法。

Newton-Raphson算法把除法转化为求除数的倒数，然后通过与被除数相乘得到商，而Goldschmidt算法通过规格化除数和被除数然后直接求得商。这两个算法每次迭代都涉及到两次乘法操作，但是在Newton-Raphson算法中这两次乘法操作不能并行操作，只能进行串行操作，而在Goldschmidt算法中这两个乘法可以并行处理，所以采用Goldschmidt算法的处理速度会比较快。

Goldschmidt算法在进行除法运算时，Q=N/D，其中N为被除数，D为除数，Q为商，通过构造一组因子F₁，F₂，…，F_M，令分子和分母分别与这一组因子相乘，当分母趋近于1时，分子即为所要求的除法结果。分子和分母与这组因子相乘是通过M次迭代实现的。

传统的Goldschmidt除法运算的迭代公式可以写成如下形式：

$\left\{\begin{array}{c}{F}_{i+1}=2-D_{i+1}\\ N_{i+1}=N_i{F}_i\\ D_{i+1}=D_i{F}_i\end{array}\right.---\left(1\right)$

其中，i＝1,2,…M-1，M代表迭代次数，N_i和D_i分别代表第i次迭代之后的被除数和除数，F_i代表因子。这里N₀和D₀为初始要计算的被除数和除数，而F₀是由除数的倒数估计出来的，但是估计的精度有限，设初始估计值的误差为s，则F₀=(1-s)/D，经过推导可以知道传统的Goldschmidt除法算法的计算误差为

由上面分析可以看出，传统Goldschmidt算法迭代结果的误差取决于初始值的估计误差s，即初始值估计精度的高低决定了迭代结果误差的大小。

针对这一问题，D.DasSarma和D.W.Matula在FaithfulBipartiteROMReciprocalTables一文，以及M.J.Schulte和J.E.Stine在SymmetricBipartiteTablesforAccurateFunctionApproximation一文中提出了多种查找表法来减少初始估计值误差以提高迭代结果的精度。

这种方法的原理是在ROM中存储1/D的初始值，为了提高运算的精度，就需要增加查找表的位数，设查找表位数为n，则查找表的资源为n·2ⁿ，即所需要的资源将随着查找表的位数以指数方式增加。

综上所述，在精度要求比较高的情况下，查找表的存储资源将随着位数急剧增大，为了提高除法器的性能，减小存储资源成为了一个亟待解决的问题。

发明内容

有鉴于此，本发明为了解决ROM存储资源大的问题，提出一种改进的Goldschmidt除法实现方法，这种除法算法不需要初始估计值，从而能够节省大量的存储资源。

为了解决上述技术问题，本发明是这样实现的：

一种基于除数映射的Goldschmidt除法实现方法，该方法包括如下步骤：

步骤1、把浮点形式的被除数N_f和除数D_f规格化为f×2^e的形式，其中，f为尾数部分且f∈[1,2)，e为指数部分；规格化后的被除数和除数记为N和D；则浮点数据相除转化为尾数部分相除和指数部分相减；

步骤2、在尾数部分相除时，根据给定最小相对误差E以及在时延满足要求的情况下所需资源最少的迭代次数M的值，采用式（I）求出分界值p；

$p=\exp \left(\frac{\ln E}{2^M}\right)+1---\left(I\right)$

步骤3、判断规格化后的除数D和分界值p的大小；若规格化后的除数D落在[1,p]区间内，则直接令N₀=N，D₀=D，执行步骤4；若规格化的除数D落在[p,2)区间内，则通过乘以映射系数将D映射到[1,p]区间，同时将被除数N也乘以所述映射系数，将映射后的D和N赋值给D₀和N₀，然后再执行步骤4；

步骤4、通过式（IV）进行M次迭代得到尾数部分的除法运算结果；

$\left\{\begin{array}{c}{F}_i=2-D_i\\ N_{i+1}=N_i{F}_i\\ D_{i+1}=D_i{F}_i\end{array}\right.---\left(\mathrm{IV}\right)$

其中，i=0,2,…,M-1，N_i和D_i分别代表第i次迭代之后的被除数和除数；

步骤5、将步骤4得到的尾数部分相除结果与指数部分相减组合起来，得到最终的除法运算结果。

优选地，步骤3中，所述通过乘以映射系数将D映射到[1,p]区间，同时将被除数N也乘以所述映射系数的步骤为：

步骤31、根据式（II）求出最小映射分段数T的值，则[p,2)被分为了T段；

T＝ceil(log_p2-1)（II）

其中，ceil代表向上取整；

步骤32、通过式（III）求出各映射范围的边界以及映射系数，找到除数D所在的映射范围，提取对应的映射系数；将除数D和被除数N均与提取的对应映射系数相乘；

$\left\{\begin{array}{c}{\mathrm{Range}}_{(j,\max )}=\frac{2}{p^{(T-1-j)}}\\ {\mathrm{Range}}_{(j,\min )}=\frac{2}{p^{T-j}}\\ c\left(j\right)=\frac{p^{T-j}}{2}\end{array}\right.---\left(\mathrm{III}\right)$

Range_（j,min）和Range_（j,max）分别为第j个映射范围的下限和上限，c(j)为第j个映射范围的映射系数；j=0,1,…,T-1。

优选地，预先对映射系数进行CSD编码，并存储；将除数D和被除数N与对应映射系数相乘时，采用CSD编码来进行固定乘数乘法。

有益效果：

（1）本发明提出的这种改进的Goldschmidt除法算法，不需要对F初始值的估计，能够节省大量的存储资源。

（2）本方法迭代结果的误差也只是和除数与迭代次数有关，使得运算结果误差可控。

（3）本方法还利用固定数最小资源乘法算法来减少除数映射算法的资源占用。

（4）本方法易于硬件的并行实现。

附图说明

图1（a）和图1（b）为除数映射过程示意图。

图2为改进Goldschmidt算法实现流程图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

首先将除法运算的迭代公式进行修改，写成如下形式：

$\left\{\begin{array}{cc}{F}_i=2-D_i& \left(2\right)\\ N_{i+1}=N_i{F}_i& \left(3\right)\\ D_{i+1}=D_i{F}_i& \left(4\right)\end{array}\right.$

其中，i＝1,2,…M-1，M代表迭代次数，N_i和D_i分别代表第i次迭代之后的被除数和除数。这里N₀=N代表最初的被除数，D₀=D代表最初的除数，而且因子的初始值F₀不需要估计，可以直接由公式（2）计算得到，这样就可以节省存储初始估计值F₀的ROM空间。

将式（2）带入式（4），并进行迭代推导，可以得到

$\begin{array}{c}{D}_{i+1}-1=D_i(2-D_i)-1\\ =-{(D_i-1)}^2\\ =-{(D-1)}^{2^{i+1}}\end{array}---\left(5\right)$

根据式（3）和式（4），可以推导出D_i和N_i的迭代关系

$N_{i+1}=\frac{N}{D}{D}_{i+1}---\left(6\right)$

把式（5）代入式（6）可得

$N_{i+1}=\frac{N}{D}[1-{(D-1)}^{2^{i+1}}]---\left(7\right)$

经过M次迭代之后，分母D_M≈1，故

$Q_M=\frac{N_M}{D_M}\≈D_M=\frac{N}{D}[1-{(D-1)}^{2^M}]---\left(8\right)$

其中，Q_i为第i次迭代后之后的商。

定义相对误差

${\mathrm{\ϵ}}_{Q_i}=\mathrm{abs}\left(\frac{Q_i-(N/D)}{N/D}\right)---\left(9\right)$

其中，则在M次迭代下的除法结果相对误差为

${\mathrm{\ϵ}}_{Q_M}={(D-1)}^{2^M}---\left(10\right)$

从公式（10）可以看出，改进的除法运算结果的相对误差只与除数的大小以及迭代次数有关。

对于浮点***来说，把N和D用科学计数法表示，即把被除数和除数都规格化为f×2^e的形式，例如2.5被规格化后为1.25×2¹，其中f为尾数部分且f∈[1,2)，e为指数部分，浮点数据的除法运算即尾数相除、指数相减，所以这里只需要考虑尾数部分的除法运算即可。则经过规格化之后D和N的值范围是[1，2)，在这个范围内，当迭代次数M固定时，相对误差函数随着D的增大为单调递增函数，除数D越接近于2则相对误差越大。

假设在迭代次数固定为M时，给定要求的参数E使得相对误差不大于E，即

${(D-1)}^{2^M}\≤E---\left(11\right)$

可以推出

$D\≤\exp \left(\frac{\ln E}{2^M}\right)+1=p---\left(12\right)$

那么p就是满足相对误差要求时除数最大值，且，p∈[1，2)。也就是说，只有当除数D∈[l，p)时，才能满足相对误差的要求。

但并不是所有除数D都能满足D∈[1，p)的要求，因此为了让所有范围内的除数都能够满足误差要求，对于D∈(p，2)的数据，采取将其映射到[1,p]范围内的方法来满足精度要求，如图1所示。

当p>3/2时，只需要一个统一的映射系数则可把（p,2）范围内的数映射到[1,p]区间内。对于一次映射，当D∈(p，2)时，定义一个映射系数γ∈(0.5，1)，有D′＝D×γ∈[1,p]，D’为映射后的D。则用N×γ和D×γ分别作为映射后的被除数和除数，则除数被映射到[1,p]范围内。此时，经过M次迭代后的相对误差将由公式（11）保证满足精度要求。

但是当p<3/2时，前面提及的统一映射系数并不能将所有（3/2,2）范围内的数映射到[1,p]区间内，因此需要将（3/2,2）分为T段，每段设定一个映射系数c，根据D的所在区段选取适当的映射系数进行映射。假设将(p,2)被分为T段，其中相邻的两段可以有交集，记第j段的映射范围为(Range_(j，min),Range_(j，max)]，j＝0,1,...,T-1，第T-1段的映射范围为(Range_(T-1,min),2)，Range表示映射范围的边界，则有如下特性

$\left\{\begin{array}{c}1\&le;{\mathrm{Range}}_{(0,\min )}\&le;p\\ {\mathrm{Range}}_{(j-1,\max )}\&GreaterEqual;{\mathrm{Range}}_{(j,\min )}\\ {\mathrm{Range}}_{(T-1,\max )}=2\end{array}\right.---\left(13\right)$

记C＝{c(0),c(1),...,c(T-1)}为每一段映射范围的映射系数，经过映射之后有如下关系

$\left\{\begin{array}{c}c\left(j\right){\mathrm{Range}}_{(j,\max )}\&le;p\\ c\left(j\right){\mathrm{Range}}_{(j,\min )}\&GreaterEqual;1\end{array}\right.---\left(14\right)$

当所有的映射范围都不重叠的时候，所需要的映射分段数T最小，故令

$\left\{\begin{array}{c}c\left(j\right){\mathrm{Range}}_{(j,\max )}=p\\ c\left(j\right){\mathrm{Range}}_{(j,\min )}=1\end{array}\right.---\left(15\right)$

则式（13）可化为

$\left\{\begin{array}{c}{\mathrm{Range}}_{(0,\min )}=p\\ {\mathrm{Range}}_{(j-1,\max )}={\mathrm{Range}}_{(j,\min )}\\ {\mathrm{Range}}_{(T-1,\max )}=2\end{array}\right.---\left(16\right)$

结合公式（15）和（16），可以求得各映射范围以及映射系数，如下所示

$\left\{\begin{array}{c}{\mathrm{Range}}_{(j,\max )}=\frac{2}{p^{(T-1-j)}}\\ {\mathrm{Range}}_{(j,\min )}=\frac{2}{p^{T-j}}\\ c\left(j\right)=\frac{p^{T-j}}{2}\end{array}\right.---\left(17\right)$

再由前面所述的公式（13）可以得到，

$1\&le;{\mathrm{Range}}_{(0,\min )}=\frac{2}{p^T}\&le;p---\left(18\right)$

则有

log_p2-1≤T≤log_p2（19）

因此，映射分段数T的取值由（19）决定。考虑到T为整数，同时为了选取最小的T，这里取

T＝ceil(log_p2-1)（20）

其中ceil代表向上取整。

通过上述除数算法的推导可以看到，在固定的M次迭代条件下，只要在迭代之前对除数和被除数在不同的范围内乘以不同的映射系数，就可以通过一次映射运算使所有范围除数的运算结果都达到精度要求。

对于要进行M次迭代的改进Goldschmidt除法算法来说，根据式（2，3，4）可知，每次迭代需要2M个乘法和M个加法操作（在这里认为加减法占用的资源相同），而在硬件实现浮点运算***中，特别是高精度计算，相对于乘法器来说，加法器可以忽略不计，因此在这里只考虑2M个乘法器。

同样还需要考虑除数映射算法中对各个映射范围的映射系数的乘法，由于分子和分母都需要与映射系数相乘，故需要两次乘法。但考虑到该情况下是对固定系数的乘法，通过对映射系数进行CSD（CanonicalSigned-DigitCode，规范符号代码）编码可以减少乘法硬件资源占用量，定义β(β∈[0，1]）为固定乘数乘法对一般乘法占用资源比，表示一个固定乘数乘法单元占用硬件资源为一个一般乘法单元占用资源的β倍。

定义本算法在M次迭代下，相对误差小于E时所需要的乘法器个数为R(M)，则

R(M)＝2M+2β（21）

当采用CSD编码来进行固定乘数乘法时，β的平均值为0.3。Macleod在Useofminimum-addermultiplierblocksinFIRdigitalfilters一文中提出一种固定乘数最小资源乘法算法，比CSD乘法降低25%的资源，此时β的平均值为0.25。

由（21）式可以看出，当给定精度要求时，可以求得函数R（M）的最小值，也就是说在给定误差精度以后，本发明所推导的算法可以求得满足该误差精度所需要的最小迭代次数，以及占用资源最小的情况下所需要的所有参数。

基于上述推导，本发明改进Goldschmidt除法运算的核心思想是，首先，把浮点形式的被除数和除数规格化为f×2^e的形式，规格化后的被除数和除数记为N和D；求出分界值p，若D落在[1,p]区间内，则直接利用公式（2）～（4）进行M次迭代；若D落在[p,2)区间内，则通过乘以映射系数将D映射到[1,p]区间，同时将被除数N也乘以所述映射系数，再进行M次迭代。M次迭代得到f部分的相除结果，最后将f部分的相除结果与2^e部分的相减组合起来，得到最终的除法运算结果。

参见图2，其具体实现过程如下：

步骤一、把浮点形式的被除数N_f和除数D_f规格化为f×2^e的形式，则浮点数据相除转化为尾数部分相除和指数部分相减，则需要进行除法运算的被除数N和除数D的取值范围是[1,2)。

$\begin{array}{c}{N}_f=N\&times;2^{e1}\\ D_f=D\&times;2^{e2}\end{array}$

步骤二、给定最小相对误差E以及在时延满足要求的情况下所需资源最少的迭代次数M的值。

步骤三、根据式（12）求出分界值p，则区间[1,2)被分为了[1,p]和[p,2)，然后判断规格化的除数D和分界值p的大小。若规格化后的除数落在[1,p]区间内，则令N₀=N，D₀=D，跳转到步骤六；若规格化的除数落在[p,2)区间内，则顺序向下执行。

步骤四、根据式（20）求出最小映射分段数T的值，则[p,2)被分为了T段。

步骤五、通过式（17）求出各映射范围的边界以及映射系数，找到除数D所在的映射范围，提取对应的映射系数；将除数D和被除数N均与提取的对应映射系数相乘。将映射后的除数D和N赋值给D₀和N₀，转入步骤六。

步骤六、通过式（2）～（4）进行M次迭代得到除法运算结果。

步骤七、将步骤六得到的尾数部分相除结果与指数部分相减组合起来，得到最终的除法运算结果。

至此，本流程结束。

以上流程是基于所有的映射范围都不重叠，也就是式（16）的情况，此时所需要的映射分段数T最小。在实际中，也可以不限定这么严格，只需要满足式（13）即可，此时映射系数选取满足式（14）的数值即可。

本发明推导出了该改进算法的迭代结果相对误差，并得到了所需硬件资源与精度要求之间的关系，使得在满足精度要求的范围内，得到占用资源最小的除法器结构。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种基于除数映射的Goldschmidt除法实现方法，其特征在于，该方法包括如下步骤：

步骤1、把浮点形式的被除数N_f和除数D_f规格化为f×2^e的形式，其中，f为尾数部分且f∈[1,2)，e为指数部分；规格化后的被除数和除数记为N和D；则浮点数据相除转化为尾数部分相除和指数部分相减；

步骤2、在尾数部分相除时，根据给定最小相对误差E以及在时延满足要求的情况下所需资源最少的迭代次数M的值，采用式（I）求出分界值p；

$p=\exp \left(\frac{\ln E}{2^M}\right)+1---\left(I\right)$

步骤3、判断规格化后的除数D和分界值p的大小；若规格化后的除数D落在[1,p]区间内，则直接令N₀=N，D₀=D，执行步骤4；若规格化的除数D落在[p,2)区间内，则通过乘以映射系数将D映射到[1,p]区间，同时将被除数N也乘以所述映射系数，将映射后的D和N赋值给D₀和N₀，然后再执行步骤4；

步骤4、通过式（IV）进行M次迭代得到尾数部分的除法运算结果；

$\left\{\begin{array}{c}{F}_i=2-D_i\\ N_{i+1}=N_i{F}_i\\ D_{i+1}=D_i{F}_i\end{array}\right.---\left(IV\right)$

其中，i=0,2,…,M-1，N_i和D_i分别代表第i次迭代之后的被除数和除数；

步骤5、将步骤4得到的尾数部分相除结果与指数部分相减组合起来，得到最终的除法运算结果。

2.如权利要求1所述的方法，其特征在于，步骤3中，所述通过乘以映射系数将D映射到[1,p]区间，同时将规格化后的被除数N也乘以所述映射系数的步骤为：

步骤31、根据式（II）求出最小映射分段数T的值，则[p,2)被分为了T段；

T＝ceil(log_p2-1)（II）

其中，ceil代表向上取整；

步骤32、通过式（III）求出各映射范围的边界以及映射系数，找到规格化后的除数D所在的映射范围，提取对应的映射系数；将规格化后的除数D和规格化后的被除数N均与提取的对应映射系数相乘；

$\left\{\begin{array}{c}{\mathrm{Range}}_{(j,max)}=\frac{2}{p^{(T-1-j)}}\\ {\mathrm{Range}}_{(j,min)}=\frac{2}{p^{T-j}}\\ c\left(j\right)=\frac{p^{T-j}}{2}\end{array}\right.---\left(III\right)$

Range_（j,min）和Range_（j,max）分别为第j个映射范围的下限和上限，c(j)为第j个映射范围的映射系数；j=0,1,…,T-1。

3.如权利要求1或2所述的方法，其特征在于，预先对映射系数进行CSD编码，并存储；将规格化后的除数D和规格化后的被除数N与对应映射系数相乘时，采用CSD编码来进行固定乘数乘法。