#include <u.h>
#include <libc.h>

u64int mulmod(u64int a, u64int b, u64int p) {
	u64int r = 0;
	while(b) {
		if(b & 1) r = (r + a) % p;
		a = (a << 1) % p;
		b >>= 1;
	}
	return r;
}

u64int modpot(u64int a, u64int b, u64int p) {
	u64int f = 1;
	while(b) {
		if(b & 1) f = mulmod(f, a, p);
		a = mulmod(a, a, p);
		b >>= 1;
	}
	return f;
}

int miller(u64int n, u64int a) {
	u64int s = 0, t = n - 1, j = 0;
	while(!(t & 1)) {
		t >>= 1;
		s++;
	}
	u64int x = modpot(a, t, n);
	if(x == 1) return 1;
	while(1) {
		if(x+1 == n) return 1;
		x = mulmod(x, x, n);
		j++;
		if((j > s) || (x == 1)) return 0;
	}

}

int main(int argc, char** argv) {
	srand(time(0));
	int i;
	for(i=0;i<1000;i++) {
		u64int n = ((u64int) rand() << 32) | rand(), a;
		u64int topbound = ceil(4 * log(log(n)));
		int p = 1;
		for(a=1;a<=topbound;a++) {
			p = miller(n, a);
			if(!p) break;
		}
//		if(p) print("%lld\n", n);
	}
	return 0;
}