//------------------------------------------------------------------------------
//
// modexpng_util.c
// ------------------------------------------------------
// Helper precomputation routines for the "modexpng" core
//
// Authors: Pavel Shatov
//
// Copyright (c) 2019, NORDUnet A/S
//
// 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 "modexpng_util.h"


//
// internal buffers
//
static uint32_t MOD_FACTOR_N[BUF_NUM_WORDS];
static uint32_t MOD_NN[BUF_NUM_WORDS+1];
static uint32_t MOD_T[BUF_NUM_WORDS+1];

static void _add32(uint32_t, uint32_t, uint32_t, uint32_t *, uint32_t *);
static void _sub32(uint32_t, uint32_t, uint32_t, uint32_t *, uint32_t *);
static void _mul32(uint32_t, uint32_t, uint32_t, uint32_t, uint32_t *, uint32_t *);


//
// calculation of the Montgomery factor
//
void _calc_montgomery_factor(uint32_t num_words, const uint32_t *N, uint32_t *N_FACTOR)
{
		// counters
		uint32_t i, j;		

		// flag
		uint32_t flag_keep;

		// carry and borrow
		uint32_t cry_in, cry_out;
		uint32_t brw_in, brw_out;

		// initially set N_FACTOR = 1
		for (i=0; i<num_words; i++)
				N_FACTOR[i] = i ? 0 : 1;
		
		// do the math
		for (i=0; i<2*(num_words * UINT32_BITS + UINT16_BITS); i++)
		{	
				// clear carry and borrow
				cry_in = 0, brw_in = 0;
		
				// calculate N_FACTOR = N_FACTOR << 1, MOD_FACTOR_N = N_FACTOR - N
				for (j=0; j<num_words; j++)
				{
						cry_out = N_FACTOR[j] >> (UINT32_BITS - 1);		// | N_FACTOR <<= 1
						N_FACTOR[j] <<= 1; N_FACTOR[j] |= cry_in;			// |

						_sub32(N_FACTOR[j], N[j], brw_in, &MOD_FACTOR_N[j], &brw_out);	// MOD_FACTOR_N = N_FACTOR - N
	
						// propagate carry & borrow
						cry_in = cry_out, brw_in = brw_out;
				}

				// obtain flag
				flag_keep = brw_out && !cry_out;

				// now select the right value
				for (j=0; j<num_words; j++)
						N_FACTOR[j] = flag_keep ? N_FACTOR[j] : MOD_FACTOR_N[j];
		}
}


//
// calculation of the modulus-dependent speed-up coefficient
//
void _calc_modulus_coeff(uint32_t num_words, const uint32_t *N, uint32_t *N_COEFF)
{
		// counters
		uint32_t i, j, k, jk;

		// indices
		uint32_t word_index, bit_index;

		// flag
		uint32_t flag_update;
	
		// carries
		uint32_t cry_in, cry_out;
	
		// temporary variables
		uint32_t mod_p, add_s, b_word;
	
		// initially set N_COEFF to 1
		for (i=0; i<=num_words; i++)
				N_COEFF[i] = i ? 0 : 1;

		// also set NN to ~N+1
		// note that since N must be odd, ~N is even, so adding 1 to it doesn't need
		// any carry propagation
		for (i=0; i<num_words; i++) MOD_NN[i] = ~N[i];
		MOD_NN[0] += 1;
		MOD_NN[num_words] = 0xffffffff;
        
		// do the math
		for (i=1; i<(num_words * UINT32_BITS + UINT16_BITS); i++)
		{
				word_index = i / UINT32_BITS;
			  bit_index = i & (UINT32_BITS - 1);

				// clear T
				for (j=0; j<=num_words; j++) MOD_T[j] = 0;
				
				// T = N_COEFF * NN mod 2 ** (modulus_length + 16)
				/*
			   * Note, that we only need the lower half of the product T, so in
			   * the outer loop we always scan entire N_COEFF, but the inner
				 * loop only scans entire NN during the first iteration, and then
			   * keeps skipping one more word every iteration, during the last
				 * iteration we only scan one word of NN.
			   *
			   */
				for (j=0; j<=num_words; j++)
				{   cry_in = 0;
						for (k=0; k<=(num_words-j); k++)
						{		jk = j + k;
								_mul32(N_COEFF[j], MOD_NN[k], MOD_T[jk], cry_in, &mod_p, &cry_out);
								MOD_T[jk] = mod_p;
								cry_in = cry_out;
								if (word_index == jk)
										flag_update = MOD_T[jk] & (1 << bit_index) ? 1 : 0;
						}
				}
				if (flag_update)
				{		cry_in = 0;
						for (j=0; j<=num_words; j++)
						{		b_word = (j == word_index) ? (1 << bit_index) : 0;
								_add32(b_word, N_COEFF[j], cry_in, &add_s, &cry_out);
								N_COEFF[j] = add_s;
								cry_in = cry_out;
						}
				}
		}
}


//
// low-level addition w/ carry
//
static void _add32(uint32_t a, uint32_t b, uint32_t c_in, uint32_t *s, uint32_t *c_out)
{
		uint64_t t;															// intermediate var
	
		t = (uint64_t)a + (uint64_t)b;					// obtain "wide" difference
		t += (uint64_t)(c_in & 1);							// take borrow into account
	
		*s = (uint32_t)t;												// return the lower part of result
		*c_out = (uint32_t)(t >> UINT32_BITS);	// return the higher part of result, ...
		*c_out &= (uint32_t)1;									// ...but truncate it to 1 bit
}


//
// low-level subtraction w/ borrow
//
static void _sub32(uint32_t a, uint32_t b, uint32_t b_in, uint32_t *d, uint32_t *b_out)
{
		uint64_t t;															// intermediate var
	
		t = (uint64_t)a - (uint64_t)b;					// obtain "wide" difference
		t -= (uint64_t)(b_in & 1);							// take borrow into account
	
		*d = (uint32_t)t;												// return the lower part of result
		*b_out = (uint32_t)(t >> UINT32_BITS);	// return the higher part of result, ...
		*b_out &= (uint32_t)1;									// ...but truncate it to 1 bit
}


//
// low-level multiplication w/ carry and pre-adder
//
static void _mul32(uint32_t a, uint32_t b, uint32_t t, uint32_t c_in, uint32_t *p, uint32_t *c_out)
{
		uint64_t r;															// intermediate result
	
		r = (uint64_t)a * (uint64_t)b;					// obtain wide product
		r += (uint64_t)t;												// handle pre-addition
		r += (uint64_t)c_in;										// take carry into account
	
		*p = (uint32_t)r;												// return the lower part of result
		*c_out = (uint32_t)(r >> UINT32_BITS);	// return the higher part of result, ...
}


//
// end-of-file
//