平方根计算的硬件设计

本文最后更新于：2024年12月20日 下午

用verilog实现开方运算是一个不断逼近迭代的过程，常见的迭代的算法有牛顿迭代法、二分迭代法、逐次逼近法、CORDIC。

将一个32位正整数开方后将会得到一个16位正数，32位能够表示的最大正整数为2^32-1=4294967295，按理来说，这个也应该是需要迭代次数最多的一个数。

牛顿迭代法

牛顿迭代法的核心思想是Taylor展开式的高阶近似，最开始是用来迭代函数零点的。

  • 令方程f(x)=0的的近似根为xk，也就是初始迭代点，xk任取即可，只是根据取的点的好坏最终会影响迭代次数。
  • 接着取函数f(x)在xk处的一阶Taylor展开式p(x)来代替f(x)：

$ p(x)=f(x_k)+f^{'}(x_k)(x-x_k)f(x) $

  • 那么待求转化为：

$ f(x_k)+f^{'}(x_k)(x-x_k)=0 $

  • 求得该线性方程的根记为xk+1，这就得到了f(x)=0的一个新的近似根，可以表示为：

$ x_{k+1}=x_k- $

  • 上式本质上就是牛顿迭代公式，称为Newton迭代格式，与之对应的Newton迭代函数是：

$ (x)=x- $

按照这个迭代函数求函数f(x)零点的过程就称为牛顿迭代法，它的几何解释如下图，就是一个不断作切线逼近的过程。

将牛顿迭代函数稍加变形，就可以得到开方过程的迭代方法：

开方的函数表示是

$ y=f(x)= 若x=x_0，求y_0，y_0=f(x_0)= $

其反函数为：

$ y^2=x, x=y^2 $

换元可得：

$ y=x^2 $

那么待求问题转换为，已知y=y0，求x0，可以运用牛顿迭代法求解：

令

$ F(x)=x^2-y_0 $

问题简化为求F(x)的零点，将F(x)代入牛顿迭代函数可得：

$ (x)=(x+) $

选取任意初始迭代值xk，进行n次迭代得xk+n至F(xk+n)逼近零即可。当然这个公式也可以推广到开任意次方，只需要求出对应的迭代函数即可。

我起初尝试使用牛顿迭代法进行逼近，但是在使用C语言测试的时候，似乎很难满足题目中y = floor(sqrt(x))的要求：

牛顿迭代法需要17次迭代达到一个很好的精确度（当初始迭代点设置为4时），但题目要求该数开方后应该是65535，那我们就需要额外再迭代一次获得65535.9999的结果。此外，如果想满足题目要求，还必须要添加额外的逻辑进行向上取整，而且牛顿迭代中还涉及到除法运算，这将会使硬件电路变得更加复杂。最后，牛顿迭代法的迭代次数非常依赖于初始迭代值，因此，该电路的运算速度是不确定的。

二分法

二分法的核心思想是“对半分”，也就是不断的对分缩小猜数的区间，直到最后猜的数字满足精度要求。

具体来说，首先设置两个初始比较值a0与b0（其中a0<b0，两个比较值代表着区间对半缩小的两个方向，一个方向为左区间，另一个方向是右区间），接着进行k次迭代，每次迭代都需要求出an与·bn的平均值cn，比较cn^2与待开方数xn大小，若大于xn，则将cn赋给bn，反之则将cn赋给an。

以'd95为例，初始左操作数a='d0，右操作数b='d15：

  1. ('d0+'d15)>>1='d7,('d7)^2='d49<'d95,令a='d7;
  2. ('d7+'d15)>>1='d11,('d11)^2='d121>'d95,令b='d11;
  3. ('d7+'d11)>>1='d9,('d9)^2='d81<'d95,令a='d9;
  4. ('d9+'d11)>>1='d10,('d10)^2='d100>'d95,令b='d10,因为此时a='d9,b='d10，两者相差1'b1，因此无需再做下一次运算，输出a即结束。若中途出现a==b也即结束运算，输出a。

可以直观的感受到，此方法的迭代次数也非常依赖于初始值的选取，但相较于牛顿迭代法，其优点在于不需要实现除法部件（求平均值可以用移位代替）。

如果将a设置为0，b设置为输入的待开方数的一半的话，需要31次迭代才能计算出结果，性能不是很好。

逐次逼近法

相较于上面两种迭代方法，逐次逼近拥有确定的迭代次数，因为对于一个32位输入16位输出的模型来说，逐次逼近就是对这16位的输出从高位到低位依次比较，总共需要16次迭代，“逐次”也因此得名。

首先最高位MSB置为1，其余低位为0，将其平方后与din比较，如果大了表示这个最高位应该为0，如果小了表示这个最高位应该为1，接着判断次高位，重复置1、平方、比较、确定数值的过程，依次计算下一位直到最低位LSB结束得到结果。以下是算法过程：

  1. 从高到低将上一次置位的下一低位置1，该位以下的低位置零，若没有上一位则置最高位，若没有以下的低位则运算结束，转第二步；
  2. 将第一步得到的数平方，与原数比较，若比原数大则上一步里置1的那一位确定置0，若比原数小则上一步里置1的那一位确定置1，转第一步。

权衡了一下，我最终还是采用了逐次逼近的方法设计开方运算的硬件电路。

由于需要输出16位数据，所以需要循环重复16次比较操作，我们可以把它设计为16级流水线结构，这样能得到最高效的性能提升。

代码写得稍微冗杂了些，不过可读性还阔以，由于16级流水线的结构十分类似，所以说也可以采用下面参考文章中使用的generate块进行优化。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
module sqrt_32 (
    input   wire            clk,
    input   wire            rst_n,
    input   wire            vld_in,
    input   wire    [31:0]  x,
    output  wire            vld_out,
    output  wire    [31:0]  y
);

reg [15:0]  temp_y  [0:15] ;        //y的暂存值
reg         set_bit [0:15] ;        //用来判断当前位是否为1
reg [15:0]  state_cnt;              //状态计数器，用来计算vld_out
reg [31:0]  temp_x  [0:15] ;        //寄存x的值
//debug数个小时得出的结论：
//要想实现流水线，必须将输入数据x也是一级一级地通过寄存器传递下去，而不是直连到每一个流水级



always @(posedge clk or negedge rst_n) begin
    if(!rst_n)begin
        temp_y[15] <= 16'b1000_0000_0000_0000;
        set_bit[15] <= 1'b0;
    end
    else if(vld_in)begin
        set_bit[15] <= ((temp_y[15] * temp_y[15]) > x) ? 1'b0 : 1'b1;
    end
    else begin
        temp_y[15] <= temp_y[15];
        set_bit[15] <= set_bit[15];
    end

end

always @(posedge clk or negedge rst_n) begin
    if(!rst_n)begin
        temp_y[14] <= 16'b0100_0000_0000_0000;
        set_bit[14] <= 1'b0;
    end
    else if(vld_in)begin
        temp_y[14] <= {set_bit[15], 15'b100_0000_0000_0000};
        set_bit[14] <= (({set_bit[15], 15'b100_0000_0000_0000} * {set_bit[15], 15'b100_0000_0000_0000}) > temp_x[15]) ? 1'b0 : 1'b1;
    end
    else begin
        temp_y[14] <= temp_y[14];
        set_bit[14] <= set_bit[14];
    end
end

always @(posedge clk or negedge rst_n) begin
    if(!rst_n)begin
        temp_y[13] <= 16'b0010_0000_0000_0000;
        set_bit[13] <= 1'b0;
    end
    else if(vld_in)begin
        temp_y[13] <= {temp_y[14][15], set_bit[14], 14'b10_0000_0000_0000};
        set_bit[13] <= (({temp_y[14][15], set_bit[14], 14'b10_0000_0000_0000} * {temp_y[14][15], set_bit[14], 14'b10_0000_0000_0000}) > temp_x[14]) ? 1'b0 : 1'b1;
    end
    else begin
        temp_y[13] <= temp_y[13];
        set_bit[13] <= set_bit[13];
    end
end

always @(posedge clk or negedge rst_n) begin
    if(!rst_n)begin
        temp_y[12] <= 16'b0001_0000_0000_0000;
        set_bit[12] <= 1'b0;
    end
    else if(vld_in)begin
        temp_y[12] <= {temp_y[13][15:14], set_bit[13], 13'b1_0000_0000_0000};
        set_bit[12] <= (({temp_y[13][15:14], set_bit[13], 13'b1_0000_0000_0000} * {temp_y[13][15:14], set_bit[13], 13'b1_0000_0000_0000}) > temp_x[13]) ? 1'b0 : 1'b1;
    end
    else begin
        temp_y[12] <= temp_y[12];
        set_bit[12] <= set_bit[12];
    end
end

always @(posedge clk or negedge rst_n) begin
    if(!rst_n)begin
        temp_y[11] <= 16'b0000_1000_0000_0000;
        set_bit[11] <= 1'b0;
    end
    else if(vld_in)begin
        temp_y[11] <= {temp_y[12][15:13], set_bit[12], 12'b1000_0000_0000};
        set_bit[11] <= (({temp_y[12][15:13], set_bit[12], 12'b1000_0000_0000} * {temp_y[12][15:13], set_bit[12], 12'b1000_0000_0000}) > temp_x[12]) ? 1'b0 : 1'b1;
    end
    else begin
        temp_y[11] <= temp_y[11];
        set_bit[11] <= set_bit[11];
    end
end

always @(posedge clk or negedge rst_n) begin
    if(!rst_n)begin
        temp_y[10] <= 16'b0000_0100_0000_0000;
        set_bit[10] <= 1'b0;
    end
    else if(vld_in)begin
        temp_y[10] <= {temp_y[11][15:12], set_bit[11], 11'b100_0000_0000};
        set_bit[10] <= (({temp_y[11][15:12], set_bit[11], 11'b100_0000_0000} * {temp_y[11][15:12], set_bit[11], 11'b100_0000_0000}) > temp_x[11]) ? 1'b0 : 1'b1;
    end
    else begin
        temp_y[10] <= temp_y[10];
        set_bit[10] <= set_bit[10];
    end
end

always @(posedge clk or negedge rst_n) begin
    if(!rst_n)begin
        temp_y[9] <= 16'b0000_0010_0000_0000;
        set_bit[9] <= 1'b0;
    end
    else if(vld_in)begin
        temp_y[9] <= {temp_y[10][15:11], set_bit[10], 10'b10_0000_0000};
        set_bit[9] <= (({temp_y[10][15:11], set_bit[10], 10'b10_0000_0000} * {temp_y[10][15:11], set_bit[10], 10'b10_0000_0000}) > temp_x[10]) ? 1'b0 : 1'b1;
    end
    else begin
        temp_y[9] <= temp_y[9];
        set_bit[9] <= set_bit[9];
    end
end

always @(posedge clk or negedge rst_n) begin
    if(!rst_n)begin
        temp_y[8] <= 16'b0000_0001_0000_0000;
        set_bit[8] <= 1'b0;
    end
    else if(vld_in)begin
        temp_y[8] <= {temp_y[9][15:10], set_bit[9], 9'b1_0000_0000};
        set_bit[8] <= (({temp_y[9][15:10], set_bit[9], 9'b1_0000_0000} * {temp_y[9][15:10], set_bit[9], 9'b1_0000_0000}) > temp_x[9]) ? 1'b0 : 1'b1;
    end
    else begin
        temp_y[8] <= temp_y[8];
        set_bit[8] <= set_bit[8];
    end
end

always @(posedge clk or negedge rst_n) begin
    if(!rst_n)begin
        temp_y[7] <= 16'b0000_0000_1000_0000;
        set_bit[7] <= 1'b0;
    end
    else if(vld_in)begin
        temp_y[7] <= {temp_y[8][15:9], set_bit[8], 8'b1000_0000};
        set_bit[7] <= (({temp_y[8][15:9], set_bit[8], 8'b1000_0000} * {temp_y[8][15:9], set_bit[8], 8'b1000_0000}) > temp_x[8]) ? 1'b0 : 1'b1;
    end
    else begin
        temp_y[7] <= temp_y[7];
        set_bit[7] <= set_bit[7];
    end
end

always @(posedge clk or negedge rst_n) begin
    if(!rst_n)begin
        temp_y[6] <= 16'b0000_0000_0100_0000;
        set_bit[6] <= 1'b0;
    end
    else if(vld_in)begin
        temp_y[6] <= {temp_y[7][15:8], set_bit[7], 7'b100_0000};
        set_bit[6] <= (({temp_y[7][15:8], set_bit[7], 7'b100_0000} * {temp_y[7][15:8], set_bit[7], 7'b100_0000}) > temp_x[7]) ? 1'b0 : 1'b1;
    end
    else begin
        temp_y[6] <= temp_y[6];
        set_bit[6] <= set_bit[6];
    end
end

always @(posedge clk or negedge rst_n) begin
    if(!rst_n)begin
        temp_y[5] <= 16'b0000_0000_0010_0000;
        set_bit[5] <= 1'b0;
    end
    else if(vld_in)begin
        temp_y[5] <= {temp_y[6][15:7], set_bit[6], 6'b10_0000};
        set_bit[5] <= (({temp_y[6][15:7], set_bit[6], 6'b10_0000} * {temp_y[6][15:7], set_bit[6], 6'b10_0000}) > temp_x[6]) ? 1'b0 : 1'b1;
    end
    else begin
        temp_y[5] <= temp_y[5];
        set_bit[5] <= set_bit[5];
    end
end

always @(posedge clk or negedge rst_n) begin
    if(!rst_n)begin
        temp_y[4] <= 16'b0000_0000_0001_0000;
        set_bit[4] <= 1'b0;
    end
    else if(vld_in)begin
        temp_y[4] <= {temp_y[5][15:6], set_bit[5], 5'b1_0000};
        set_bit[4] <= (({temp_y[5][15:6], set_bit[5], 5'b1_0000} * {temp_y[5][15:6], set_bit[5], 5'b1_0000}) > temp_x[5]) ? 1'b0 : 1'b1;
    end
    else begin
        temp_y[4] <= temp_y[4];
        set_bit[4] <= set_bit[4];
    end
end

always @(posedge clk or negedge rst_n) begin
    if(!rst_n)begin
        temp_y[3] <= 16'b0000_0000_0000_1000;
        set_bit[3] <= 1'b0;
    end
    else if(vld_in)begin
        temp_y[3] <= {temp_y[4][15:5], set_bit[4], 4'b1000};
        set_bit[3] <= (({temp_y[4][15:5], set_bit[4], 4'b1000} * {temp_y[4][15:5], set_bit[4], 4'b1000}) > temp_x[4]) ? 1'b0 : 1'b1;
    end
    else begin
        temp_y[3] <= temp_y[3];
        set_bit[3] <= set_bit[3];
    end
end

always @(posedge clk or negedge rst_n) begin
    if(!rst_n)begin
        temp_y[2] <= 16'b0000_0000_0000_0100;
        set_bit[2] <= 1'b0;
    end
    else if(vld_in)begin
        temp_y[2] <= {temp_y[3][15:4], set_bit[3], 3'b100};
        set_bit[2] <= (({temp_y[3][15:4], set_bit[3], 3'b100} * {temp_y[3][15:4], set_bit[3], 3'b100}) > temp_x[3]) ? 1'b0 : 1'b1;
    end
    else begin
        temp_y[2] <= temp_y[2];
        set_bit[2] <= set_bit[2];
    end
end

always @(posedge clk or negedge rst_n) begin
    if(!rst_n)begin
        temp_y[1] <= 16'b0000_0000_0000_0010;
        set_bit[1] <= 1'b0;
    end
    else if(vld_in)begin
        temp_y[1] <= {temp_y[2][15:3], set_bit[2], 2'b10};
        set_bit[1] <= (({temp_y[2][15:3], set_bit[2], 2'b10} * {temp_y[2][15:3], set_bit[2], 2'b10}) > temp_x[2]) ? 1'b0 : 1'b1;
    end
    else begin
        temp_y[1] <= temp_y[1];
        set_bit[1] <= set_bit[1];
    end
end

always @(posedge clk or negedge rst_n) begin
    if(!rst_n)begin
        temp_y[0] <= 16'b0000_0000_0000_0001;
        set_bit[0] <= 1'b0;
    end
    else if(vld_in)begin
        temp_y[0] <= {temp_y[1][15:2], set_bit[1], 1'b1};
        set_bit[0] <= (({temp_y[1][15:2], set_bit[1], 1'b1} * {temp_y[1][15:2], set_bit[1], 1'b1}) > temp_x[1]) ? 1'b0 : 1'b1;
    end
    else begin
        temp_y[0] <= temp_y[0];
        set_bit[0] <= set_bit[0];
    end
end


always @(posedge clk or negedge rst_n) begin
    if(!rst_n) begin
        state_cnt <= 16'd0;
        temp_x[15] <= 32'b0;
        temp_x[14] <= 32'b0;
        temp_x[13] <= 32'b0;
        temp_x[12] <= 32'b0;
        temp_x[11] <= 32'b0;
        temp_x[10] <= 32'b0;
        temp_x[9] <= 32'b0;
        temp_x[8] <= 32'b0;
        temp_x[7] <= 32'b0;
        temp_x[6] <= 32'b0;
        temp_x[5] <= 32'b0;
        temp_x[4] <= 32'b0;
        temp_x[3] <= 32'b0;
        temp_x[2] <= 32'b0;
        temp_x[1] <= 32'b0;
        temp_x[0] <= 32'b0;
    end
    else begin
        state_cnt[15] <= vld_in;
        state_cnt[14] <= state_cnt[15];
        state_cnt[13] <= state_cnt[14];
        state_cnt[12] <= state_cnt[13];
        state_cnt[11] <= state_cnt[12];
        state_cnt[10] <= state_cnt[11];
        state_cnt[9] <= state_cnt[10];
        state_cnt[8] <= state_cnt[9];
        state_cnt[7] <= state_cnt[8];
        state_cnt[6] <= state_cnt[7];
        state_cnt[5] <= state_cnt[6];
        state_cnt[4] <= state_cnt[5];
        state_cnt[3] <= state_cnt[4];
        state_cnt[2] <= state_cnt[3];
        state_cnt[1] <= state_cnt[2];
        state_cnt[0] <= state_cnt[1];
        temp_x[15] <= x;
        temp_x[14] <= temp_x[15];
        temp_x[13] <= temp_x[14];
        temp_x[12] <= temp_x[13];
        temp_x[11] <= temp_x[12];
        temp_x[10] <= temp_x[11];
        temp_x[9] <= temp_x[10];
        temp_x[8] <= temp_x[9];
        temp_x[7] <= temp_x[8];
        temp_x[6] <= temp_x[7];
        temp_x[5] <= temp_x[6];
        temp_x[4] <= temp_x[5];
        temp_x[3] <= temp_x[4];
        temp_x[2] <= temp_x[3];
        temp_x[1] <= temp_x[2];
        temp_x[0] <= temp_x[1];
    end
end

assign y = {temp_y[0][15:1], set_bit[0]};
assign vld_out = state_cnt[0];

endmodule

最后仿真结果也确实能够实现数据的流水输出，并且延迟是16个时钟周期。

vld_out仅仅是将vld_in打了16拍，也符合逻辑。

对于vld_in信号无效时，则不会产生对应的数据输出，能够节省功耗。

打开Schematic原理图查看，确实如我们所预料的那样综合生成出了16级流水线结构：

在zynq7020的芯片（xc7z020clg400-1）上消耗的硬件资源如下：

至于功耗的话，由于我一如既往的没有分配管脚，所以它也一如既往的功耗爆表了，所以功耗分析就当看个笑话了。

三种常见平方根算法的电路设计及Verilog实现与仿真_verilog开根号-CSDN博客

Newton牛顿法（一）| 基本思想+迭代公式_idl 函数newton 的用法-CSDN博客

verilog 整数开方算法实现（逐次逼近法）_verilog开根号运算-CSDN博客

中国科学院大学数字集成电路作业开源——时序逻辑与存储器章节 - sasasatori - 博客园

学习笔记 > 数字IC
#硬件算法 #sqrt
平方根计算的硬件设计
http://example.com/2024/12/20/平方根计算的硬件设计/
作者
叶逸昇
发布于
2024年12月20日
许可协议