modular_multiplicate.c

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <assert.h>

#define Assert assert
#define uint64 uint64_t

/*
 * Calculate (x * y) % m, where x and y in [0, 2^64), m in [1, 2^64).
 *
 * If x or y is greater than 2^32, improved interleaved modular
 * multiplication algorithm is used to avoid overflow.
 */
static uint64 modular_multiplicate(uint64 x, uint64 y, const uint64 m)
{
	int		i, bits;
	uint64		r = 0;

	Assert(1 <= m);

	/* Because of (x * y) % m = (x % m * y % m) % m */
	if (x >= m)
		x %= m;
	if (y >= m)
		y %= m;

	/* Return the trivial result. */
	if (x == 0 || y == 0 || m == 1)
		return 0;

	/* Return the result if (x * y) can be multiplicated without overflow. */
	if ((x | y) < (0xffffffff))
		return (x * y) % m;

	/* To reduce the for loop in the algorithm below. */
	if (x < y)
	{
		uint64 tmp = x;
		x = y;
		y = tmp;
	}

	/* Interleaved modular multiplication algorithm [1]
	 *
	 * This algorithm is usually used in the field of digital circuit
	 * design.
	 *
	 * Input: X, Y, M; 0 <= X, Y <= M;
	 * Output: R = X *  Y mod M;
	 * bits: number of bits of Y
	 * Y[i]: i th bit of Y
	 *
	 * 1. R = 0;
	 * 2. for (i = bits - 1; i >= 0; i--) {
	 * 3. 	R = 2 * R;
	 * 4. 	if (Y[i] == 0x1)
	 * 5. 		R += X;
	 * 6. 	if (R >= M) R -= M;
	 * 7.	if (R >= M) R -= M;
	 *   }
	 *
	 * In Steps 3 and 5, overflow should be avoided.
	 * Steps 6 and 7 can be instead of a modular operation (R %= M).
	 *
	 * Reference
	 * [1] D.N. Amanor, et al, "Efficient hardware architecture for
	 *    modular multiplication on FPGAs", in Field Programmable
	 *    Logic and Apllications, 2005. International Conference on,
	 *    Aug 2005, pp. 539-542.
	 */

	bits = 64;
	while (bits > 0 && (y >> (64 - bits) | 0x1) == 0)
		bits--;

	for (i = bits - 1; i >= 0; i--)
	{
		if (r > 0x7fffffffffffffff)
			/* To avoid overflow, transform from (2 * r) to
			 * (2 * r) % m, and further transform to
			 * mathematically equivalent form shown below:
			 */
			r = m - ((m - r) << 1);
		else
			r <<= 1;

		if ((y >> i) & 0x1)
		{
			/* Calculate (r + x) without overflow using same
			 * transformations described in the above comment.
			 */
			if (m > 0x7fffffffffffffff)
				r = ((m - r) > x) ? r + x : r + x - m;
			else
				r = (r > m) ? r - m + x : r + x;
		}

		r %= m;
	}

	return r;
}

int main(int argc, char **argv) {

	if (argc != 4) {
		printf("Syntax Error:\n\tUsage:%s A B N\n", argv[0]);
		return -1;
	}

	uint64_t a  = strtouq(argv[1], NULL, 10);
	uint64_t b  = strtouq(argv[2], NULL, 10);
	uint64_t n  = strtouq(argv[3], NULL, 10);

	uint64_t r =  modular_multiplicate(a, b, n);
	printf("(%llu * %llu) %% %llu = %llu\n", a, b, n, modular_multiplicate(a, b, n));
       	return 0;
}