Merge branch 'Oliver_FPA'

scripts/fpu_test_gen.py

@@ -18,6 +18,7 @@ def dtype_size(dtype):
     else:
         raise ValueError(f"Unknown dtype {dtype}")
 
+
 def generate_numbers(cases, dtype=np.float16):
     dsize = dtype_size(dtype)
     x = np.frombuffer(os.urandom(cases * dsize), dtype=dtype)
@@ -38,7 +39,7 @@ def generate_fp_vector(cases, filename, dtype=np.float16, big_endian=False, comp
                 t(y.tobytes()[i * dsize:i * dsize + dsize]).hex(),
                 t(sum.tobytes()[i * dsize:i * dsize + dsize]).hex(),
                 t(mul.tobytes()[i * dsize:i * dsize + dsize]).hex(),
-            ]) + '\n')
+            ]) + f'  // {x[i]:10.6e} {y[i]:10.6e} {sum[i]:10.6e} {mul[i]:10.6e}\n')
     if comp_file is not None:
         gt = x > y
         lt = x < y

simulation/modelsim/wave_floating_tb.do

@@ -0,0 +1,28 @@
+onerror {resume}
+quietly WaveActivateNextPane {} 0
+add wave -noupdate -label CLK /floating_tb/clk
+add wave -noupdate -label RESET /floating_tb/reset
+add wave -noupdate -label {INPUT A} -radix hexadecimal /floating_tb/input_a
+add wave -noupdate -label {INPUT B} -radix hexadecimal /floating_tb/input_b
+add wave -noupdate -label ADD -radix hexadecimal /floating_tb/result_add
+add wave -noupdate -label {ADD expected} -radix hexadecimal /floating_tb/expected_add
+add wave -noupdate -label MULT -radix hexadecimal /floating_tb/result_mult
+add wave -noupdate -label {MULT expected} -radix hexadecimal /floating_tb/expected_mult
+TreeUpdate [SetDefaultTree]
+WaveRestoreCursors {{Cursor 1} {81 ps} 0}
+quietly wave cursor active 1
+configure wave -namecolwidth 150
+configure wave -valuecolwidth 100
+configure wave -justifyvalue left
+configure wave -signalnamewidth 0
+configure wave -snapdistance 10
+configure wave -datasetprefix 0
+configure wave -rowmargin 4
+configure wave -childrowmargin 2
+configure wave -gridoffset 0
+configure wave -gridperiod 1
+configure wave -griddelta 40
+configure wave -timeline 0
+configure wave -timelineunits ns
+update
+WaveRestoreZoom {0 ps} {124 ps}

src/FPA_module_test.sv

@@ -156,69 +156,156 @@ module floating_product #(parameter N=16, M=4)(input_1, input_2, product, clk, r
 	input logic [N-1:0] input_1, input_2;
 	input logic clk, reset;
 	output logic [N-1:0] product;
-
+	
+	logic [2*N-1:0] D0 [7:0];
+	logic [2*N-1:0] Q0 [7:0];
+	logic [2*N-1:0] Q1 [7:0];
+	logic [2*N-1:0] Q2 [7:0];
+	
+	always_comb
+	begin
+		D0[0][N-1:0] = input_1;
+		D0[1][N-1:0] = input_2;
+		D0[2] = 0; // product (output)
+		D0[3] = 0; // multiple for multiplying the mantissa
+		D0[4] = 0; // flag for return NaN
+		D0[5] = 0; // flag for return infinity
+		D0[6] = 0; // flag for return zero
+		D0[7][N-3-M:0] = input_1[N-3-M:0]; // for storing the mantissa of input 1
+	end
 	// sign_x = x[N-1]
 	// exponent_x = x[N-2:N-2-M]
 	// mantissa_x = x[N-3-M:0]
-
-	logic [N-2:N-2-M] sum;
-	logic [2*(N-3-M):0] mult;
-	logic [2*(N-3-M):0] D0 [4:0];
-	logic [2*(N-3-M):0] Q0 [4:0];
-	logic [2*(N-3-M):0] Q1 [4:0];
-	logic [2*(N-3-M):0] Q2 [4:0];
 	
-	// First pipeline stage
+	pipe #(.N(2*N-1), .K(7)) pipe0(.clk(clk), .reset(reset), .D(D0), .Q(Q0));
+	
 	always_comb
+	begin
+		// if input_1 or input_2 is NaN then return NaN
+		if ((Q0[0][N-2:N-2-M] == (1<<M) && Q0[0][N-3-M:0] != 0) || (Q1[0][N-2:N-2-M] == (1<<M) && Q1[0][N-3-M:0] != 0))
 		begin
-			D0[0] = input_1;
-			D0[1] = input_2;
-			D0[2] = 0; // product
-			D0[3] = 0; // sum
-			D0[4] = 0; // mult
+			Q0[2][N-1] = 1;
+			Q0[2][N-2:N-2-M] = (1 << (M+1)) - 1;
+			Q0[2][N-3-M] = 1;
+			Q0[2][N-4-M:0] = 0;
+			Q0[4][0] = 1;
 		end
-		
-	pipe#(.N(2*(N-3-M)), .K(4)) pipe0(.clk(clk), .reset(reset), .D(D0), .Q(Q0));
-
-	// We have assigned an {M+1} bit exponent so we must have a 2^{M} offset
-	assign Q0[3] = Q0[0][N-2:N-2-M] + Q0[1][N-2:N-2-M];
-	assign Q0[2][N-2:N-2-M] = Q0[3] - (1'b1 << M) + 2;
+		// if input 1 is infinity then return infinity
+		else if (Q0[0][N-2:N-2-M] == (1<<M))
+		begin
+			Q0[2][N-1] = Q0[0][N-1] ^ Q0[1][N-1];
+			Q0[2][N-2:N-2-M] = (1 << (M+1)) - 1;
+			Q0[2][N-3-M:0] = 0;
+			Q0[5] = 1;
+			// if input 2 is zero then return NaN
+			if (($signed(Q0[1][N-2:N-2-M]) == (-1*((1<<M)-1))) && (Q0[1][N-3-M:0] == 0))
+			begin
+				Q0[2][N-1] = 1;
+				Q0[2][N-2:N-2-M] = (1 << (M+1)) - 1;
+				Q0[2][N-3-M] = 1;
+				Q0[2][N-4-M:0] = 0;
+				Q0[4][0] = 1;
+			end
+		end
+		// if input 2 is infinity then return infinity
+		else if (Q0[1][N-2:N-2-M] == (1<<M))
+		begin
+			Q0[2][N-1] = Q0[0][N-1] ^ Q0[1][N-1];
+			Q0[2][N-2:N-2-M] = (1 << (M+1)) - 1;
+			Q0[2][N-3-M:0] = 0;
+			Q0[5][0] = 1;
+			// if input 1 is zero then return NaN
+			if (($signed(Q0[0][N-2:N-2-M]) == (-1*((1<<M)-1))) && (Q0[0][N-3-M:0] == 0))
+			begin
+				Q0[2][N-1] = 1;
+				Q0[2][N-2:N-2-M] = (1 << (M+1)) - 1;
+				Q0[2][N-3-M] = 1;
+				Q0[2][N-4-M:0] = 0;
+				Q0[4][0] = 1;
+			end
+		end
+		// if input 1 is zero then return zero
+		else if (($signed(Q0[0][N-2:N-2-M]) == (-1*((1<<M)-1))) && (Q0[0][N-3-M:0] == 0))
+		begin
+			Q0[2][N-1] = Q0[0][N-1] ^ Q0[1][N-1];
+			Q0[2][N-2:N-2-M] = 0;
+			Q0[2][N-3-M:0] = 0;
+			Q0[6][0] = 1;
+		end
+		// if input 2 is zero then return zero
+		else if (($signed(Q0[1][N-2:N-2-M]) == (-1*((1<<M)-1))) && (Q0[1][N-3-M:0] == 0))
+		begin
+			Q0[2][N-1] = Q0[0][N-1] ^ Q0[1][N-1];
+			Q0[2][N-2:N-2-M] = 0;
+			Q0[2][N-3-M:0] = 0;
+			Q0[6][0] = 1;
+		end
+	end
+	
+	pipe #(.N(2*N-1), .K(7)) pipe1(.clk(clk), .reset(reset), .D(Q0), .Q(Q1));
 	
-	// Second pipeline stage
-	pipe#(.N(2*(N-3-M)), .K(4)) pipe1(.clk(clk), .reset(reset), .D(Q0), .Q(Q1));
-
 	always_comb
+	begin
+	// If none of the return flags have been triggered
+		if ((Q1[4][0] && Q1[5][0] && Q1[6][0]) != 1)
 		begin
-				// Setting the mantissa of the output
-				Q1[4] = Q1[0][N-3-M:0] * Q1[1][N-3-M:0];
-				if (Q1[4][N-3-M]) Q1[2][N-3-M:0] = Q1[4][2*(N-3-M):2*(N-3-M)-9];
-				else Q1[2][N-3-M:0] = Q1[4][2*(N-3-M):2*(N-3-M)-9] << 1;
-				Q1[2][N-1] = Q1[0][N-1] ^ Q1[1][N-1];
+			// if msb of input 1 mantissa is not 1 then shift left by 1 and reduce exponent by 1
+			if (Q1[0][N-3-M] != 1)
+			begin
+				Q1[0][N-3-M:0] = Q1[0][N-3-M:0] << 1;
+				Q1[0][N-2:N-2-M] = Q1[0][N-2:N-2-M] - 1;
+			end
+			// if msb of input 2 mantissa is not 1 then shift left by 1 and reduce exponent by 1
+			if (Q1[1][N-3-M] != 1)
+			begin
+				Q1[1][N-3-M:0] = Q1[1][N-3-M:0] << 1;
+				Q1[1][N-2:N-2-M] = Q1[1][N-2:N-2-M] - 1;
+			end
+			Q1[2][N-1] = Q1[2][N-1] ^ Q1[2][N-1]; // ouput sign = input_1 sign xor input_2 sign
+			Q1[2][N-2:N-2-M] = Q1[0][N-2:N-2-M] + Q1[1][N-2:N-2-M] - (1<<M); // out exp = in1 exp + in2 exp - 2**M
+			// multiplying mantissa
+			for (int i = 0; i<N-2-M; i++)
+			begin
+				// multiplying each digit of input 2 by all of input 1 shifting the bits left each time and adding result to the array
+				if (Q1[1][i] == 1)
+				begin
+					Q1[3] = Q1[3] +(Q1[7]<<i);
+				end
+				else
+				begin
+					Q1[3] = Q1[3] + 0;
+				end
+			end
+			// Assigning the top set of bits to the mantissa of the output 
+			Q1[2][N-3-M:0] = Q1[3][2*N-1:2*N-1-N-3-M];
 		end
+	end
 	
-	//Final Pipeline stage
-	pipe#(.N(2*N-3-M), .K(4)) pipe2(.clk(clk), .reset(reset), .D(Q1), .Q(Q2));
+	pipe #(.N(2*N-1), .K(7)) pipe2(.clk(clk), .reset(reset), .D(Q1), .Q(Q2));
 	assign product = Q2[2][N-1:0];
+	
 endmodule : floating_product
 
 
 
-module pipe #(parameter N=16)(clk, reset, Q, D);
+module pipe #(parameter N, K)(clk, reset, Q, D);
 	input logic clk, reset;
-	input logic [N-1:0] D;
-	output reg [N-1:0] Q;
-	reg [N-1:0] in_pipe;
-	
-	always @(posedge clk or negedge reset)
-		begin
-			if(reset) in_pipe = 0;
-			else in_pipe = D;
-		end
+	input logic [N:0] D [K:0];
+	output reg [N:0] Q [K:0];
+	reg [N:0] in_pipe [K:0];
 	
-	always @(posedge clk or negedge reset)
+	always_ff @(posedge clk or negedge reset)
 		begin
-			if(reset) Q = 0;
-			else Q = in_pipe;
+			if(reset) 
+				begin
+					in_pipe <= '{default:0};
+					Q <= '{default:0};
+				end
+			else 
+				begin
+					in_pipe <= D;
+					Q <= in_pipe;
+				end
 		end
 endmodule : pipe
 
@@ -228,22 +315,37 @@ module floating_tb;
 	reg reset, clk;
 	logic [15:0] input_a, input_b, result_add, result_mult;
 	logic [4:0] diff;
+	logic [15:0] expected_add, expected_mult;
 
-	floating_add adder1(.input_1(input_a), .input_2(input_b), .sum(result_add), .diff(diff));
-
-	floating_product multiplier1(.input_1(input_a), .input_2(input_b), .product(result_mult));
+	floating_add adder1(.input_1(input_a), .input_2(input_b), .sum(result_add), .diff(diff), .clk(clk), .reset(reset));
 
+	floating_product multiplier1(.input_1(input_a), .input_2(input_b), .product(result_mult), .clk(clk), .reset(reset));
+	
+	initial forever #5 clk = ~clk;
+	localparam PIPELINES = 3;
 
 	reg [15:0] test_mem [29:0][3:0];
 
-	initial $readmemh("../../scripts/fp16_test.hex", test_mem);
-
+	initial $readmemh("scripts/fp16_test.hex", test_mem);
 
 	initial begin
         static int num_err = 0;
         static int num_tests = $size(test_mem) * 2;
 
-        for (int i=0; i < $size(test_mem); i++) begin
+		clk = 0;
+		reset = 1;
+		  
+		#15;
+		reset = 0;
+
+		expected_add = 0;
+		expected_mult = 0;
+
+        for (int i=0; i < $size(test_mem)+PIPELINES; i++) begin
+			if(i >= PIPELINES) begin
+				expected_add = test_mem[i-PIPELINES][2];
+				expected_mult = test_mem[i-PIPELINES][3];
+			end
             input_a = test_mem[i][0];
             input_b = test_mem[i][1];
 
@@ -260,7 +362,10 @@ module floating_tb;
                 num_err = num_err + 1;
             end
 
-        end
+		end
+		expected_add = 0;
+		expected_mult = 0;
+		#50;
         $display("Passed %d of %d tests", num_tests-num_err, num_tests);
         $finish();
 	end

src/fpu16/fp_adder.sv

@@ -0,0 +1,168 @@
+typedef enum logic [1:0]{
+    greater_a,
+    greater_b,
+    equal_ab
+
+} grater_state;
+
+module fp_adder#(parameter N=16, M=5)(input_a, input_b, output_z, clk, reset);
+    localparam K=N-M-1;  // Size of mantissa
+
+    input reg [N-1:0] input_a, input_b;
+    input logic clk, reset;
+    output reg [N-1:0] output_z;
+
+    reg [K-1:0] a_m0, b_m0; // mantissa
+    reg [K-1:0] a_m1, b_m1, z_m2;
+
+    reg [K*2-1:0] z_m1a, z_m1b, z_m1z;  // Double mantissa
+    reg z_m1s;
+
+    reg [M-1:0] a_e0, b_e0; // exponent
+    reg [M-1:0] z_e1, z_e2;
+
+    reg a_s0, b_s0; // sign
+    reg a_s1, b_s1, z_s1, z_s2;
+
+    grater_state greater; // 01 for a, 10 for b, 11 for both and 00 for neither
+    reg [M:0] abs; // For the absolute difference between exponents
+
+    always_comb begin
+        output_z = {z_s2, z_e2, z_m2};
+        z_m1a = {a_m1, {K{1'b0}}};
+        z_m1b = {b_m1, {K{1'b0}}};
+
+        // If a has the bigger exponent
+        if (greater == greater_a)
+            begin
+                // If the signs are the same then add
+                if (a_s1 == b_s1) {z_m1s, z_m1z} = z_m1a+(z_m1b >> abs - 2);
+                    // If they are different then subtract
+                else {z_m1s, z_m1z} = z_m1a-(z_m1b >> abs - 2);
+            end
+            // If b has the bigger exponent
+        else if (greater == greater_b)
+            begin
+                // If the signs are the same then add
+                if (a_s1 == b_s1) {z_m1s, z_m1z} = z_m1b+(z_m1a >> abs - 2);
+                    // If they are different then subtract
+                else {z_m1s, z_m1z} = z_m1b-(z_m1a >> abs - 2);
+            end
+            // If the exponents are equal
+        else
+            begin
+                // If the signs are the same then add
+                if (a_s1 == b_s1) {z_m1s, z_m1z} = (z_m1a + z_m1b) >> 1;
+                    // If the signs are different then subtract
+                else
+                    begin
+                        // First checking which has the bigger mantissa
+                        if (a_m1 > b_m1) {z_m1s, z_m1z} = z_m1a-z_m1b;
+                        else if (b_m1 > a_m1) {z_m1s, z_m1z} = z_m1b-z_m1a;
+                            // If the mantissa are the same as well then the result should be 0
+                        else {z_m1s, z_m1z} = 0;
+                    end
+            end
+    end
+
+    always_ff @(posedge clk)
+        begin
+            if (~reset)
+                begin
+                    // Unpacking the inputs
+                    a_m0 <= input_a[K-1:0];
+                    a_e0 <= input_a[N-2:K];
+                    a_s0 <= input_a[N-1];
+
+                    b_m0 <= input_b[K-1:0];
+                    b_e0 <= input_b[N-2:K];
+                    b_s0 <= input_b[N-1];
+
+                    // Second stage
+                    a_m1 <= a_m0;
+                    a_s1 <= a_s0;
+                    b_m1 <= b_m0;
+                    b_s1 <= b_s0;
+
+                    z_e2 <= z_e1;
+                    z_s2 <= z_s1;
+
+                    // If input_a has the bigger exponent then flag it with greater and find the absolute difference
+                    if (a_e0 > b_e0)
+                        begin
+                            greater <= greater_a;
+                            abs <= a_e0-b_e0;
+                            z_s1 <= a_s0;
+                            z_e1 <= a_e0;
+                        end
+
+                        // If input_a has the bigger exponent then flag it with greater and find the absolute difference
+                    else if (b_e0 > a_e0)
+                        begin
+                            greater <= greater_b;
+                            abs <= b_e0-a_e0;
+                            z_s1 <= b_s0;
+                            z_e1 <= b_e0;
+                        end
+
+                        // If the inputs have equal exponent
+                    else
+                        begin
+                            greater <= equal_ab;
+                            abs <= -1;
+                            z_e1 <= a_e0;
+                            // Assigning the overall sign based on the difference between the mantissa
+                            if (a_m0 > b_m0) z_s1 <= a_s0;
+                            else if (b_m0 > a_m0) z_s1 <= b_s0;
+                            else z_s1 <= 0;
+                        end
+
+                    // Condition for overflow is that it sets the output to the larger input
+                    if (abs > K) // Shifting by N-1-M would give 0
+                        begin
+                            z_m2 <= (greater == greater_a) ? a_m1 : b_m1;
+
+                            // Input a is larger and is translated to the output
+                            // if (greater == greater_a) z_m0 <= a_m1;
+
+                            // Input b is larger and is translated to the output
+                            // else if (greater == greater_b) z_m0 <= b_m1;
+
+                            // Shouldn't happen as abs should be 0 for this to occur
+                            // else begin
+                            // 	if (a_m1 >= b_m1) z_m0 <= a_m1; // Equal exponents but a has the larger mantissa
+                            // 	else if (b_m1 > a_m1) z_m0 <= b_m1; // Equal exponents but b has the larger mantissa
+                            // end
+                        end
+
+                    else
+                        begin
+							z_m2 <= z_m1z[K*2-1:K];
+                        end
+                end
+            else
+                begin
+                    a_m0 <= 0;
+                    a_e0 <= 0;
+                    a_s0 <= 0;
+
+                    b_m0 <= 0;
+                    b_e0 <= 0;
+                    b_s0 <= 0;
+
+                    a_m1 <= 0;
+                    b_m1 <= 0;
+                    z_e2 <= 0;
+                    z_s2 <= 0;
+
+                    z_s1 <= 0;
+                    z_e1 <= 0;
+                    z_s2 <= 0;
+                    z_e2 <= 0;
+                    z_m2 <= 0;
+
+                    greater <= equal_ab;
+                    abs <= 0;
+                end
+        end
+endmodule : fp_adder

src/fpu16/fp_mult.sv

@@ -0,0 +1,250 @@
+`include "lead_zero.sv"
+
+
+module lead_zeros_enc #(K=10)(in, out);
+    localparam L = $clog2(K) +1;
+
+    input wire [K-1:0] in;
+    output logic [L-1:0] out;
+    always_comb begin
+        if(K == 10) begin
+            casez (in)
+                10'b00000_00000 :out = 4'd10;
+                10'b00000_0000? :out = 4'd9;
+                10'b00000_000?? :out = 4'd8;
+                10'b00000_00??? :out = 4'd7;
+                10'b00000_0???? :out = 4'd6;
+                10'b00000_????? :out = 4'd5;
+                10'b0000?_????? :out = 4'd4;
+                10'b000??_????? :out = 4'd3;
+                10'b00???_????? :out = 4'd2;
+                10'b0????_????? :out = 4'd1;
+                default : out = 4'd0;
+            endcase
+        end
+        else if (K == 23) begin
+            casez (in)
+                23'b000_00000_00000_00000_00000 :out = 5'd23;
+                23'b000_00000_00000_00000_0000? :out = 5'd22;
+                23'b000_00000_00000_00000_000?? :out = 5'd21;
+                23'b000_00000_00000_00000_00??? :out = 5'd20;
+                23'b000_00000_00000_00000_0???? :out = 5'd19;
+                23'b000_00000_00000_00000_????? :out = 5'd18;
+                23'b000_00000_00000_0000?_????? :out = 5'd17;
+                23'b000_00000_00000_000??_????? :out = 5'd16;
+                23'b000_00000_00000_00???_????? :out = 5'd15;
+                23'b000_00000_00000_0????_????? :out = 5'd14;
+                23'b000_00000_00000_?????_????? :out = 5'd13;
+                23'b000_00000_0000?_?????_????? :out = 5'd12;
+                23'b000_00000_000??_?????_????? :out = 5'd11;
+                23'b000_00000_00???_?????_????? :out = 5'd10;
+                23'b000_00000_0????_?????_????? :out = 5'd9;
+                23'b000_00000_?????_?????_????? :out = 5'd8;
+                23'b000_0000?_?????_?????_????? :out = 5'd7;
+                23'b000_000??_?????_?????_????? :out = 5'd6;
+                23'b000_00???_?????_?????_????? :out = 5'd5;
+                23'b000_0????_?????_?????_????? :out = 5'd4;
+                23'b000_?????_?????_?????_????? :out = 5'd3;
+                23'b00?_?????_?????_?????_????? :out = 5'd2;
+                23'b0??_?????_?????_?????_????? :out = 5'd1;
+                default : out = 5'd0;
+            endcase
+        end
+    end
+endmodule : lead_zeros_enc
+
+module fp_mult #(parameter N=16, M=5)(input_a, input_b, output_z, clk, reset);
+    localparam K=N-M-1;  // Size of mantissa
+    localparam L=$clog2(10);  // Size of shifting representation
+
+    input logic [N-1:0] input_a, input_b;
+    input logic clk, reset;
+    output logic [N-1:0] output_z;
+
+    // ====================
+    // Stage 0 store input
+    // ====================
+
+    reg [M-1:0] a_e0, b_e0;
+    reg [K-1:0] a_m0, b_m0;
+    reg  a_s0, b_s0;
+
+    always_ff @(posedge clk) begin
+        {a_s0, a_e0, a_m0} <= input_a;
+        {b_s0, b_e0, b_m0} <= input_b;
+    end
+
+    // ====================
+    // Stage 1 denormalise0
+    // ====================
+
+    reg [M-1:0] a_e1, b_e1;
+    reg [K-1:0] a_m1, b_m1;
+    reg z_s1;
+
+    // leading zeros
+    wire [L-1:0] a_z1w, b_z1w;
+    reg [L-1:0] a_z1, b_z1;
+
+    lead_zeros_enc #(.K(K)) lz_a1(a_m1, a_z1w);
+    lead_zeros_enc #(.K(K)) lz_b1(b_m1, b_z1w);
+
+    always_ff @(posedge clk) begin
+        z_s1 <= a_s0 ^ b_s0;  // We don't need to track and propagate signs
+        {a_e1, a_m1} <= {a_e0, a_m0};
+        {b_e1, b_m1} <= {b_e0, b_m0};
+        a_z1 <= a_z1w;
+        b_z1 <= b_z1w;
+    end
+
+    // ====================
+    // Stage 2 denormalise1
+    // ====================
+
+    reg [M-1:0] a_e2, b_e2;
+    reg [K-1:0] a_m2, b_m2;
+    reg z_s2;
+
+    always_ff @(posedge clk) begin
+        a_e2 <= a_e1 - a_z1;
+        a_m2 <= a_m1 << a_z1;
+
+        b_e2 <= b_e1 - b_z1;
+        b_m2 <= b_m1 << b_z1;
+
+        z_s2 <= z_s1;
+    end
+
+    // ====================
+    // Stage 3 product
+    // ====================
+
+    reg z_s3;
+    reg [M-1:0] z_e3;
+    reg [K*2-1:0] z_p3;  // product is double mantissa
+
+    always_ff @(posedge clk) begin
+        z_e3 <= a_e2 + b_e2 + 1;
+        z_p3 <= a_m2 * b_m2;
+        z_s3 <= z_s2;
+    end
+
+    // ====================
+    // Stage 4 unpack/normalise0
+    // ====================
+
+    reg z_s4;
+    reg [M-1:0] z_e4;
+    reg [K-1:0] z_m4;
+    reg [1:0] z_bits4; // guard and round_bit
+    reg z_sticky4;
+
+
+    wire [K-1:0] p3;
+    // leading zeros
+    wire [3:0] z_z4w;
+    reg [3:0] z_z4;
+
+    assign p3 = z_p3[K*2-1:K];
+    lead_zeros_enc #(.K(K)) lz_p3(p3, z_z4w);
+
+    always_ff @(posedge clk) begin
+        z_m4 <= p3;
+        z_bits4 <= z_p3[K:K-2];
+        z_sticky4 <= z_p3[K-2:0] != 0;
+        z_z4 <= z_z4w;
+
+        z_e4 <= z_e3;
+        z_s4 <= z_s3;
+    end
+
+    // ====================
+    // Stage 5 normalise1
+    // ====================
+
+    reg z_s5;
+    reg [M-1:0] z_e5;
+    reg [K-1:0] z_m5;
+
+    reg [1:0] z_bits5; // guard and round_bit
+    reg z_sticky5;
+
+    always_ff @(posedge clk) begin
+        z_e5 <= z_e4 - z_z4;
+        z_m5 <= z_m4 << z_z4;
+
+        case (z_z4)
+            4'd0: begin
+                z_bits5 <= z_bits4;
+            end
+            4'd1: begin
+                z_m5[0] <= z_bits4[0];
+                z_bits5 <= {z_bits4[1], 1'b0};
+            end
+            4'd2: begin
+                z_m5[1-:2] <= z_bits4; // ?
+                z_bits5 <= 2'b00;
+            end
+            default : begin
+                z_m5[z_z4-1-:2] <= z_bits4; // ?
+                z_bits5 <= 2'b00;
+            end
+        endcase
+
+        z_sticky5 <= z_sticky4;
+        z_s5 <= z_s4;
+    end
+
+    // ====================
+    // Stage 6 normalise2
+    // ====================
+    reg z_s6;
+    reg [M-1:0] z_e6;
+    reg [K-1:0] z_m6;
+
+    reg [1:0] z_bits6;
+    reg z_sticky6;
+
+
+    always_ff @(posedge clk) begin
+        // TODO: Fix this stage
+        // if ($signed(z_e5) < -126) begin
+        //     shift6 <= $signed(z_e5) -126;
+        // end else begin
+        //     shift6 <= 0;
+        // end
+
+        z_s6 <= z_s5;
+        z_e6 <= z_e5;
+        z_m6 <= z_m5;
+        z_bits6 <= z_bits5;
+        z_sticky6 <= z_sticky5;
+    end
+
+    // ====================
+    // Stage 7 round
+    // ====================
+    reg z_s7;
+    reg [M-1:0] z_e7;
+    reg [K-1:0] z_m7;
+
+    always_ff @(posedge clk) begin
+        if (z_bits6[0] && (z_bits6[1] | z_sticky6 | z_m6[0])) begin
+            z_m7 <= z_m6 + 1;
+            // TODO: Check for all 1s
+        end else begin
+            z_m7 <= z_m6;
+        end
+        z_e7 <= z_e6;
+        z_s7 <= z_s6;
+    end
+
+    // ====================
+    // Stage 8 pack
+    // ====================
+    // This stage is skipped as it checks for overflow
+
+
+    always_comb output_z = {z_s7, z_e7, z_m7};
+
+endmodule : fp_mult

src/fpu16/fp_product.sv

@@ -0,0 +1,222 @@
+typedef enum logic [1:0]{
+    flag_none,
+    flag_nan,
+    flag_inf,
+    flag_zero
+} product_flags;
+
+module fp_product #(parameter N=16, M=5)(input_a, input_b, output_z, clk, reset);
+    localparam K=N-M-1;  // Size of mantissa
+
+	input logic [N-1:0] input_a, input_b;
+	input logic clk, reset;
+	output logic [N-1:0] output_z;
+	
+	reg [K-1:0] a_m0, b_m0, z_m1; // mantissa
+	reg [K-1:0] a_m1, b_m1, z_m2; // mantissa
+	reg [K-1:0] a_m2, b_m2, z_m3; // mantissa
+	reg [K-1:0] z_m4;
+	
+	reg [M:0] a_e0, b_e0, z_e1; // exponent
+	reg [M:0] a_e1, b_e1, z_e2; // exponent
+	reg [M:0] a_e2, b_e2, z_e3; // exponent
+	reg [M:0] z_e4;
+	
+	reg a_s0, b_s0, z_s1; // sign
+	reg a_s1, b_s1, z_s2; // sign
+	reg a_s2, b_s2, z_s3; // sign
+	reg z_s4; // sign
+
+    product_flags flags1; // 00 for no flag, 01 for NaN, 10 for infinity, 11 for zero
+	reg [2*K-1:0] z_p3; // For storing the product of the two mantissa
+
+	always_comb begin
+		// Packing the output
+		output_z = {z_s4, z_e4 + (1<<M), z_m4};
+	end
+
+	always_ff @(posedge clk)
+	begin
+		if (~reset)
+		begin
+			// Unpacking the inputs
+			a_s0 <= input_a[N-1];
+			a_e0 <= input_a[N-2:N-2-M] - (1<<M);
+			a_m0 <= input_a[N-3-M:0];
+			
+			b_s0 <= input_b[N-1];
+			b_e0 <= input_b[N-2:N-2-M] - (1<<M);
+			b_m0 <= input_b[N-3-M:0];
+			
+			// Untouched pipelined registers
+			a_s1 <= a_s0;
+			a_s2 <= a_s1;
+			a_e1 <= a_e0;
+			a_m1 <= a_m0;
+			b_s1 <= b_s0;
+			b_s2 <= b_s1;
+			b_e1 <= b_e0;
+			b_m1 <= b_m0;
+			z_s2 <= z_s1;
+			z_s4 <= z_s3;
+			z_e4 <= z_e3;
+			z_m2 <= z_m1;
+			z_m3 <= z_m2;
+
+			z_s3 <= a_s2 ^ b_s2; //signs xor together
+			z_e3 <= a_e2 + b_e2 - (1<<M); // exponents added together subtracting one offset
+			z_p3 <= a_m2 * b_m2; // mantissa multiplied together and the most significant bits are stored in the output mantissa
+			z_m4 <= z_p3[2*K-1:K];
+
+			// If input a or input b is NaN then return NaN
+			if (((a_e0 ==(1<<M)) && (a_m0 != 0)) || ((b_e0 == (1<<M)) && (b_m0 != 0)))
+			begin
+				z_s1 <= 1;
+				z_e1 <= (1 << (M+1)) - 1;
+				z_m1[K-1] <= 1;
+				z_m1[K-2:0] <= 0;
+				flags1 <= flag_nan;
+			end
+			// If a is infinity then return infinity
+			else if (a_e0 == (1<<M))
+			begin
+				// Unless b is zero, then you return NaN instead
+				if (($signed(b_e0) == (-1*((1<<M)-1))) && b_m0 == 0)
+				begin
+					z_s1 <= 1;
+					z_e1 <= (1 << (M+1)) - 1;
+					z_m1[K-1] <= 1;
+					z_m1[K-2:0] <= 0;
+					flags1 <= flag_nan;
+				end
+				else
+				// Returning infinity
+				begin
+					z_s1 <= a_s0 ^ b_s0;
+					z_e1 <= (1 << (M+1)) - 1;
+					z_m1 <= 0;
+					flags1 <= flag_inf;
+				end
+			end
+			// If b is infinity then return infinity
+			else if (b_s0 == (1<<M))
+			begin
+				//Unless a is zero, then return NaN instead
+				if (($signed(a_e0) == (-1*((1<<M)-1))) && a_m0 == 0)
+				begin
+					z_s1 <= 1;
+					z_e1 <= (1 << (M+1)) - 1;
+					z_m1[K-1] <= 1;
+					z_m1[K-2:0] <= 0;
+					flags1 <= flag_nan;
+				end
+				else
+				// Returning infinity
+				begin
+					z_s1 <= a_s0 ^ b_s0;
+					z_e1 <= (1 << (M+1)) - 1;
+					z_m1 <= 0;
+					flags1 <= flag_inf;
+				end
+			end
+			// If either input is zero then return zero
+			else if ((($signed(a_e0) == (-1*((1<<M)-1))) && (a_m0 == 0)) || (($signed(b_e0) == (-1*((1<<M)-1))) && (b_m0 == 0)))
+			begin
+				z_s1 <= a_s0 ^ b_s0;
+				z_e1 <= 0;
+				z_m1 <= 0;
+				flags1 <= flag_zero;
+			end
+			// If b is zero then return zero
+			
+			// If none of the return flags have been set
+			if (flags1 == flag_none)
+			begin
+				// If msb of a_m is 0 then shift left and reduce exponent by 1
+				if (a_m1[K-1] == 0)
+				begin
+					a_m2 <= a_m1 << 1;
+					a_e2 <= a_e1 - 1;
+				end
+				//Just for completion of logic
+				else
+				begin
+					a_m2 <= a_m1;
+					a_e2 <= a_e1;
+				end
+				// If msb of b_m is 0 then shift left and reduce exponent by 1
+				if (b_m1[K-1] == 0)
+				begin
+					b_m2 <= b_m1 << 1;
+					b_e2 <= b_e1 - 1;
+				end
+				//Just for completion of logic
+				else
+				begin
+					b_m2 <= b_m1;
+					b_e2 <= b_e1;
+				end
+			end
+			
+			else
+			begin
+				a_m2 <= a_m1;
+				a_e2 <= a_e1;
+				b_m2 <= b_m1;
+				b_e2 <= b_e1;
+				z_s3 <= z_s2;
+				z_e2 <= z_e1;
+				z_e3 <= z_e2;
+				z_m4 <= z_m3;
+				z_p3 <= 0;
+			end
+		end
+		else
+		begin
+			a_s0 <= 0;
+			a_e0 <= 0;
+			a_m0 <= 0;
+			
+			a_s1 <= 0;
+			a_e1 <= 0;
+			a_m1 <= 0;
+			
+			a_s2 <= 0;
+			a_e2 <= 0;
+			a_m2 <= 0;
+			
+			b_s0 <= 0;
+			b_e0 <= 0;
+			b_m0 <= 0;
+			
+			b_s1 <= 0;
+			b_e1 <= 0;
+			b_m1 <= 0;
+			
+			b_s2 <= 0;
+			b_e2 <= 0;
+			b_m2 <= 0;
+			
+			z_s1 <= 0;
+			z_e1 <= 0;
+			z_m1 <= 0;
+			
+			z_s2 <= 0;
+			z_e2 <= 0;
+			z_m2 <= 0;
+			
+			z_s3 <= 0;
+			z_e3 <= 0;
+			z_m3 <= 0;
+			
+			z_s4 <= 0;
+			z_e4 <= 0;
+			z_m4 <= 0;
+			
+			z_p3 <= 0;
+			
+			flags1 <= flag_none;
+		end
+	end
+
+endmodule : fp_product

src/fpu16/fpu16.sv

@@ -0,0 +1,91 @@
+`include "fp_adder.sv"
+`include "fp_product.sv"
+`include "fp_mult.sv"
+
+
+module fpu16_tb;
+	localparam N=32;
+	localparam M=8;
+
+	reg reset, clk;
+	logic [N-1:0] input_a, input_b, result_add, result_mult;
+	logic [31:0] expected_add, expected_mult;
+
+	fp_adder #(.N(N), .M(M)) adder1(
+		.clk(clk),
+		.reset(reset),
+		.input_a(input_a),
+		.input_b(input_b),
+		.output_z(result_add)
+	);
+
+	fp_product #(.N(N), .M(M)) multiplier1(
+		.clk(clk),
+		.reset(reset),
+		.input_a(input_a),
+		.input_b(input_b),
+		.output_z(result_mult)
+	);
+	
+	initial forever #5 clk = ~clk;
+	localparam PIPELINES_ADD = 2;
+	localparam PIPELINES_MUL = 7;
+
+	reg [N-1:0] test_mem [29:0][3:0];
+
+	initial $readmemh("scripts/fp32_test.hex", test_mem);
+
+	reg [N-2-M:0] exp0_m;
+	reg [M-1:0] exp0_e;
+	reg exp0_s;
+
+	reg [N-2-M:0] exp1_m;
+	reg [M-1:0] exp1_e;
+	reg exp1_s;
+
+	always_comb begin
+		{exp0_s, exp0_e, exp0_m} = expected_add;
+		{exp1_s, exp1_e, exp1_m} = expected_mult;
+	end
+
+	initial begin
+        static int num_err = 0;
+        static int num_tests = $size(test_mem) * 2;
+
+		clk = 0;
+		reset = 1;
+		  
+		#15;
+		reset = 0;
+
+		expected_add = 0;
+		expected_mult = 0;
+
+        for (int i=0; i < $size(test_mem)+PIPELINES_MUL; i++) begin
+			if(i >= PIPELINES_ADD) expected_add = test_mem[i-PIPELINES_ADD][2];
+			if(i >= PIPELINES_MUL) expected_mult = test_mem[i-PIPELINES_MUL][3];
+
+            input_a = test_mem[i][0];
+            input_b = test_mem[i][1];
+
+            #10;
+            if(result_add != expected_add) begin
+                if(num_err < 20)
+                    $display("FAIL ADD: %H + %H = %H, expected %H", test_mem[i-PIPELINES_ADD][0], test_mem[i-PIPELINES_ADD][1], result_add, test_mem[i][2]);
+                num_err = num_err + 1;
+            end
+
+			if(result_mult != expected_mult) begin
+                if(num_err < 20)
+                    $display("FAIL MULTIPLY: %H x %H = %H, expected %H", test_mem[i-PIPELINES_MUL][0], test_mem[i-PIPELINES_MUL][1], result_mult, test_mem[i][3]);
+                num_err = num_err + 1;
+            end
+
+		end
+		expected_add = 0;
+		expected_mult = 0;
+		#50;
+        $display("Passed %d of %d tests", num_tests-num_err, num_tests);
+        $finish();
+	end
+endmodule : fpu16_tb

src/fpu16/lead_zero.sv

@@ -0,0 +1,58 @@
+module count_lead_zero #(
+    parameter W_IN = 64, // Must be power of 2, >=2
+    parameter W_OUT = $clog2(W_IN) // Let this default
+) (
+    input wire  [W_IN-1:0] in,
+    output wire [W_OUT-1:0] out
+);
+
+    generate
+        if (W_IN == 2) begin: base
+            assign out = !in[1];
+        end else begin: recurse
+            wire [W_OUT-2:0] half_count;
+            wire [W_IN / 2-1:0] lhs = in[W_IN / 2 +: W_IN / 2];
+            wire [W_IN / 2-1:0] rhs = in[0        +: W_IN / 2];
+            wire left_empty = ~|lhs;
+
+            count_lead_zero #(
+                .W_IN (W_IN / 2)
+            ) inner (
+                .in  (left_empty ? rhs : lhs),
+                .out (half_count)
+            );
+
+            assign out = {left_empty, half_count};
+        end
+    endgenerate
+
+endmodule : count_lead_zero
+
+
+module count_lead_zero_tb;
+
+    reg [15:0] x0;
+    wire [3:0] y0;
+    count_lead_zero #(16) clz16(x0, y0);
+
+
+    task test_clz;
+        input b;
+        #5;
+        if(y0 != b) $display("FAIL %B: got %d expected %d", x0, y0, b);
+        else $display("OK %B: got %d", x0, y0);
+    endtask : test_clz
+
+    initial begin
+        x0 = 0;
+        test_clz(16);
+        x0 = 1;
+        test_clz(15);
+        for(int i=14; i>=0; i--) begin
+            x0 <<= 1;
+            x0[0] = 1;
+            test_clz(i);
+        end
+    end
+
+endmodule : count_lead_zero_tb