Initial commit of faster modular exponentiation model based on systolic architecture.

1 files changed, 293 insertions, 0 deletions

diff --git a/modexp_fpga_model.cpp b/modexp_fpga_model.cpp
new file mode 100644
index 0000000..b19cdeb
--- /dev/null
+++ b/modexp_fpga_model.cpp
@@ -0,0 +1,293 @@
+//

+// modexp_fpga_model.cpp

+// -----------------------------------------------

+// Model of fast modular exponentiation on an FPGA

+//

+// Authors: Pavel Shatov

+// Copyright (c) 2017, NORDUnet A/S

+// All rights reserved.

+//

+// Redistribution and use in source and binary forms, with or without

+// modification, are permitted provided that the following conditions are

+// met:

+// - Redistributions of source code must retain the above copyright notice,

+//   this list of conditions and the following disclaimer.

+//

+// - Redistributions in binary form must reproduce the above copyright

+//   notice, this list of conditions and the following disclaimer in the

+//   documentation and/or other materials provided with the distribution.

+//

+// - Neither the name of the NORDUnet nor the names of its contributors may

+//   be used to endorse or promote products derived from this software

+//   without specific prior written permission.

+//

+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS

+// IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED

+// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A

+// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT

+// HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,

+// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED

+// TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR

+// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF

+// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING

+// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS

+// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

+//

+

+

+//----------------------------------------------------------------

+// Headers

+//----------------------------------------------------------------

+#include <stdio.h>

+#include <stdlib.h>

+#include <string.h>

+#include "modexp_fpga_model.h"

+#include "modexp_fpga_model_montgomery.h"

+#include "test/modexp_fpga_model_vectors.h"

+

+

+//----------------------------------------------------------------

+// Defined values

+//----------------------------------------------------------------

+#define OPERAND_WIDTH_384	384

+#define OPERAND_WIDTH_512	512

+

+#define OPERAND_NUM_WORDS_384	(OPERAND_WIDTH_384 / (CHAR_BIT * sizeof(FPGA_WORD)))

+#define OPERAND_NUM_WORDS_512	(OPERAND_WIDTH_512 / (CHAR_BIT * sizeof(FPGA_WORD)))

+

+

+//----------------------------------------------------------------

+// Test vectors

+//----------------------------------------------------------------

+static const FPGA_WORD N_384_ROM[] = N_384;		// 384-bit

+static const FPGA_WORD M_384_ROM[] = M_384;		//

+static const FPGA_WORD D_384_ROM[] = D_384;		//

+static const FPGA_WORD S_384_ROM[] = S_384;		//

+

+static const FPGA_WORD N_512_ROM[] = N_512;		// 512-bit

+static const FPGA_WORD M_512_ROM[] = M_512;		//

+static const FPGA_WORD D_512_ROM[] = D_512;		//

+static const FPGA_WORD S_512_ROM[] = S_512;		//

+

+

+//----------------------------------------------------------------

+// Prototypes

+//----------------------------------------------------------------

+void print_fpga_buffer		(const char *str, const FPGA_WORD *buf, size_t len);

+bool compare_fpga_buffers	(const FPGA_WORD *src, const FPGA_WORD *dst, size_t len);

+void load_value_from_rom	(const FPGA_WORD *src, FPGA_WORD *dst, size_t len);

+

+void modexp			(const FPGA_WORD *M, const FPGA_WORD *D,

+					 const FPGA_WORD *N,       FPGA_WORD *R, size_t len);

+

+bool test_modexp	(const FPGA_WORD *n_rom, const FPGA_WORD *m_rom,

+					 const FPGA_WORD *d_rom, const FPGA_WORD *s_rom, size_t len);

+

+

+//----------------------------------------------------------------

+int main()

+//----------------------------------------------------------------

+{

+	bool ok;

+

+	printf("Trying to sign 384-bit message...\n\n");

+	ok = test_modexp(N_384_ROM, M_384_ROM, D_384_ROM, S_384_ROM, OPERAND_NUM_WORDS_384);

+	if (!ok) return EXIT_FAILURE;

+

+	printf("Trying to sign 512-bit message...\n\n");

+	ok = test_modexp(N_512_ROM, M_512_ROM, D_512_ROM, S_512_ROM, OPERAND_NUM_WORDS_512);

+	if (!ok) return EXIT_FAILURE;

+

+	return EXIT_SUCCESS;

+}

+

+

+//----------------------------------------------------------------

+// Modular exponentiation routine

+//----------------------------------------------------------------

+void modexp(	const FPGA_WORD *M,

+				const FPGA_WORD *D,

+				const FPGA_WORD *N,

+				      FPGA_WORD *R,

+					  size_t     len)

+//----------------------------------------------------------------

+//

+// R = A ** B mod N

+//

+//----------------------------------------------------------------

+{

+		// temporary buffers

+	FPGA_WORD FACTOR [MAX_OPERAND_WORDS];

+	FPGA_WORD N_COEFF[MAX_OPERAND_WORDS];

+	FPGA_WORD M_FACTOR[MAX_OPERAND_WORDS];

+

+		// pre-calculate modulus-dependant coefficients

+	montgomery_calc_factor(N, FACTOR, len);

+	montgomery_calc_n_coeff(N, N_COEFF, len);

+		

+		// bring M into Montgomery domain

+	montgomery_multiply(M, FACTOR, N, N_COEFF, M_FACTOR, len);

+

+		/*

+		 * Montgomery multiplication adds an extra factor of 2 ^ -w to every product.

+		 * We pre-calculate a special factor of 2 ^ 2w and multiply the message

+		 * by this factor using our Montgomery multiplier. This way we get the message

+		 * with the an extra factor of just 2 ^ w:

+		 * (m) * (2 ^ 2w) * (2 ^ -w) = m * 2 ^ w 

+		 *

+		 * Now we feed this message with that extra factor to the binary exponentiation

+		 * routine. The current power of m will always keep that additional factor:

+		 * (p * 2 ^ w) * (p * 2 ^ w) * (2 ^ -w) = p ^ 2 * 2 ^ w

+		 *

+		 * The result starts at 1, i.e. without any extra factors. If at any particular

+		 * iteration it gets multiplied with the current power of m, the product will

+		 * not carry any extra factors, because the power's factor gets eliminated

+		 * by the extra factor of Montgomery multiplication:

+		 * (r) * (p * 2 ^ w) * (2 ^ -w) = r * p

+		 *

+		 * This way we don't need any extra post-processing to convert the final result

+		 * from Montgomery domain. 

+		 *

+		 */

+

+		// exponentiate

+	montgomery_exponentiate(M_FACTOR, D, N, N_COEFF, R, len);

+}

+

+

+//----------------------------------------------------------------

+// Copies words from src into dst reversing their order

+//----------------------------------------------------------------

+void load_value_from_rom(const FPGA_WORD *src, FPGA_WORD *dst, size_t len)

+//----------------------------------------------------------------

+//

+// This routine copies src into dst word-by-word reversing their order

+// in the process. This reversal is necessary because of the way C

+// arrays are initialized. This model requires the least significant

+// word of operand to be stored at array offset 0, while C places

+// the most significant word there instead.

+//

+// Suppose that the operand is 0xFEDCBA9876543210, now the following line

+// uint32_t X[2] = {0xFEDCBA98, 0x76543210}

+// will place the most significant word 0xFEDCBA98 at index [0].

+//

+//----------------------------------------------------------------

+{

+	size_t w;

+

+	for (w=0; w<len; w++)

+		dst[w] = src[len - (w + 1)];

+}

+

+

+//----------------------------------------------------------------

+// Compare two operands

+//----------------------------------------------------------------

+bool compare_fpga_buffers(const FPGA_WORD *src, const FPGA_WORD *dst, size_t len)

+//----------------------------------------------------------------

+//

+// This routine compares two multi-word intergers, it is used to compare

+// the calculated value against the reference one.

+//

+//----------------------------------------------------------------

+{

+	size_t w;	// word counter

+

+		// print all the values

+	print_fpga_buffer("  Expected:   M = ", src, len);

+	print_fpga_buffer("  Calculated: R = ", dst, len);

+

+		// compare values

+	for (w=0; w<len; w++)

+	{

+			// compare

+		if (src[w] != dst[w]) return false;

+	}

+

+		// values are the same

+	return true;

+}

+

+

+//----------------------------------------------------------------

+// Prints large multi-word integer

+//----------------------------------------------------------------

+void print_fpga_buffer(const char *str, const FPGA_WORD *buf, size_t len)

+//----------------------------------------------------------------

+{

+	size_t w, s;	// word counter

+

+		// print header

+	printf("%s", str);

+

+		// print all bytes

+	for (w=len; w>0; w--)

+	{	

+		printf("%08x", buf[w-1]);

+

+			// insert space after every group of 4 bytes or print new line

+		if (w > 1)

+		{

+			if (((len - w) % 4) == 3)

+			{	printf("\n");

+				s = strlen(str);

+				while (s)

+				{	printf(" ");

+					s--;

+				}

+			}

+			else printf(" ");

+		}

+	}

+

+		// print footer

+	printf("\n\n");

+}

+

+

+//----------------------------------------------------------------

+// Test the modular multiplication model

+//----------------------------------------------------------------

+bool test_modexp(const FPGA_WORD *n_rom, const FPGA_WORD *m_rom, const FPGA_WORD *d_rom, const FPGA_WORD *s_rom, size_t len)

+//----------------------------------------------------------------

+//

+// This routine uses the Montgomery exponentiation routine to

+// calculate r = m ** d mod m, and then compares it to the

+// reference value s.

+//

+//----------------------------------------------------------------

+{

+	bool ok;	// flag

+

+		// buffers

+	FPGA_WORD N[MAX_OPERAND_WORDS];

+	FPGA_WORD M[MAX_OPERAND_WORDS];

+	FPGA_WORD D[MAX_OPERAND_WORDS];

+	FPGA_WORD S[MAX_OPERAND_WORDS];

+	FPGA_WORD R[MAX_OPERAND_WORDS];

+

+		// fill buffers with test vector

+	load_value_from_rom(n_rom, N, len);

+	load_value_from_rom(m_rom, M, len);

+	load_value_from_rom(d_rom, D, len);

+	load_value_from_rom(s_rom, S, len);

+

+		// calculate power

+	modexp(M, D, N, R, len);

+

+		// check result

+	ok = compare_fpga_buffers(S, R, len);

+	if (!ok)

+	{	printf("\n    ERROR\n\n");

+		return false;

+	}

+

+		// everything went just fine

+	printf("\n    OK\n\n");

+	return true;

+}

+

+

+//----------------------------------------------------------------

+// End of file

+//----------------------------------------------------------------