// clang -mavx2 -march=native -std=c99 -O3 -o shortshuffle shortshuffle.c -Wall -Wextra

#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <math.h>
#include <string.h>
#include <assert.h>

#include "pcg.h"


// map random value to [0,range) with slight bias, redraws to avoid bias if needed
static inline uint32_t pcg32_random_bounded_divisionless(uint32_t range) {
    uint64_t random32bit, multiresult;
    uint32_t leftover;
    uint32_t threshold;
    random32bit =  pcg32_random();
    multiresult = random32bit * range;
    leftover = (uint32_t) multiresult;
    if(leftover < range ) {
        threshold = -range % range ;
        while (leftover < threshold) {
            random32bit =  pcg32_random();
            multiresult = random32bit * range;
            leftover = (uint32_t) multiresult;
        }
    }
    return multiresult >> 32; // [0, range)
}









#define RDTSC_START(cycles)                                                   \
    do {                                                                      \
        register unsigned cyc_high, cyc_low;                                  \
        __asm volatile(                                                       \
            "cpuid\n\t"                                                       \
            "rdtsc\n\t"                                                       \
            "mov %%edx, %0\n\t"                                               \
            "mov %%eax, %1\n\t"                                               \
            : "=r"(cyc_high), "=r"(cyc_low)::"%rax", "%rbx", "%rcx", "%rdx"); \
        (cycles) = ((uint64_t)cyc_high << 32) | cyc_low;                      \
    } while (0)

#define RDTSC_FINAL(cycles)                                                   \
    do {                                                                      \
        register unsigned cyc_high, cyc_low;                                  \
        __asm volatile(                                                       \
            "rdtscp\n\t"                                                      \
            "mov %%edx, %0\n\t"                                               \
            "mov %%eax, %1\n\t"                                               \
            "cpuid\n\t"                                                       \
            : "=r"(cyc_high), "=r"(cyc_low)::"%rax", "%rbx", "%rcx", "%rdx"); \
        (cycles) = ((uint64_t)cyc_high << 32) | cyc_low;                      \
    } while (0)


/**
Sorting stuff

*/

int qsort_compare_uint32_t(const void *a, const void *b) {
    return ( *(uint32_t *)a - *(uint32_t *)b );
}
uint32_t *create_sorted_array(size_t length) {
    uint32_t *array = malloc(length * sizeof(uint32_t));
    for (size_t i = 0; i < length; i++) array[i] = (uint32_t) pcg32_random();
    qsort(array, length, sizeof(*array), qsort_compare_uint32_t);
    return array;
}

uint32_t *create_random_array(size_t count) {
    uint32_t *targets = malloc(count * sizeof(uint32_t));
    for (size_t i = 0; i < count; i++) targets[i] = (uint32_t) (pcg32_random() & 0x7FFFFFF);
    return targets;
}




// flushes the array from cache
void array_cache_flush(uint32_t* B, int32_t length) {
    const int32_t CACHELINESIZE = 64;// 64 bytes per cache line
    for(int32_t  k = 0; k < length; k += CACHELINESIZE/sizeof(uint32_t)) {
        __builtin_ia32_clflush(B + k);
    }
}

// tries to put the array in cache
void array_cache_prefetch(uint32_t* B, int32_t length) {
    const int32_t CACHELINESIZE = 64;// 64 bytes per cache line
    for(int32_t  k = 0; k < length; k += CACHELINESIZE/sizeof(uint32_t)) {
        __builtin_prefetch(B + k);
    }
}



/*
 * 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 {                                                              \
            printf("%-60s: ", #test);                                        \
            fflush(NULL);                                                 \
            uint64_t cycles_start, cycles_final, cycles_diff;             \
            uint64_t min_diff = (uint64_t)-1;                             \
            for (int i = 0; i < repeat; i++) {                            \
                pre;                                                       \
                __asm volatile("" ::: /* pretend to clobber */ "memory"); \
                RDTSC_START(cycles_start);                                \
                test;                     \
                RDTSC_FINAL(cycles_final);                                \
                cycles_diff = (cycles_final - cycles_start);              \
                if (cycles_diff < min_diff) min_diff = cycles_diff;       \
            }                                                             \
            uint64_t S = size;                                            \
            float cycle_per_op = (min_diff) / (double)S;                  \
            printf(" %.2f cycles per input key ", cycle_per_op);           \
            printf("\n");                                                 \
            fflush(NULL);                                                 \
 } while (0)



int sortAndCompare(uint32_t * shuf, uint32_t * orig, uint32_t size) {
    qsort(shuf, size, sizeof(uint32_t), qsort_compare_uint32_t);
    qsort(orig, size, sizeof(uint32_t), qsort_compare_uint32_t);
    for(uint32_t k = 0; k < size; ++k)
        if(shuf[k] != orig[k]) {
            printf("[bug]\n");
            return -1;
        }
    return 0;
}



// good old Fisher-Yates shuffle, shuffling an array of integers, without division
void  shuffle_pcg_divisionless(uint32_t *storage, uint32_t size) {
    uint32_t i;
    for (i=size; i>1; i--) {
        uint32_t nextpos = pcg32_random_bounded_divisionless(i);
        uint32_t tmp = storage[i-1];// likely in cache
        uint32_t val = storage[nextpos]; // could be costly
        storage[i - 1] = val;
        storage[nextpos] = tmp; // you might have to read this store later
    }
}

// precomputed factorials... 2!, 3!..., 12!
static uint32_t precomputedFactorial [] = {2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600};


// magic numbers to divide by... 2!, 3!..., 12!
static uint32_t MagicNumberMultipliers [] = {0x80000000, 0xaaaaaaab, 0xaaaaaaab, 0x88888889, 0xb60b60b7, 0xd00d00d1, 0xd00d00d1, 0xb8ef1d2b, 0x24fc9f6f, 0x035cc8ad, 0x011eed8f};
static uint32_t MagicNumberShifts [] = {0, 2, 4, 6, 9, 12, 15, 18, 19, 19, 21};


// divide val by N! for N=2,..., 12!, assumes that val <(1<<31)
static inline uint32_t magicDivide(uint32_t val, uint32_t N) {
    uint32_t highbits = (uint32_t)(((uint64_t) val * MagicNumberMultipliers[N - 2]) >> 32);
    uint32_t answer =  highbits >> MagicNumberShifts[N - 2];
    return answer;
}

// optimized for short arrays, uses division
void  shuffle_pcg_short(uint32_t *storage, uint32_t size) {
    uint32_t i;
    for (i=size; i>12; i--) {
        uint32_t nextpos = pcg32_random_bounded_divisionless(i);
        uint32_t tmp = storage[i-1];
        uint32_t val = storage[nextpos];
        storage[i - 1] = val;
        storage[nextpos] = tmp;
    }
    if (i < 1) return; // nothing to do
    // at this point i <= 12, i >=2
    uint32_t key = pcg32_random_bounded_divisionless(precomputedFactorial[i - 2]);
    for (; i >= 3; i--) {
        uint32_t fa = precomputedFactorial[i - 3];
        uint32_t divresult = key / fa;
        key = key % fa;
        // we swap
        uint32_t tmp = storage[i-1];
        uint32_t val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
    }
    // here i = 2, key is zero or one
    uint32_t tmp = storage[key];
    storage[key] = storage[1];
    storage[1] = tmp;
}


// optimized for short arrays, does not use division
void  shuffle_pcg_divisionless_short(uint32_t *storage, uint32_t size) {
    uint32_t i = size;
   for (; i>12; i--) {
        uint32_t nextpos = pcg32_random_bounded_divisionless(i);
        uint32_t tmp = storage[i-1];
        uint32_t val = storage[nextpos];
        storage[i - 1] = val;
        storage[nextpos] = tmp;
    }
    if (i < 1) return; // nothing to do
    // at this point i <= 12, i >=2
    uint32_t key = pcg32_random_bounded_divisionless(precomputedFactorial[i - 2]);
    for (; i >= 3; i--) {
        uint32_t fa = precomputedFactorial[i - 3];
        uint32_t divresult = magicDivide(key, i - 1);
        key = key - divresult * fa;       // value in [0,(i-1)!)
        // we swap
        uint32_t tmp = storage[i-1];
        uint32_t val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
    }
    // here i = 2, key is zero or one
    uint32_t tmp = storage[key];
    storage[key] = storage[1];
    storage[1] = tmp;
}


static inline uint32_t multiplyAndShift(uint32_t val, uint32_t MagicNumberMultipliers, int MagicNumberShifts) {
    uint32_t highbits = (uint32_t)(((uint64_t) val * MagicNumberMultipliers) >> 32);
    uint32_t answer =  highbits >> MagicNumberShifts;
    return answer;
}

// optimized for short arrays, does not use division
void  shuffle_pcg_divisionless_duff_short(uint32_t *storage, uint32_t size) {
    uint32_t i = size;
   for (; i>12; i--) {
        uint32_t nextpos = pcg32_random_bounded_divisionless(i);
        uint32_t tmp = storage[i-1];
        uint32_t val = storage[nextpos];
        storage[i - 1] = val;
        storage[nextpos] = tmp;
    }
    if (i < 1) return; // nothing to do
    // at this point i <= 12, i >=2
    uint32_t key = pcg32_random_bounded_divisionless(precomputedFactorial[i - 2]);
    uint32_t divresult;
    uint32_t tmp, val;
    switch (i) {
      case 12:
      // we divide by 11!
        divresult = multiplyAndShift(key, 0x035cc8ad ,19);
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
        key -= divresult * 39916800;       // value in [0,(i-1)!)
      case 11:
        // we divide by 10!
        divresult = multiplyAndShift(key, 0x24fc9f6f,19);
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
        key -= divresult * 3628800;       // value in [0,(i-1)!)
      case 10:
        // we divide by 9!
        divresult = multiplyAndShift(key, 0xb8ef1d2b,  18);
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
        key -= divresult * 362880;       // value in [0,(i-1)!)
      case 9:
        // we divide by 8!
        divresult = multiplyAndShift(key, 0xd00d00d1,15);
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
        key -= divresult * 40320;       // value in [0,(i-1)!)
      case 8:
        // we divide by 7!
        divresult = multiplyAndShift(key, 0xd00d00d1, 12);
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
        key -= divresult * 5040;       // value in [0,(i-1)!)
      case 7:
        // we divide by 6!
        divresult = multiplyAndShift(key, 0xb60b60b7,  9);
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
        key -= divresult * 720;       // value in [0,(i-1)!)
      case 6:
        // we divide by 5!
        divresult = multiplyAndShift(key, 0x88888889,6);
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
        key -= divresult * 120;       // value in [0,(i-1)!)
      case 5:
        // we divide by 4!
        divresult = multiplyAndShift(key, 0xaaaaaaab,4);
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
        key -= divresult * 24;       // value in [0,(i-1)!)
      case 4:
        // we divide by 3!
        divresult = multiplyAndShift(key, 0xaaaaaaab,  2);
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
        key -= divresult * 6;       // value in [0,(i-1)!)

      case 3:
        // we divide by 2!
        divresult = key >> 1;
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
        key = key & 1;
      case 2:
        tmp = storage[i-1];
        val = storage[key];
        storage[i - 1] = val;
        storage[key] = tmp;
        i--;
      // no division needed
      case 1:;
      // nothing to do
      default:;
      // nothing
    }
}


// optimized for short arrays, does not use division
void  shuffle_pcg_divisionless_duffcompiler_short(uint32_t *storage, uint32_t size) {
    uint32_t i = size;
   for (; i>12; i--) {
        uint32_t nextpos = pcg32_random_bounded_divisionless(i);
        uint32_t tmp = storage[i-1];
        uint32_t val = storage[nextpos];
        storage[i - 1] = val;
        storage[nextpos] = tmp;
    }
    if (i < 1) return; // nothing to do
    // at this point i <= 12, i >=2
    uint32_t key = pcg32_random_bounded_divisionless(precomputedFactorial[i - 2]);
    uint32_t divresult;
    uint32_t tmp, val;
    switch (i) {
      case 12:
      // we divide by 11!
        divresult = key / 39916800;
        key = key % 39916800;
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
      case 11:
        // we divide by 10!
        divresult = key / 39916800;
        key = key % 39916800;
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
      case 10:
        // we divide by 9!
        divresult = key / 362880;
        key = key % 362880;
        divresult = multiplyAndShift(key, 0xb8ef1d2b,  18);
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
        key -= divresult * 362880;       // value in [0,(i-1)!)
      case 9:
        // we divide by 8!
        divresult = key / 40320;
        key = key % 40320;
        divresult = multiplyAndShift(key, 0xd00d00d1,15);
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
        key -= divresult * 40320;       // value in [0,(i-1)!)
      case 8:
        // we divide by 7!
        divresult = key / 5040;
        key = key % 5040;
        divresult = multiplyAndShift(key, 0xd00d00d1, 12);
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
        key -= divresult * 5040;       // value in [0,(i-1)!)
      case 7:
        // we divide by 6!
        divresult = key / 720;
        key = key % 720;
        divresult = multiplyAndShift(key, 0xb60b60b7,  9);
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
      case 6:
        // we divide by 5!
        divresult = key / 120;
        key = key % 120;
        divresult = multiplyAndShift(key, 0x88888889,6);
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
      case 5:
        // we divide by 4!
        divresult = key / 24;
        key = key % 24;
        divresult = multiplyAndShift(key, 0xaaaaaaab,4);
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
      case 4:
        // we divide by 3!
        divresult = key / 6;
        key = key % 6;
        divresult = multiplyAndShift(key, 0xaaaaaaab,  2);
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;

      case 3:
        // we divide by 2!
        divresult = key / 2;
        key = key % 2;
        divresult = key >> 1;
        tmp = storage[i-1];
        val = storage[divresult];
        storage[i - 1] = val;
        storage[divresult] = tmp;
        i--;
      case 2:
        tmp = storage[i-1];
        val = storage[key];
        storage[i - 1] = val;
        storage[key] = tmp;
        i--;
      // no division needed
      case 1:;
      // nothing to do
      default:;
      // nothing
    }
}


void demo(int size) {
    printf("Shuffling arrays of size %d \n",size);
    printf("Time reported in number of cycles per array element.\n");
    printf("Tests assume that array is in cache as much as possible.\n");
    int repeat = 5000;
    uint32_t * testvalues = create_random_array(size);
    uint32_t * pristinecopy = malloc(size * sizeof(uint32_t));
    memcpy(pristinecopy,testvalues,sizeof(uint32_t) * size);
    if(sortAndCompare(testvalues, pristinecopy, size)!=0) return;


    BEST_TIME(shuffle_pcg_divisionless(testvalues,size), array_cache_prefetch(testvalues,size), repeat, size);
    if(sortAndCompare(testvalues, pristinecopy, size)!=0) return;

    BEST_TIME(shuffle_pcg_short(testvalues,size), array_cache_prefetch(testvalues,size), repeat, size);
    if(sortAndCompare(testvalues, pristinecopy, size)!=0) return;

    BEST_TIME(shuffle_pcg_divisionless_short(testvalues,size), array_cache_prefetch(testvalues,size), repeat, size);
    if(sortAndCompare(testvalues, pristinecopy, size)!=0) return;

    BEST_TIME(shuffle_pcg_divisionless_duff_short(testvalues,size), array_cache_prefetch(testvalues,size), repeat, size);
    if(sortAndCompare(testvalues, pristinecopy, size)!=0) return;

    BEST_TIME(shuffle_pcg_divisionless_duffcompiler_short(testvalues,size), array_cache_prefetch(testvalues,size), repeat, size);
    if(sortAndCompare(testvalues, pristinecopy, size)!=0) return;

    free(testvalues);
    free(pristinecopy);
    printf("\n");
}

int main() {
    demo(10);
    demo(100);
    demo(1000);
    return 0;
}
