// gcc -march=native -O3 -o simdunion16 simdunion16.c
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <stdbool.h>
#include <string.h>
#include <inttypes.h>
#include <assert.h>
#include "benchmark.h"


#include <x86intrin.h>

// print array (duh)
static void array_print(const uint16_t * x, uint32_t len) {
    for(uint32_t i = 0; i < len; i++) {
        printf("%u,",x[i]);
    }
    printf("\n");
}

// check that array is sorted (duh)
static bool array_is_sorted(const uint16_t * x, uint32_t len) {
    for(uint32_t i = 1; i < len; i++) {
        if(x[i] < x[i-1]) return false;
    }
    return true;
}


// check that arrays are equal
static bool arrays_are_equal(const uint16_t * x,const uint16_t * y, uint32_t len) {
    for(uint32_t i = 0; i < len; i++) {
        if(x[i] != y[i]) return false;
    }
    return true;
}

// check that array contains strictly sorted 16-bit values
static bool array_is_strictly_sorted(const uint16_t * x, uint32_t len) {
    for(uint32_t i = 1; i < len; i++) {
        if(x[i] <= x[i-1]) return false;
    }
    return true;
}


// check that vector contains sorted 16-bit values
static bool is_sorted(__m128i x) {
    return ( (uint16_t)_mm_extract_epi16(x,0)<=(uint16_t)_mm_extract_epi16(x,1))
           && ( (uint16_t)_mm_extract_epi16(x,1)<=(uint16_t)_mm_extract_epi16(x,2))
           && ( (uint16_t)_mm_extract_epi16(x,2)<=(uint16_t)_mm_extract_epi16(x,3))
           && ( (uint16_t)_mm_extract_epi16(x,3)<=(uint16_t)_mm_extract_epi16(x,4))
           && ( (uint16_t)_mm_extract_epi16(x,4)<=(uint16_t)_mm_extract_epi16(x,5))
           && ( (uint16_t)_mm_extract_epi16(x,5)<=(uint16_t)_mm_extract_epi16(x,6))
           && ( (uint16_t)_mm_extract_epi16(x,6)<=(uint16_t)_mm_extract_epi16(x,7));
}

// print vector as if it were an array of 16-bit unsigned integers
static void print16(__m128i x) {
    printf("%u %u %u %u %u %u %u %u \n", _mm_extract_epi16(x,0),_mm_extract_epi16(x,1),_mm_extract_epi16(x,2),_mm_extract_epi16(x,3)
           , _mm_extract_epi16(x,4),_mm_extract_epi16(x,5),_mm_extract_epi16(x,6),_mm_extract_epi16(x,7));
}

// Assuming that vInput1 and vInput2 are sorted, produces a sorted output going from vecMin all the way to vecMax
// developed originally for merge sort using SIMD instructions.
// Standard merge. See, e.g., Inoue and Taura, SIMD- and Cache-Friendly Algorithm for Sorting an Array of Structures
static void sse_merge(__m128i *vInput1, __m128i *vInput2, // input 1 & 2
                      __m128i *vecMin, __m128i *vecMax) { // output
    __m128i vecTmp;
    vecTmp = _mm_min_epu16(*vInput1, *vInput2);
    *vecMax = _mm_max_epu16(*vInput1, *vInput2);
    vecTmp = _mm_alignr_epi8(vecTmp, vecTmp, 2);
    *vecMin = _mm_min_epu16(vecTmp, *vecMax);
    *vecMax = _mm_max_epu16(vecTmp, *vecMax);
    vecTmp = _mm_alignr_epi8(*vecMin, *vecMin, 2);
    *vecMin = _mm_min_epu16(vecTmp, *vecMax);
    *vecMax = _mm_max_epu16(vecTmp, *vecMax);
    vecTmp = _mm_alignr_epi8(*vecMin, *vecMin, 2);
    *vecMin = _mm_min_epu16(vecTmp, *vecMax);
    *vecMax = _mm_max_epu16(vecTmp, *vecMax);
    vecTmp = _mm_alignr_epi8(*vecMin, *vecMin, 2);
    *vecMin = _mm_min_epu16(vecTmp, *vecMax);
    *vecMax = _mm_max_epu16(vecTmp, *vecMax);
    vecTmp = _mm_alignr_epi8(*vecMin, *vecMin, 2);
    *vecMin = _mm_min_epu16(vecTmp, *vecMax);
    *vecMax = _mm_max_epu16(vecTmp, *vecMax);
    vecTmp = _mm_alignr_epi8(*vecMin, *vecMin, 2);
    *vecMin = _mm_min_epu16(vecTmp, *vecMax);
    *vecMax = _mm_max_epu16(vecTmp, *vecMax);
    vecTmp = _mm_alignr_epi8(*vecMin, *vecMin, 2);
    *vecMin = _mm_min_epu16(vecTmp, *vecMax);
    *vecMax = _mm_max_epu16(vecTmp, *vecMax);
    *vecMin = _mm_alignr_epi8(*vecMin, *vecMin, 2);
}



// used by store_unique, generated by simdunion.py
static uint8_t uniqshuf[]= {0x0,0x1,0x2,0x3,0x4,0x5,0x6,0x7,0x8,0x9,0xa,0xb,0xc,0xd,0xe,0xf,
                            0x2,0x3,0x4,0x5,0x6,0x7,0x8,0x9,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x6,0x7,0x8,0x9,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,
                            0x4,0x5,0x6,0x7,0x8,0x9,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x6,0x7,0x8,0x9,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,
                            0x2,0x3,0x6,0x7,0x8,0x9,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x6,0x7,0x8,0x9,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x6,0x7,0x8,0x9,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x8,0x9,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x8,0x9,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x8,0x9,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x8,0x9,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x8,0x9,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x8,0x9,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x8,0x9,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x8,0x9,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x6,0x7,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x6,0x7,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x6,0x7,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x6,0x7,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x6,0x7,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x6,0x7,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x6,0x7,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x6,0x7,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0xa,0xb,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x6,0x7,0x8,0x9,0xc,0xd,0xe,0xf,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x6,0x7,0x8,0x9,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x6,0x7,0x8,0x9,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x6,0x7,0x8,0x9,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x6,0x7,0x8,0x9,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x6,0x7,0x8,0x9,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x6,0x7,0x8,0x9,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x6,0x7,0x8,0x9,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x8,0x9,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x8,0x9,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x8,0x9,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x8,0x9,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x8,0x9,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x8,0x9,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x8,0x9,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x8,0x9,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x6,0x7,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x6,0x7,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x6,0x7,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x6,0x7,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x6,0x7,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x6,0x7,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x6,0x7,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x6,0x7,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0xc,0xd,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x6,0x7,0x8,0x9,0xa,0xb,0xe,0xf,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x6,0x7,0x8,0x9,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x6,0x7,0x8,0x9,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x6,0x7,0x8,0x9,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x6,0x7,0x8,0x9,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x6,0x7,0x8,0x9,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x6,0x7,0x8,0x9,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x6,0x7,0x8,0x9,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x8,0x9,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x8,0x9,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x8,0x9,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x8,0x9,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x8,0x9,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x8,0x9,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x8,0x9,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x8,0x9,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x6,0x7,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x6,0x7,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x6,0x7,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x6,0x7,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x6,0x7,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x6,0x7,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x6,0x7,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x6,0x7,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0xa,0xb,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x6,0x7,0x8,0x9,0xe,0xf,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x6,0x7,0x8,0x9,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x6,0x7,0x8,0x9,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x6,0x7,0x8,0x9,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x6,0x7,0x8,0x9,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x6,0x7,0x8,0x9,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x6,0x7,0x8,0x9,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x6,0x7,0x8,0x9,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x8,0x9,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x8,0x9,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x8,0x9,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x8,0x9,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x8,0x9,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x8,0x9,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x8,0x9,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x8,0x9,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x6,0x7,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x6,0x7,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x6,0x7,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x6,0x7,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x6,0x7,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x6,0x7,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x6,0x7,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x6,0x7,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0xe,0xf,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x6,0x7,0x8,0x9,0xa,0xb,0xc,0xd,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x6,0x7,0x8,0x9,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x6,0x7,0x8,0x9,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x6,0x7,0x8,0x9,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x6,0x7,0x8,0x9,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x6,0x7,0x8,0x9,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x6,0x7,0x8,0x9,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x6,0x7,0x8,0x9,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x8,0x9,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x8,0x9,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x8,0x9,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x8,0x9,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x8,0x9,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x8,0x9,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x8,0x9,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x8,0x9,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x6,0x7,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x6,0x7,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x6,0x7,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x6,0x7,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x6,0x7,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x6,0x7,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x6,0x7,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x6,0x7,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0xa,0xb,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x6,0x7,0x8,0x9,0xc,0xd,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x6,0x7,0x8,0x9,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x6,0x7,0x8,0x9,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x6,0x7,0x8,0x9,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x6,0x7,0x8,0x9,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x6,0x7,0x8,0x9,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x6,0x7,0x8,0x9,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x6,0x7,0x8,0x9,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x8,0x9,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x8,0x9,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x8,0x9,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x8,0x9,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x8,0x9,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x8,0x9,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x8,0x9,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x8,0x9,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x6,0x7,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x6,0x7,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x6,0x7,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x6,0x7,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x6,0x7,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x6,0x7,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x6,0x7,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x6,0x7,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0xc,0xd,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x6,0x7,0x8,0x9,0xa,0xb,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x6,0x7,0x8,0x9,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x6,0x7,0x8,0x9,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x6,0x7,0x8,0x9,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x6,0x7,0x8,0x9,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x6,0x7,0x8,0x9,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x6,0x7,0x8,0x9,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x6,0x7,0x8,0x9,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x8,0x9,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x8,0x9,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x8,0x9,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x8,0x9,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x8,0x9,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x8,0x9,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x8,0x9,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x8,0x9,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x6,0x7,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x6,0x7,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x6,0x7,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x6,0x7,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x6,0x7,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x6,0x7,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x6,0x7,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x6,0x7,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0xa,0xb,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x6,0x7,0x8,0x9,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x6,0x7,0x8,0x9,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x6,0x7,0x8,0x9,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x6,0x7,0x8,0x9,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x6,0x7,0x8,0x9,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x6,0x7,0x8,0x9,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x6,0x7,0x8,0x9,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x6,0x7,0x8,0x9,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x8,0x9,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x8,0x9,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x8,0x9,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x8,0x9,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x8,0x9,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x8,0x9,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x8,0x9,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x8,0x9,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0x6,0x7,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0x6,0x7,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0x6,0x7,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0x6,0x7,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x6,0x7,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x6,0x7,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x6,0x7,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x6,0x7,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0x4,0x5,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0x4,0x5,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x4,0x5,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x4,0x5,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0x2,0x3,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x2,0x3,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0x0,0x1,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
                            0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF
                           };

// write vector new, while omitting repeated values assuming that previously written vector was "old"
static inline int store_unique(__m128i old,__m128i new, uint16_t *  output) {
    __m128i vecTmp = _mm_alignr_epi8(new, old, 16 - 2);
    int M = _mm_movemask_epi8(_mm_cmpeq_epi16(vecTmp,new));
    M = _pext_u32(M,0x5555);
    int numberofnewvalues = 8 - _mm_popcnt_u32(M);
    __m128i key =  _mm_lddqu_si128((const __m128i* )uniqshuf + M);
    __m128i val =  _mm_shuffle_epi8(new,key);
    _mm_storeu_si128((__m128i* )output,val);
    return numberofnewvalues;
}

static inline int store_unique_pack(__m128i old,__m128i new, uint16_t *  output) {
    __m128i vecTmp = _mm_alignr_epi8(new, old, 16 - 2);
    int M = _mm_movemask_epi8(_mm_packs_epi16(_mm_cmpeq_epi16(vecTmp,new),_mm_setzero_si128()));
    int numberofnewvalues = 8 - _mm_popcnt_u32(M);
    __m128i key =  _mm_lddqu_si128((const __m128i* )uniqshuf + M);
    __m128i val =  _mm_shuffle_epi8(new,key);
    _mm_storeu_si128((__m128i* )output,val);
    return numberofnewvalues;
}

// working in-place, this function overwrites the repeated values
static uint32_t unique(uint16_t *  out, uint32_t len) {
    uint32_t pos = 1;
    for(uint32_t i = 1; i < len; ++i) {
        if(out[i] != out[i-1]) {
            out[pos++] = out[i];
        }
    }
    return pos;
}


// the simple union algo. first compute a union with repeated values, then uses the "unique" function (as a second pass) to prune them out.
static uint32_t sse_unite_simple(uint16_t * __restrict__ array1, uint32_t length1 ,
                                 uint16_t *  __restrict__ array2, uint32_t length2, uint16_t *  __restrict__  output) {
    __m128i vA, vB, V, vecMin,vecMax;
    uint16_t * initoutput = output;
    assert(length1 % 8 == 0); // because I am lazy
    assert(length2 % 8 == 0); // because I am lazy
    assert(length1 != 0); // because I am lazy
    assert(length2 != 0); // because I am lazy
    uint32_t len1 = length1 / 8;
    uint32_t len2 = length2 / 8;
    uint32_t pos1 = 0;
    uint32_t pos2 = 0;
    // we start the machine
    vA = _mm_lddqu_si128((const __m128i* )array1 + pos1);
    pos1++;
    vB = _mm_lddqu_si128((const __m128i* )array2 + pos2);
    pos2++;
    uint16_t curA, curB;
    curA = (pos1 < len1) ? array1[8*pos1] : 0;
    curB = (pos2 < len2) ? array2[8*pos2] : 0;
    sse_merge(&vA,&vB,&vecMin,&vecMax);
    _mm_storeu_si128((__m128i* )output,vecMin);
    output+= 8;
    while((pos1 < len1) && (pos2<len2)) {
        // code is not optimal, too much branching?
        if(curA <= curB) {
            V = _mm_lddqu_si128((const __m128i* )array1 + pos1);
            pos1++;
            curA = (pos1 < len1) ? array1[8*pos1] : 0;
        } else {
            V = _mm_lddqu_si128((const __m128i* )array2 + pos2);
            pos2++;
            curB = (pos2 < len2) ? array2[8*pos2] : 0;
        }
        sse_merge(&V,&vecMax,&vecMin,&vecMax);
        _mm_storeu_si128((__m128i* )output,vecMin);
        output+= 8;
    }
    //
    while(pos1<len1) {
        V = _mm_lddqu_si128((const __m128i* )array1 + pos1);
        pos1++;
        sse_merge(&V,&vecMax,&vecMin,&vecMax);
        _mm_storeu_si128((__m128i* )output,vecMin);
        output+= 8;
    }
    while(pos2<len2) {
        V = _mm_lddqu_si128((const __m128i* )array2 + pos2);
        pos2++;
        sse_merge(&V,&vecMax,&vecMin,&vecMax);
        _mm_storeu_si128((__m128i* )output,vecMin);
        output+=8;
    }
    _mm_storeu_si128((__m128i* )output,vecMax);
    output+= 8;
    uint32_t len = output-initoutput;
    return unique(initoutput, len);
}

// a one-pass SSE union algorithm
static uint32_t sse_unite(uint16_t * __restrict__ array1, uint32_t length1 ,
                          uint16_t *  __restrict__ array2, uint32_t length2, uint16_t *  __restrict__  output) {
    __m128i vA, vB, V, vecMin,vecMax;
    __m128i laststore;
    uint16_t * initoutput = output;
    assert(length1 % 8 == 0); // because I am lazy
    assert(length2 % 8 == 0); // because I am lazy
    assert(length1 != 0); // because I am lazy
    assert(length2 != 0); // because I am lazy
    uint32_t len1 = length1 / 8;
    uint32_t len2 = length2 / 8;
    uint32_t pos1 = 0;
    uint32_t pos2 = 0;
    // we start the machine
    vA = _mm_lddqu_si128((const __m128i* )array1 + pos1);
    pos1++;
    vB = _mm_lddqu_si128((const __m128i* )array2 + pos2);
    pos2++;
    uint16_t curA, curB;
    curA = (pos1 < len1) ? array1[8*pos1] : 0;
    curB = (pos2 < len2) ? array2[8*pos2] : 0;
    sse_merge(&vA,&vB,&vecMin,&vecMax);
    laststore = _mm_set1_epi32(array1[0] < array2[0] ? array1[0] -1 : array2[0] -1);
    output+= store_unique(laststore,vecMin, output);
    laststore=vecMin;
    while((pos1 < len1) && (pos2<len2)) {
        // code is not optimal, too much branching?
        if(curA <= curB) {
            V = _mm_lddqu_si128((const __m128i* )array1 + pos1);
            pos1++;
            curA = (pos1 < len1) ? array1[8*pos1] : 0;
        } else {
            V = _mm_lddqu_si128((const __m128i* )array2 + pos2);
            pos2++;
            curB = (pos2 < len2) ? array2[8*pos2] : 0;
        }
        sse_merge(&V,&vecMax,&vecMin,&vecMax);
        output+= store_unique(laststore,vecMin, output);
        laststore = vecMin;
    }
    while(pos1<len1) {
        V = _mm_lddqu_si128((const __m128i* )array1 + pos1);
        pos1++;
        sse_merge(&V,&vecMax,&vecMin,&vecMax);
        output+= store_unique(laststore,vecMin, output);
        laststore = vecMin;
    }
    while(pos2<len2) {
        V = _mm_lddqu_si128((const __m128i* )array2 + pos2);
        pos2++;
        sse_merge(&V,&vecMax,&vecMin,&vecMax);
        output+= store_unique(laststore,vecMin, output);
        laststore = vecMin;
    }
    output+= store_unique(laststore,vecMax,output);

    uint32_t len = output-initoutput;

    return len;
}


// standard scalar algorithm
uint32_t union_uint16(const uint16_t *set_1, uint32_t size_1, const uint16_t *set_2,
                      uint32_t size_2, uint16_t *buffer) {
    uint32_t pos = 0, idx_1 = 0, idx_2 = 0;

    if (0 == size_2) {
        memcpy(buffer, set_1, size_1 * sizeof(uint16_t));
        return size_1;
    }
    if (0 == size_1) {
        memcpy(buffer, set_2, size_2 * sizeof(uint16_t));
        return size_2;
    }

    uint16_t val_1 = set_1[idx_1], val_2 = set_2[idx_2];

    while (true) {
        if (val_1 < val_2) {
            buffer[pos++] = val_1;
            ++idx_1;
            if (idx_1 >= size_1) break;
            val_1 = set_1[idx_1];
        } else if (val_2 < val_1) {
            buffer[pos++] = val_2;
            ++idx_2;
            if (idx_2 >= size_2) break;
            val_2 = set_2[idx_2];
        } else {
            buffer[pos++] = val_1;
            ++idx_1;
            ++idx_2;
            if (idx_1 >= size_1 || idx_2 >= size_2) break;
            val_1 = set_1[idx_1];
            val_2 = set_2[idx_2];
        }
    }
    if (idx_1 < size_1) {
        const uint32_t n_elems = size_1 - idx_1;
        memcpy(buffer + pos, set_1 + idx_1, n_elems * sizeof(uint16_t));
        pos += n_elems;
    } else if (idx_2 < size_2) {
        const uint32_t n_elems = size_2 - idx_2;
        memcpy(buffer + pos, set_2 + idx_2, n_elems * sizeof(uint16_t));
        pos += n_elems;
    }
    return pos;
}

static int uint16_compare (const void * a, const void * b) {
    return ( *(uint16_t*)a - *(uint16_t*)b );
}

// a one-pass SSE union algorithm
static uint32_t sse_unite_opti(uint16_t * __restrict__ array1, uint32_t length1 ,
                               uint16_t *  __restrict__ array2, uint32_t length2, uint16_t *  __restrict__  output) {
    if ((length1 < 8) || (length2 < 8)) {
       return union_uint16(array1, length1, array2, length2,output);
    }
    __m128i vA, vB, V, vecMin,vecMax;
    __m128i laststore;
    uint16_t * initoutput = output;
    uint32_t len1 = length1 / 8;
    uint32_t len2 = length2 / 8;
    uint32_t pos1 = 0;
    uint32_t pos2 = 0;
    // we start the machine
    vA = _mm_lddqu_si128((const __m128i* )array1 + pos1);
    pos1++;
    vB = _mm_lddqu_si128((const __m128i* )array2 + pos2);
    pos2++;
    sse_merge(&vA,&vB,&vecMin,&vecMax);
    laststore = _mm_set1_epi16(-1);
    output+= store_unique(laststore,vecMin, output);
    laststore=vecMin;
    if((pos1 < len1) && (pos2<len2)) {
        uint16_t curA, curB;
        curA = array1[8*pos1];
        curB = array2[8*pos2];
        while(true) {
            if(curA <= curB) {
                V = _mm_lddqu_si128((const __m128i* )array1 + pos1);
                pos1++;
                if(pos1 < len1 ) {
                    curA =  array1[8*pos1];
                } else {
                    break;
                }
            } else {
                V = _mm_lddqu_si128((const __m128i* )array2 + pos2);
                pos2++;
                if(pos2 < len2 ) {
                    curB =  array2[8*pos2];
                } else {
                    break;
                }
            }
            sse_merge(&V,&vecMax,&vecMin,&vecMax);
            output+= store_unique(laststore,vecMin, output);
            laststore = vecMin;
        }
        sse_merge(&V,&vecMax,&vecMin,&vecMax);
        output+= store_unique(laststore,vecMin, output);
        laststore = vecMin;
    }
    // we finish the rest off using a scalar algorithm
    // copy the small end on a tmp buffer
    uint32_t len = (uint32_t) (output-initoutput);
    uint16_t buffer[16];
    uint32_t leftoversize = store_unique(laststore,vecMax, buffer);
    if(pos1 == len1) {
        memcpy(buffer + leftoversize, array1 + 8 * pos1, (length1 - 8 * len1) * sizeof(uint16_t));
        leftoversize += length1 - 8 * len1;
        qsort (buffer, leftoversize, sizeof(uint16_t), uint16_compare);


        leftoversize= unique(buffer, leftoversize);
        len += union_uint16(buffer,leftoversize, array2 + 8 * pos2,
                            length2 - 8 * pos2, output);
    } else {
        memcpy(buffer + leftoversize, array2 + 8 * pos2, (length2 - 8 * len2) * sizeof(uint16_t));
        leftoversize += length2 - 8 * len2;
        qsort (buffer, leftoversize, sizeof(uint16_t), uint16_compare);
        leftoversize= unique(buffer, leftoversize);
        len += union_uint16(buffer,leftoversize, array1 + 8 * pos1,
                            length1 - 8 * pos1, output);
    }
    return len;

}


// a one-pass SSE union algorithm
static uint32_t sse_unite_opti_pack(uint16_t * __restrict__ array1, uint32_t length1 ,
                               uint16_t *  __restrict__ array2, uint32_t length2, uint16_t *  __restrict__  output) {
    if ((length1 < 8) || (length2 < 8)) {
       return union_uint16(array1, length1, array2, length2,output);
    }
    __m128i vA, vB, V, vecMin,vecMax;
    __m128i laststore;
    uint16_t * initoutput = output;
    uint32_t len1 = length1 / 8;
    uint32_t len2 = length2 / 8;
    uint32_t pos1 = 0;
    uint32_t pos2 = 0;
    // we start the machine
    vA = _mm_lddqu_si128((const __m128i* )array1 + pos1);
    pos1++;
    vB = _mm_lddqu_si128((const __m128i* )array2 + pos2);
    pos2++;
    sse_merge(&vA,&vB,&vecMin,&vecMax);
    laststore = _mm_set1_epi16(-1);
    output+= store_unique_pack(laststore,vecMin, output);
    laststore=vecMin;
    if((pos1 < len1) && (pos2<len2)) {
        uint16_t curA, curB;
        curA = array1[8*pos1];
        curB = array2[8*pos2];
        while(true) {
            if(curA <= curB) {
                V = _mm_lddqu_si128((const __m128i* )array1 + pos1);
                pos1++;
                if(pos1 < len1 ) {
                    curA =  array1[8*pos1];
                } else {
                    break;
                }
            } else {
                V = _mm_lddqu_si128((const __m128i* )array2 + pos2);
                pos2++;
                if(pos2 < len2 ) {
                    curB =  array2[8*pos2];
                } else {
                    break;
                }
            }
            sse_merge(&V,&vecMax,&vecMin,&vecMax);
            output+= store_unique_pack(laststore,vecMin, output);
            laststore = vecMin;
        }
        sse_merge(&V,&vecMax,&vecMin,&vecMax);
        output+= store_unique_pack(laststore,vecMin, output);
        laststore = vecMin;
    }
    // we finish the rest off using a scalar algorithm
    // copy the small end on a tmp buffer
    uint32_t len = (uint32_t) (output-initoutput);
    uint16_t buffer[16];
    uint32_t leftoversize = store_unique_pack(laststore,vecMax, buffer);
    if(pos1 == len1) {
        memcpy(buffer + leftoversize, array1 + 8 * pos1, (length1 - 8 * len1) * sizeof(uint16_t));
        leftoversize += length1 - 8 * len1;
        qsort (buffer, leftoversize, sizeof(uint16_t), uint16_compare);


        leftoversize= unique(buffer, leftoversize);
        len += union_uint16(buffer,leftoversize, array2 + 8 * pos2,
                            length2 - 8 * pos2, output);
    } else {
        memcpy(buffer + leftoversize, array2 + 8 * pos2, (length2 - 8 * len2) * sizeof(uint16_t));
        leftoversize += length2 - 8 * len2;
        qsort (buffer, leftoversize, sizeof(uint16_t), uint16_compare);
        leftoversize= unique(buffer, leftoversize);
        len += union_uint16(buffer,leftoversize, array1 + 8 * pos1,
                            length1 - 8 * pos1, output);
    }
    return len;

}





/******
* Everything below this point is testing and benchmarking code.
******/




void generic_unittesting(uint32_t len1, uint32_t len2) {
    uint16_t * array1= malloc(len1 * sizeof(uint16_t));
    uint16_t * array2= malloc(len2 * sizeof(uint16_t));
    uint16_t last;
    last = 0;
    for(int k = 0; k < len1; ++k) {
        last = last + (rand()&0x0F) + 1;
        array1[k] =last;
    }
    if(!array_is_strictly_sorted(array1,len1)) return;
    last = 0;
    for(int k = 0; k < len2; ++k) {
        last = last + (rand()&0x0F) + 1;
        array2[k] = last;
    }
    if(!array_is_strictly_sorted(array2,len2)) return;

    assert(array_is_strictly_sorted(array1,len1));
    assert(array_is_strictly_sorted(array2,len2));

    uint16_t * output = malloc((len1 + len2) * sizeof(uint16_t));
    uint16_t * outputopti = malloc((len1 + len2) * sizeof(uint16_t));
    uint32_t lenopti = sse_unite_opti(array1, len1,array2, len2,outputopti);
    for(int k = 0 ; k < lenopti; ++k) {
        if(k>0) {
            if(outputopti[k]<=outputopti[k-1]) {
                printf("bug opti SSE at k =%d out of %d:  %u %u",k,lenopti,outputopti[k-1],outputopti[k]);
                abort();
            }
        }
    }

    lenopti = sse_unite_opti_pack(array1, len1,array2, len2,outputopti);
    for(int k = 0 ; k < lenopti; ++k) {
        if(k>0) {
            if(outputopti[k]<=outputopti[k-1]) {
                printf("bug opti SSE at k =%d out of %d:  %u %u",k,lenopti,outputopti[k-1],outputopti[k]);
                abort();
            }
        }
    }

    uint16_t * output2 = malloc((len1 + len2) * sizeof(uint16_t));

    uint32_t lenslow =  union_uint16(array1, len1,array2, len2,output2);
    if(lenopti != lenslow) {
        printf("bad opti length = %d %d \n", lenopti,lenslow);
        array_print(outputopti,lenopti);
        array_print(output2,lenslow);
        abort();
    }
    for(int k = 0 ; k < lenopti; ++k) {
        if(outputopti[k]!=output2[k]) {
            printf("bug opti SSE at k =%d out of %d:  %u %u",k,lenopti,outputopti[k],output2[k]);
            abort();
        }
    }
    free(output2);
    free(outputopti);
    free(array1);
    free(array2);
}


void unittesting(uint32_t len1, uint32_t len2) {
    assert(len1 % 8 == 0);
    uint16_t * array1= malloc(len1 * sizeof(uint16_t));
    assert(len2 % 8 == 0);
    uint16_t * array2= malloc(len2 * sizeof(uint16_t));
    uint16_t last;
    last = 0;
    for(int k = 0; k < len1; ++k) {
        last = last + (rand()&0x0F) + 1;
        array1[k] =last;
    }
    assert(array_is_strictly_sorted(array1,len1));

    last = 0;
    for(int k = 0; k < len2; ++k) {
        last = last + (rand()&0x0F) + 1;
        array2[k] = last;
    }
    assert(array_is_strictly_sorted(array2,len2));

    uint16_t * output = malloc((len1 + len2) * sizeof(uint16_t));
    uint32_t len = sse_unite(array1, len1,array2, len2,output);
    for(int k = 0 ; k < len; ++k) {
        if(k>0) {
            if(output[k]<=output[k-1]) {
                printf("bug SSE at k =%d out of %d:  %u %u",k,len,output[k-1],output[k]);
                abort();
            }
        }
    }
    uint16_t * outputs = malloc((len1 + len2) * sizeof(uint16_t));
    uint32_t lens = sse_unite_simple(array1, len1,array2, len2,outputs);
    for(int k = 0 ; k < lens; ++k) {
        if(k>0) {
            if(outputs[k]<=outputs[k-1]) {
                printf("bug SSE simple at k =%d out of %d:  %u %u",k,lens,outputs[k-1],outputs[k]);
                abort();
            }
        }
    }
    uint16_t * outputopti = malloc((len1 + len2) * sizeof(uint16_t));
    uint32_t lenopti = sse_unite_opti(array1, len1,array2, len2,outputopti);
    for(int k = 0 ; k < lenopti; ++k) {
        if(k>0) {
            if(outputopti[k]<=outputopti[k-1]) {
                printf("bug opti SSE at k =%d out of %d:  %u %u",k,lenopti,outputopti[k-1],outputopti[k]);
                abort();
            }
        }
    }
    lenopti = sse_unite_opti_pack(array1, len1,array2, len2,outputopti);
    for(int k = 0 ; k < lenopti; ++k) {
      if(k>0) {
         if(outputopti[k]<=outputopti[k-1]) {
            printf("bug opti SSE at k =%d out of %d:  %u %u",k,lenopti,outputopti[k-1],outputopti[k]);
            abort();
         }
      }
    }

    uint16_t * output2 = malloc((len1 + len2) * sizeof(uint16_t));

    uint32_t lenslow =  union_uint16(array1, len1,array2, len2,output2);
    if(len != lenslow) {
        printf("bad length = %d %d \n", len,lenslow);
        array_print(output,len);
        array_print(output2,lenslow);
        abort();
    }
    if(lens != lenslow) {
        printf("bad s length = %d %d \n", lens,lenslow);
        array_print(outputs,lens);
        array_print(output2,lenslow);
        abort();
    }
    if(lenopti != lenslow) {
        printf("bad opti length = %d %d \n", lenopti,lenslow);
        array_print(outputopti,lenopti);
        array_print(output2,lenslow);
        abort();
    }
    for(int k = 0 ; k < len; ++k) {
        if(output[k]!=output2[k]) {
            printf("bug SSE at k =%d out of %d:  %u %u",k,len,output[k],output2[k]);
            abort();
        }
        if(outputs[k]!=output2[k]) {
            printf("bug simple SSE at k =%d out of %d:  %u %u",k,len,outputs[k],output2[k]);
            abort();
        }
        if(outputopti[k]!=output2[k]) {
            printf("bug opti SSE at k =%d out of %d:  %u %u",k,len,outputopti[k],output2[k]);
            abort();
        }
    }
    free(output2);
    free(outputs);
    free(outputopti);
    free(output);
    free(array1);
    free(array2);
}


void benchmark(uint32_t len1, uint32_t len2) {
    printf("== benchmark sizes = %d %d \n", len1, len2);
    assert(len1 % 8 == 0);
    assert(len2 % 8 == 0);
    uint16_t * array1= malloc(len1 * sizeof(uint16_t));
    uint16_t * array2= malloc(len2 * sizeof(uint16_t));
    uint16_t last;
    last = 0;
    for(int k = 0; k < len1; ++k) {
        last = last + (rand()&0x0F) + 1;
        array1[k] =last;
    }
    assert(array_is_strictly_sorted(array1,len1));

    last = 0;
    for(int k = 0; k < len2; ++k) {
        last = last + (rand()&0x0F) + 1;
        array2[k] = last;
    }
    assert(array_is_strictly_sorted(array2,len2));

    uint16_t * output1 = malloc((len1+len2) * sizeof(uint16_t));
    uint16_t * output2 = malloc((len1+len2) * sizeof(uint16_t));
    uint16_t * output3 = malloc((len1+len2) * sizeof(uint16_t));
    uint16_t * output4 = malloc((len1+len2) * sizeof(uint16_t));

    uint64_t cycles_start = 0, cycles_final = 0;
    RDTSC_START(cycles_start);
    uint32_t lensse = sse_unite(array1, len1,array2, len2,output1);
    RDTSC_STOP(cycles_final);
    printf("SSE took %" PRIu64 " cycles\n",
           cycles_final - cycles_start);
    RDTSC_START(cycles_start);
    uint32_t optilensse = sse_unite_opti(array1, len1,array2, len2,output2);
    RDTSC_STOP(cycles_final);
    printf("opti SSE took %" PRIu64 " cycles\n",
           cycles_final - cycles_start);
    RDTSC_START(cycles_start);
    uint32_t optilenssepack = sse_unite_opti_pack(array1, len1,array2, len2,output2);
    RDTSC_STOP(cycles_final);
    printf("opti pack SSE took %" PRIu64 " cycles\n",
           cycles_final - cycles_start);

    RDTSC_START(cycles_start);
    uint32_t slensse = sse_unite_simple(array1, len1,array2, len2,output3);
    RDTSC_STOP(cycles_final);
    printf("simple SSE took %" PRIu64 " cycles\n",
           cycles_final - cycles_start);

    RDTSC_START(cycles_start);
    uint32_t lenslow =  union_uint16(array1, len1,array2, len2,output4);
    RDTSC_STOP(cycles_final);
    printf("scalar took %" PRIu64 " cycles\n",
           cycles_final - cycles_start);
    printf("lensse = %d slensse = %d lenslow = %d \n", lensse, slensse,lenslow);
    assert(lensse == lenslow);
    assert(slensse == lenslow);
    assert(optilensse == lenslow);

    assert(arrays_are_equal(output1,output4,lenslow));
    assert(arrays_are_equal(output2,output4,lenslow));
    assert(arrays_are_equal(output3,output4,lenslow));

    free(output1);
    free(output2);
    free(output3);
    free(output4);
    free(array1);
    free(array2);
}



void randomtest() {
  const int repeat = 50;
  int lenA = (rand() % 4096);
  int lenB = (rand() % 4096);
  uint16_t *A = (uint16_t *)malloc(lenA * sizeof(uint16_t));
  uint16_t *B = (uint16_t *)malloc(lenB * sizeof(uint16_t));
  if (lenA > 0) {
    A[0] = rand() % 10;
    for (int i = 1; i < lenA; i++)
      A[i] = A[i - 1] + 1 + (rand() % (50000 / lenA));
  }
  if (lenB > 0) {
    B[0] = rand() % 10;
    for (int i = 1; i < lenB; i++)
      B[i] = B[i - 1] + 1 + (rand() % (50000 / lenB));
  }
  uint16_t *buffer1 = (uint16_t *)malloc((lenB + lenA) * sizeof(uint16_t));
  uint16_t *buffer2 = (uint16_t *)malloc((lenB + lenA) * sizeof(uint16_t));

  // first side
  int f1 = union_uint16(A, lenA, B, lenB, buffer1);

  int g1 = sse_unite_opti_pack(A, lenA, B, lenB, buffer2);
  if (f1 != g1) {
    printf("A:");
    for (int k = 0; k < lenA; ++k)
      printf(" %d ", A[k]);
    printf("\n");
    printf("B:");
    for (int k = 0; k < lenB; ++k)
      printf(" %d ", B[k]);
    printf("\n");

    printf("xor:");
    for (int k = 0; k < f1; ++k)
      printf(" %d ", buffer1[k]);
    printf("\n");
    printf("simd xor:");
    for (int k = 0; k < g1; ++k)
      printf(" %d ", buffer2[k]);
    printf("\n");
  }

  assert(f1 == g1);
  for (int k = 0; k < g1; ++k)
    assert(buffer1[k] == buffer2[k]);
  BEST_TIME(union_uint16(A, lenA, B, lenB, buffer1), f1, , repeat, lenA + lenB,
            true);
  BEST_TIME(sse_unite_opti_pack(A, lenA, B, lenB, buffer2), f1, , repeat, lenA + lenB,
            true);

  // second side
  int f2 = union_uint16(B, lenB, A, lenA, buffer1);
  int g2 = sse_unite_opti_pack(B, lenB, A, lenA, buffer2);
  assert(f2 == g2);
  for (int k = 0; k < g2; ++k)
    assert(buffer1[k] == buffer2[k]);
  BEST_TIME(union_uint16(B, lenB, A, lenA, buffer1), f2, , repeat, lenA + lenB,
            true);
  BEST_TIME(sse_unite_opti_pack(B, lenB, A, lenA, buffer2), f2, , repeat, lenA + lenB,
            true);

  // cleaning
  free(A);
  free(B);
  free(buffer1);
  free(buffer2);
}

int main() {
    for(int k = 8; k <32; k+=8) {
        printf("k=%d\n",k);
        for(int j = 0; j< 1000; ++j) {
            unittesting(k, k);
        }
        for(int j = 0; j< 1000; ++j) {
            unittesting(k, 2*k);
        }
        for(int j = 0; j< 1000; ++j) {
            unittesting(2*k, k);
        }
    }
    for(int k = 0; k <512; k++) {
        for(int l = 0; l<=k; l++) {
            generic_unittesting(k, l);
            generic_unittesting(l,k);

        }
    }


    benchmark(512,512);
    benchmark(1024,512);
    benchmark(512,1024);
    benchmark(2048,512);
    for(int k = 0; k < 100; ++k) { randomtest(); printf("\n");}

}
