牛顿-拉夫逊迭代法原理与除法器的软件与硬件实现
牛顿-拉夫逊迭代法
牛顿迭代法是一种快速求解方程根的方法，它的迭代原理如下：
选择一个初始点（x0, f(x0) ），通过该点做切线。切线与X轴的交点记为x1，做为下一次迭代的点，再在（x1, f(x1) ）处做切线，依次类推，直道求解满足一定的精度即可。
我们以求解根的该曲线的近似解为例，其动态求解过程如下：
在进行叙述之前，需要对求近似解的原因做几点说明：
- 很多方程是无法准确求解的，如我们知道 ex=2 的解为，x= ln2. 但是你无法表达出ln2
- 计算机的硬件资源有限, 无法准确将数字完全的表达出来，而且有的时候完全没有必要
- 计算机是用二进制表达数据，有些数据用二进制可以准确表示，但是有的数据二进制格式只能近似表达
那么，在牛顿法在数学上该如何表达呢？
过曲线上一点做切线，求切线与X轴的交点，需要用到导数
假设切线的方程为：
y
=
k
x
+
b
y=kx+b
y=kx+b
那么过( xn, f(xn）)点的切线为
y
−
f
(
x
n
)
=
f
′
(
x
n
)
(
x
−
x
n
)
y-f(x_n)=f'(x_n)(x-x_n)
y−f(xn)=f′(xn)(x−xn)
令 y=0, 我们有
x
n
+
1
∗
f
′
(
x
n
)
=
f
′
(
x
n
)
∗
x
n
−
f
(
x
n
)
x_{n+1}*f'(x_n)= f'(x_n)*x_n-f(x_n)
xn+1∗f′(xn)=f′(xn)∗xn−f(xn)
有
x
n
+
1
=
x
n
−
f
(
x
n
)
f
′
(
x
n
)
x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}
xn+1=xn−f′(xn)f(xn)
停止迭代方法有很多。常用的方法有设定迭代次数和相对迭代相对误差小于一定精度，根据自己的情况选择合适的方法。
后续我们将以设计一个定点除法器为例，手把手带大家动手实践如何使用牛顿迭代法。不过在进入设计之前，很有必要了解计算机中数字的表示方法。
浮点和定点问题
计算机的首要工作是计算，计算需要则需要处理数据。那么计算机是如何表示数据呢？通常来说，在计算机领域编程中，你通常会发现数字的类型通常有两种：整数型和浮点型。这恰恰对应了两种数据的表示方式。
浮点
浮点方法表示数据常用在科学计数法，尤其是在当今AI的训练中。例如：
1.6E+19 =1.6 x 1019
9E-15 =9x10-15
由于二进制的特质，任何数据在都可以表示成如下形式：
1.
x
x
x
x
∗
2
y
y
y
y
y
y
1.xxxx *2^{yyyyyy}
1.xxxx∗2yyyyyy
而浮点之所以叫浮点，在运算前后数据的1后面的小数点的位置是会改变的
单精度浮点数与10进制数的转换关系如下：
31位：符号位，1代表负数，0代表正数
30~23位 ：指数位
22~0位： 小数位
bias 为偏置= 2指数位-1-1， 在单精度浮点数下，bias为 2^7-1=127
例如
符号位:1，该数据为负数
指数位:1000001（2）=129(10)
小数为：01000000… =0.25
所以该数据为
(
−
1
)
S
∗
(
1
+
0.25
)
∗
2
129
−
127
=
−
5.0
(-1)^S *(1+0.25)*2^{129-127}=-5.0
(−1)S∗(1+0.25)∗2129−127=−5.0
定点
定点格式的表示
在计算机中int类型是一种特殊的定点数据的表示方式，它直接忽略了小数位。然而在实际一定精度运算中，小数位是不能忽略的，那么该如何处理呢？接下来，定点格式提供解决该问题
具体定点数代表的数值以 FN 代为缩写，可由以下公式 4-5 给出，以字长 （Word length）表示定点数的整体长度，简写为wl；符号位的数值以s表示；以 Xi表示除了符号位以外每个位置上的数值；以分数长（Fractional length）表示定 点数中分数位的长度，简写为fl；以整数长（Integer length）表示定点数中整数 位的长度，简写为il
一般来说，定点数在嵌入式系统中更常见，因为使用定点方式进行算术时，尤其是进行乘除法时，资源消耗和延迟远比采用浮点方法更好。需要说明的是，其实在硬件并不能识别小数点，它默认数据都为整数，那为什么能够表示小数呢？这是设计者凭空想象出来的。计算机处理的结果经过移位操作，就可以得到设计者规定的定点模式。
定点格式的转换
- 数据的扩位与截位
一般来说，数据的扩位不会造成精度的影响； 数据的截位会造成误差
扩展
无符号数据（以8位整数位，8位小数位，即S0I8F8, Q8.8格式为例）
Q8.8–>Q16.16 ， 在原整数位最高位加入8bit，在小数位最低位再加8bit
有符号数据（S1I8F8为例）
Q8.8–>Q16.16–>如果原来的数为正数，在原最高整数位补8个零，小数位再补8bit；为负数的话，则需要在原最高整数位补8个1，小数再补8bit
截位
无符号
Q16.16–>Q8.8–> 截掉整数最高8位和低小数最低8位，整数需要判断是否进行溢出处理
有符号
Q16.16->Q8.8–> 小数位直接截掉即可，整数位在截取的时候需要判断是否溢出，如果溢出进行饱和处理，不然确保符号正确的情况下，截掉整数最高8位
比较
浮点
-
指数部分动态，范围很大
-
精度灵活，开发便捷（无需处理溢出和缩放问题）
-
运算起来复杂（需要对齐指数、规格化、舍入等）
-
需要使用浮点运算单元，能耗较大
定点
-
整数部分和小数部分固定，表示范围小
-
运算起来简单（移位，加减）
-
无需硬件支持，适合低功耗设备
-
存在确定性误差，容易溢出
Python
案例：使用牛顿–拉夫逊迭代法，对除法硬件进行描述
输入(dividend, divisor)数据均为S0I8F0, 输出(quotient)为S0I8F8
思路分析：
q
u
o
t
i
e
n
t
（商）
=
d
i
v
i
d
e
n
d
(
被除数
)
/
d
i
v
i
s
o
r
(
除数
)
quotient（商）=dividend(被除数)/divisor(除数)
quotient（商）=dividend(被除数)/divisor(除数)
可以将除法转换为乘法
q
u
o
t
i
e
n
t
=
d
i
v
i
d
e
n
d
∗
1
d
i
v
i
s
o
r
quotient= dividend*\frac{1}{divisor}
quotient=dividend∗divisor1
所以，我们可以使用牛顿迭代法，快速求出1/divisor的近似值，然后将其与dividend相乘，即可以获得结果
因为
$$
1/\frac{1}{divisor}=divisor
令
令
令
f(x)= \frac{1}{x}-divisor
$$
我们求f(x)=0的根即可
根据我们前面推到的公式有：
x
n
e
x
t
=
x
c
u
r
r
e
n
t
−
f
(
x
c
u
r
r
n
e
t
)
f
′
(
x
c
u
r
r
e
n
t
)
=
x
c
u
r
r
e
n
t
−
x
c
u
r
r
e
n
t
−
1
−
d
i
v
i
s
o
r
x
c
u
r
r
e
n
t
−
2
=
x
c
u
r
r
e
n
t
+
x
c
u
r
r
e
n
t
2
(
x
c
u
r
r
e
n
t
−
1
−
d
i
v
i
s
o
r
)
x_{next}=x_{current}-\frac{f(x_{currnet})}{f'(x_{current})} =x_{current}-\frac{x^{-1}_{current}-divisor}{x^{-2}_{current}} =x_{current}+x^2_{current}(x^{-1}_{current}-divisor)
xnext=xcurrent−f′(xcurrent)f(xcurrnet)=xcurrent−xcurrent−2xcurrent−1−divisor=xcurrent+xcurrent2(xcurrent−1−divisor)
$$
=2x_{current}-divisorx^2_{current}=x_{current}(2-divisorx_{current})
$$
代码
import math
#四舍五入
def round_half_up(x):
return int(math.floor(x + 0.5))
#牛顿迭代法
def newton_raphson_division(dividend, divisor, iterations=4):
if divisor == 0:
return None
x = 1.0 / max(divisor, 1e-9) #给定初始值
for _ in range(iterations):
x *= 2 - divisor * x
quotient = dividend * x
return round_half_up(quotient * (1 << 8))
#满精度除法，用于比较
def full_precision_divide(dividend, divisor):
if divisor == 0:
return None
quotient = (dividend / divisor) * (1 << 8)
return round_half_up(quotient)
def main():
methods = [
#我实验了很多除法器，此处只使用牛顿迭代法
#('direct_division', dm.direct_division),
#('restoring_remainder', dm.restoring_remainder_division),
#('non_restoring_division', dm.non_restoring_division),
#('lut_division', dm.lut_division),
#('binary_bit_method', dm.binary_bit_division),
('newton_raphson_division', newton_raphson_division),
]
error_counts = {name: 0 for name, _ in methods}
total_cases = 0
for dividend in range(256):
for divisor in range(1, 256):
total_cases += 1
expected =full_precision_divide(dividend, divisor)
for name, method in methods:
result = method(dividend, divisor)
if result != expected:
error_counts[name] += 1
print(f"Mismatch: {name} ({dividend}/{divisor}) Got: {result} Expected: {expected}")
print("\nError Summary:")
for name, count in error_counts.items():
print(f"{name}: {count} errors ({(count/total_cases)*100:.2f}%)")
if __name__ == "__main__":
main()
结果
错误率为0%
Verilog
软件语言和硬件语言区别在两个例子中，体现的就会非常明显：软件（如python）直接调用的专用浮点处理模块，因而实现起来很简答，因为复杂的部分，做硬件的人已经解决掉了。而硬件语言（Verilog）则需要考虑硬件的实际情况，对位宽精度进行处理。此外在硬件中，是没有小数的概念的，小数点的位置完全由人为控制。
除法器
`timescale 1ns / 1ps
module Newton_Raphson_Division (
input [7:0] dividend, // 被除数（Q8.0）
input [7:0] divisor, // 除数（Q8.0）
output reg [15:0] quotient // 商（Q8.8）
);
// 参数定义
localparam Q_FORMAT = 16; // 迭代中间值用Q16.16格式
localparam ITERATIONS = 4; // 迭代次数
// 牛顿迭代核心计算
reg [31:0] x; // Q16.16格式的倒数
reg [47:0]product; //Q24.24格式
reg [47:0]term;
reg [79:0]x_temp;
integer i;
always @(*) begin
// 初始值：1/divisor的近似值（查表简化版）
x = 32'hffff_ffff / {divisor, 16'h0000}; // 近似计算
// 四次牛顿迭代
for (i=0; i<ITERATIONS; i=i+1) begin
product = (divisor<<8)* x; // (Q8.0-->Q8.8）* Q16.16 = Q24.24
term = 48'h200_0000-product; //2 - D*x, 定点格式需要对齐 Q24.24
x_temp = x * term; // x = x * term = Q40.40
if(x_temp[79:56]!=0) //x_temp的高整数位溢出
x={16'hffff,x_temp[39:24]}; //截取Q16.16
else
x=x_temp[55:24];
end
end
// 最终乘法与格式转换
reg [39:0] final_product;
always @(*) begin
if(divisor==0)quotient=16'hffff;
else if (divisor==1)quotient=dividend<<8;
else begin
final_product = dividend * x; // 实际乘法,Q8.0 * Q16.16 = Q24.16
//quotient = final_product[31:16] + |final_product[15:0]; // 四舍五入到Q8.8
quotient=final_product[23:8]+(final_product[7]?1:0);
end
end
endmodule
Testbench
`timescale 1ns / 1ps
module div_test();
reg [7:0] dividend;
reg [7:0]divisor;
wire [15:0]quotient;
initial begin
#5;
dividend=8'hFF;
divisor =8'h1;
#10;
dividend=8'h03;
divisor =8'h0A;
#10;
dividend=8'h08;
divisor =8'h03;
#10;
dividend = 8'd145;
divisor = 8'd252;
#10;
dividend = 8'd255;
divisor = 8'd239;
end
Newton_Raphson_Division uut(
.dividend(dividend), // 被除数（S0I8F0）
.divisor(divisor), // 除数（S0I8F0）
.quotient(quotient) // 商（S0I8F8）
);
endmodule
结果
- 被除数为ff（255） ，除数为01（1） , 商为ff00，高两位16进制数代表整数，低两位16进制数，结果为255.00
- 被除数为03（3） ，除数为0a（10）, 商为004d, 即0.30078125 ， 真实值为0.3333333
- 被除数为08（8） ，除数为03（03）, 商为02ab, 即2.66796875 ， 真实值为2.6666666
- 被除数为91（145），除数为fc(252) , 商为0093， 即0.57421875 ， 真实值为0.5753968
- 被除数为ff(255) , 除数为ef(239) , 商为0111， 即1. 06640625， 真实值为1.0669456
可以看出，定点计算存在一定误差，但在资源有限下的情况下，采用定点运算可以达到不错精度的效果
，
- 除数为01(1), 商为ff00，高两位16进制数代表整数，低两位16进制数，结果为255.00
- 被除数为03(3), 除数为0a（10）, 商为004d, 即0.30078125 ， 真实值为0.3333333
- 被除数为08(8),除数为03（03）, 商为02ab, 即2.66796875 ， 真实值为2.6666666
- 被除数为91(145)，除数为fc(252) , 商为0093， 即0.57421875 ， 真实值为0.5753968
- 被除数为ff(255), 除数为ef(239) , 商为0111， 即1. 06640625， 真实值为1.0669456
可以看出，定点计算存在一定误差，但在资源有限下的情况下，采用定点运算可以达到不错精度的效果