FPGA实现S15.16定点牛顿-拉夫逊除法近0输入的初值优化咨询

现有实现基础 #

当前已在FPGA平台完成S15.16格式牛顿-拉夫逊除法算法开发，整体支持输入取值范围为[-16:16]。
针对输入幅值处于[1:16]区间的使用场景，现有方案采用各取值区间中点取倒数的方式设置迭代初值$a_0$，对应公式为：
$$a_0 = \frac{2}{Min+Max}$$
典型初值配置示例：

• 1 ≤ x < 2 区间取 $a_0 = 2/3$
• 6 ≤ x < 7 区间取 $a_0 = 2/13$
• 15 ≤ x < 16 区间取 $a_0 = 2/31$
  该初值方案仅需3次迭代即可实现10e-9的MSE，精度表现符合预期。

[0:1]输入区间的待解问题 #

目前算法的短板集中在x ∈ [0:1]的输入区间：
参考公开方案的处理思路，针对该区间的通用流程为：通过比特移位将除数D缩放至0.5 ≤ D ≤ 1的归一化区间，采用线性近似式 $X_0 = T_1 + T_2 D \approx \frac{1}{D}$ 作为迭代初值。
为最小化[0.5,1]区间内近似值的最大绝对误差，当前已选用针对该区间优化的线性近似公式：
$$X_0 = \frac{48}{17} - \frac{32}{17}D$$
该公式在[0.5:1]区间可达到较好的近似精度，但仍存在两个核心问题待明确：

• 当输入值逐步趋近于0（如0.1、0.01、0.001、0.00001等不同量级）时，是否需要为[0.01:0.1]、[0.001:0.01]等逐段细分区间分别配置适配初值？
• 针对幅值极小的输入值，是否更适合采用Goldschmidt除法算法替代现有牛顿-拉夫逊实现方案？

参考仿真代码 #

以下为已实现的S15.16定点牛顿-拉夫逊除法仿真代码，供方案分析参考：

i = 0
# 16 fractional bits
SHIFT = 2 ** 16 
# Lut with "optimal?" values to init NRD algorithm
LUT = np.round(1 / np.arange(0.5, 16, 1) * SHIFT).astype(np.int64)
LUT_f = 1 / np.arange(0.5, 16, 1)
# Function to simulates the NRD algorithm in S15.16 over a FPGA
def FIXED_RECIPROCAL(x):
    # Smart adressing to the initial iteration value
    adress = x >> 16
    # Get the value from LUT
    a0 = LUT[adress]
    # Algorithm with only 3 iterations
    for i in range(3):
        s1 = (a0*a0) >> 16
        s2 = (x*s1) >> 16
        a0 = (a0 << 1) - s2
    # Return rescaled value (Only for analysis purposes)
    return(a0 / SHIFT)

# ----- TEST ----- #
t = np.arange(1, 16, 0.1)
teor = 1 / t
t_fixed = (t * SHIFT).astype(np.int32)
prac = np.zeros(len(t))
for value in t_fixed:
    prac[i] = FIXED_RECIPROCAL(value)
    i = i + 1

# Get and print Errors
errors = abs(prac - teor)
mse = ((prac - teor)**2).mean(axis=None)

print("Max Error : %s" % ('{:.3E}'.format(np.max(errors))))
print("MSE:      : %s" % ('{:.3E}'.format(mse)))

# Print the obtained values:
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.plot(t, teor, label='Theorical division')
plt.plot(t, prac, '.', label='Newton-Raphson Division')
plt.legend(fontsize=16)
plt.title('Float 32 theorical division Vs. S15.16 Newton-Raphson division', fontsize=22)
plt.xlabel('x', fontsize=20)
plt.ylabel('1 / x', fontsize=16)
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)

内容的提问来源于stack exchange，提问作者Diego Ruiz