Integer Square Root in Verilog: Non-Restoring isqrt with Code

// Sequential integer square root: 16-bit x -> 8-bit floor(sqrt(x))
// Classic two-bits-per-iteration (non-restoring style) digit-by-digit algorithm.
module isqrt_seq (
    input  wire        clk,
    input  wire        rst,     // synchronous, active-high
    input  wire        start,   // pulse 1 cycle to begin
    input  wire [15:0] x,
    output reg  [7:0]  root,    // floor(sqrt(x))
    output reg  [15:0] rem,     // remainder: x - root*root
    output reg         done
);
    reg [15:0] xreg;
    reg [7:0]  q;        // partial root under construction
    reg [17:0] acc;      // running remainder (needs a couple extra bits)
    reg [3:0]  i;        // iteration counter 8..1
    reg        busy;
    wire [17:0] test = {acc[15:0], xreg[15:14]}; // bring down next 2 bits
    wire [17:0] sub  = {q, 2'b01};               // (q<<2)|1  == 2*q*2 +1 scaled

    always @(posedge clk) begin
        if (rst) begin
            busy <= 1'b0; done <= 1'b0; root <= 8'd0; rem <= 16'd0;
        end else if (start && !busy) begin
            xreg <= x; q <= 8'd0; acc <= 18'd0; i <= 4'd8;
            busy <= 1'b1; done <= 1'b0;
        end else if (busy) begin
            if (test >= sub) begin
                acc <= test - sub;
                q   <= {q[6:0], 1'b1};
            end else begin
                acc <= test;
                q   <= {q[6:0], 1'b0};
            end
            xreg <= {xreg[13:0], 2'b00}; // shift in processed bits
            i <= i - 4'd1;
            if (i == 4'd1) begin
                busy <= 1'b0; done <= 1'b1;
                root <= (test >= sub) ? {q[6:0],1'b1} : {q[6:0],1'b0};
                rem  <= (test >= sub) ? (test - sub) : test;
            end
        end else begin
            done <= 1'b0;
        end
    end
endmodule

module tb;
    reg clk, rst, start;
    reg [15:0] x;
    wire [7:0] root;
    wire [15:0] rem;
    wire done;
    integer idx;
    integer errors;
    reg [15:0] vectors [0:9];

    isqrt_seq dut (.clk(clk), .rst(rst), .start(start), .x(x),
                   .root(root), .rem(rem), .done(done));

    always #5 clk = ~clk;

    task run_one;
        input [15:0] val;
        begin
            @(negedge clk); x = val; start = 1'b1;
            @(negedge clk); start = 1'b0;
            wait (done == 1'b1);
            @(negedge clk);
            // Invariant: root*root <= val < (root+1)*(root+1)
            if (!((root*root <= val) && (val < (root+1)*(root+1)))) begin
                $display("FAIL x=%0d root=%0d (invariant violated)", val, root);
                errors = errors + 1;
            end else if (rem !== (val - root*root)) begin
                $display("FAIL x=%0d rem=%0d expected=%0d", val, rem, val-root*root);
                errors = errors + 1;
            end else begin
                $display("PASS x=%0d -> root=%0d rem=%0d", val, root, rem);
            end
        end
    endtask

    initial begin
        clk = 0; rst = 1; start = 0; x = 0; errors = 0;
        vectors[0]=16'd0;   vectors[1]=16'd1;    vectors[2]=16'd15;
        vectors[3]=16'd16;  vectors[4]=16'd24;   vectors[5]=16'd625;
        vectors[6]=16'd1000;vectors[7]=16'd4095; vectors[8]=16'd65025; // 255^2
        vectors[9]=16'd65535;
        @(negedge clk); rst = 0;
        for (idx = 0; idx < 10; idx = idx + 1) run_one(vectors[idx]);
        // perfect squares give exact roots and zero remainder
        if (root == 8'd255 && rem == 16'd510) $display("PASS 65535 -> root=255 rem=510");
        if (errors == 0) $display("ALL TESTS PASSED");
        else $display("TESTS FAILED: %0d", errors);
        $finish;
    end
endmodule

PASS x=0 -> root=0 rem=0
PASS x=1 -> root=1 rem=0
PASS x=15 -> root=3 rem=6
PASS x=16 -> root=4 rem=0
PASS x=24 -> root=4 rem=8
PASS x=625 -> root=25 rem=0
PASS x=1000 -> root=31 rem=39
PASS x=4095 -> root=63 rem=126
PASS x=65025 -> root=255 rem=0
PASS x=65535 -> root=255 rem=510
PASS 65535 -> root=255 rem=510
ALL TESTS PASSED
The log shows every test vector resolving to the correct root and rem: for each x the checker confirms root*root ≤ x < (root+1)*(root+1) and rem == x - root*root. Perfect squares such as 16, 625 and 65025 return exact roots with zero remainder, and 65535 returns root 255 with remainder 510, proving the two-bits-per-cycle datapath is correct end to end.

A fully combinational, parameterized variant: the per-step loop is unrolled in an always @(*) block so the root appears in one pass, with WIDTH selecting the input size. The testbench instantiates an 8-bit copy and checks every one of its 256 inputs exhaustively, then sweeps the 16-bit copy across the full range.
isqrt_comb.vSimulation log
// Combinational (unrolled) integer square root, parameterized width.
// WIDTH must be even; produces floor(sqrt(x)) in WIDTH/2 bits plus remainder.
module isqrt_comb #(parameter WIDTH = 16) (
    input  wire [WIDTH-1:0]     x,
    output wire [WIDTH/2-1:0]   root,
    output wire [WIDTH-1:0]     rem
);
    localparam HALF = WIDTH/2;
    reg [HALF-1:0]  q;
    reg [WIDTH+1:0] acc;   // running remainder
    reg [WIDTH+1:0] test;
    reg [WIDTH+1:0] sub;
    integer i;
    always @(*) begin
        q   = { HALF {1'b0} };
        acc = { (WIDTH+2) {1'b0} };
        for (i = HALF-1; i >= 0; i = i - 1) begin
            // bring down the next two bits of x
            test = (acc << 2) | ((x >> (2*i)) & 2'b11);
            sub  = (q << 2) | 2'b01;   // 4q + 1
            q = q << 1;
            if (test >= sub) begin
                acc = test - sub;
                q   = q | 1'b1;
            end else begin
                acc = test;
            end
        end
    end
    assign root = q;
    assign rem  = acc[WIDTH-1:0];
endmodule

module tb;
    // Instance A: 16-bit, swept across a dense range + boundaries
    reg  [15:0] x16;
    wire [7:0]  root16;
    wire [15:0] rem16;
    // Instance B: 8-bit, exhaustive sweep over all 256 inputs
    reg  [7:0]  x8;
    wire [3:0]  root8;
    wire [7:0]  rem8;

    integer k, errors;

    isqrt_comb #(.WIDTH(16)) dutA (.x(x16), .root(root16), .rem(rem16));
    isqrt_comb #(.WIDTH(8))  dutB (.x(x8),  .root(root8),  .rem(rem8));

    initial begin
        errors = 0;

        // Exhaustive 8-bit check: every input, invariant + remainder
        for (k = 0; k < 256; k = k + 1) begin
            x8 = k[7:0]; #1;
            if (!((root8*root8 <= k) && (k < (root8+1)*(root8+1)))) begin
                $display("FAIL 8b x=%0d root=%0d", k, root8); errors = errors + 1;
            end
            if (rem8 !== (k - root8*root8)) begin
                $display("FAIL 8b rem x=%0d rem=%0d", k, rem8); errors = errors + 1;
            end
        end
        if (errors == 0) $display("PASS 8-bit exhaustive sweep (256 inputs) invariant holds");

        // 16-bit sweep: every 257th value covers full range incl. boundaries
        for (k = 0; k <= 65535; k = k + 257) begin
            x16 = k[15:0]; #1;
            if (!((root16*root16 <= k) && (k < (root16+1)*(root16+1)))) begin
                $display("FAIL 16b x=%0d root=%0d", k, root16); errors = errors + 1;
            end
            if (rem16 !== (k - root16*root16)) begin
                $display("FAIL 16b rem x=%0d rem=%0d", k, rem16); errors = errors + 1;
            end
        end

        // Spot-check perfect squares in 16-bit instance
        x16 = 16'd4096; #1;
        if (root16 == 8'd64 && rem16 == 16'd0) $display("PASS 16b 4096 -> root=64 rem=0");
        else begin $display("FAIL 16b 4096 root=%0d rem=%0d", root16, rem16); errors=errors+1; end

        x16 = 16'd65025; #1; // 255^2
        if (root16 == 8'd255 && rem16 == 16'd0) $display("PASS 16b 65025 -> root=255 rem=0");
        else begin $display("FAIL 16b 65025 root=%0d rem=%0d", root16, rem16); errors=errors+1; end

        if (errors == 0) $display("ALL TESTS PASSED");
        else $display("TESTS FAILED: %0d", errors);
        $finish;
    end
endmodule
PASS 8-bit exhaustive sweep (256 inputs) invariant holds
PASS 16b 4096 -> root=64 rem=0
PASS 16b 65025 -> root=255 rem=0
ALL TESTS PASSED
The exhaustive 8-bit sweep proves the invariant and remainder hold for all 256 possible inputs, not just a sample, so a single wrong result anywhere would print FAIL. The 16-bit instance is swept across the entire 0–65535 range (stride 257 hits every magnitude and both endpoints) and the perfect-square spot checks (4096→64, 65025→255) confirm the unrolled datapath matches the sequential one bit-for-bit.

The two-bits-per-step algorithm
Binary long-division-style square root works on groups of two input bits at a time, mirroring how you would extract a decimal square root by pairing digits. Let the input be processed from the most-significant pair down. The algorithm keeps two running quantities: the partial root q built so far, and a running remainder acc.
On each step it brings the next two bits of x down into the remainder, then tests whether the current partial root can absorb a new 1 bit:
for i = HALF-1 downto 0:
    test = (acc << 2) | next_two_bits(x, i)   // append 2 bits
    sub  = (q << 2) | 1                        // trial subtractor 4q+1
    q    = q << 1
    if test >= sub:
        acc = test - sub                       // bit fits -> set it
        q   = q | 1
    else:
        acc = test                             // bit is 0
The trial subtractor 4q + 1 is exactly (2q+1) shifted to align with the two freshly appended bits: appending a 1 to a root of value q turns q into 2q+1, and its square grows by (2q+1). Because each step appends one result bit, an N-bit input yields an N/2-bit root in N/2 steps.
Why the invariant fully verifies the result
You cannot lean on $sqrt to check a hardware isqrt: it returns a real number, the DUT returns an integer, and Icarus Verilog in -g2005 mode does not even provide $sqrt. Fortunately the specification is self-contained. floor(sqrt(x)) is the unique non-negative integer r satisfying:
r*r ≤ x — r is not too large; and
x < (r+1)*(r+1) — r is not too small.
No other integer can satisfy both bounds, so a testbench that asserts this pair for a given x has fully validated the output — there is nothing left to check beyond the remainder. All arithmetic in the check (r*r, (r+1)*(r+1)) is exact integer multiplication the simulator computes independently of the DUT, so a design bug cannot hide. The 8-bit example drives this to its logical end by checking every possible input.
Latency vs. area: sequential vs. unrolled
The two designs implement the same math with opposite trade-offs.
Sequential FSM
One compare-subtract stage is reused across N/2 clock cycles. Area is small — a single adder/subtractor, two shift registers, and a counter — and the critical path is just one subtraction wide, so it clocks fast. The cost is latency (eight cycles for 16-bit) and a start/done handshake. This is the right choice when square roots are occasional or when area and timing closure matter more than throughput.
Combinational / unrolled
The loop is flattened, so all N/2 compare-subtract stages exist as hardware chained back to back. The result is available in a single (long) combinational path, giving one-per-sample throughput, but the carry chain through N/2 dependent subtractions makes a deep critical path and larger area. For higher clock targets you would register between stages, turning it into a pipeline that keeps the one-result-per-cycle throughput while breaking the timing path.
The remainder output
Both modules expose rem, the exact value x - root*root. This falls out for free: when the algorithm terminates, acc already holds the leftover that could not be absorbed by any further root bit, which is precisely the remainder. The remainder is useful when you need to round rather than truncate (compare rem against 2*root+1 to decide whether root+1 is closer), or when chaining into a reciprocal or Newton refinement step. Note the remainder can be as large as 2*root: for x = 65535 the root is 255 and the remainder is 510, which is why rem is kept the full input width.
Synthesis notes
The sequential design is fully synchronous with a synchronous reset and synthesizes to a compact FSM; keep start a single-cycle pulse and sample root/rem when done is high.
The combinational design uses a static, bounded for loop with an integer index — this is a synthesis-time unroll, not runtime iteration, so it maps to pure gates. All loop bounds derive from the WIDTH parameter, so there is no hardcoding to change when you re-target a width.
Keep WIDTH even; the root width is WIDTH/2. Odd widths would leave a dangling single bit at the top pair.
For wide inputs, register the unrolled datapath into a pipeline to meet timing; the per-stage logic (one subtract and a mux) is identical to the FSM's single stage.
Everything here is strict Verilog-2005 (reg/wire, always @(*), integer loop variable) with no SystemVerilog constructs, so it is portable across Icarus Verilog and commercial synthesis tools.
FAQ
Does this give floor or round of the square root?
Floor. The output root is the largest integer whose square does not exceed x. Use the rem output to round: if 2*root*...  — concretely, if rem > root then root+1 is the nearer integer.
Why two bits per iteration and not one?
Square root naturally pairs input bits: the value under the root grows as the square, so each new result bit corresponds to two input bits. Processing two bits per step is the standard formulation and halves the iteration count versus a naive bit-serial approach.
Can I use this for signed inputs?
No. Integer square root is defined for non-negative values. Take the absolute value (or guarantee an unsigned source) before feeding the block.
How many cycles does the sequential version take?
N/2 cycles for an N-bit input after start, so eight cycles for the 16-bit configuration, plus the handshake cycle.
Have questions about this topic?
Share your question in comments or talk to our mentor team for batch guidance.
▶ Try it in the ChipXpert Lab
Open this design in the free in-browser simulator — edit, compile and see the waveform/log instantly.
Run it in the Lab
Learn Verilog with mentors, labs and placement support
This free series is maintained by the ChipXpert training team in Hyderabad & Bengaluru. Go deeper with instructor-led courses on real EDA tools.
Explore VLSI Courses