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

#define MIN(x,y)	(((x)<(y))?(x):(y))
#define MAX(x,y)	(((x)>(y))?(x):(y))

double
elapsed_time(struct timeval tp[2])
{
  return tp[1].tv_sec-tp[0].tv_sec+1e-6*(tp[1].tv_usec-tp[0].tv_usec);
}

int systhr_create(void * (*start_func)(void *), void *arg){
  int status = 0;
  pthread_t tid;
  pthread_attr_t attr;
  pthread_attr_init(&attr);
  status = pthread_attr_setscope(&attr, PTHREAD_SCOPE_SYSTEM);
  if(status == 0) status = pthread_create(&tid, &attr, start_func, arg);
  if(status != 0) status = pthread_create(&tid, 0,     start_func, arg);
  return status;
}

pthread_mutex_t mut = PTHREAD_MUTEX_INITIALIZER;
pthread_cond_t cond = PTHREAD_COND_INITIALIZER;
int progress = 0;

void
barrier_sync(int P, int myid)
{
  int prg;
  pthread_mutex_lock(&mut);
  if((prg = progress + 1) < P)
    {
      progress = prg;
      pthread_cond_wait(&cond, &mut);
    }
  else
    {
      progress = 0;
      pthread_cond_broadcast(&cond);
    }
  pthread_mutex_unlock(&mut);
}

void
sieve_w (int P, int myid, int n, char *tab){
  int n2=n*n,i,j,k,i2,i0,i1,k0,k1;
  /* initialize the table. */
  i0=n2*myid/P; i1=n2*(myid+1)/P;
  for(i=i0; i<i1; i++) tab[i] = 0;
  barrier_sync(P, myid);
  /* main loop {i=2,i2=4}=>{i=4,i2=16}=>{i=16,i2=256}=>... */
  for(i=2,i2=MIN(4,n); i<n; i=i2,i2=MIN(i*i,n))
    /* sieve (write) tab[i*i .. (n2-1)] by reading tab[i .. (i*i-1)]. */
    {
      k0 = i*i + (n2-i*i)*myid/P;
      k1 = i*i + (n2-i*i)*(myid+1)/P;
      for(j=i; j<i2; j++)
	if(tab[j] == 0)
	  if(j*j>k1)
	    break;
	  else
	    for(k=MAX(j*j,k0),k+=j-1,k-=k%j; k<k1; k+=j)
	      tab[k] = 1;
      barrier_sync(P, myid);
    }
}

struct workload {
  int a;
  int P;
  int myid;
  int n;
  char *tab;
};

void *exec_workload(void *w0)
{
  struct workload *w = w0;
  switch(w->a){
  case 1: sieve_w(w->P, w->myid, w->n, w->tab);
  }
}

void
print_prims(int n, char *tab){
  int i;
  for(i=2; i<n*n; i++)
    if(tab[i] == 0)
      printf(" %d", i);
  putchar('\n');
}

void
p_sieve(int P, int n, char *tab)
{
  int i;
  struct workload *w = (struct workload *)malloc(sizeof(struct workload)*P);
  progress = 0;
  for(i=0;i<P;i++)
    {
      w[i].a = 1;
      w[i].P = P;
      w[i].myid = i;
      w[i].n = n;
      w[i].tab = tab;
    }
  for(i=1;i<P;i++) systhr_create(exec_workload, w+i);
  exec_workload(w);
}

main(int argc, char *argv[])
{
  char *tab;
  struct timeval tp[2];
  int n = 10, d = 0, P = 1;
  if(argc > 1) n = atol(argv[1]);
  if(argc > 2) P = atol(argv[2]);
  if(argc > 3) d = 1;
  tab = (char *) malloc(sizeof(char)*n*n);
  gettimeofday(tp, 0);
  p_sieve(P, n, tab);
  gettimeofday(tp+1, 0);
  if(d) print_prims(n, tab);
  printf("%lf\n", elapsed_time(tp));
}