//
// modexp_fpga_model_montgomery.cpp
// -------------------------------------------------------------
// Montgomery modular multiplication and exponentiation routines
//
// 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 "modexp_fpga_model.h"
#include "modexp_fpga_model_pe.h"
#include "modexp_fpga_model_montgomery.h"


//----------------------------------------------------------------
// 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)
//----------------------------------------------------------------
//
// R = A * B * 2^-len mod N
//
// The high-level algorithm is:
//
// 1. AB =  A * B
// 2. Q  = AB * N_COEFF
// 3. QN =  Q * N
// 4. S  = AB + QN
// 5. SN =  S - N
// 6. R  = (SN < 0) ? S : SN
// 7. R  = R >> len
//
//----------------------------------------------------------------
{
	size_t i, j, k;													// counters

	bool select_s;													// flag

	FPGA_WORD t_ab[MAX_SYSTOLIC_CYCLES][SYSTOLIC_NUM_WORDS];		// accumulators
	FPGA_WORD t_q [MAX_SYSTOLIC_CYCLES][SYSTOLIC_NUM_WORDS];		//
	FPGA_WORD t_qn[MAX_SYSTOLIC_CYCLES][SYSTOLIC_NUM_WORDS];		//

	FPGA_WORD s_ab[MAX_SYSTOLIC_CYCLES][SYSTOLIC_NUM_WORDS];		// intermediate products
	FPGA_WORD s_q [MAX_SYSTOLIC_CYCLES][SYSTOLIC_NUM_WORDS];		//
	FPGA_WORD s_qn[MAX_SYSTOLIC_CYCLES][SYSTOLIC_NUM_WORDS];		//

	FPGA_WORD c_in_ab[MAX_SYSTOLIC_CYCLES][SYSTOLIC_NUM_WORDS];		// input carries
	FPGA_WORD c_in_q [MAX_SYSTOLIC_CYCLES][SYSTOLIC_NUM_WORDS];		//
	FPGA_WORD c_in_qn[MAX_SYSTOLIC_CYCLES][SYSTOLIC_NUM_WORDS];		//
	FPGA_WORD c_out_ab[MAX_SYSTOLIC_CYCLES][SYSTOLIC_NUM_WORDS];	// output carries
	FPGA_WORD c_out_q [MAX_SYSTOLIC_CYCLES][SYSTOLIC_NUM_WORDS];	//
	FPGA_WORD c_out_qn[MAX_SYSTOLIC_CYCLES][SYSTOLIC_NUM_WORDS];	//

	FPGA_WORD c_in_s;												// 1-bit carry and borrow
	FPGA_WORD b_in_sn;												//
	FPGA_WORD c_out_s;												//
	FPGA_WORD b_out_sn;												//

	FPGA_WORD AB[2 * MAX_OPERAND_WORDS];							// final products
	FPGA_WORD Q [2 * MAX_OPERAND_WORDS];							//
	FPGA_WORD QN[2 * MAX_OPERAND_WORDS];							//

	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
	size_t num_systolic_cycles = len / SYSTOLIC_NUM_WORDS;

		// initialize arrays of accumulators and carries to zeroes
	for (i=0; i<num_systolic_cycles; i++)
		for (j=0; j<SYSTOLIC_NUM_WORDS; j++)
			c_in_ab[i][j] = 0, c_in_q [i][j] = 0, c_in_qn[i][j] = 0,
			t_ab[i][j]    = 0, t_q [i][j]    = 0, t_qn[i][j]    = 0;

		// initialize 1-bit carry and borrow to zeroes too
	c_in_s = 0, b_in_sn = 0;

		// simultaneously calculate AB, Q, QN, S, SN
	for (i = 0; i < (2 * len); i++)
	{
			// multiply entire B by current word of A to get AB
			// multiply entire N_COEFF by current word of AB to get Q
			// multiply entire N by current word of Q to get QN
		for (k = 0; k < num_systolic_cycles; k++)
		{
				// simulate how a systolic array would work
			for (j = 0; j < 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];

					// current word of A
				FPGA_WORD Aj_ab = (i < len) ? A[i] : 0;

					// AB = A * B		
				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];

					// current word of AB
				FPGA_WORD Aj_q = (i < len) ? AB[i] : 0;

					// Q = AB * N		
				pe_mul(Aj_q, N_COEFFj, t_q[k][j], c_in_q[k][j], &s_q[k][j], &c_out_q[k][j]);

					// store current word of Q
				if ((k == 0) && (j == 0)) Q[i] = s_q[0][0];

					// current word of Q
				FPGA_WORD Aj_qn = (i < len) ? Q[i] : 0;

					// QN = Q * N
				pe_mul(Aj_qn, Nj, t_qn[k][j], c_in_qn[k][j], &s_qn[k][j], &c_out_qn[k][j]);

					// store next word of QN
				if ((k == 0) && (j == 0)) QN[i] = s_qn[0][0];
			}

				// propagate carries
			for (j=0; j<SYSTOLIC_NUM_WORDS; j++)
				c_in_ab[k][j] = c_out_ab[k][j],
				c_in_q [k][j] = c_out_q [k][j],
				c_in_qn[k][j] = c_out_qn[k][j];

				// update accumulators
			for (j=1; j<SYSTOLIC_NUM_WORDS; j++)
			{
				t_ab[k][j-1] = s_ab[k][j];
				t_q [k][j-1] = s_q [k][j];
				t_qn[k][j-1] = s_qn[k][j];
			}
			
				// update accumulators
			if (k > 0)
				t_ab[k-1][SYSTOLIC_NUM_WORDS-1] = s_ab[k][0],
				t_q [k-1][SYSTOLIC_NUM_WORDS-1] = s_q [k][0],
				t_qn[k-1][SYSTOLIC_NUM_WORDS-1] = s_qn[k][0];

		}
	
			// now it's time to simultaneously add and subtract

			// current operand words
		FPGA_WORD QNi = QN[i];
		FPGA_WORD Ni  = (i < len) ? 0 : N[i-len];

			// add, subtract
		pe_add(AB[i], QNi, c_in_s,  &S[i],  &c_out_s);
		pe_sub(S [i], Ni,  b_in_sn, &SN[i], &b_out_sn);

			// propagate carry and borrow
		c_in_s  = c_out_s;
		b_in_sn = b_out_sn;
	}

		// flag select the right result
	select_s = b_out_sn && !c_out_s;

		// copy product into output buffer
	for (i=0; i<len; i++)
		R[i] = select_s ? S[i+len] : SN[i+len];
}


//----------------------------------------------------------------
// Classic binary exponentiation
//----------------------------------------------------------------
void montgomery_exponentiate(const FPGA_WORD *A, const FPGA_WORD *B, const FPGA_WORD *N, const FPGA_WORD *N_COEFF, FPGA_WORD *R, size_t len)
//----------------------------------------------------------------
//
// R = A ** B mod N
//
//----------------------------------------------------------------
{
	size_t word_cnt,   bit_cnt;			// counters
	size_t word_index, bit_index;		// indices

	bool flag_update_r;					// flag

	FPGA_WORD P[MAX_OPERAND_WORDS];		// power of A
	FPGA_WORD mask;						// mask		

		// R = 1, P = 1
	for (word_cnt=0; word_cnt<len; word_cnt++)
		R[word_cnt] = (word_cnt > 0) ? 0 : 1,
		P[word_cnt] = A[word_cnt];

	FPGA_WORD M_PP[MAX_OPERAND_WORDS];	// intermediate buffer for next power
	FPGA_WORD M_RP[MAX_OPERAND_WORDS];	// intermediate buffer for next result

		// 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
		
		word_index = bit_cnt / (CHAR_BIT * sizeof(FPGA_WORD));
		bit_index = bit_cnt & ((CHAR_BIT * sizeof(FPGA_WORD)) - 1);

		mask = 1 << bit_index;	// current bit of exponent

			// whether we need to update R (non-zero exponent bit)
		flag_update_r = (B[word_index] & mask) == mask;

			// always update P
		for (word_cnt=0; word_cnt<len; word_cnt++)
			P[word_cnt] = M_PP[word_cnt];

			// only update R when necessary
		if (flag_update_r)
		{
			for (word_cnt=0; word_cnt<len; word_cnt++)
				R[word_cnt] = M_RP[word_cnt];
		}
	}
}


//----------------------------------------------------------------
// Montgomery factor calculation
//----------------------------------------------------------------
void montgomery_calc_factor(const FPGA_WORD *N, FPGA_WORD *FACTOR, size_t len)
//----------------------------------------------------------------
//
// FACTOR = 2 ** (2*len) mod N
//
// This routine calculates the factor that is necessary to bring
// numbers into Montgomery domain. The high-level algorithm is:
//
// 1. f = 1
// 2. for i=0 to 2*len-1
// 3.   f1 = f << 1
// 4.   f2 = f1 - n
// 5.   f = (f2 < 0) ? f1 : f2
//
//----------------------------------------------------------------
{
	size_t i, j;		// counters

	bool flag_keep_f;	// flag
	
		// temporary buffer
	FPGA_WORD FACTOR_N[MAX_OPERAND_WORDS];

		// carry and borrow
	FPGA_WORD carry_in, carry_out;
	FPGA_WORD borrow_in, borrow_out;

		// FACTOR = 1
	for (i=0; i<len; i++)
		FACTOR[i] = (i > 0) ? 0 : 1;
		
		// do the math
	for (i=0; i<(2 * len * CHAR_BIT * sizeof(FPGA_WORD)); i++)
	{
			// clear carry and borrow
		carry_in = 0, borrow_in = 0;
		
			// calculate f1 = f << 1, f2 = f1 - n
		for (j=0; j<len; j++)
		{
			carry_out = FACTOR[j] >> (sizeof(FPGA_WORD) * CHAR_BIT - 1);		// | M <<= 1
			FACTOR[j] <<= 1, FACTOR[j] |= carry_in;								// |

			pe_sub(FACTOR[j], N[j], borrow_in, &FACTOR_N[j], &borrow_out);		// MN = M - N
	
			carry_in = carry_out, borrow_in = borrow_out;						// propagate carry & borrow
		}

			// obtain flag
		flag_keep_f = (borrow_out && !carry_out);

			// now select the right value
		for (j=0; j<len; j++)
			FACTOR[j] = flag_keep_f ? FACTOR[j] : FACTOR_N[j];
	}
}


//----------------------------------------------------------------
// Montgomery modulus-dependent coefficient calculation
//----------------------------------------------------------------
void montgomery_calc_n_coeff(const FPGA_WORD *N, FPGA_WORD *N_COEFF, size_t len)
//----------------------------------------------------------------
//
// N_COEFF = -N ** -1 mod 2 ** len
//
// This routine calculates the coefficient that is used during the reduction
// phase of Montgomery multiplication to zero out the lower half of product.
//
// The high-level algorithm is:
//
// 1. R = 1
// 2. NN = ~N + 1
// 3. for i=1 to len-1
// 4.   T = R * NN mod 2 ** len
// 5.   if T[i] then
// 6.     R = R + (1 << i)
//
//----------------------------------------------------------------
{
	size_t i, j, k;							// counters

	FPGA_WORD NN[MAX_OPERAND_WORDS];		// temporary buffers
	FPGA_WORD T [MAX_OPERAND_WORDS];		//
	FPGA_WORD R [MAX_OPERAND_WORDS];		//
	FPGA_WORD R1[MAX_OPERAND_WORDS];		//

	bool flag_update_r;						// flag

	FPGA_WORD nw, pwr;						//
	FPGA_WORD sum_c_in, sum_c_out;			//
	FPGA_WORD carry_in, carry_out;			//
	FPGA_WORD mul_s, mul_c_in, mul_c_out;	//

		// NN = -N mod 2 ** len = ~N + 1 mod 2 ** len
	carry_in = 0;
	for (i=0; i<len; i++)
	{	nw = (i > 0) ? 0 : 1;								// nw = 1
		pe_add(~N[i], nw, carry_in, &NN[i], &carry_out);	// NN = ~N + nw
		carry_in = carry_out;								// propagate carry
	}

		// R = 1
	for (i=0; i<len; i++)
		R[i] = (i > 0) ? 0 : 1;

		// calculate T = R * NN
	for (k=1; k<(len * sizeof(FPGA_WORD) * CHAR_BIT); k++)
	{
			// T = 0
		for (i=0; i<len; i++) T[i] = 0;

			// T = NN * R
		for (i=0; i<len; i++)
		{
			mul_c_in = 0;

			for (j=0; j<(len-i); j++)
			{
			
				pe_mul(R[j], NN[i], T[i+j], mul_c_in, &mul_s, &mul_c_out);
				
				T[i+j] = mul_s;
				
				mul_c_in = mul_c_out;
			}
		}

			// get word and index indices
		size_t word_index = k / (CHAR_BIT * sizeof(FPGA_WORD));
		size_t bit_index = k & ((CHAR_BIT * sizeof(FPGA_WORD)) - 1);

			// update bit mask
		FPGA_WORD bit_mask = (1 << bit_index);

		sum_c_in = 0;			// clear carry
		flag_update_r = false;	// reset flag

			// calculate R1 = R + 1 << (2 * len)
		for (i=0; i<len; i++)
		{
			if (i == word_index) flag_update_r = (T[i] & bit_mask) == bit_mask;

			pwr = (i == word_index) ? bit_mask : 0;
			pe_add(R[i], pwr, sum_c_in,  &R1[i], &sum_c_out);
			carry_in = carry_out;
		}

			// update r
		for (i=0; i<len; i++)
			R[i] = flag_update_r ? R1[i] : R[i];
	}

		// store output
	for (i=0; i<len; i++)
		N_COEFF[i] = R[i];
}


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