A Verilog Function for finding SQUARE ROOT

   

UPDATE: A CLOCKED SQUARE ROOT CALCULATOR  IN VERILOG IS AVAILABLE HERE.

    In Verilog, there are no built-in operator to find the square root of a number. There are algorithms for doing this, but you have to write the code from scratch.

    Here, I want to share a Verilog function for finding the square root of a binary number. The function is based on "Non-Restoring Square Root algorithm". You can learn more about the algorithm from this paper. The function takes a 32 bit input number and returns a 16 bit square root. The block diagram of the algorithm is given below:

     Here D is the input number. R is the remainder of the operation for non-perfect squares. Q contains the square root of 'D'.

The Verilog function along with the testbench code is given below:

module testbench;

reg [15:0] sqr;

//Verilog function to find square root of a 32 bit number.
//The output is 16 bit.
function [15:0] sqrt;
    input [31:0] num;  //declare input
    //intermediate signals.
    reg [31:0] a;
    reg [15:0] q;
    reg [17:0] left,right,r;    
    integer i;
begin
    //initialize all the variables.
    a = num;
    q = 0;
    i = 0;
    left = 0;   //input to adder/sub
    right = 0;  //input to adder/sub
    r = 0;  //remainder
    //run the calculations for 16 iterations.
    for(i=0;i<16;i=i+1) begin 
        right = {q,r[17],1'b1};
        left = {r[15:0],a[31:30]};
        a = {a[29:0],2'b00};    //left shift by 2 bits.
        if (r[17] == 1) //add if r is negative
            r = left + right;
        else    //subtract if r is positive
            r = left - right;
        q = {q[14:0],!r[17]};       
    end
    sqrt = q;   //final assignment of output.
end
endfunction //end of Function

//simulation-Apply inputs.
    initial begin
        sqr = sqrt(32'd4000000);    #100;
        sqr = sqrt(32'd96100);  #100;
        sqr = sqrt(32'd25); #100;
        sqr = sqrt(32'd100000000);  #100;
        sqr = sqrt(32'd33); #100;
        sqr = sqrt(32'd3300);   #100;
        sqr = sqrt(32'd330000); #100;
        sqr = sqrt(32'd3300000000); #100;
    end
      
endmodule

Simulation waveform:

The code was synthesised and simulated using Xilinx ISE 14.6.

From the waveform, you can see that, for perfect squares we get the exact square root as result. But for other inputs, there is an error of +/- 1. 

Another way to increase the precision is adding multiples of two zeros to the right of input.

Notice the last 4 inputs in the testbench.

Actual Square root of 33 = 5.7445, But we get 5.

Further for input, 3300 we get 57 as output. Adding 2 zeros we get one fractional digit in the output.

Further for input, 330000 we get 574 as output. Adding 4 zeros we get two fractional digits in the output.

Further for input, 3300000000 we get 57445 as output. Adding 8 zeros we get four fractional digits in the output.

Synthesisable Verilog code for Division of two binary numbers

For doing division, Verilog has an operator, '/' defined. But this operator has some limitation when it comes to certain synthesis tools such as Xilinx XST. The '/' operator is synthesisable only when the second operand is a power of 2. 

For division of generic numbers Xilinx returns the following error message during synthesis:

Can not simplify operator DIV

 Few year back, I wrote a VHDL function for implementing division which was synthesisable. The design was based on Restoring Division algorithm. In this article, I have converted the same design into Verilog.

Division code:

The size of operands to the division module are defined through a parameter named WIDTH. This way you can use the same code for implementing 8 or 16 or 32 or any sized division. Note that, both the input operands and output have the same size.

module division(A,B,Res);

    //the size of input and output ports of the division module is generic.
    parameter WIDTH = 8;
    //input and output ports.
    input [WIDTH-1:0] A;
    input [WIDTH-1:0] B;
    output [WIDTH-1:0] Res;
    //internal variables    
    reg [WIDTH-1:0] Res = 0;
    reg [WIDTH-1:0] a1,b1;
    reg [WIDTH:0] p1;   
    integer i;

    always@ (A or B)
    begin
        //initialize the variables.
        a1 = A;
        b1 = B;
        p1= 0;
        for(i=0;i < WIDTH;i=i+1)    begin //start the for loop
            p1 = {p1[WIDTH-2:0],a1[WIDTH-1]};
            a1[WIDTH-1:1] = a1[WIDTH-2:0];
            p1 = p1-b1;
            if(p1[WIDTH-1] == 1)    begin
                a1[0] = 0;
                p1 = p1 + b1;   end
            else
                a1[0] = 1;
        end
        Res = a1;   
    end 

endmodule

Testbench code:

module tb_division;
    parameter WIDTH = 8;
    // Inputs
    reg [WIDTH-1:0] A;
    reg [WIDTH-1:0] B;
    // Outputs
    wire [WIDTH-1:0] Res;

    // Instantiate the division module (UUT)
    division #(WIDTH) uut (
        .A(A), 
        .B(B), 
        .Res(Res)
    );

    initial begin
        // Initialize Inputs and wait for 100 ns
        A = 0;  B = 0;  #100;  //Undefined inputs
        //Apply each set of inputs and wait for 100 ns.
        A = 100;    B = 10; #100;
        A = 200;    B = 40; #100;
        A = 90; B = 9;  #100;
        A = 70; B = 10; #100;
        A = 16; B = 3;  #100;
        A = 255;    B = 5;  #100;
    end
      
endmodule

Simulation Waveform:

The codes were simulated in Xilinx ISE 13.1 to get the following waveform.

Synthesis Results:

The design was synthesised for Virtex 4 fpga and a maximum combinational path delay of 20 ns was obtained for 8 bit division.