// gcc -O3 -march=native -std=c99 -o bitsetaggregation bitsetaggregation.c

#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
void *aligned_malloc(size_t alignment, size_t size) {
    void *mem;
    if (posix_memalign(&mem, alignment, size)) exit(1);
    return mem;
}


void __attribute__ ((noinline)) intersection(const uint64_t * restrict array_1,
        const uint64_t * restrict array_2,
		  size_t length, uint64_t * restrict out) {
  for (size_t i =  0; i < length; i ++) {
   const uint64_t word_1 = (array_1[i]) & (array_2[i]);
   out[i] = word_1;
  }
}


void __attribute__ ((noinline)) interleave(const uint64_t * restrict array_1,
        const uint64_t * restrict array_2,
		  size_t length, uint64_t * restrict out) {
  for (size_t i =  0; i + 1< length; i +=2) {
    out[i] = array_1[i];
    out[i + 1] = array_2[i+1];
  }
}


void  __attribute__ ((noinline)) copy(const uint64_t * restrict array_1,
		  size_t length, uint64_t * restrict out) {
  for (size_t i =  0; i < length; i ++) {
   const uint64_t word_1 = array_1[i];
   out[i] = word_1;
  }
}



#define RDTSC_START(cycles)                                             \
    do {                                                                \
        uint32_t cyc_high, cyc_low;                                     \
        __asm volatile("cpuid\n"                                        \
                       "rdtsc\n"                                        \
                       "mov %%edx, %0\n"                                \
                       "mov %%eax, %1" :                                \
                       "=r" (cyc_high),                                 \
                       "=r"(cyc_low) :                                  \
                       : /* no read only */                             \
                       "%rax", "%rbx", "%rcx", "%rdx" /* clobbers */    \
                       );                                               \
        (cycles) = ((uint64_t)cyc_high << 32) | cyc_low;                \
    } while (0)

#define RDTSC_STOP(cycles)                                              \
    do {                                                                \
        uint32_t cyc_high, cyc_low;                                     \
        __asm volatile("rdtscp\n"                                       \
                       "mov %%edx, %0\n"                                \
                       "mov %%eax, %1\n"                                \
                       "cpuid" :                                        \
                       "=r"(cyc_high),                                  \
                       "=r"(cyc_low) :                                  \
                       /* no read only registers */ :                   \
                       "%rax", "%rbx", "%rcx", "%rdx" /* clobbers */    \
                       );                                               \
        (cycles) = ((uint64_t)cyc_high << 32) | cyc_low;                \
    } while (0)

static __attribute__ ((noinline))
uint64_t rdtsc_overhead_func(uint64_t dummy) {
    return dummy;
}

uint64_t global_rdtsc_overhead = (uint64_t) UINT64_MAX;

#define RDTSC_SET_OVERHEAD(test, repeat)			      \
  do {								      \
    uint64_t cycles_start, cycles_final, cycles_diff;		      \
    uint64_t min_diff = UINT64_MAX;				      \
    for (int i = 0; i < repeat; i++) {			      \
      __asm volatile("" ::: /* pretend to clobber */ "memory");	      \
      RDTSC_START(cycles_start);				      \
      test;							      \
      RDTSC_STOP(cycles_final);                                       \
      cycles_diff = (cycles_final - cycles_start);		      \
      if (cycles_diff < min_diff) min_diff = cycles_diff;	      \
    }								      \
    global_rdtsc_overhead = min_diff;				      \
    printf("rdtsc_overhead set to %d\n", (int)global_rdtsc_overhead);     \
  } while (0)							      \


/*
 * Prints the best number of operations per cycle where
 * test is the function call, answer is the expected answer generated by
 * test, repeat is the number of times we should repeat and size is the
 * number of operations represented by test.
 */
#define BEST_TIME(test, pre, repeat, size)                         \
        do {                                                              \
            if (global_rdtsc_overhead == UINT64_MAX) {                    \
               RDTSC_SET_OVERHEAD(rdtsc_overhead_func(1), repeat);        \
            }                                                             \
            printf("%-60s\t: ", #test);                                        \
            fflush(NULL);                                                 \
            uint64_t cycles_start, cycles_final, cycles_diff;             \
            uint64_t min_diff = (uint64_t)-1;                             \
            uint64_t sum_diff = 0;                                        \
            for (int i = 0; i < repeat; i++) {                            \
                pre;                                                      \
                __asm volatile("" ::: /* pretend to clobber */ "memory"); \
                RDTSC_START(cycles_start);                                \
                test;                    \
                RDTSC_STOP(cycles_final);                                \
                cycles_diff = (cycles_final - cycles_start - global_rdtsc_overhead);           \
                if (cycles_diff < min_diff) min_diff = cycles_diff;       \
                sum_diff += cycles_diff;                                  \
            }                                                             \
            uint64_t S = size;                                            \
            float cycle_per_op = (min_diff) / (double)S;                  \
            float avg_cycle_per_op = (sum_diff) / ((double)S * repeat);   \
            printf(" %.2f cycles per operation (best) ", cycle_per_op);   \
            printf("\t%.2f cycles per operation (avg) ", avg_cycle_per_op);   \
            printf("\n");                                                 \
            fflush(NULL);                                                 \
 } while (0)



void demo(int size) {
    printf("size = %d words or %lu bytes \n",size,  size*sizeof(uint64_t));
    int repeat = 500;
    uint64_t * dataA = aligned_malloc(32,size * sizeof(uint64_t));
    uint64_t * dataB = aligned_malloc(32,size * sizeof(uint64_t));
    uint64_t * out = aligned_malloc(32,size * sizeof(uint64_t));

    for(int k = 0; k < size; ++k) {
    	dataA[k]  = -k;
    	dataB[k] = k;
    }
    BEST_TIME(interleave(dataA, dataB, size, out),, repeat, size);
    BEST_TIME(intersection(dataA, dataB, size, out),, repeat, size);
    BEST_TIME(copy(dataA, size, out),, repeat, size);
    BEST_TIME(memcpy(out, dataA, sizeof(uint64_t)*size),, repeat, size);

    free(dataA);
    free(dataB);
    free(out);
    printf("\n");
}

int main() {
    for(int w = 8; w <= 8192; w *= 2) {
      demo(w);
    }
    return 0;
}
