Program calculating LTPs (Left-Truncatable Primes)

ltp-2.cpp

/*
 * Oryginal by "ProgramistaBezDoswiadczenia!"
 * At this point highly modified.
 * See http://hoppke.jogger.pl/2010/01/13/cwiczenie-programistyczne/
 */

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

#define COUNT 2166
#define THREADS 2


static int MaR(unsigned long long x, unsigned long long n) {
	if (x >= n)
		return 0;

	unsigned long long d = 1, y;
	unsigned long t = 0, l;

	for (l = n - 1; !(l & 1); l >>= 1) {
		++t;
	}

	for (; l > 0; l >>= 1) {
		if (l & 1)
			d = (d*x) % n;
		y = (x*x) % n;
		if (y == 1 && x != 1 && x != n-1)
			return 1;
		x = y;
	}

	for (x = d; t; --t) {
		y = (x*x) % n;
		if (y == 1 && x != 1 && x != n-1)
			return 1;
		x = y;
	}

	return x != 1;
}

static int isPrime(unsigned long x) {
	/* http://primes.utm.edu/prove/prove2_3.html */
	if (x == 1) {
		return 0;
	} else if (x < 1373653) {
		return !(MaR(2, x) || MaR(3, x));
	} else if (x < 9080191) {
		return !(MaR(31, x) || MaR(73, x));
	} else {
		return !(MaR(2, x) || MaR(7, x) || MaR(61, x));
	}
}


static unsigned long prev[545], *prev_end;
static unsigned long curr[9][100], *curr_end[9];

static unsigned long q;
static unsigned next_step, busy = THREADS;

static pthread_mutex_t mutex   = PTHREAD_MUTEX_INITIALIZER;
static pthread_cond_t  ready   = PTHREAD_COND_INITIALIZER;
static pthread_cond_t  waiting = PTHREAD_COND_INITIALIZER;
static pthread_t       threads[THREADS];

static void worker_step(unsigned step) {
	const unsigned long front = q * (step + 1);
	unsigned long *out = curr[step];

	const unsigned long *it = prev, *const end = prev_end;
	for (; it != end; ++it) {
		const unsigned long x = front + *it;
		if (isPrime(x)) {
			*out++ = x;
		}
	}

	curr_end[step] = out;
}

static void worker_do(void) {
	while (next_step != 9) {
		unsigned step = next_step;
		++next_step;
		pthread_mutex_unlock(&mutex);
		worker_step(step);
		pthread_mutex_lock(&mutex);
	}
}

static void *worker(void *_ignore) {
	pthread_mutex_lock(&mutex);
	for(;;){
		--busy;
		if (!busy) {
			pthread_cond_signal(&waiting);
		}
		pthread_cond_wait(&ready, &mutex);
		if (!prev_end) {
			break;
		}

		worker_do();
	}
	pthread_mutex_unlock(&mutex);
	return 0;
}

static void finish(void) {
	unsigned i;

	prev_end = 0;
	pthread_cond_broadcast(&ready);
	pthread_mutex_unlock(&mutex);

	for (i = 0; i < THREADS; ++i) {
		pthread_join(threads[i], 0);
	}

	exit(EXIT_SUCCESS);
}


static unsigned n;
static int printAll;

static void findNth() {
	unsigned i;

	/* Let workers do the job */
	pthread_cond_broadcast(&ready);
	next_step = 0;
	worker_do();
	if (busy) {
		pthread_cond_wait(&waiting, &mutex);
	}
	busy = THREADS;

	/* Join lists */
	prev_end = prev;
	for (i = 0; i < 9; ++i) {
		unsigned long *it = curr[i], *end = curr_end[i];
		for (; it != end; ++it) {
			*prev_end++ = *it;
		}
	}

	/* Print */
	if (printAll) {
		unsigned long *it = prev, *end = prev_end;
		for (; it != end; ++it) {
			printf("%lu\n", *it);
			--n;
			if (!n) finish();
		}
	} else if (n > (prev_end - prev)) {
		n -= (prev_end - prev);
	} else {
		printf("%lu\n", prev[n - 1]);
		finish();
	}
}


int main() {
	{
		int _n;
		if (scanf("%d", &_n) != 1 || _n == 0 || _n > COUNT || _n < -COUNT) {
			return EXIT_FAILURE;
		}
		printAll = _n < 0;
		n = _n < 0 ? -_n : _n;
	}

	pthread_mutex_lock(&mutex);

	{
		unsigned i;
		for (i = 0; i < THREADS; ++i) {
			pthread_create(threads + i, 0, worker, 0);
		}
	}

	prev[0] = 0;
	prev_end = prev + 1;
	pthread_cond_wait(&waiting, &mutex);
	busy = THREADS;
	for (q = 1; ; q *= 10) {
		findNth();
	}

	pthread_mutex_unlock(&mutex);
	assert(0);
	return EXIT_FAILURE;
}