// REFRACT.v — single-file version: REFRACT (EP1 top), ray_top and every module it uses.
// Merged from ray_top (was ray_top.v), ray_z0.v, z_surface.v, test.v, unit_normal_vector.v, normal_vector.v, refract_vector.v, sqrt_seq.v, div_seq.v.
// Those files still exist for the other testbenches; this file does not
// follow later edits to them.

// ==========================================================================
// REFRACT (EP1 top level)
// ==========================================================================
// EP1 interface: the host holds RST high with RI (2..15) set, then releases it;
// RI stays the same until the end. REFRACT computes the crossing with z = 0
// (zx, zy) of the refracted ray from every grid point X, Y = 0..15 and writes it
// into the host's SRAM, 512 words of unsigned Q4.12:
//   ray i = 16*Y + X: zx_i at 2i, zy_i at 2i + 1 (zx0 0, zy0 1, zx1 2, zy1 3, ...)
// The rays are written in order i = 0..255, so SRAM_A goes 0, 1, 2, ... 511.
// Each ray is two writes on back-to-back clocks with SRAM_WE high, zx first;
// points are ~148 clocks apart. DONE goes high the clock after the last write
// and stays high until RST. SRAM_Q is not used.
// n = 1/RI (Q1.16, rounded) comes from a table, so no divider is needed.
module REFRACT(
    input  wire        CLK,
    input  wire        RST,                     // active high, asynchronous
    input  wire [3:0]  RI,                      // refractive index 2..15, n = 1/RI
    output reg  [8:0]  SRAM_A,
    output reg  [15:0] SRAM_D,                  // q4.12
    input  wire [15:0] SRAM_Q,                  // not used
    output reg         SRAM_WE,
    output reg         DONE
);
    wire        rst_n = ~RST;                   // the modules below use an active-low reset
    // q1.16, round(65536 / RI); a function + assign rather than always @(*),
    // which a simulator need not run at time 0 when RI is already constant
    function [16:0] n_of(input [3:0] ri);
        case (ri)
            4'd2:    n_of = 17'd32768;
            4'd3:    n_of = 17'd21845;
            4'd4:    n_of = 17'd16384;
            4'd5:    n_of = 17'd13107;
            4'd6:    n_of = 17'd10923;
            4'd7:    n_of = 17'd9362;
            4'd8:    n_of = 17'd8192;
            4'd9:    n_of = 17'd7282;
            4'd10:   n_of = 17'd6554;
            4'd11:   n_of = 17'd5958;
            4'd12:   n_of = 17'd5461;
            4'd13:   n_of = 17'd5041;
            4'd14:   n_of = 17'd4681;
            4'd15:   n_of = 17'd4369;
            default: n_of = 17'd65536;          // RI = 0, 1: n = 1, no refraction
        endcase
    endfunction
    wire [16:0] n = n_of(RI);

    reg         start, started;                 // one start pulse after reset
    wire [3:0]  out_x, out_y;
    wire [15:0] X, Y;
    wire        valid, done, busy;
    reg  [15:0] zy;                             // Y, written the clock after zx
    reg         zy_next;                        // write zy on this clock
    reg         last;                           // the point being written is (15, 15)

    ray_top ray_top_inst (.CLOCK2_50(CLK), .reset(rst_n), .start(start), .n(n),
        .out_x(out_x), .out_y(out_y), .X(X), .Y(Y), .valid(valid), .done(done), .busy(busy));

    always @(posedge CLK or posedge RST) begin
        if (RST) begin
            start   <= 1'b0;
            started <= 1'b0;
            SRAM_A  <= 9'd0;
            SRAM_D  <= 16'd0;
            SRAM_WE <= 1'b0;
            DONE    <= 1'b0;
            zy      <= 16'd0;
            zy_next <= 1'b0;
            last    <= 1'b0;
        end else begin
            start   <= !started;
            started <= 1'b1;
            SRAM_WE <= 1'b0;
            zy_next <= 1'b0;
            if (valid) begin                    // zx; done comes with the last valid
                SRAM_A  <= {out_y, out_x, 1'b0};  // 2i, i = 16*Y + X
                SRAM_D  <= X;
                SRAM_WE <= 1'b1;
                zy      <= Y;
                zy_next <= 1'b1;
                last    <= done;
            end else if (zy_next) begin         // zy
                SRAM_A  <= {SRAM_A[8:1], 1'b1};
                SRAM_D  <= zy;
                SRAM_WE <= 1'b1;
            end
            if (SRAM_WE && SRAM_A[0] && last)   // the SRAM takes the last zy on this edge
                DONE <= 1'b1;
        end
    end

endmodule

// ==========================================================================
// ray_top (was ray_top.v)
// ==========================================================================
// For every grid point x, y = 0..15: the line through P = (x, y, z) on the
// surface z = 6 - 2((x-8)/8)^8 - 2((y-8)/8)^8, in the direction of the
// refracted ray T = (Tx, Ty, Tz), and where it crosses z = 0:
//   z_surface : x, y -> z                       (q4.24)
//   test      : x, y, n -> unit normal N -> T   (q8.84)
//   ray_z0    : x, y, z, T -> X, Y at z = 0     (q4.12)
//
// A start pulse sweeps y = 0..15 (outer), x = 0..15 (inner). For each point:
//   1. pulse test's start: unit normal N (sqrt 44 + divide 36 clocks), then
//      T from N (sqrt 28 clocks); wait for T_done
//   2. pulse ray_z0's start, wait for its done (29 clocks)
// then out_x, out_y, X, Y are valid for the one clock valid is high. done
// pulses after the last point (15, 15). A start while busy is ignored.
// A point takes ~148 clocks, the sweep 37893 (RI does not change it).
module ray_top (
    input  wire        CLOCK2_50,
    input  wire        reset,                   // active low
    input  wire        start,
    input  wire [16:0] n,                       // q1.16, n = n1/n2, constant during a sweep
    output reg  [3:0]  out_x,                   // q4, the point X, Y belong to
    output reg  [3:0]  out_y,                   // q4
    output reg  [15:0] X,                       // q4.12, crossing with z = 0
    output reg  [15:0] Y,                       // q4.12
    output reg         valid,                   // one clock per point
    output reg         done,                    // one clock after the last point
    output wire        busy
);
    reg  [3:0] x, y;                            // point being computed
    wire [27:0] z;
    wire signed [91:0] Tx, Ty, Tz;
    wire n_done, T_done;
    reg  n_start, ray_start;
    wire [15:0] ray_X, ray_Y;
    wire ray_done;

    z_surface z_surface_inst (.x(x), .y(y), .z(z));

    test test_inst (.CLOCK2_50(CLOCK2_50), .reset(reset), .x(x), .y(y), .z(4'd0), .n(n),
        .start(n_start), .Tx(Tx), .Ty(Ty), .Tz(Tz), .n_done(n_done), .T_done(T_done));

    ray_z0 ray_z0_inst (.CLOCK2_50(CLOCK2_50), .reset(reset), .start(ray_start),
        .x(x), .y(y), .z(z), .Tx(Tx), .Ty(Ty), .Tz(Tz), .X(ray_X), .Y(ray_Y), .done(ray_done));

    localparam IDLE = 3'd0, START_T = 3'd1, WAIT_T = 3'd2, START = 3'd3, WAIT_RAY = 3'd4;
    reg [2:0] state;
    assign busy = state != IDLE;

    always @(posedge CLOCK2_50 or negedge reset) begin
        if (!reset) begin
            state     <= IDLE;
            x         <= 4'd0;
            y         <= 4'd0;
            n_start   <= 1'b0;
            ray_start <= 1'b0;
            out_x     <= 4'd0;
            out_y     <= 4'd0;
            X         <= 16'd0;
            Y         <= 16'd0;
            valid     <= 1'b0;
            done      <= 1'b0;
        end else begin
            n_start   <= 1'b0;
            ray_start <= 1'b0;
            valid     <= 1'b0;
            done      <= 1'b0;
            case (state)
                IDLE:
                    if (start) begin
                        x     <= 4'd0;
                        y     <= 4'd0;
                        state <= START_T;
                    end
                START_T: begin                  // test latches x, y on the next edge
                    n_start <= 1'b1;
                    state   <= WAIT_T;
                end
                WAIT_T:                         // T for this x, y
                    if (T_done) state <= START;
                START: begin                    // ray_z0 latches x, y, z, T on the next edge
                    ray_start <= 1'b1;
                    state     <= WAIT_RAY;
                end
                WAIT_RAY:
                    if (ray_done) begin
                        out_x <= x;
                        out_y <= y;
                        X     <= ray_X;
                        Y     <= ray_Y;
                        valid <= 1'b1;
                        if (x == 4'd15 && y == 4'd15) begin
                            done  <= 1'b1;
                            state <= IDLE;
                        end else begin
                            if (x == 4'd15) y <= y + 4'd1;
                            x     <= x + 4'd1;
                            state <= START_T;
                        end
                    end
                default: state <= IDLE;
            endcase
        end
    end

endmodule

// ==========================================================================
// ray_z0.v
// ==========================================================================
// Where the line through P = (x, y, z) with direction T = (Tx, Ty, Tz) crosses
// the plane z = 0:
//   P + t*T, t = -z/Tz  ->  X = x - z*Tx/Tz,  Y = y - z*Ty/Tz
// x, y are the grid point (Q4), z its height from z_surface (Q4.24), T the
// refracted ray from refract_vector (Q8.84). X, Y are unsigned Q4.12, rounded
// to nearest and clamped to 0..15.9998.
//
// z*|Tx| / |Tz| and z*|Ty| / |Tz| are done by two div_seq in parallel, 29
// clocks. x, y, z, T are latched on a start pulse (ignored while busy); done
// goes high for one clock with X, Y. Tz must not be 0 (a ray parallel to the
// plane never crosses it); |Tz| >= 1/3 for this surface.
module ray_z0 (
    input  wire               CLOCK2_50,
    input  wire               reset,          // active low
    input  wire               start,
    input  wire        [3:0]  x,              // q4
    input  wire        [3:0]  y,              // q4
    input  wire        [27:0] z,              // q4.24, 2..6 on this surface
    input  wire signed [91:0] Tx,             // q8.84, |T| = 1
    input  wire signed [91:0] Ty,             // q8.84
    input  wire signed [91:0] Tz,             // q8.84
    output reg         [15:0] X,              // q4.12
    output reg         [15:0] Y,              // q4.12
    output reg                done
);
    // |T| <= 1, so q8.84 bits [85:54] are T as q2.30 (the bits above are sign copies)
    wire signed [31:0] tx = Tx[85:54];
    wire signed [31:0] ty = Ty[85:54];
    wire signed [31:0] tz = Tz[85:54];
    wire        [31:0] ax = tx[31] ? -tx : tx;                  // |Tx| q2.30
    wire        [31:0] ay = ty[31] ? -ty : ty;
    wire        [31:0] az = tz[31] ? -tz : tz;

    // z*|T| is q4.24 * q2.30 -> q.54; / |Tz| (q2.30) -> quotient q5.24, < 32
    wire [59:0] num_x = z * ax;
    wire [59:0] num_y = z * ay;
    wire [28:0] qx, qy;                                         // z*|Tx|/|Tz|, z*|Ty|/|Tz|
    wire        div_done;

    div_seq #(.WN(60), .WD(32), .WQ(29)) div_x (.clk(CLOCK2_50), .rst(reset), .start(start),
        .dividend(num_x), .divisor(az), .quotient(qx), .done(div_done));
    div_seq #(.WN(60), .WD(32), .WQ(29)) div_y (.clk(CLOCK2_50), .rst(reset), .start(start),
        .dividend(num_y), .divisor(az), .quotient(qy), .done());

    // x, y and the signs of Tx/Tz, Ty/Tz that the running divide belongs to
    reg  [3:0] x_l, y_l;
    reg        sx, sy;                                          // 1: Tx/Tz < 0, so X = x + q
    reg        busy;
    always @(posedge CLOCK2_50 or negedge reset) begin
        if (!reset) begin
            busy <= 1'b0;
            x_l  <= 4'd0;
            y_l  <= 4'd0;
            sx   <= 1'b0;
            sy   <= 1'b0;
        end else if (start && (!busy || div_done)) begin        // div_seq latches on this edge too
            busy <= 1'b1;
            x_l  <= x;
            y_l  <= y;
            sx   <= tx[31] ^ tz[31];
            sy   <= ty[31] ^ tz[31];
        end else if (div_done)
            busy <= 1'b0;
    end

    // X = x -/+ q in q5.24 (signed), then round to q4.12 and clamp to 0..0xFFFF
    wire signed [31:0] xs = $signed({4'd0, x_l, 24'd0});
    wire signed [31:0] ys = $signed({4'd0, y_l, 24'd0});
    wire signed [31:0] x24 = sx ? xs + $signed({3'd0, qx}) : xs - $signed({3'd0, qx});
    wire signed [31:0] y24 = sy ? ys + $signed({3'd0, qy}) : ys - $signed({3'd0, qy});
    wire signed [31:0] x12 = (x24 + 32'sd2048) >>> 12;
    wire signed [31:0] y12 = (y24 + 32'sd2048) >>> 12;

    always @(posedge CLOCK2_50 or negedge reset) begin
        if (!reset) begin
            X    <= 16'd0;
            Y    <= 16'd0;
            done <= 1'b0;
        end else begin
            done <= div_done;
            if (div_done) begin
                X <= x12 < 0 ? 16'd0 : (x12 > 32'sd65535 ? 16'hFFFF : x12[15:0]);
                Y <= y12 < 0 ? 16'd0 : (y12 > 32'sd65535 ? 16'hFFFF : y12[15:0]);
            end
        end
    end

endmodule

// ==========================================================================
// z_surface.v
// ==========================================================================
//Z=6−2((x-8)/8)^8-2((y-8)/8)^8
module z_surface (
    input wire [3:0] x, //q4
    input wire [3:0] y, //q4
    output reg [27:0] z //q4.24
);
    reg signed [39:0] temp, dx, dy;
    always @(*) begin
        // ** with a non-constant base can't be synthesized, so multiply out (x-8)^8
        dx = {36'b0,x} - 40'sd8;
        dy = {36'b0,y} - 40'sd8;
        temp = (40'sd6 << 24)-((dx*dx*dx*dx*dx*dx*dx*dx)<<1) - ((dy*dy*dy*dy*dy*dy*dy*dy)<<1);
        z=temp[27:0];
    end

endmodule

// ==========================================================================
// test.v
// ==========================================================================
module test (
    input wire CLOCK2_50,
    input wire reset, //active low
    input wire [3:0] x, //q4
    input wire [3:0] y, //q4
    input wire [3:0] z,  //q4
    input wire [16:0] n, //q1.16, n = 1/RI
    input wire start,    //latch x, y and compute N, then T

    output wire signed [91:0] Tx, //q8.84
    output wire signed [91:0] Ty, //q8.84
    output wire signed [91:0] Tz, //q8.84
    output wire n_done,           //unit normal done: nx, ny, nz update with it
    output wire T_done            //one clock, Tx, Ty, Tz valid from this clock on
);
    wire signed [1:0] Ix = 2'sd0;
    wire signed [1:0] Iy = 2'sd0;
    wire signed [1:0] Iz = -2'sd1;

    wire [35:0] nx,ny,nz;
    wire done0;
    assign n_done = done0;
    unit_normal_vector unit_normal_vector_inst(
        .CLOCK2_50(CLOCK2_50),
        .x(x), //q4
        .y(y), //q4
        // .z(z), //q4
        .reset(reset), //active low
        .start(start),
        .nx(nx), //q2.34
        .ny(ny), //q2.34
        .nz(nz),  //q2.34
        .done(done0)
    );

    refract_vector refract_vector_inst(
        .done0(done0),
        .Ix(Ix), //q2
        .Iy(Iy), //q2
        .Iz(Iz), //q2
        .Nx(nx), //q2.34
        .Ny(ny), //q2.34
        .Nz(nz), //q2.34
        .n(n), //q1.16
        .Tx(Tx), //q8.84
        .Ty(Ty), //q8.84
        .Tz(Tz),  //q8.84
        .done(T_done),
        .reset(reset), //active low
        .CLOCK2_50(CLOCK2_50)
    );
    
endmodule

// ==========================================================================
// unit_normal_vector.v
// ==========================================================================
//N= (gx,gy,gz)/sqrt(gx^2+gy^2+gz^2)
module unit_normal_vector (
    input  wire        CLOCK2_50,
    input  wire [3:0] x, //q4
    input  wire [3:0] y, //q4
    input wire reset, //active low
    input wire start, //latch x, y and begin: sqrt W/2 = 44 clocks, then divide 36 clocks
    output reg [35:0] nx, //q2.34
    output reg [35:0] ny, //q2.34
    output reg [35:0] nz,  //q2.34
    output wire done
);
    wire signed [35:0] x1,y1,z1;
    wire [71:0] sum;
    normal_vector normal_vector_inst (
        .x(x),      //q4
        .y(y),      //q4
        .x1(x1),    //q2.34
        .y1(y1),    //q2.34
        .z1(z1),    //q2.34
        .sum(sum)   //q4.68
    );

    wire [43:0] root;
    wire sqrt_done;
    sqrt_seq #(.W(88)) sqrt_seq_inst (
        .clk    (CLOCK2_50),
        .rst    (reset),
        .start  (start),
        .radicand({sum,16'd0}), //q4.84
        .root   (root), //q2.42
        .done   (sqrt_done)
    );
    reg [35:0] x1_l, y1_l, z1_l;           // the x1, y1, z1 that the running sqrt belongs to

    // Q2.34 << 42 divided by Q2.42 root -> Q2.34; adding half the divisor
    // to |numerator| rounds instead of truncating. A combinational 88-bit divide
    // can't be synthesized (lpm_divide is limited to 64 bits), so div_seq does
    // it in 36 clocks, started on sqrt_done
    wire [87:0] half  = {45'd0, root[43:1]};
    wire [87:0] num_x = {10'd0, (x1_l[35] ? -x1_l : x1_l), 42'd0} + half;   // |x1_l| << 42 + half
    wire [87:0] num_y = {10'd0, (y1_l[35] ? -y1_l : y1_l), 42'd0} + half;
    wire [87:0] num_z = {10'd0, (z1_l[35] ? -z1_l : z1_l), 42'd0} + half;
    wire [35:0] qx, qy, qz;                // |nx|, |ny|, |nz|
    reg         sx, sy, sz;                // signs of the running divide

    div_seq #(.WN(88), .WD(44), .WQ(36)) div_x (.clk(CLOCK2_50), .rst(reset), .start(sqrt_done),
        .dividend(num_x), .divisor(root), .quotient(qx), .done(done));
    div_seq #(.WN(88), .WD(44), .WQ(36)) div_y (.clk(CLOCK2_50), .rst(reset), .start(sqrt_done),
        .dividend(num_y), .divisor(root), .quotient(qy), .done());
    div_seq #(.WN(88), .WD(44), .WQ(36)) div_z (.clk(CLOCK2_50), .rst(reset), .start(sqrt_done),
        .dividend(num_z), .divisor(root), .quotient(qz), .done());

    always @(posedge CLOCK2_50) begin
        if (start) begin
            x1_l <= x1;                     // sqrt latches sum on this same edge
            y1_l <= y1;
            z1_l <= z1;
        end
        if (sqrt_done) begin
            sx <= x1_l[35];                 // the divides latch num_x.. on this same edge
            sy <= y1_l[35];
            sz <= z1_l[35];
        end
        if (done) begin                     // done = divide done: nx, ny, nz update with it
            nx <= sx ? -qx : qx;
            ny <= sy ? -qy : qy;
            nz <= sz ? -qz : qz;
        end
    end


    endmodule

// ==========================================================================
// normal_vector.v
// ==========================================================================
//gx=2*(x-8)^7/8^7 gy=2*(y-8)^7/8^7 gz=1

module normal_vector (
    input wire [3:0] x, //q4
    input wire [3:0] y, //q4
    output reg [35:0] x1, //q2.34
    output reg [35:0] y1, //q2.34
    output reg [35:0] z1, //q2.34
    output reg [71:0] sum  //q4.68, x1^2+y1^2+z1^2 <= 9
);
    reg signed [39:0] dx, dy;   // x-8, y-8
    always @(*) begin
        // reg signed [39:0] temp;
        // reg [27:0] x1,y1,z1;
        // reg [63:0] len;
        // reg [63:0] sum; //q8.48
        // 2(x-8)^7/8^7 = (x-8)^7 / 2^20 -> Q2.34 is (x-8)^7 << 14 (-2.0 at x=0 fits)
        z1 = 40'sd1 << 34;
        // ** with a non-constant base can't be synthesized, so multiply out (x-8)^7
        dx = {36'b0,x} - 40'sd8;
        dy = {36'b0,y} - 40'sd8;
        x1 = (dx*dx*dx*dx*dx*dx*dx)<<14;
        y1 = (dy*dy*dy*dy*dy*dy*dy)<<14;
        sum = (($signed(x1)*$signed(x1) + $signed(y1)*$signed(y1) + $signed(z1)*$signed(z1)));
        // sqrt_seq #(.W(32)) sqrt_seq_inst (
        //     .clk    (1'b0),
        //     .rst    (1'b0),
        //     .start  (1'b1),
        //     .radicand(sum[31:0]),
        //     .root   (len),
        //     .done   ()
        // );
        // nx = ((64'd0+x1) << 24)/len;
        // ny = ((64'd0+y1) << 24)/len;
        // nz = ((64'd0+z1) << 24)/len;
    end
endmodule

// ==========================================================================
// refract_vector.v
// ==========================================================================
//T=n*I+(n*(-N*I)-(1-n^2*(1-(N*I)^2))^0.5)*N
//I is unit input vector, N is unit normal vector
//n is refractive index ratio 1/1,1/2~1/16
//T is unit output vector
module refract_vector(
    input wire CLOCK2_50,
    input wire done0,
    input wire signed[1:0] Ix, //q2
    input wire signed[1:0] Iy, //q2
    input wire signed[1:0] Iz, //q2
    input wire signed[35:0] Nx, //q2.34
    input wire signed[35:0] Ny, //q2.34
    input wire signed[35:0] Nz, //q2.34
    input wire [16:0] n,   //q1.16
    input wire reset, //active low
    output wire signed [91:0] Tx, //q8.84
    output wire signed [91:0] Ty, //q8.84
    output wire signed [91:0] Tz, //q8.84
    output wire done              //one clock, Tx, Ty, Tz valid from this clock on
);
    // done0 pulses on the edge that updates N, so the sqrt starts one clock
    // later, from the new N; T is latched when it is done (W/2 = 28 clocks)
    reg sqrt_start, done_r;
    wire [27:0] sqrt_one_m_n_sq_times_one_m_NI_sq; //q3.25
    wire done1;                   //sqrt done
    always @(posedge CLOCK2_50 or negedge reset) begin
        if (!reset) begin
            sqrt_start <= 1'b0;
            done_r     <= 1'b0;
        end else begin
            sqrt_start <= done0;
            done_r     <= done1;
        end
    end

    reg signed [37:0] NI; //q4.34
    reg signed [41:0] one_m_NI_sq; //q8.34
    reg signed [55:0] one_m_n_sq_times_one_m_NI_sq; //q6.50
    reg signed [55:0] n_sq_times_one_m_NI_sq; //q6.50
    
    always@(*) begin
        //N*I
        NI = $signed(Nx)*$signed(Ix) + $signed(Ny)*$signed(Iy) + $signed(Nz)*$signed(Iz); //q4.34
        //1-(N*I)^2
        one_m_NI_sq = (8'sd1<<34) - (($signed({{38{NI[37]}},NI})* NI)>>>34); //q8.34
        //n^2*(1-(N*I)^2)
        n_sq_times_one_m_NI_sq = ((($signed({1'b0,n})*$signed({1'b0,n}))) * $signed({{36{one_m_NI_sq[41]}},one_m_NI_sq})) >>> 16; //q2.32*q8.34=.66 -> q6.50
        //1-n^2*(1-(N*I)^2)
        one_m_n_sq_times_one_m_NI_sq = ((12'sd1<<50) - n_sq_times_one_m_NI_sq); //q6.50

    end
    //sqrt(1-n^2*(1-(N*I)^2))
    sqrt_seq #(.W(56)) sqrt_seq_inst (
            .clk    (CLOCK2_50),
            .rst    (reset),
            .start  (sqrt_start),
            .radicand(one_m_n_sq_times_one_m_NI_sq), //q6.50
            .root   (sqrt_one_m_n_sq_times_one_m_NI_sq), //q3.25
            .done   (done1)
        );

    // vectors as arrays: index 0 = x, 1 = y, 2 = z
    wire signed [1:0]  I [0:2];                         //q2
    wire signed [35:0] N [0:2];                         //q2.34
    reg  signed [91:0] T [0:2];                         //q2.34
    assign I[0] = Ix;  assign I[1] = Iy;  assign I[2] = Iz;
    assign N[0] = Nx;  assign N[1] = Ny;  assign N[2] = Nz;

    // wire signed [27:0] sqrt_one_m_n_sq_times_one_m_NI_sq;   //q3.25
    // reg signed [37:0] NI; //q4.34
    //T=n*I+(n*(-N*I)-(1-n^2*(1-(N*I)^2))^0.5)*N
    // T = n*I + (n*(-N*I)-sqrt_one_m_n_sq_times_one_m_NI_sq)*N
    reg signed [19:0] nI;                               //q4.16
    reg signed [91:0] vN;                               //q8.84
    reg signed [55:0] n_multiply_minusNI;               //q6.50
    integer k;
    always @(posedge CLOCK2_50) begin
        if (done1) begin                                // root is for the current N
            for (k = 0; k < 3; k = k + 1) begin
                nI = $signed({1'b0, n}) * I[k];           // q2.16 * q2   -> q4.16
                n_multiply_minusNI = $signed({1'b0, n}) * (-NI); // q2.16 * q4.34 -> q6.50
                vN = (n_multiply_minusNI - $signed({19'd0, sqrt_one_m_n_sq_times_one_m_NI_sq, 25'd0}))*N[k];// (q6.50-q3.25<<25) * q2.34 -> q8.84
                T[k] <= (nI <<< 68) + (vN);    // q4.16<<68 -> q8.84
            end
        end
    end
    assign Tx = T[0];  assign Ty = T[1];  assign Tz = T[2];
    assign done = done_r;

    
endmodule

// ==========================================================================
// sqrt_seq.v
// ==========================================================================
// Sequential integer square root: root = floor(sqrt(radicand)).
//
// Digit-by-digit (restoring) method, one result bit per clock, so a W-bit
// radicand takes W/2 cycles. radicand is latched on a start pulse; done goes
// high for one clock with the result, and root holds until the next result.
// A start while busy is ignored. W must be even.
module sqrt_seq #(
    parameter W = 16
)(
    input  wire           clk,
    input  wire           rst,         // active-low, asynchronous
    input  wire           start,
    input  wire [W-1:0]   radicand,
    output reg  [W/2-1:0] root,
    output reg            done
);
    localparam N = W / 2;

    reg [W-1:0]           x;           // radicand bits not yet consumed, MSB first
    reg [N+1:0]           rem;         // partial remainder
    reg [N-1:0]           q;           // partial root
    reg [$clog2(N+1)-1:0] cnt;         // bit pairs left
    reg                   busy;

    // Bring down the next two radicand bits, then try subtracting (q << 2) | 1
    wire [N+1:0] rem_next = {rem[N-1:0], x[W-1:W-2]};
    wire [N+1:0] trial    = {q, 2'b01};
    wire         fits     = (rem_next >= trial);

    always @(posedge clk or negedge rst) begin
        if (!rst) begin
            x    <= 0;
            rem  <= 0;
            q    <= 0;
            cnt  <= 0;
            busy <= 1'b0;
            root <= 0;
            done <= 1'b0;
        end else if (!busy) begin
            done <= 1'b0;
            if (start) begin
                x    <= radicand;
                rem  <= 0;
                q    <= 0;
                cnt  <= N;
                busy <= 1'b1;
            end
        end else begin
            x   <= x << 2;
            rem <= fits ? rem_next - trial : rem_next;
            q   <= {q[N-2:0], fits};
            cnt <= cnt - 1'b1;
            if (cnt == 1) begin
                busy <= 1'b0;
                root <= {q[N-2:0], fits};
                done <= 1'b1;
            end
        end
    end

endmodule

// ==========================================================================
// div_seq.v
// ==========================================================================
// Sequential unsigned divide: quotient = floor(dividend / divisor).
//
// Restoring method, one quotient bit per clock. The quotient must fit in WQ
// bits (dividend < divisor << WQ), so only the low WQ dividend bits need a
// step and a divide takes WQ cycles. dividend and divisor are latched on a
// start pulse; done goes high for one clock with the result. A start while
// busy is ignored.
module div_seq #(
    parameter WN = 88,                 // dividend width
    parameter WD = 44,                 // divisor width
    parameter WQ = 36                  // quotient width
)(
    input  wire          clk,
    input  wire          rst,          // active-low, asynchronous
    input  wire          start,
    input  wire [WN-1:0] dividend,
    input  wire [WD-1:0] divisor,
    output reg  [WQ-1:0] quotient,
    output reg           done
);
    reg [WQ-1:0]           lo;         // dividend bits not yet brought down, MSB first
    reg [WD:0]             rem;        // partial remainder, < d before each step
    reg [WD-1:0]           d;
    reg [WQ-1:0]           q;          // partial quotient
    reg [$clog2(WQ+1)-1:0] cnt;        // quotient bits left
    reg                    busy;

    // Bring down the next dividend bit, then try subtracting the divisor
    wire [WD:0] rem_next = {rem[WD-1:0], lo[WQ-1]};
    wire        fits     = (rem_next >= {1'b0, d});

    always @(posedge clk or negedge rst) begin
        if (!rst) begin
            lo       <= 0;
            rem      <= 0;
            d        <= 0;
            q        <= 0;
            cnt      <= 0;
            busy     <= 1'b0;
            quotient <= 0;
            done     <= 1'b0;
        end else begin
            done <= 1'b0;
            if (!busy) begin
                if (start) begin
                    rem  <= dividend[WN-1:WQ];   // < divisor, so it fits in WD bits
                    lo   <= dividend[WQ-1:0];
                    d    <= divisor;
                    q    <= 0;
                    cnt  <= WQ;
                    busy <= 1'b1;
                end
            end else begin
                lo  <= lo << 1;
                rem <= fits ? rem_next - {1'b0, d} : rem_next;
                q   <= {q[WQ-2:0], fits};
                cnt <= cnt - 1'b1;
                if (cnt == 1) begin
                    busy     <= 1'b0;
                    quotient <= {q[WQ-2:0], fits};
                    done     <= 1'b1;
                end
            end
        end
    end

endmodule