Improved the model:

 * added CRT support
 * fixed bug in systolic array when operand width is not a multiple of array width

1 files changed, 14 insertions, 8 deletions

diff --git a/modexp_fpga_model_montgomery.cpp b/modexp_fpga_model_montgomery.cpp
index d1cca60..34ef2b6 100644
--- a/modexp_fpga_model_montgomery.cpp
+++ b/modexp_fpga_model_montgomery.cpp
@@ -46,7 +46,7 @@
 //----------------------------------------------------------------

 // Montgomery modular multiplier

 //----------------------------------------------------------------

-void montgomery_multiply(const FPGA_WORD *A, const FPGA_WORD *B, const FPGA_WORD *N, const FPGA_WORD *N_COEFF, FPGA_WORD *R, size_t len)

+void montgomery_multiply(const FPGA_WORD *A, const FPGA_WORD *B, const FPGA_WORD *N, const FPGA_WORD *N_COEFF, FPGA_WORD *R, size_t len, bool reduce_only)

 //----------------------------------------------------------------

 //

 // R = A * B * 2^-len mod N

@@ -94,8 +94,11 @@ void montgomery_multiply(const FPGA_WORD *A, const FPGA_WORD *B, const FPGA_WORD
 	FPGA_WORD S [2 * MAX_OPERAND_WORDS];							// final sum

 	FPGA_WORD SN[2 * MAX_OPERAND_WORDS];							// final difference

 

-		// number of systolic cycles needed to multiply entire B by one word of A

+		// number of full systolic cycles needed to multiply entire B by one word of A

 	size_t num_systolic_cycles = len / SYSTOLIC_NUM_WORDS;

+	

+		// adjust number of cycles

+	if ((num_systolic_cycles * SYSTOLIC_NUM_WORDS) < len) num_systolic_cycles++;

 

 		// initialize arrays of accumulators and carries to zeroes

 	for (i=0; i<num_systolic_cycles; i++)

@@ -117,10 +120,12 @@ void montgomery_multiply(const FPGA_WORD *A, const FPGA_WORD *B, const FPGA_WORD
 				// simulate how a systolic array would work

 			for (j = 0; j < SYSTOLIC_NUM_WORDS; j++)

 			{

+				size_t j_index = k * SYSTOLIC_NUM_WORDS + j;

+

 					// current words of B, N_COEFF, N

-				FPGA_WORD Bj       = B      [k * SYSTOLIC_NUM_WORDS + j];

-				FPGA_WORD N_COEFFj = N_COEFF[k * SYSTOLIC_NUM_WORDS + j];

-				FPGA_WORD Nj       = N      [k * SYSTOLIC_NUM_WORDS + j];

+				FPGA_WORD Bj       = (j_index < len) ? B      [k * SYSTOLIC_NUM_WORDS + j] : 0;

+				FPGA_WORD N_COEFFj = (j_index < len) ? N_COEFF[k * SYSTOLIC_NUM_WORDS + j] : 0;

+				FPGA_WORD Nj       = (j_index < len) ? N      [k * SYSTOLIC_NUM_WORDS + j] : 0;

 

 					// current word of A

 				FPGA_WORD Aj_ab = (i < len) ? A[i] : 0;

@@ -129,7 +134,7 @@ void montgomery_multiply(const FPGA_WORD *A, const FPGA_WORD *B, const FPGA_WORD
 				pe_mul(Aj_ab, Bj, t_ab[k][j], c_in_ab[k][j], &s_ab[k][j], &c_out_ab[k][j]);

 

 					// store current word of AB

-				if ((k == 0) && (j == 0)) AB[i] = s_ab[0][0];

+				if ((k == 0) && (j == 0)) AB[i] = reduce_only ? A[i] : s_ab[0][0];

 

 					// current word of AB

 				FPGA_WORD Aj_q = (i < len) ? AB[i] : 0;

@@ -225,8 +230,8 @@ void montgomery_exponentiate(const FPGA_WORD *A, const FPGA_WORD *B, const FPGA_
 		// scan all bits of the exponent

 	for (bit_cnt=0; bit_cnt<(len * CHAR_BIT * sizeof(FPGA_WORD)); bit_cnt++)

 	{

-		montgomery_multiply(P, P, N, N_COEFF, M_PP, len);	// M_PP = P * P

-		montgomery_multiply(R, P, N, N_COEFF, M_RP, len);	// M_RP = R * P

+		montgomery_multiply(P, P, N, N_COEFF, M_PP, len, false);	// M_PP = P * P

+		montgomery_multiply(R, P, N, N_COEFF, M_RP, len, false);	// M_RP = R * P

 		

 		word_index = bit_cnt / (CHAR_BIT * sizeof(FPGA_WORD));

 		bit_index = bit_cnt & ((CHAR_BIT * sizeof(FPGA_WORD)) - 1);

@@ -308,6 +313,7 @@ void montgomery_calc_factor(const FPGA_WORD *N, FPGA_WORD *FACTOR, size_t len)
 		for (j=0; j<len; j++)

 			FACTOR[j] = flag_keep_f ? FACTOR[j] : FACTOR_N[j];

 	}

+

 }