#include <stdio.h>
#include <pthread.h> 
#include <omp.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>

#define N_THREADS 8
#define N_DIMENSIONS 200000000

pthread_barrier_t b;
long *p_s;
long *q_s;
double res[N_THREADS];
int active_threads[N_THREADS];

void setup() {
  p_s = (long *) malloc(N_DIMENSIONS * sizeof(long));
  q_s = (long *) malloc(N_DIMENSIONS * sizeof(long));
  for (long i = 0; i < N_DIMENSIONS; i++) {
    p_s[i] = rand() % 10;
    q_s[i] = rand() % 10;
  }
}

void *calculate(void* arg) {
  int *myrank = (int *)arg;
  long mystart = (long)(N_DIMENSIONS/ (double)N_THREADS * *myrank);
  long myend = (long)(N_DIMENSIONS / (double)N_THREADS * (*myrank + 1));
  double sum = 0;


  for (long i = mystart; i < myend; i++) {
    sum +=  pow(p_s[i] - q_s[i], 2);
  }

  res[*myrank] = sum;
  printf("Thread2 %d finished\n", *myrank);
  fflush(stdout);
  pthread_barrier_wait( &b );

  int n_threads = N_THREADS;

  if (n_threads != 1) {

    n_threads = N_THREADS / 2;
    if (*myrank >= n_threads) {
      active_threads[*myrank] = -1;
    }
  }

  while(1) {
    if(active_threads[*myrank] != -1) {
      res[*myrank] = res[*myrank * 2] + res[*myrank *2 + 1];
    }
    pthread_barrier_wait( &b );

    n_threads = n_threads / 2;
    if (n_threads == 1 || n_threads == 0) {
      break;
    }
    if (*myrank >= n_threads) {
      //printf("Thread %d is on pause\n", *myrank);
      active_threads[*myrank] = -1;
    }
    pthread_barrier_wait( &b );
  }
  
  return NULL;
}


int main() {
  pthread_t t[N_THREADS];
  setup();
  double dt = omp_get_wtime(); // start timing
  int p[N_THREADS];

  if(N_THREADS %2 != 0 && N_THREADS >= 2) {
    printf("Number of threads must be even!\n");
    return 1;
  }
  pthread_barrier_init( &b, NULL, N_THREADS );

  for(int i = 0; i < N_THREADS; i++ ) {
    p[i] = i;
    pthread_create(&t[i], NULL, calculate, (void *)&p[i]);
  }

  for(int i=0; i<N_THREADS; i++)
      pthread_join(t[i], NULL);

  dt = omp_get_wtime() - dt; // finish timing
 
  printf("Execution time: %lf, caluclation result: %lf ", dt, sqrt(res[0]));
  return 0;

  // 8 THREADS: Execution time: 0.807376, caluclation result: 49748.358094
  // 4 THREADS: Execution time: 1.614672, caluclation result: 49748.358094
  // 2 THREADS: Execution time: 3.223820, caluclation result: 57445.749660
  // 1 THREADS: Execution time: 4.698835, caluclation result: 57445.749660 
}