/*
** gcc -Wall -Wextra -pedantic -pedantic-errors
** -std=gnu99 -lpthread -lgmp -lprimesieve
*/

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

#include <gmp.h>
#include <primesieve.h>

#define FILENAME "/dev/stdout"
#define SIZE 1000000000000
#define NTHREADS 4
#define PRECISION 512
#define LENGTH 1000
#define RESOLUTION 1000

mpf_t Sum;
pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
pthread_barrier_t barrier;

unsigned long int primes[NTHREADS] = { 0 };

void *job (void *arg)
{
	unsigned long int index, n, p, size, end;
	primesieve_iterator it;
	time_t t;
	struct tm *time_p;
	FILE *f;
	char buf[BUFSIZ];
	mpf_t N, Localsum;

	pthread_mutex_lock(&mutex);

	index = *((unsigned long int*)arg);
	size = SIZE/NTHREADS;
	n = size*index;

	if (index == NTHREADS - 1)
	{
		end = SIZE - 1;
	}
	else
	{
		end = n + size - 1;
	}

	mpf_set_default_prec((mp_bitcnt_t)PRECISION);
	mpf_inits(N, Localsum, NULL);

	f = fopen(FILENAME, "w");
	setvbuf(f, buf, _IOLBF, BUFSIZ);

	fprintf(f, "%ld: from %ld to %ld\n", index, n + 1, end + 1);

	time(&t);
	time_p = localtime(&t);
	fprintf(f, "[%04d-%02d-%02d %02d:%02d:%02d] %ld: skipping...\n",
		time_p->tm_year + 1900, time_p->tm_mon + 1, time_p->tm_mday,
		time_p->tm_hour, time_p->tm_min, time_p->tm_sec, index);

	if (primes[index])
	{
		primesieve_init(&it);
		p = primes[index];
		primesieve_skipto(&it, p - 1, p + size);
	}
	else
	{
		primesieve_init(&it);
		p = primesieve_nth_prime(n + 1, 0);
		printf("%ld\n", p);
		primesieve_skipto(&it, p - 1, p + size);
	}

	pthread_mutex_unlock(&mutex);

	pthread_barrier_wait(&barrier);

	time(&t);
	time_p = localtime(&t);
	fprintf(f, "[%04d-%02d-%02d %02d:%02d:%02d] %ld: calculating...\n",
		time_p->tm_year + 1900, time_p->tm_mon + 1, time_p->tm_mday,
		time_p->tm_hour, time_p->tm_min, time_p->tm_sec, index);

	while (n <= end)
	{
		mpf_set_ui(N, ++n);
		mpf_div_ui(N, N, (unsigned long int)primesieve_next_prime(&it));
		mpf_sub(Localsum, Localsum, N);

		time(&t);
		time_p = localtime(&t);
		fprintf(f, "[%04d-%02d-%02d %02d:%02d:%02d] %ld: %ld\n",
			time_p->tm_year + 1900, time_p->tm_mon + 1,
			time_p->tm_mday, time_p->tm_hour,
			time_p->tm_min, time_p->tm_sec, index, n);

		mpf_set_ui(N, ++n);
		mpf_div_ui(N, N, (unsigned long int)primesieve_next_prime(&it));
		mpf_add(Localsum, Localsum, N);

		time(&t);
		time_p = localtime(&t);
		fprintf(f, "[%04d-%02d-%02d %02d:%02d:%02d] %ld: %ld\n",
			time_p->tm_year + 1900, time_p->tm_mon + 1,
			time_p->tm_mday, time_p->tm_hour,
			time_p->tm_min, time_p->tm_sec, index, n);

		for (unsigned long int i = 0; i < size/RESOLUTION - 1; i++)
		{
			mpf_set_ui(N, ++n);
			mpf_div_ui(N, N,
				(unsigned long int)primesieve_next_prime(&it));
			mpf_sub(Localsum, Localsum, N);

			mpf_set_ui(N, ++n);
			mpf_div_ui(N, N,
				(unsigned long int)primesieve_next_prime(&it));
			mpf_add(Localsum, Localsum, N);
		}
	}

	primesieve_free_iterator(&it);

	time(&t);
	time_p = localtime(&t);
	fprintf(f, "[%04d-%02d-%02d %02d:%02d:%02d] %ld: done! partial result: ",
		time_p->tm_year + 1900, time_p->tm_mon + 1, time_p->tm_mday,
		time_p->tm_hour, time_p->tm_min, time_p->tm_sec, index);
	mpf_out_str(f, 10, (size_t)LENGTH, Localsum);
	fprintf(f, "\n");

	fclose(f);

	mpf_add(Sum, Sum, Localsum);

	mpf_clears(N, Localsum, NULL);

	return NULL;
}

int main (int argc, char *argv[])
{
	pthread_t threads[NTHREADS];
	unsigned long int args[NTHREADS];
	int i;
	FILE *f;
	char buf[BUFSIZ];
	time_t t;
	struct tm *time_p;

	mpf_set_default_prec((mp_bitcnt_t)PRECISION);
	mpf_init(Sum);

	pthread_barrier_init(&barrier, NULL, NTHREADS);

	f = fopen(FILENAME, "w");
	setvbuf(f, buf, _IOLBF, BUFSIZ);

	time(&t);
	time_p = localtime(&t);

	fprintf(f, "[%04d-%02d-%02d %02d:%02d:%02d] M: started\n",
		time_p->tm_year + 1900, time_p->tm_mon + 1, time_p->tm_mday,
		time_p->tm_hour, time_p->tm_min, time_p->tm_sec);

	fprintf(f, "precision: %d bits\n\n", (int)mpf_get_prec(Sum));

	for (i = 0; i < argc - 1; i++)
	{
		primes[i] = strtoul(argv[i + 1], NULL, 10);
		printf("%ld\n", primes[i]);
	}

	for (i = 0; i < NTHREADS; i++)
	{
		args[i] = i;
		pthread_create(&threads[i], NULL, job, &args[i]);
	}

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

	fprintf(f, "\nM: total result: ");
	mpf_out_str(f, 10, (size_t)LENGTH, Sum);
	fprintf(f, "\n");

	fclose(f);

	mpf_clear(Sum);

	return EXIT_SUCCESS;
}