#include <windows.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <omp.h>
#include <time.h>
#include <windows.h>

# ifndef __INLINE
#   ifdef __GNUC__
#     define __INLINE inline __attribute__((always_inline, __target__("default")))
#   endif
# endif

template <typename T>
__INLINE
static
void volatilize(T &r) noexcept
{
  asm volatile ("" : "=m" (r) : "m" (r) :);
}

__INLINE
static
void cpu_pause()
{
  asm volatile ("pause" : : :);
}

__INLINE
static
uint64_t xchg(uint64_t &recv, uint64_t src)
{
  uint64_t register ret;
  asm volatile ("\
      xchg %[ret], %[recv];\n"
    : [ret] "=r" (ret), [recv] "=m" (recv)
    : "0" (src)
    :
  );
  return ret; // old value
}

__INLINE
static
uint64_t tsc_get_ticks() noexcept
{
  register size_t tsc_lo;
  register size_t tsc_hi;
  register size_t tsc_id;
  register uint64_t ret;
  asm volatile ("\
    mov $0, %[tsc_hi];\n\
    mov $0, %[tsc_lo];\n\
    mov $0, %[tsc_id];\n\
    rdtscp;\n"
    : [tsc_lo] "=a" (tsc_lo), [tsc_hi] "=d" (tsc_hi), [tsc_id] "=c" (tsc_id)
    :
    :
  );
  return (tsc_hi << 32) | tsc_lo;
}


struct TSC
{
  uint64_t  value;
  timespec  tmsp;
};

static
double tsc_get_ticks_per_nanosec
    (
      uint32_t pre_delay,
      uint32_t delay,
      TSC *tsc = NULL
    )
{
  alignas(16)
  timespec tmsp_before;
  alignas(8)
  timespec tmsp_after;

  register double ticks_per_nanosec = 0;
  register uint64_t
      ticks_before,
      ticks_after,
      nsec_delta;

  Sleep(pre_delay);
  // -- barrier of insructions
  tsc_get_ticks();
  clock_gettime(CLOCK_REALTIME, &tmsp_before);
  ticks_before = tsc_get_ticks();
  Sleep(delay);
  // -- barrier of instructions
  tsc_get_ticks();
  clock_gettime(CLOCK_REALTIME, &tmsp_after);
  ticks_after = tsc_get_ticks();
  if (NULL != tsc)
  {
    tsc->value = ticks_after;
    tsc->tmsp = tmsp_after;
  }
  if (ticks_after < ticks_before)
    ticks_after += ~(--ticks_before);
  else
    ticks_after -= ticks_before;

  tmsp_after.tv_sec -= tmsp_before.tv_sec;
  tmsp_after.tv_nsec -= tmsp_before.tv_nsec;
  if (0 > tmsp_after.tv_nsec)
  {
    tmsp_after.tv_sec--;
    tmsp_after.tv_nsec += 1000000000L;
  }
  nsec_delta = tmsp_after.tv_sec * 1000000000 + tmsp_after.tv_nsec;
  ticks_per_nanosec = (double)ticks_after / nsec_delta;

  return ticks_per_nanosec;
}

class TSC_ELAPSED
{
  double ticks_per_nsec;
  uint64_t before;
  uint64_t after;

public:
  __INLINE
  TSC_ELAPSED() : before(0), after(0)
  {
    this->ticks_per_nsec = tsc_get_ticks_per_nanosec(16, 7000);
  };

  __INLINE
  void start()
  {
    Sleep(16);
    this->before = tsc_get_ticks();
  }

  __INLINE
  void stop()
  {
    this->after = tsc_get_ticks();
  }

  __INLINE
  double get_ticks_per_nsec()
  {
    return this->ticks_per_nsec;
  }

  __INLINE
  double evaluate_nsec()
  {
    register double ret = 1;
    TSC_ELAPSED & _ = *this;
    if (_.after < _.before)
      ret += _.after + (~_.before);
    else
      ret = _.after - _.before;
    return ret /= _.ticks_per_nsec;
  }
};

constexpr size_t dim = 16384;
//constexpr size_t dim = 5;
typedef int TArr[dim][dim];
static TArr arr;

int64_t diamond1()
{
  const size_t Half = dim / 2;
  int64_t res = 0;

  for (size_t i=0; i < Half; ++i)
    for (size_t j = Half-i; j <= Half+i; ++j)
    {
      res += arr[i][j];
      res += arr[dim-i-1][j];
    }

  for (int i=0; i < dim; ++i)
    res += arr[Half][i];

  return res;
}

int64_t diamond2()
{
  const size_t Half = dim / 2;
  int64_t res = 0;

#pragma omp parallel for reduction (+: res)
  for (size_t i=0; i < Half; ++i)
    for (size_t j = Half-i; j <= Half+i; ++j)
    {
      res += arr[i][j];
      res += arr[dim-i-1][j];
    }

#pragma omp parallel for reduction (+: res)
  for (int i=0; i < dim; ++i)
    res += arr[Half][i];

  return res;
}

int main(int argc, char** argv)
{
  printf("\r\nTSC calibrating, please wait...\r\n");
  TSC_ELAPSED elapsed;
  volatilize(elapsed);
  printf("TSC ticks per nanosec: %lf\r\n", elapsed.get_ticks_per_nsec());

  if (NULL == arr)
    return 1;

  printf("\r\nSize of matrix: %lli\r\n", sizeof(arr));

  constexpr size_t dim_plain = dim * dim;
  int *e = (int*)arr;
  for (size_t i = 0; i < dim_plain; i++) e[i] = i + 1;

  int64_t val;

#pragma omp flush(arr)
  elapsed.start();
  val = diamond1();
  elapsed.stop();
  printf
      (
        "\r\nDIAMOND_1\r\nSUM: %lli\r\nElapsed time(nsec): %lf\r\n",
        val,
        elapsed.evaluate_nsec()
      );

#pragma omp flush(arr)
  elapsed.start();
  val = diamond2();
  elapsed.stop();
  printf
      (
        "\r\nOpenMP warmingup\r\nSUM: %lli\r\nElapsed time(nsec): %lf\r\n",
        val,
        elapsed.evaluate_nsec()
      );

#pragma omp flush(arr)
  elapsed.start();
  volatilize(val);
  val = diamond2();
  elapsed.stop();
  printf
      (
        "\r\nDIAMOND_2\r\nSUM: %lli\r\nElapsed time(nsec): %lf\r\n",
        val,
        elapsed.evaluate_nsec()
      );

  return (EXIT_SUCCESS);
}