#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;

pthread_mutex_t mut2 = PTHREAD_MUTEX_INITIALIZER;
pthread_cond_t cond2 = PTHREAD_COND_INITIALIZER;
int progress2 = 4;

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);
}

int wait_progress2(int i){
  int prg;
  pthread_mutex_lock(&mut2);
  while(! (i<(prg=progress2)))
    pthread_cond_wait(&cond2, &mut2);
  pthread_mutex_unlock(&mut2);
  return prg;
}

int wait_and_signal_progress2(int newp, int minprg){
  int prg;
  pthread_mutex_lock(&mut2);
  if((prg=progress2) < newp){
    while(! (progress2>=minprg))
      pthread_cond_wait(&cond2, &mut2);
    prg = progress2 = newp;
    pthread_cond_broadcast(&cond2);
  }
  pthread_mutex_unlock(&mut2);
  return prg;
}

void
sieve_w (int P, int myid, int n, char *tab){
  int n2=n*n,i,k,i2,newp,maxnewp,i0,i1,k0,prg;
  i0=n2*myid/P; i1=n2*(myid+1)/P;
  /* initialize the table. */
  for(i=i0; i<i1; i++) tab[i] = 0;
  maxnewp = MIN(i1, n); prg = 4;
  /* main loop */
  for(i=2; i<n; i++){
    if(! (i<prg))
      prg = wait_progress2(i);
    if(tab[i] == 0){
      i2=i*i;
      newp = MIN(i2, maxnewp);
      if(newp > i0)
	prg = wait_and_signal_progress2(newp, i0);
      if(i2 >= i1) break;
      for(k=MAX(i2,i0),k+=i-1,k-=k%i; k<i1; k+=i)
	tab[k] = 1;
    }
  }
  if(maxnewp > i0)
    prg = wait_and_signal_progress2(maxnewp, i0);
  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));
}