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

uint64_t nextnum = 3;
/* Fifty million takes a while
 * For testing, you might want to reduce this quite a bit
 */ 
#define max 50000000
uint64_t pcount = 1;
uint64_t primes[max] = {2};
pthread_mutex_t report_mutex = PTHREAD_MUTEX_INITIALIZER;
pthread_mutex_t num_mutex = PTHREAD_MUTEX_INITIALIZER;

/* Only needed if using quicksort */
int intcmp(const uint64_t *a, const uint64_t *b){
	return *a - *b;
}

char check_prime(uint64_t num){
	for(uint64_t idx = 0; idx < pcount && primes[idx] * primes[idx] <= num; idx++)
		if(!(num % primes[idx]))
			return 0;
	return 1;
}

void* prime_thread(void*){
	uint64_t num;
	uint64_t idx;
	while(pcount < max){
		pthread_mutex_lock(&num_mutex);
		num = nextnum;
		nextnum += 2;
		pthread_mutex_unlock(&num_mutex);

		if(check_prime(num)){
			pthread_mutex_lock(&report_mutex);
			if(pcount < max)
				primes[pcount++] = num;
			else 
				if(primes[pcount - 1] > num)
					primes[pcount - 1] = num;

			/* Backbubble process */
			for(uint64_t i = pcount - 1; i > 1; i--) {
				if(primes[i-1] > primes[i]) {
					register uint64_t t = primes[i-1];
					primes[i-1] = primes[i];
					primes[i] = t;
				} else break;
			}
			pthread_mutex_unlock(&report_mutex);
			if(primes[pcount] * primes[pcount] >= (uint64_t)max * 20)
				break;
		}
	}
#define BLOCKSIZE 20000 // Note:  Must be even
//	printf("Moving to phase 2, pcount = %lu, primes[pcount-1] = %lu\n", pcount, primes[pcount-1]);
	uint64_t local_results[BLOCKSIZE];
	uint64_t local_end = 0;
	while(pcount < max){
		pthread_mutex_lock(&num_mutex);
		num = nextnum;
		nextnum += BLOCKSIZE;
		pthread_mutex_unlock(&num_mutex);
		local_end = 0;
		for(int i = 0; i < BLOCKSIZE/2; i++) {
			if(check_prime(num))
				local_results[local_end++] = num;
			num += 2;
		}
		
		pthread_mutex_lock(&report_mutex);
		for(int i = 0; i < local_end; i++){
			if(pcount < max)
				primes[pcount++] = local_results[i];	
		}
		pthread_mutex_unlock(&report_mutex);
	}
}

/* Helper function for test code near the end of main */
char is_in_list(uint64_t num){
	for(int64_t i = 0; i < max; i++)
		if(primes[i] == num)
			return 1;
	return 0;
}

int main(int argc, char ** argv){
	int threads = 1;
	if(argc > 1)
		threads = atoi(argv[1]);
	pthread_t handles[threads];

	for(int i = 0; i < threads; i++)	
		pthread_create(handles + i, 0, prime_thread, 0);

	for(int i = 0; i < threads; i++)
		pthread_join(handles[i], 0);

	/* Duplicate prime detection assumes the list is sorted
 	 * Uncomment this if your list isn't already sorted at this point
 	 */
	printf("Sorting List\n");
	qsort(primes, max, sizeof(uint64_t), (int (*)(const void*, const void*))intcmp);

	printf("Middle 25 primes:\n");
	for(uint64_t i = max/2; i < max/2 + 25; i++)
		printf("%lu\n", primes[i]);

	printf("Last 25 primes:\n");
	for(uint64_t i = max - 25; i < max; i++)
		printf("%lu\n", primes[i]);

	/* These tests take a while, so they'll throw off any timing we do */

	for(uint64_t i = 0; i < pcount - 1; i++)
		if(primes[i] == primes[i+1])
			printf("Error:  Duplicate prime recorded at index %d and %d (%lu, %lu)\n", i, i + 1, primes[i], primes[i+1]);

	printf("Done checking for duplicate primes, checking for out of order primes\n");

	for(uint64_t i = 0; i < pcount - 1; i++)
		if(primes[i] > primes[i+1]) 
			printf("Note:  Out of order primes %lu and %lu\n", primes[i], primes[i+1]);

	printf("Done checking for out-of-order primes, checking for non-primes\n");
/*
	for(uint64_t i = 0; i < pcount; i++) 
		if(!check_prime(primes[i])) 
			printf("Error:  %lu at index %lu is not prime!\n", primes[i], i);

	printf("Done checking for non-primes, checking for missed primes\n");

	for(uint64_t i = 2; i < primes[pcount-1]; i++)
		if(check_prime(i))
			if(!is_in_list(i))
				printf("Error:  %lu is prime, but not in the list\n", i);
*/
	return 0;
}