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

/* vypocet a^b mod n */
unsigned long long int power(unsigned long long int a, unsigned long long int b, unsigned long long int n)
{
	unsigned int retval = 1;
	
	/* priklad: a^61 = 1 * a * a^4 * a^8 * a^ 16 * a^32 */
	while (b != 0) {
		if (b & 0x01) retval = (retval * a) % n;
		a = (a * a) % n;  
		b = b / 2;
	}
	
	return retval;
}

int millrab(unsigned long long int n)
{
	unsigned long long int s, t, i, j, b, l;
	
	/* n - 1 = 2^s * t */
	s = 0; t = n - 1;
	while (!(t & 0x01)) {
		s++; t = t / 2;
	}
	
	for (i = 0; i < 20; i++) {
		b = 1 + (int) ((double) (n - 1) * (double) rand() / (1.0f + (double) RAND_MAX));
		l = power(b, t, n);
		if (l == 1 || l == (n - 1)) continue;
		for (j = 1; (j < s) && (l != (n - 1)); j++) {
			l = (l * l) % n;
		}
		if (l != (n - 1)) return 0; /* nie je prvocislo */
	}
	
	return 1; /* je (pravdepodobne) prvocislo */
}


int main(int argc, char * argv[])
{
	unsigned long int n;

	if (argc == 2) {
		n = strtol(argv[1], NULL, 10);
	} else {
		n = 65521;
	}
	
	if (millrab(n)) {
		fprintf(stdout, "Cislo %d preslo Miller-Rabinovym testom. (Pravdepodobne je prvocislo.)\n\n", n);
	} else {
		fprintf(stdout, "Cislo %d nepreslo Miller-Rabinovym testom. (Nie je prvocislo.)\n\n", n);
	}

	return 0;
}
