一、平方根及三種常見平方根算法簡介

數(shù)學是物理的基礎(chǔ)，是廣大世界的基本組成部分，而數(shù)學運算是數(shù)學理論的核心部分，數(shù)學運算有加減乘除乘方等基本運算，拓展的運算里有一項是開方運算，開方運算在數(shù)字計算、圖形顯示等領(lǐng)域具有重要的地位，所以如何在硬件上實現(xiàn)該運算可以提高計算單元的性能，加快計算速度。

本文實現(xiàn)的算法包括二分迭代法、牛頓迭代法、逐次逼近法，前兩種方法來源于數(shù)值計算方法，第三種方法類似于逐次漸進型ADC的原理，以下分別介紹這三種算法。本篇文章約定被開方數(shù)為16位無符號數(shù)，輸出開方結(jié)果為8位無符號數(shù)，采用多時鐘周期計算結(jié)果。

（一）、二分迭代法

二分法本質(zhì)上是一種區(qū)間迭代算法，通過不斷縮小隔根區(qū)間長度使區(qū)間中點不斷逼近根值。對于求一個連續(xù)函數(shù)f(x)在[a,b]區(qū)間上等于0的根，首先判斷是否f(a)*f(b)<0，若小于0則包含根，若大于0那么判斷是否單調(diào)，若單調(diào)則無根若不單調(diào)則含多根。以下簡介它的算法過程：

第一，計算y1=f(a)，y2=f(b);

第二，計算x0=0.5*(a+b)，y0=f(x0)，若|y0|<ε0，則輸出x0結(jié)束否則轉(zhuǎn)第三步；

第三，若y0*y1<0,則置b=x0，y2=y0，否則置a=x0，y1=y0，轉(zhuǎn)第四步；

第四，若|b-a|>ε1，則轉(zhuǎn)第二步，否則輸出x0結(jié)束。

（二）、牛頓迭代法

牛頓迭代法是一種局部收斂的算法，它的作法是在方程f(x)=0的根的鄰近點x0用f(x)的泰勒展開式的前兩項代替f(x)，得f(x0)+f'(x0)(x-x0)=0，如果f'(x0)!=0，可得到根的近似值x1=x0-f(x0)/f'(x0)。重復以上過程得到迭代公式。以下是算法過程：

第一，輸入根的初試近似值x0以及允許誤差ε，置n=0；

第二，計算；

第三，判斷，若|xn+1-xn|>=ε，則置n=n+1，轉(zhuǎn)第二步，否則輸出xn+1結(jié)束。

（三）、逐次逼近法

在平時的生活中我們會用到天平來稱物體的重量，那么稱重的過程是怎么樣的呢？

首先我們先放一個最大重量的砝碼，如果太重了表明這個砝碼應(yīng)該不用加，如果太輕了表明這個砝碼應(yīng)該加著，接著我們放一個第二重量的砝碼，重復以上的判斷過程，接著第三個，第四個...直到最輕的砝碼判斷完畢或者天平達到平衡結(jié)束稱重。

逐次逼近法也是如此，對于一個定寬的數(shù)據(jù)din，假設(shè)為16bits，開方后得到的數(shù)result應(yīng)該是8bits的，那么我們可以對8bits的數(shù)進行稱重過程的放置判斷操作，首先最高位MSB置為1，其余低位為0判斷平方后與din比較，如果大了表示這個最高位應(yīng)該為0，如果小了表示這個最高位應(yīng)該為1，接著判斷次高位，重復置1、平方、比較、確定數(shù)值的過程，依次計算下一位直到最低位LSB結(jié)束得到結(jié)果。以下是算法過程：

第一，從高到低將上一次置位的下一低位置1，該位以下的低位置零，若沒有上一位則置最高位，若沒有以下的低位則運算結(jié)束，轉(zhuǎn)第二步；

第二，將第一步得到的數(shù)平方，與原數(shù)比較，若比原數(shù)大則上一步里置1的那一位確定置0，若比原數(shù)小則上一步里置1的那一位確定置1，轉(zhuǎn)第一步。

二、Verilog狀態(tài)機的描述

狀態(tài)機來源于數(shù)字電路里的時序邏輯電路，設(shè)計狀態(tài)機的方法就是設(shè)計數(shù)字電路時序邏輯電路的方法，包括狀態(tài)編碼，狀態(tài)轉(zhuǎn)換圖的繪制，狀態(tài)轉(zhuǎn)換表的建立，狀態(tài)方程、輸出方程、驅(qū)動方程的書寫，競爭冒險的檢查和自啟動的檢查。

狀態(tài)機有摩爾Moore型和米利Mealy型，其中Moore型輸出僅依賴于狀態(tài)，Mealy型輸出與狀態(tài)和輸入有關(guān)。采用Verilog描述的狀態(tài)機有兩段式和三段式，兩段式性能最好但時序差，三段式時序性能兼顧。

兩段式描述分為狀態(tài)時序段state timing和狀態(tài)跳轉(zhuǎn)段state jump。狀態(tài)時序段采用時序邏輯always@(posedge clk)非阻塞賦值語句描述現(xiàn)態(tài)cur_state到次態(tài)nxt_state的轉(zhuǎn)換，狀態(tài)跳轉(zhuǎn)段采用組合邏輯always@(*)阻塞賦值語句配合case、if-else根據(jù)現(xiàn)態(tài)描述次態(tài)和輸出的邏輯。

三段式描述分為狀態(tài)時序段和狀態(tài)跳轉(zhuǎn)段和輸出信號段。狀態(tài)時序段和狀態(tài)跳轉(zhuǎn)段與二段式描述一致，只是將輸出信號的輸出邏輯的描述單獨拿出來在輸出信號段采用時序邏輯always@(posedge clk)根據(jù)次態(tài)nxt_state生成輸出信號。

三、算法的Verilog實現(xiàn)

在使用Verilog描述電路前必須做到心中有電路，這是采用HDL設(shè)計數(shù)字電路的前提。數(shù)字電路根據(jù)功能大體上可以分為數(shù)據(jù)通路和控制通路，至于先設(shè)計哪一部分取決于總體電路的功能是偏向數(shù)據(jù)處理運算還是偏向控制。根據(jù)以上的說明將對以下三種算法進行電路設(shè)計并用Verilog描述。

（一）、二分迭代法

由于在無符號二進制數(shù)運算里沒有乘積小于零的判斷結(jié)果，因此每次求出平均值后作平方與被開方數(shù)比較，若小于等于被開方數(shù)則將平均值賦給左操作數(shù)，若大于等于平均值則將平均值賦給右操作數(shù)。

以'd95為例，初始左操作數(shù)a='d0，右操作數(shù)b='d15。

第一次，('d0+'d15)>>1='d7,('d7)^2='d49<'d95,令a='d7;

第二次，('d7+'d15)>>1='d11,('d11)^2='d121>'d95,令b='d11;

第三次，('d7+'d11)>>1='d9,('d9)^2='d81<'d95,令a='d9;

第四次，('d9+'d11)>>1='d10,('d10)^2='d100>'d95,令b='d10,因為此時a='d9,b='d10，兩者相差1'b1，因此無需再做下一次運算，輸出a即結(jié)束。若中途出現(xiàn)a==b也即結(jié)束運算，輸出a。

首先分析運算過程考慮器件復用，決定采用時序電路多周期計算?？梢詫⒖刂仆返臓顟B(tài)劃分為三個狀態(tài)：IDLE，CALC，DONE。IDLE表示空閑，只有輸入一個使能信號calc_en才能啟動計算，即進入CALC狀態(tài)，這個狀態(tài)主要用于輸出數(shù)據(jù)通路使用的觸發(fā)器flip-flop的使能端，用以存儲計算中產(chǎn)生的信號，待計算達到一定的程度時輸入信號calc_end有效后狀態(tài)進入DONE，即完成狀態(tài)，此時輸出done信號表示計算結(jié)束。為了比較各個算法實現(xiàn)的電路的性能，在CALC狀態(tài)還應(yīng)該輸出一個計算器的計數(shù)使能，用于對計算所用時鐘周期的計數(shù)。通過以上分析可以得到以下的狀態(tài)轉(zhuǎn)換圖和輸出信號表。

以下是數(shù)據(jù)通路電路圖

通過result乘方與din比較決定是否刷新左右操作數(shù)的選擇信號selLeft和selRight。

result的輸出邏輯

那么問題來了，怎么判斷計算結(jié)束了呢？我們通過上文二分法的例子發(fā)現(xiàn)計算完成時左操作數(shù)與右操作數(shù)不是相等就是差1，于是可以有以下的邏輯輸出calc_end，再輸入給狀態(tài)機使狀態(tài)跳轉(zhuǎn)。

經(jīng)過以上分析已經(jīng)可以用Verilog描述電路了，模塊名為mysqrt1，文件名一致。

/******************************************************************************
* mysqrt1.v
*******************************************************************************/

module mysqrt1(
  input        clk,
  input        calc_en,
  input        rst_n,
  input [15:0] din,
  output [7:0] result_o,
  output [3:0] cnt_clk,
  output       done_o
);

encode state
  localparam IDLE = 2'b00;
  localparam CALC = 2'b01;
  localparam DONE = 2'b10;

middle signals
  reg  [1:0]  cur_state,nxt_state;//state
  reg         ff_en;//enable flip-flop
  reg         cnt_en;//enable counter
  reg         done;//finish
  reg  [8:0]  opLeft,opRight;//for operation "opLeft"+"opRight"
  wire [7:0]  result;//store result
  wire [8:0]  adder_tmp;//temp adder output data
  wire [8:0]  opLeft_tmp;//temp opLeft data
  wire [8:0]  opRight_tmp;//temp opRight data
  wire        opOr1,opOr2;//gate Or inputs
  wire [15:0] multi_tmp;//temp multi out data
  wire        calc_end;//end calculate
  wire        co;//from counter
  wire        selLeft,selRight;//select store which to opLeft and opRight


assign output signals
  assign result_o = result;
  assign done_o   = done;

state timing
  always@(posedge clk,negedge rst_n) begin
    if(!rst_n)
  cur_state <= IDLE;
    else
  cur_state <= nxt_state;
  end

state jump->nxt_state
  always@(*) begin
    case(cur_state)
  IDLE:nxt_state = calc_en  ? CALC : IDLE;
  CALC:nxt_state = calc_end ? DONE : CALC;
  DONE:nxt_state = IDLE;
  default:nxt_state = IDLE;
    endcase
  end

control signals logic to data path
  always@(posedge clk,negedge rst_n) begin
    if(!rst_n) begin
  ff_en  <= 1'b0;
  cnt_en <= 1'b0;
  done   <= 1'b0;
    end
    else
  case(nxt_state)
IDLE:begin
  ff_en  <= 1'b0;
  cnt_en <= 1'b0;
  done   <= 1'b0;
end
CALC:begin
  ff_en  <= 1'b1;
  cnt_en <= 1'b1;
  done   <= 1'b0;
end
DONE:begin
  ff_en  <= 1'b0;
  cnt_en <= 1'b0;
  done   <= 1'b1;
end
default:begin
  ff_en  <= 1'b0;
  cnt_en <= 1'b0;
  done   <= 1'b0;
end
  endcase
  end

data path
//counter
  cnt_ceil cnt_ceil_x(
    .clk   (clk),
    .en    (cnt_en),
    .rst_n (rst_n),
    .ceil  (4'b1111),
    .cnt   (cnt_clk),
    .co    (co)
  );
//"selLeft" and "selRight" logic
  assign multi_tmp = result * result;
  assign selLeft   = (multi_tmp <= din);
  assign selRight  = (multi_tmp >= din);
//"calc_end" logic
  assign opOr1    = ((1'b1+opLeft)==opRight);
  assign opOr2    = (opLeft==opRight);
  assign calc_end = opOr1 || opOr2;
//"result" logic
  assign opLeft_tmp  = selLeft  ? {1'b0,result} : opLeft;
  assign opRight_tmp = selRight ? {1'b0,result} : opRight;
  assign adder_tmp   = opLeft + opRight;
  assign result      = adder_tmp[8:1];
  always@(posedge clk,negedge rst_n) begin
    if(!rst_n) begin
      opLeft  <= 9'b0_0000_0000;//set left to minimal
      opRight <= 9'b0_1111_1111;//set right to maximal
    end
    else
  if(ff_en) begin
        opLeft  <= opLeft_tmp;
        opRight <= opRight_tmp;
      end
  end
  
endmodule

計數(shù)器模塊cnt_ceil代碼如下。

/******************************************************************************
* cnt_ceil.v
*******************************************************************************/

module cnt_ceil(
  input        clk,
  input        en,
  input        rst_n,
  input  [3:0] ceil,
  output [3:0] cnt,
  output       co
);
  reg [3:0] cnt_temp;
  always @(posedge clk,negedge rst_n) begin
if(!rst_n)
  cnt_temp <= 4'b0000;
else 
  if(en)
    if(cnt_temp>=ceil)
          cnt_temp <= 4'b0000;
        else
  cnt_temp <= cnt_temp+1'b1;
  end
  assign cnt = cnt_temp;
  assign co  = en && (cnt_temp==ceil);
endmodule

（二）、牛頓迭代法

電路分為控制通路和數(shù)據(jù)通路，控制通路與二分法一致，不再贅述。以下分析數(shù)據(jù)通路。

計算的核心是公式x=(a/x+x)/2，使用除法器和加法器可以構(gòu)成計算核心。如下圖所示。

問題是怎么判斷計算結(jié)束了呢？

舉個例子，被開方數(shù)a='d5343，初始迭代數(shù)x='d255。

第一次，(5343/255+255)/2=137；第二次，(5343/137+137)/2=88；

第三次，(5343/88+88)/2=74；第四次，(5343/74+74)/2=73；第五次，(5343/73+73)/2=73

可以發(fā)現(xiàn)經(jīng)過迭代后最后中間數(shù)穩(wěn)定不變即可判斷計算結(jié)束，還發(fā)現(xiàn)最終的結(jié)果與上一次迭代的結(jié)果僅差1，那么也可由此判斷計算已經(jīng)結(jié)束無需作下一次迭代。于是可以畫出以下的電路，輸出calc_end，再輸入給狀態(tài)機使狀態(tài)跳轉(zhuǎn)。

經(jīng)過以上分析，可以用Verilog描述電路，定義模塊名為mysqrt2，文件名一致。

/******************************************************************************
* mysqrt2.v
*******************************************************************************/

module mysqrt2(
  input        clk,
  input        calc_en,
  input        rst_n,
  input [15:0] din,
  output [7:0] result_o,
  output [3:0] cnt_clk,
  output       done_o
);

encode state
  localparam IDLE = 2'b00;
  localparam CALC = 2'b01;
  localparam DONE = 2'b11;

middle signals
  reg  [1:0] cur_state,nxt_state;//state
  reg  [7:0] result;//result
  reg        done;//finish
  reg        ff_en;//enable flip-flop store
  reg        cnt_en;//enable counter
  wire       calc_end;//end calculate
  wire [7:0] div_tmp;//temp divide data
  wire [8:0] opAdder1,opAdder2;//to adder
  wire [8:0] adder_tmp;//output from adder
  wire [7:0] r_tmp;//temp result
  wire       opOr1,opOr2;//to gate Or
  wire       co;//counter output
  

assign output signals
  assign result_o = result;
  assign done_o   = done;

state timing
  always@(posedge clk,negedge rst_n) begin
    if(!rst_n)
  cur_state <= IDLE;
    else
  cur_state <= nxt_state;
  end

state jump->nxt_state
  always@(*) begin
    case(cur_state)
  IDLE:nxt_state = calc_en  ? CALC : IDLE;
  CALC:nxt_state = calc_end ? DONE : CALC;
  DONE:nxt_state = IDLE;
  default:nxt_state = IDLE;
    endcase
  end

control signals logic to data path
  always@(posedge clk,negedge rst_n) begin
    if(!rst_n) begin
  ff_en  <= 1'b0;
  cnt_en <= 1'b0;
  done   <= 1'b0;
    end
    else
  case(nxt_state)
        IDLE:begin
      ff_en  <= 1'b0;
      cnt_en <= 1'b0;
      done   <= 1'b0;
    end
        CALC:begin
      ff_en  <= 1'b1;
      cnt_en <= 1'b1;
      done   <= 1'b0;
    end
        DONE:begin
      ff_en  <= 1'b0;
      cnt_en <= 1'b0;
      done   <= 1'b1;
    end
        default:begin
      ff_en  <= 1'b0;
      cnt_en <= 1'b0;
      done   <= 1'b0;
    end
  endcase
  end

data path
//counter
  cnt_ceil cnt_ceil_1(
    .clk   (clk),
    .en    (cnt_en),
    .rst_n (rst_n),
    .ceil  (4'b1111),
    .cnt   (cnt_clk),
    .co    (co)
  );
//"calc_end" logic
  assign opOr1    = ((r_tmp+1'b1)==result);
  assign opOr2    = (r_tmp==result);
  assign calc_end = opOr1 || opOr2;
//"result" logic
  assign div_tmp   = din / result;
  assign opAdder1  = {1'b0,div_tmp};
  assign opAdder2  = {1'b0,result};
  assign adder_tmp = opAdder1 + opAdder2;
  assign r_tmp     = adder_tmp[8:1];
  always@(posedge clk,negedge rst_n) begin
    if(!rst_n)
  result <= 8'b1111_1111;//set to maximal---'d255
    else
  if(ff_en)
result <= r_tmp;
  end
endmodule

（三）、逐次逼近法

控制通路的狀態(tài)機與二分法一致，不再贅述。以下分析數(shù)據(jù)通路。

數(shù)據(jù)通路的關(guān)鍵是如何依次對result每一位判斷是否為1，這里引入中間信號index[2:0]，用來記錄當前應(yīng)該對哪一位數(shù)操作，這里的index可以為計數(shù)器輸出的低三位，再由index和乘法器比較器輸出中間信號judge，用來判斷該位是否為1，由index和常數(shù)進行比較產(chǎn)生selx信號，作為MUX的選擇信號，最后用judge和selx驅(qū)動MUX輸出result的每一位。

定義模塊名為mysqrt3，Verilog文件名一致。

/******************************************************************************
*mysqrt3.v
*******************************************************************************/
module mysqrt3(
  input         clk,
  input         calc_en,
  input         rst_n,
  input  [15:0] din,
  output [7:0]  result_o,
  output        done_o
);
//
//encode states
  localparam IDLE = 2'b00;
  localparam CALC = 2'b01;
  localparam DONE = 2'b10;
/
//middle signals
  reg  [1:0] cur_state,nxt_state;//state
  reg        done;//done
  reg  [7:0] result;//store result
  reg        cnt_en;//enable counter
  wire [3:0] cnt;//counter output
  wire [2:0] index;//calculated index from 'd0 to 'd7
  wire       judge;//judge if result[index] is 1'b1
  wire       co;//counter output co
  wire       sel0;//if index=='d7
  wire       sel1;//if index=='d6
  wire       sel2;//if index=='d5
  wire       sel3;//if index=='d4
  wire       sel4;//if index=='d3
  wire       sel5;//if index=='d2
  wire       sel6;//if index=='d1
  wire       sel7;//if index=='d0
  reg        ff_en;//enable store result
  wire       calc_end;//calculate end signal
  wire       r0_tmp,r1_tmp,r2_tmp,r3_tmp,r4_tmp,r5_tmp,r6_tmp,r7_tmp;//temp data
  wire       j_tmp;//temp data
  reg  [7:0] op_tmp;//temp data
  wire [15:0] multi_tmp;//temp data
/
//assign output signals
  assign result_o = result;
  assign done_o   = done;
/
//state timing
  always @(posedge clk,negedge rst_n) begin
    if(!rst_n)
  cur_state <= IDLE;
    else
  cur_state <= nxt_state;
  end
/
//state jump
  always @(*) begin
case(cur_state)
  IDLE:nxt_state = calc_en ? CALC : IDLE;
  CALC:nxt_state = calc_end ? DONE : CALC;
  DONE:nxt_state = IDLE;
  default:nxt_state = IDLE;
endcase
  end
/
//control signals logic to data path  
  always @(posedge clk,negedge rst_n) begin
    if(!rst_n) begin
  ff_en  <= 1'b0;
  cnt_en <= 1'b0;
  done   <= 1'b0;
    end
    else
  case(nxt_state)
IDLE: begin
  ff_en  <= 1'b0;
          cnt_en <= 1'b0;
  done   <= 1'b0;
end
CALC: begin
  ff_en  <= 1'b1;
  cnt_en <= 1'b1;
  done   <= 1'b0;
end
DONE: begin
  ff_en  <= 1'b0;
  cnt_en <= 1'b0;
  done   <= 1'b1;
end
default: begin
  ff_en  <= 1'b0;
      cnt_en <= 1'b0;
  done   <= 1'b0;
end

  endcase
  end
/
data path
//"index" logic
  cnt_ceil cnt_ceil_0(
    .clk   (clk),
    .en    (cnt_en),
    .rst_n (rst_n),
    .ceil  (4'd7),
    .cnt   (cnt),
    .co    (co)
  );
  assign index = cnt[2:0];
//"judge" logic
  always @(*) begin
    case(index)
  3'd0:op_tmp = 8'b1000_0000;
  3'd1:op_tmp = {result[7],7'b100_0000};
  3'd2:op_tmp = {result[7:6],6'b10_0000};
  3'd3:op_tmp = {result[7:5],5'b1_0000};
  3'd4:op_tmp = {result[7:4],4'b1000};
  3'd5:op_tmp = {result[7:3],3'b100};
  3'd6:op_tmp = {result[7:2],2'b10};
  3'd7:op_tmp = {result[7:1],1'b1};
    endcase
  end
  assign multi_tmp = op_tmp * op_tmp;
  assign judge     = (multi_tmp <= din);
//"selx" logic
  assign sel7 = (index==3'd0);
  assign sel6 = (index==3'd1);
  assign sel5 = (index==3'd2);
  assign sel4 = (index==3'd3);
  assign sel3 = (index==3'd4);
  assign sel2 = (index==3'd5);
  assign sel1 = (index==3'd6);
  assign sel0 = (index==3'd7);
//"result" logic
  assign j_tmp  = judge ? 1'b1  : 1'b0;
  assign r7_tmp = sel7  ? j_tmp : result[7];
  assign r6_tmp = sel6  ? j_tmp : result[6];
  assign r5_tmp = sel5  ? j_tmp : result[5];
  assign r4_tmp = sel4  ? j_tmp : result[4];
  assign r3_tmp = sel3  ? j_tmp : result[3];
  assign r2_tmp = sel2  ? j_tmp : result[2];
  assign r1_tmp = sel1  ? j_tmp : result[1];
  assign r0_tmp = sel0  ? j_tmp : result[0];
  always @(posedge clk,negedge rst_n) begin
    if(!rst_n) begin
  result <= 8'b0000_0000;
    end
    else
  if(ff_en) begin
    result[7] <= r7_tmp;
    result[6] <= r6_tmp;
    result[5] <= r5_tmp;
    result[4] <= r4_tmp;
    result[3] <= r3_tmp;
    result[2] <= r2_tmp;
    result[1] <= r1_tmp;
    result[0] <= r0_tmp;
  end
  end
//"calc_end" logic
  assign calc_end = co;

endmodule

四、進行仿真

（一）統(tǒng)一的testbench

用Verilog編寫一個統(tǒng)一的testbench，在其中分別例化三個算法實現(xiàn)模塊。

result1對應(yīng)mysqrt1的結(jié)果，result2對應(yīng)mysqrt2的結(jié)果，result3對應(yīng)mysqrt3的結(jié)果；

done1對應(yīng)mysqrt1的完成，done2對應(yīng)mysqrt2的完成，done3對應(yīng)mysqrt3的完成；

rst_n1對應(yīng)mysqrt1的異步重置，rst_n2對應(yīng)mysqrt2的異步重置，rst_n3對應(yīng)mysqrt3的異步重置；

在每個模塊的每次計算完畢后使用rst_nx異步重置內(nèi)部寄存器數(shù)據(jù)并輸入新的din進行下一次運算。

定義模塊名為mysqrt_tb，文件名一致，程序如下。

//testbenchf for mysqrt1 and mysqrt2 and mysqrt3
//mysqrt_tb.v
`timescale 1ns/1ns
module mysqrt_tb();
  reg        clk;
  reg        ensqrt1,ensqrt2,ensqrt3;
  reg        rst_n1,rst_n2,rst_n3;
  reg [15:0] din1,din2,din3;
  wire [7:0] result1,result2,result3;
  wire       done1,done2,done3;
  wire [3:0] cnt_clk1;
  wire [3:0] cnt_clk2;
//initialize
  initial begin
    clk     = 1'b0;
    ensqrt1 = 1'b1;
    ensqrt2 = 1'b1;
    ensqrt3 = 1'b1;
    rst_n1  = 1'b1;
    rst_n2  = 1'b1;
    rst_n3  = 1'b1;
    din1    = 'd95;
    din2    = 'd95;
    din3    = 'd95;
  end
  initial begin//clk
    forever #1 clk = ~clk;
  end
  initial begin//the first rst_n
    #1 rst_n1 = 1'b0; rst_n2 = 1'b0; rst_n3 = 1'b0;
    #1 rst_n1 = 1'b1; rst_n2 = 1'b1; rst_n3 = 1'b1;
  end
//din1 and rst_n1
  always@(posedge clk,posedge done1) begin
    if(done1) begin//when done1 is pulled high
      #2 rst_n1 <= 1'b0;
         din1   <= {$random}%16'b1111_1111_1111_1111;
      #1 rst_n1 <= 1'b1;
    end
  end
//din2 and rst_n2
  always@(posedge clk,posedge done2) begin
    if(done2) begin//when done2 is pulled high
      #2 rst_n2 <= 1'b0;
         din2   <= {$random}%16'b1111_1111_1111_1111;
      #1 rst_n2 <= 1'b1;
    end
  end
//din3 and rst_n3
  always@(posedge clk,posedge done3) begin
    if(done3) begin//when done3 is pulled high
      #2 rst_n3 <= 1'b0;
         din3   <= {$random}%16'b1111_1111_1111_1111;
      #1 rst_n3 <= 1'b1;
    end
  end
//instances
  mysqrt1 mysqrt1_0(
    .clk      (clk),
    .calc_en  (ensqrt1),
    .rst_n    (rst_n1),
    .din      (din1),
    .result_o (result1),
    .cnt_clk  (cnt_clk1),
    .done_o   (done1)
  );
  mysqrt2 mysqrt2_0(
    .clk      (clk),
    .calc_en  (ensqrt2),
    .rst_n    (rst_n2),
    .din      (din2),
    .result_o (result2),
    .cnt_clk  (cnt_clk2),
    .done_o   (done2)
  );
  mysqrt3 mysqrt3_0(
    .clk      (clk),
    .calc_en  (ensqrt3),
    .rst_n    (rst_n3),
    .din      (din3),
    .result_o (result3),
    .done_o   (done3)
  );

endmodule

（二）、波形結(jié)果

附上使用Modelsim仿真的波形結(jié)果

附上使用Verilator和GTKWave的仿真結(jié)果

使用Verilator仿真基于Verilog的testbench可以參考我寫的另一篇博客：使用Verilator仿真基于Verilog的testbench并用GTKWave查看波形。

五、分析與反思

二分法

計算性能取決于起始區(qū)間的位置，由于設(shè)計中沒有設(shè)定讀入起始區(qū)間的信號，而且被開方數(shù)遍布于整個16位空間，只能將區(qū)間設(shè)為最大，即左為0，右為’b1111_1111=’d255，這就使得每次計算都需要8個周期。那么為什么是8周期呢？首先尋找一個4位數(shù)用最大區(qū)間需要找4次，這好比在一個邊長為4的正方形里找一個邊長為1的正方形x，每次劃分后剩下區(qū)域都為先前的一半，經(jīng)過4次迭代才找到x，所以找8位數(shù)需要8周期，再經(jīng)過一周期的拉高done信號的周期，總共9周期。有一個問題是數(shù)據(jù)通路中左右操作數(shù)經(jīng)過加法器和乘法器的串聯(lián)再經(jīng)過MUX回到觸發(fā)器D端，中間的延時太長了，如果在加法和乘法中間加一級寄存器則會減小延時，但是這會導致計算周期翻倍為16周期，所以這是時序和性能的權(quán)衡，這個問題和權(quán)衡在牛頓法中（除法器和加法器串聯(lián)）也是如此。

牛頓法

性能最好，但是用了一個除法器。設(shè)計中為了滿足存儲中間數(shù)據(jù)的定寬9位數(shù)不溢出，迭代初始值設(shè)為’b1111_1111=’d255，在較大的數(shù)輸入時能夠較快地迭代出結(jié)果（2-4周期），但是遇到較小的數(shù)比如’d95就需要迭代7周期。因為輸出一定是介于0到255，所以設(shè)想如果將初始值設(shè)為’d127是否能對所有數(shù)據(jù)迭代周期平均化為4周期，數(shù)學上是可行的，但是當前的電路不能滿足要求，因為中間數(shù)據(jù)(a/x+x)/2可能會超出9位造成結(jié)果不收斂無法拉高calc_end信號，這里計算一下在初始值為127時最大的中間數(shù)據(jù)，(65534/127+127)=643=’b10_1000_0011，有十位數(shù)，除2后為9位數(shù)，那么存儲除法結(jié)果應(yīng)該用10位數(shù)，計算加法應(yīng)該用10位數(shù)，而存儲加法除二后的結(jié)果應(yīng)該用9位數(shù)，這樣才不會數(shù)據(jù)溢出。其實testbench里隨機的數(shù)據(jù)最大為65534，如果為65535(‘b1111_1111_1111_1111)，計算的中間數(shù)據(jù)其實是’b10_0000_0000，這是十位數(shù)，在當前設(shè)計中也是會數(shù)據(jù)溢出的，如果不改變設(shè)計會導致錯誤，因為下一次計算會得到除數(shù)為0的情況。除了修改設(shè)計外的解決辦法是額外增加邏輯判斷輸入的數(shù)是否為65535，是則直接輸出結(jié)果為255結(jié)束，否則按初始值為255迭代計算。

逐次逼近法

性能最穩(wěn)定，計算周期為8周期，最后完成狀態(tài)輸出加1周期，使用了一個乘法器和較多的比較器，result每一位的觸發(fā)器使能在計算狀態(tài)均處于有效狀態(tài)增加了動態(tài)功耗，文中的selx信號是由比較器組輸出的，與寫操作的每一位的觸發(fā)器一一對應(yīng)，于是思考可以由selx直接作為觸發(fā)器的使能端而不需由控制通路輸出，那么設(shè)計可以改成由index作為輸入的三-八譯碼器輸出獨熱碼作為selx組合同時也作為觸發(fā)器的使能端組合，使得僅對當前操作位的觸發(fā)器寫入而對其余觸發(fā)器均寫無效以減少動態(tài)功耗和面積。還有一個可以改進的地方是當中間數(shù)據(jù)平方后等于輸入其實已經(jīng)可以結(jié)束計算了，但是本文的設(shè)計中沒有該判斷邏輯，所以無論怎樣都要計算8周期，對于很多數(shù)這是多余的，有興趣的讀者可以自己設(shè)計加上該邏輯。