Verilog generator for RNS multipliers based on difference of squares

Choose module (from 2 to 1000):
The main idea of the method is as follows: : X*Y = (1/4)*(X+Y)² - (1/4)*(X-Y)²
The method is not applicable to even moduli
|4^-1|₇ = 2
module multiplication_mod_7(inp1, inp2, out);
	input [2:0] inp1;
	input [2:0] inp2;
	output [2:0] out;
	wire [3:0] plus;
	wire [3:0] minus;
	wire [2:0] plout;
	wire [2:0] miout;
	assign plus = inp1 + inp2; // [0; 12]
	assign minus = 6 + inp1 - inp2; // [0; 12]
	lut_sqr_sum ls1(plus, plout);
	lut_sqr_sub ls2(minus, miout);
	sub_mod_7 sub1(plout, miout, out);
endmodule

module lut_sqr_sum (in, out);
	input [3:0] in;
	output reg [2:0] out;
	always @ (in)
	begin
	case (in)
	0: out = 0;
	1: out = 2;
	2: out = 1;
	3: out = 4;
	4: out = 4;
	5: out = 1;
	6: out = 2;
	7: out = 0;
	8: out = 2;
	9: out = 1;
	10: out = 4;
	11: out = 4;
	12: out = 1;
	default: out = 0;
	endcase
	end
endmodule

module lut_sqr_sub (in, out);
	input [3:0] in;
	output reg [2:0] out;
	always @ (in)
	begin
	case (in)
	0: out = 2;
	1: out = 1;
	2: out = 4;
	3: out = 4;
	4: out = 1;
	5: out = 2;
	6: out = 0;
	7: out = 2;
	8: out = 1;
	9: out = 4;
	10: out = 4;
	11: out = 1;
	12: out = 2;
	default: out = 0;
	endcase
	end
endmodule

module sub_mod_7 (in1, in2, out);
	input [2:0] in1, in2;
	output [2:0] out;
	wire [3:0] data;

	assign data = 7 + in1 - in2;
	mod_7_13 mdval(data, out);
endmodule

module mod_7_13 (in, out);
	input [3:0] in;
	output reg [2:0] out;
	always @ (in)
	begin
	// we have small max value so we can use table here
	case (in)
		0: out = 0;
		1: out = 1;
		2: out = 2;
		3: out = 3;
		4: out = 4;
		5: out = 5;
		6: out = 6;
		7: out = 0;
		8: out = 1;
		9: out = 2;
		10: out = 3;
		11: out = 4;
		12: out = 5;
		13: out = 6;
		default: out = 0;
	endcase
	end
endmodule

module atop_testbench();
	reg [2:0] inp1;
	reg [2:0] inp2;
	wire [2:0] out;

	integer i, j, l, m, k, t;
	reg dummy;
	integer fori, forj;
	multiplication_mod_7 mul1(inp1, inp2, out);
	initial
	begin
	k = 1;
	for (fori = 0; fori < 7; fori = fori + 1)
	begin
	for (forj = 0; forj < 7; forj = forj + 1)
	begin
		inp1 = fori;
		inp2 = forj;
		#1 dummy = 1;
		i = (fori*forj)%7;
		$display ("!!! OP1 = (%d) OP2 = (%d) RES = (%d) EXPECT = (%d)", fori, forj, out, i);
		l = out;
		if (l != i)
		begin
			$display ("!!! Error (%d, %d)!!!", i, l);
		end
		#1 dummy = 1;
	end
	end
	end
endmodule