#ifndef MAGYAR_SORT_H
#define MAGYAR_SORT_H
/**
 * single header lib: In-place, fast heavily modified and optimized radix sort.
 *
 * Only unsigned ints for now, but should be able to modify for int and float...
 * This is the counting variant with smart changes (not per-bit).
 *
 * LICENCE: CC3 - look it up, you need to mention me but that is all
 */

#include <cstdio>
#include <cstdint>
#include <cstring> // memset

namespace MagyarSort {

	// Only change these if you know what you are doing
	// I use these because I want to see if nibbles are
	// better or something...
	//
	// Bytes of nibbles only:
	// - DIGIT_RANGE and BITS_PER_DIGIT should correspond
	// - DIGITS should also correspond with the uint32_t
	// - and DIGIT_RANGE should be 2^n value (16 or 256)
	static constexpr int DIGITS = 8; // "helyiérték"
	static constexpr int BITS_PER_DIGIT = 4; // "bit / helyiérték"
	static constexpr int DIGIT_RANGE = 16; // "helyiérték állapottér"

	template<int DIGIT_CHOICE>
	static inline uint32_t getDigit(uint32_t num) noexcept {
		static constexpr int SHIFT = DIGIT_CHOICE * BITS_PER_DIGIT;

		uint32_t shifted = num >> SHIFT;
		return shifted & (DIGIT_RANGE - 1);
	}

	/** Recursive Functor: no class should be generated I think (compiler should be smart) */
	template<int DIGIT>
	struct OccurenceMagic : public OccurenceMagic<DIGIT - 1> {
		inline OccurenceMagic(uint32_t arr[], size_t i, size_t *radicsOut) noexcept
				: OccurenceMagic<DIGIT -1 >(arr, i, radicsOut) {
			// Parents run first so template recursion runs DIGIT=0 first...
			++radicsOut[getDigit<DIGIT>(arr[i]) + DIGIT_RANGE * DIGIT];
		}
	};
	/** Ends template recursion */
	template<>
	struct OccurenceMagic<-1> {
		inline OccurenceMagic(uint32_t arr[], size_t i, size_t *radicsOut) noexcept {}
	};

	static inline void countOccurences(uint32_t arr[], size_t size, size_t *radicsOut) noexcept {
		for(size_t i = 0; i < size; ++i) {
			// Creates no object, struct is empty
			OccurenceMagic<DIGITS - 1>(arr, i, radicsOut);
		}
	}

	/** Recursive Functor: no class should be generated I think (compiler should be smart) */
	template<int DIGIT>
	struct PrefixMagic : public PrefixMagic<DIGIT - 1> {
		inline PrefixMagic(size_t *radics, size_t *prev, int i) noexcept
				: PrefixMagic<DIGIT - 1>(radics, prev, i) {
			static constexpr int DSTART = (DIGIT * DIGIT_RANGE);
			radics[DSTART + i] += prev[DIGIT];
			prev[DIGIT] = radics[DSTART + i];
		}
	};
	/** Ends template recursion */
	template<>
	struct PrefixMagic<-1> {
		inline PrefixMagic(size_t *radics, size_t *prev, int i) noexcept {}
	};

	static inline void calcPrefixSums(size_t *radics) noexcept {
		static thread_local size_t prev[DIGITS];
		memset(prev, 0, sizeof(prev));

		for(int i = 0; i < DIGIT_RANGE; ++i) {
			// This is a template-unrolled loop too
			PrefixMagic<DIGITS - 1>(radics, prev, i);
		}
	}

	void debugIt(size_t *radics) {
		for(size_t j = 0; j < DIGITS; ++j) {
			printf("d%d: ", j);
			for(size_t i = 0; i < DIGIT_RANGE; ++i) {
				printf("%d,", radics[i + DIGIT_RANGE*j]);
			}
			printf("\n\n");
		}
	}

	/** Sort the given array (in-place sorting) with the given size */
	inline void sort(uint32_t arr[], size_t size) noexcept {
		// Holds "digit" occurences, prefix sums, whatevers
		// First "DIGIT_RANGE" elem is for MSB "DIGITS", last is for LSB
		static thread_local size_t radics[DIGITS * DIGIT_RANGE];
		memset(radics, 0, sizeof(radics));

		// Calculate occurences of digits
		countOccurences(arr, size, radics);

		debugIt(radics);

		// Calculate prefix sums
		calcPrefixSums(radics);

		debugIt(radics);
	}
};

#endif