//
// 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 P_384_ROM[]  = P_384;	// 192-bit
static const FPGA_WORD Q_384_ROM[]  = Q_384;	//
static const FPGA_WORD DP_384_ROM[] = DP_384;	//
static const FPGA_WORD DQ_384_ROM[] = DQ_384;	//
static const FPGA_WORD MP_384_ROM[] = MP_384;	//
static const FPGA_WORD MQ_384_ROM[] = MQ_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;		//

static const FPGA_WORD P_512_ROM[]  = P_512;	// 256-bit
static const FPGA_WORD Q_512_ROM[]  = Q_512;	//
static const FPGA_WORD DP_512_ROM[] = DP_512;	//
static const FPGA_WORD DQ_512_ROM[] = DQ_512;	//
static const FPGA_WORD MP_512_ROM[] = MP_512;	//
static const FPGA_WORD MQ_512_ROM[] = MQ_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);

void modexp_crt			(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);

bool test_modexp_crt	(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 exponentiate 384-bit message with 192-bit prime P and exponent dP...\n\n");
	ok = test_modexp_crt(P_384_ROM, M_384_ROM, DP_384_ROM, MP_384_ROM, OPERAND_NUM_WORDS_384 >> 1);
	if (!ok) return EXIT_FAILURE;

	printf("Trying to exponentiate 384-bit message with 192-bit prime Q and exponent dQ...\n\n");
	ok = test_modexp_crt(Q_384_ROM, M_384_ROM, DQ_384_ROM, MQ_384_ROM, OPERAND_NUM_WORDS_384 >> 1);
	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;

	printf("Trying to exponentiate 512-bit message with 256-bit prime P and exponent dP...\n\n");
	ok = test_modexp_crt(P_512_ROM, M_512_ROM, DP_512_ROM, MP_512_ROM, OPERAND_NUM_WORDS_512 >> 1);
	if (!ok) return EXIT_FAILURE;

	printf("Trying to exponentiate 512-bit message with 256-bit prime Q and exponent dQ...\n\n");
	ok = test_modexp_crt(Q_512_ROM, M_512_ROM, DQ_512_ROM, MQ_512_ROM, OPERAND_NUM_WORDS_512 >> 1);
	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, false);

		/*
		 * 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);
}


//----------------------------------------------------------------
// Modular exponentiation routine with CRT support
//----------------------------------------------------------------
void modexp_crt(	const FPGA_WORD *M,
					const FPGA_WORD *D,
					const FPGA_WORD *N,
					      FPGA_WORD *R,
						  size_t     len)
//----------------------------------------------------------------
//
// R = (A mod N) ** B mod N
//
//----------------------------------------------------------------
{
		// temporary buffers
	FPGA_WORD M0     [MAX_OPERAND_WORDS];
	FPGA_WORD M1     [MAX_OPERAND_WORDS];
	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);
	
		// reduce M to make it smaller than N
	montgomery_multiply(M, NULL, N, N_COEFF, M0, len, true);

		// bring M into Montgomery domain
	montgomery_multiply(M0, FACTOR, N, N_COEFF, M1,       len, false);
	montgomery_multiply(M1, FACTOR, N, N_COEFF, M_FACTOR, len, false);

		/*
		 * Montgomery multiplication adds an extra factor of 2 ^ -w to every product,
		 * Montgomery reduction adds that factor too. The message must be reduced before
		 * exponentiation, because in CRT mode it is twice larger, than the modulus
		 * and the exponent. After reduction the message carries an extra factor of
		 * 2 ^ -w. We pre-calculate a special factor of 2 ^ 2w and multiply the message
		 * by this factor *twice* using our Montgomery multiplier. This way we get the
		 * message with an extra factor of just 2 ^ w:
		 * 1. (m * 2 ^ -w) * (2 ^ 2w) * (2 ^ -w) = m
		 * 2. (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 exponentiation 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 model to
// calculate r = m ** d mod n, 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("    ERROR\n\n\n");
		return false;
	}

		// everything went just fine
	printf("    OK\n\n\n");
	return true;
}


//----------------------------------------------------------------
// Test the modular exponentiation model with CRT enabled
//----------------------------------------------------------------
bool test_modexp_crt(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 model to
// calculate r = (m mod n) ** d mod n, and then compares it to the
// reference value s. The difference from test_modexp() is that
// m_rom is twice larger than n_rom and d_rom.
//
//----------------------------------------------------------------
{
	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 (message is twice as large!)
	load_value_from_rom(n_rom, N, len);
	load_value_from_rom(m_rom, M, len << 1);
	load_value_from_rom(d_rom, D, len);
	load_value_from_rom(s_rom, S, len);

		// calculate power
	modexp_crt(M, D, N, R, len);

		// check result
	ok = compare_fpga_buffers(S, R, len);
	if (!ok)
	{	printf("    ERROR\n\n\n");
		return false;
	}

		// everything went just fine
	printf("    OK\n\n\n");
	return true;
}


//----------------------------------------------------------------
// End of file
//----------------------------------------------------------------