Miller-Rabin in C for uint32_t

gistfile1.txt

int MillerRabin(uint32_t a, int n)
{
	/* x^(n-1) % n == 1  */
	/* k == n-1 == d*(2^m) */

	 /*  davon :
	 *  a^d % n == 1
	 * oder
	 *  a^(d * (2^i))  %  n == n-1
	 * davon, i = 1, 2, 3 ...
	 */

	/* (x*y) % n = (x%n)*(y%n) % n */
	/* (x^k) % n = (x%n)^k % n */

	uint32_t x = a, k = n-1; /* x^k */
	uint32_t r = 1; /* Rest */
	uint32_t d, m = 0; /* k = d * (2^m) */


	while (! (k&1)) {
		x = (((uint64_t)x)*x) % n;
		k >>= 1;
		m++;
	}
	d = k;

	while (k != 1) {
		if (k&1) {
			r = (((uint64_t)r)*x) % n;
		}
		x = (((uint64_t)x)*x) % n;
		k >>= 1;
	}

	x = (((uint64_t)x)*r) % n;
	if (x == 1) {
		/* a ist wahrscheinlich Prime */
		/*  a^( d * (2^i) ) mod n == n-1,
		 *  davon :
		 *  i= 0,1,2,...
		 *  d*2^i == k == n-1 */

		/* x := a^d % n */
		x = a;
		r = 1;
		while (d != 1) {
			if (d&1) {
				r = (((uint64_t)r)*x) % n;
			}
			x = (((uint64_t)x)*x) % n;
			d >>= 1;
		}
		x = (((uint64_t)x)*r) % n;
		k = n-1;
		/* a^d % n == 1*/
		if (x == 1 || x == k) {
			return 1;
		}
		while (m--) {
			x = (((uint64_t)x)*x) % n;
			if (x == k) {
				return 1;
			}
		}
		return 0;
	}
	return 0;
}