/**
  * Tento software je copyright (c) 2020 Západočeská univerzita v Plzni.
  * Užití pro potřeby AŽD Praha je dovoleno v souladu s podmínkami
  * stanovenými ve smlouvě o dílo mezi ZČU v Plzni a AŽD Praha
  * s číslem smlouvy objednatele SOD 032/20/40.
  */


#include "txcode.h"


#include <vector>
#include <memory>

#include <iostream>
#include <iomanip>
#include <algorithm>
#include <atomic>
#include <span>

#include <tbb/parallel_for.h>
#include <tbb/global_control.h>

#include <nmmintrin.h>

#ifdef AVX512
	#include <immintrin.h>
#endif

#undef min
#undef max

const char* dsSingle_Pass = "single_pass";
const char* dsTwo_Passes = "two_passes";

void Join_Vah_Limits(TVah_Limits& src_and_result, const TVah_Limits& other) {

	// as Vah_1st, VahP_1st and VahM_1st could be -1, I assume adding 1 and then subtracting 1 from resulting count is faster, than a condition

	const auto Vah_end = std::max(src_and_result.Vah_1st + 1 + src_and_result.Vah_cnt, other.Vah_1st + 1 + other.Vah_cnt);
	src_and_result.Vah_1st = std::min(src_and_result.Vah_1st, other.Vah_1st);
	src_and_result.Vah_cnt = Vah_end - src_and_result.Vah_1st - 1;

	const auto VahP_end = std::max(src_and_result.VahP_1st + 1 + src_and_result.VahP_cnt, other.VahP_1st + 1 + other.VahP_cnt);
	src_and_result.VahP_1st = std::min(src_and_result.VahP_1st, other.VahP_1st);
	src_and_result.VahP_cnt = VahP_end - src_and_result.VahP_1st - 1;

	const auto VahM_end = std::max(src_and_result.VahM_1st + 1 + src_and_result.VahM_cnt, other.VahM_1st + 1 + other.VahM_cnt);
	src_and_result.VahM_1st = std::min(src_and_result.VahM_1st, other.VahM_1st);
	src_and_result.VahM_cnt = VahM_end - src_and_result.VahM_1st - 1;
}


CUI_64_Nvah_plus_2::CUI_64_Nvah_plus_2(const size_t first_non_zero, const size_t non_zero_count) :
		mFirst_Non_Zero(first_non_zero), mSize(non_zero_count), mNon_Zero_Count(non_zero_count) {
	
	const size_t byte_size = sizeof(UI_64) * non_zero_count;
	mData = static_cast<UI_64*>(_aligned_malloc(byte_size, SIMD_Alignment));
	mOwner = mData != nullptr;
	if (mOwner) {
		memset(mData, 0, byte_size);
	}
}

/*
CUI_64_Nvah_plus_2::CUI_64_Nvah_plus_2(UI_64* data, const size_t nvah) :
	mData(data), mNvah(nvah), mOwner(false) {
	//
}
*/

CUI_64_Nvah_plus_2::~CUI_64_Nvah_plus_2() {
	if (mOwner)
		_aligned_free(mData);
}

void CUI_64_Nvah_plus_2::add(const CUI_64_Nvah_plus_2& addend) {
	/*for (size_t copy_index = 0; copy_index < mSize; copy_index++) {
		mData[copy_index] += addend.mData[copy_index];
	}
	*/

	//tohle je sice obrovsky overhead pro mala pole, ale u velkych poli to pomuze
	//je rychleji secist + to vytizi vsechna jadra, takze nebudou pretaktovana
	//a vetrak pojede na mesi obratky => mensi hluk v mistnosti
	tbb::parallel_for(static_cast<size_t>(0), mSize, [&](const size_t copy_index) {
		mData[copy_index] += addend.mData[copy_index];
	});
}

size_t CUI_64_Nvah_plus_2::size() const {
	return mSize;
}

UI_64* CUI_64_Nvah_plus_2::data() {
	return mData;
}

void CUI_64_Nvah_plus_2::reduce() {
	//fills non-zero markers and copies the non-zero part at the offset zero
	decltype(mFirst_Non_Zero) current_non_zero = 0;
	while ((current_non_zero < mNon_Zero_Count) && (mData[current_non_zero] == 0)) {
		current_non_zero++;
	}
	
	decltype(mFirst_Non_Zero) last_non_zero = mNon_Zero_Count -1;
	while ((last_non_zero > current_non_zero) && (mData[last_non_zero] == 0)) {
		last_non_zero--;
	}

	mNon_Zero_Count = last_non_zero - current_non_zero + 1;

	const auto bytes_to_copy = mNon_Zero_Count * sizeof(UI_64);
	memmove_s(mData, mSize*sizeof(UI_64), &mData[current_non_zero], bytes_to_copy);

	mFirst_Non_Zero += current_non_zero;
}

CGVah_Intermediate::CGVah_Intermediate(const TVah_Limits& limits) {
	Vah = VahP = 0;
	Coset_Leader = Ci1 = Ci2 = 0;

	found_limits = limits;

	GVah = std::make_unique<CUI_64_Nvah_plus_2>(limits.Vah_1st, limits.Vah_cnt);
	GVahM = std::make_unique<CUI_64_Nvah_plus_2>(limits.VahM_1st, limits.VahM_cnt);
	GVahP = std::make_unique<CUI_64_Nvah_plus_2>(limits.VahP_1st, limits.VahP_cnt);
}

CGVah_Intermediate::CGVah_Intermediate(const nullptr_t null) {
	Vah = VahP = 0;
	Coset_Leader = Ci1 = Ci2 = 0;

	found_limits.VahM_1st = found_limits.VahP_1st = found_limits.Vah_1st = std::numeric_limits<decltype(found_limits.VahM_cnt)>::max();
	found_limits.VahM_cnt = found_limits.VahP_cnt = found_limits.Vah_cnt = 0;

	GVah.reset();
	GVahM.reset();
	GVahP.reset();
}

void CGVah_Intermediate::add(const CGVah_Intermediate& other) {
	GVah->add(*(other.GVah.get()));
	GVahM->add(*(other.GVahM.get()));
	GVahP->add(*(other.GVahP.get()));
}

bool CGVah_Intermediate::is_valid() const {
	return  GVah && GVahM && GVahP &&				
			(GVah->data() != nullptr) &&
			(GVahM->data() != nullptr) &&
			(GVahP->data() != nullptr);
}

template <bool find_limits_no_write>
void Calculate_Initial_Coset_Leader_Values(CGVah_Intermediate& intermediate,
					  const UI_8_256 BVah, const UI_32 N1, const UI_64 Gi, const UI_64 NMarkC[2],
					  const bool do_add,
					  TVah_Limits &limits)  {	//found limits
	

	if (find_limits_no_write) {
		limits.VahM_1st = limits.VahP_1st = limits.Vah_1st = std::numeric_limits<decltype(limits.VahM_cnt)>::max();
		limits.VahM_cnt = limits.VahP_cnt = limits.Vah_cnt = 0;
	} else
		limits = intermediate.found_limits;

	UI_64 Ci1 = intermediate.Coset_Leader;	// pocatecni nastaveni obou LFSR /* pomocna a stavove prom. LFSR1/2 */
	UI_64 Ci2 = intermediate.Coset_Leader;	// pocatecni nastaveni obou LFSR /* pomocna a stavove prom. LFSR1/2 */

	UI_32 Vah = 0;
	{
		union {
			UI_64 Ci2;
			uint8_t indexes[8];
		} Ci2_index{ Ci2 };

		for (size_t i = 0; i < 8; i++)
			Vah += BVah[Ci2_index.indexes[i]];
	}

	// vypocet vahy - pro CN<=56
	UI_32 Tmp32 = 0;
	for (UI_32 i = 0; i < N1; i++)			// posun prvniho LFSR o delku kodoveho slova
	{
		const UI_64 Tmp = Ci1 & Gi;
		Tmp32 = _mm_popcnt_u64(Tmp) & 1;
		Ci1 = (Ci1 >> 1) | (NMarkC[Tmp32]);	// uprava LFSR1
		Vah += Tmp32;
	}
	

	if (find_limits_no_write) { //limity aktualizujeme jenom tehdy, kdyz na danem indexu skutecne dojde k zapisu
		const size_t vahm = static_cast<size_t>(Vah - Tmp32);
		limits.VahM_1st = std::min(limits.VahM_1st, vahm);
		limits.VahM_cnt = std::max(limits.VahM_cnt, vahm);

		limits.Vah_1st = std::min(limits.Vah_1st, static_cast<size_t>(Vah));
		limits.Vah_cnt = std::max(limits.Vah_cnt, static_cast<size_t>(Vah));
	} else {
		const size_t vahm = static_cast<size_t>(Vah - Tmp32 - limits.VahM_1st);
		if (do_add) intermediate.GVahM->data()[vahm]++;		// inkrementace položek pole vahy
		if (do_add) intermediate.GVah->data()[Vah - limits.Vah_1st]++;				// inkrementace položek pole vahy
	}

											// prvni krok obou LFSR pred cyklem
	UI_64 Tmp = Ci1 & Gi;
	Tmp32 = _mm_popcnt_u64(Tmp) & 1;
	Ci1 = (Ci1 >> 1) | NMarkC[Tmp32];		// uprava LFSR1
	UI_32 VahP = Vah + Tmp32;

	if (find_limits_no_write) {
		limits.VahP_1st = std::min(limits.VahP_1st, static_cast<size_t>(VahP));
		limits.VahP_cnt = std::max(limits.VahP_cnt, static_cast<size_t>(VahP));
	}
	else {
		if (do_add) intermediate.GVahP->data()[VahP - limits.VahP_1st]++;				// inkrementace polozek pole vahy
	}

	{
		const UI_32 common_dec = Ci2 & 1;
		VahP -= common_dec;
		Vah -= common_dec;
	}

	Tmp = Ci2 & Gi;
	Tmp32 = _mm_popcnt_u64(Tmp) & 1;
	Ci2 = (Ci2 >> 1) | NMarkC[Tmp32];		// vypocet parity slova Tmp a uprava LFSR2
	
	intermediate.Ci1 = Ci1;
	intermediate.Ci2 = Ci2;

	intermediate.Vah = Vah;
	intermediate.VahP = VahP;
}

template <bool find_limits_no_write, bool avoid_vah_offset>
void Calculate_Code_Single_Pass(//std::span<TCoset_Leader> coset_leaders, -- cilovy prekladac neumi C++20
	
	const UI_32 mPocetCL,  const UI_64 mGi, const UI_64_2 mNMarkC,
	tbb::concurrent_queue<std::unique_ptr<CGVah_Intermediate>>& intermediates,
	const TCoset_Leader* CosetLeader, const UI_8_256 BVah, const UI_32 N1) {

	//vypocitame pracovni hodnoty
	//pouzijeme UI_64_Nvah_plus_2, protoze nam to zarovna pamet, jak potrebujeme
	constexpr size_t parity_count = 256;
	CUI_64_Nvah_plus_2 parity_holder(0, parity_count);
	CUI_64_Nvah_plus_2 mark_C1_holder(0, parity_count);
	{
		UI_64* parity = parity_holder.data();
		UI_64* mark_C1 = mark_C1_holder.data();
		for (size_t i = 0; i < parity_count; i++) {
			const UI_64 par = static_cast<UI_64>(__popcnt(static_cast<unsigned int>(i))) & static_cast<UI_64>(1);
			parity[i] = par;
			mark_C1[i] = par != 0 ? mNMarkC[1] : 0;
		}
	}

	const UI_32 Pocet_CL_by_simd = (mPocetCL / static_cast<UI_32>(simd_count)) + 
								   ((mPocetCL % static_cast<UI_32>(simd_count)) != 0 ? static_cast<UI_32>(1) : static_cast<UI_32>(0));
	
	std::atomic<uint64_t> completed_leaders{ 0 };
	std::atomic<uint64_t> last_completed_leaders{ 0 };
	const uint64_t int_one_percent = std::max(static_cast<uint64_t>(1), static_cast<uint64_t>(Pocet_CL_by_simd) / static_cast<uint64_t>(100));
	const double inv_100_percent = 1.0 / static_cast<double>(Pocet_CL_by_simd);
		
	std::cout << "Spoustim SIMD verzi pro TBB..." << std::endl;
	std::cout << "Hotovo 0%... ";

	tbb::parallel_for(tbb::blocked_range<UI_32>(0, Pocet_CL_by_simd), [&completed_leaders, &last_completed_leaders, int_one_percent, inv_100_percent,
																mPocetCL, &parity_holder, &mark_C1_holder, mGi, mNMarkC,
																Pocet_CL_by_simd, &intermediates, N1, &CosetLeader, &BVah]
																(const tbb::blocked_range<UI_32> &r) {
																																
		
		
		UI_64 *parity= parity_holder.data();
		UI_64* mark_C1 = mark_C1_holder.data();		
		//presunuto sem, aby byl aligned access		
		UI_64_SIMD Gi;
		for (size_t simd_idx = 0; simd_idx < simd_count; simd_idx++) {			
			Gi[simd_idx] = mGi;
		}

		// ziskame pracovni pole presunem pointeru
		std::unique_ptr<CGVah_Intermediate> local_intermediate;
		//v podstate delame spinlock do te doby, nez ziskame pracovni pamet
		while (!intermediates.try_pop(local_intermediate)) {
			std::this_thread::yield(); //pokud by nam TBB prece jenom planovalo vice vlaken, nez kolik mame pameti
		}

		 UI_64* __restrict local_GVah = local_intermediate->GVah ? local_intermediate->GVah->data() : nullptr;
		 UI_64* __restrict local_GVahM = local_intermediate->GVahM ? local_intermediate->GVahM->data() : nullptr;
		 UI_64* __restrict local_GVahP = local_intermediate->GVahP ? local_intermediate->GVahP->data() : nullptr;
		 auto& __restrict limits = local_intermediate->found_limits;


		 UI_64_SIMD CCi;// = Copy_UI_64_SIMD(&CosetLeader[First_Coset_Leader_Idx]);
		 UI_64_SIMD Ci1;// = CCi;
		 UI_64_SIMD Ci2;// = CCi;		 

		 UI_64_SIMD Vah;
		 UI_64_SIMD VahP;

		 UI_64_SIMD Cycle_Length;

		//provedeme vychozi, inicializacni vypocet pro jednoho simd-CosetLeader
		 const size_t rend = static_cast<size_t>(r.end());
		 for (size_t id = static_cast<size_t>(r.begin()); id != rend; id++) {
			 const size_t First_Coset_Leader_Idx = id * simd_count;

			 UI_64 least_common_cycle_count = std::numeric_limits<UI_64>::max();
			 UI_64 current_cycle_iteration = 1;


			 for (size_t simd_idx = 0; simd_idx < simd_count; simd_idx++) {
				 CGVah_Intermediate& leader_init_data = *local_intermediate.get();

				 const auto effective_leader_idx = First_Coset_Leader_Idx + simd_idx;
				 const bool do_add = effective_leader_idx < mPocetCL;

				 leader_init_data.Cycle_Length = do_add ? CosetLeader[effective_leader_idx].cycle_length : 0;
				 least_common_cycle_count = std::min(least_common_cycle_count, leader_init_data.Cycle_Length);
				 Cycle_Length[simd_idx] = leader_init_data.Cycle_Length;

				 CCi[simd_idx] = do_add ? CosetLeader[effective_leader_idx].coset_leader : CosetLeader[mPocetCL - 1].coset_leader;
				 leader_init_data.Coset_Leader = CCi[simd_idx];
				 TVah_Limits simd_limits;
				 Calculate_Initial_Coset_Leader_Values<find_limits_no_write>(leader_init_data, BVah, N1, mGi, mNMarkC, do_add, simd_limits);

				 if constexpr (find_limits_no_write) {
					 if (do_add) Join_Vah_Limits(limits, simd_limits);
				 }

				 Ci1[simd_idx] = leader_init_data.Ci1;
				 Ci2[simd_idx] = leader_init_data.Ci2;
				 Vah[simd_idx] = leader_init_data.Vah;
				 VahP[simd_idx] = leader_init_data.VahP;				 
			 }


			 UI_64_SIMD VahP_offset, VahP_cnt;
			 UI_64_SIMD Vah_offset, Vah_cnt;
			 UI_64_SIMD VahM_offset, VahM_cnt;
			 for (size_t simd_idx = 0; simd_idx < simd_count; simd_idx++) {
				 VahP_offset[simd_idx] = limits.VahP_1st;
				 Vah_offset[simd_idx] = limits.Vah_1st;
				 VahM_offset[simd_idx] = limits.VahM_1st;

				 VahP_cnt[simd_idx] = limits.VahP_cnt;
				 Vah_cnt[simd_idx] = limits.Vah_cnt;
				 VahM_cnt[simd_idx] = limits.VahM_cnt;
			 }

			 bool Ci2_NEQ_CCi_all = false;
			 UI_64 Ci2_NEQ_CCi[simd_count];

			 for (size_t simd_idx = 0; simd_idx < simd_count; simd_idx++) {
				 Ci2_NEQ_CCi[simd_idx] = (Ci2[simd_idx] != CCi[simd_idx]) && (First_Coset_Leader_Idx + simd_idx < mPocetCL) ? static_cast<UI_64>(-1) : 0;
				 if (Ci2_NEQ_CCi[simd_idx] != 0) {
					 Ci2_NEQ_CCi_all = true;
				 }
			 }

			 // a ted uz spustime vlastni, dlouhotrvajici vypocet

			 while (Ci2_NEQ_CCi_all) { // cyklus posunovani obou LFSR
			 //const size_t effective_cycle_length = local_intermediate.get()->Cycle_Length;
			 //for (size_t sub_i = 1; sub_i < effective_cycle_length; sub_i++) {				 

				 UI_64_SIMD GVahM_update;
				 UI_64_SIMD GVah_update;
				 UI_64_SIMD GVahP_update;
				 UI_64_SIMD Tmp64_1, Tmp64_2, TT2;

#ifdef __INTEL_COMPILER
#pragma ivdep
#pragma vector always
#elif _MSC_VER
#pragma loop(ivdep)
#endif
				 //vektorizovatelna cast vypoctu
				 for (size_t simd_idx = 0; simd_idx < simd_count; simd_idx++) {
					 //ted jeste nebudeme zapisovat do pameti, protoze by nam to zhorsilo autovektorizaci
					 GVahM_update[simd_idx] = Vah[simd_idx];
					 GVah_update[simd_idx] = VahP[simd_idx];
					 GVahP_update[simd_idx] = VahP[simd_idx];

					 // další krok obou LFSR
					 Tmp64_1[simd_idx] = Ci1[simd_idx] & Gi[simd_idx];
					 Ci1[simd_idx] >>= 1;

					 TT2[simd_idx] = Ci2[simd_idx] & 1;
					 VahP[simd_idx] = VahP[simd_idx] - TT2[simd_idx];
					 Vah[simd_idx] = VahP[simd_idx];

					 Tmp64_2[simd_idx] = Ci2[simd_idx] & Gi[simd_idx];
					 Ci2[simd_idx] >>= 1;

					 if constexpr (!find_limits_no_write && !avoid_vah_offset) {
						 GVahM_update[simd_idx] -= VahM_offset[simd_idx]; //limits.VahM_1st;
						 GVah_update[simd_idx] -= Vah_offset[simd_idx]; //limits.Vah_1st;
						 GVahP_update[simd_idx] -= VahP_offset[simd_idx]; //limits.VahP_1st;
					 }


					 TT2[simd_idx] = Tmp64_1[simd_idx] ^ (Tmp64_1[simd_idx] >> 32);
					 TT2[simd_idx] = TT2[simd_idx] ^ (TT2[simd_idx] >> 16);
					 TT2[simd_idx] = TT2[simd_idx] ^ (TT2[simd_idx] >> 8);
					 Tmp64_1[simd_idx] = TT2[simd_idx] & 0xFF;

					 TT2[simd_idx] = Tmp64_2[simd_idx] ^ (Tmp64_2[simd_idx] >> 32);
					 TT2[simd_idx] = TT2[simd_idx] ^ (TT2[simd_idx] >> 16);
					 TT2[simd_idx] = TT2[simd_idx] ^ (TT2[simd_idx] >> 8);
					 Tmp64_2[simd_idx] = TT2[simd_idx] & 0xFF;
				 }

#ifdef AVX512
				 * reinterpret_cast<__m512i*>(&TT2) = _mm512_i64gather_epi64(reinterpret_cast<const long long*>(parity), *reinterpret_cast<__m512i*>(&Tmp64_1), 8);
				 *reinterpret_cast<__m512i*>(&Tmp64_1) = _mm512_i64gather_epi64(reinterpret_cast<const long long*>(mark_C1), *reinterpret_cast<__m512i*>(&Tmp64_1), 8);
				 *reinterpret_cast<__m512i*>(&Tmp64_2) = _mm512_i64gather_epi64(reinterpret_cast<const long long*>(mark_C1), *reinterpret_cast<__m512i*>(&Tmp64_2), 8);
#elif __AVX2__
				 * reinterpret_cast<__m256i*>(&TT2) = _mm256_i64gather_epi64(reinterpret_cast<const long long*>(parity), *reinterpret_cast<__m256i*>(&Tmp64_1), 8);
				 *reinterpret_cast<__m256i*>(&Tmp64_1) = _mm256_i64gather_epi64(reinterpret_cast<const long long*>(mark_C1), *reinterpret_cast<__m256i*>(&Tmp64_1), 8);
				 *reinterpret_cast<__m256i*>(&Tmp64_2) = _mm256_i64gather_epi64(reinterpret_cast<const long long*>(mark_C1), *reinterpret_cast<__m256i*>(&Tmp64_2), 8);
#else
				 for (size_t simd_idx = 0; simd_idx < simd_count; simd_idx++) {
					 TT2[simd_idx] = parity[Tmp64_1[simd_idx]];
					 Tmp64_1[simd_idx] = mark_C1[Tmp64_1[simd_idx]];
					 Tmp64_2[simd_idx] = mark_C1[Tmp64_2[simd_idx]];
				 }
#endif

				 for (size_t simd_idx = 0; simd_idx < simd_count; simd_idx++) {
					 VahP[simd_idx] += TT2[simd_idx];
					 GVahP_update[simd_idx] += TT2[simd_idx];

					 Ci1[simd_idx] |= Tmp64_1[simd_idx];
					 Ci2[simd_idx] |= Tmp64_2[simd_idx];
				 }


				 //vektorizovatelna cast s manualne vektorizovanym prologem
				 if constexpr (find_limits_no_write) {
					 UI_64_SIMD m_1, m_c;
					 UI_64_SIMD v_1, v_c;
					 UI_64_SIMD p_1, p_c;

					 //Pri pouziti skalarniho if, bychom jeden limit porovnali na 2 pomale nebo 3 rychlejsi instrukce
					 //Pro AVX2, je to 8 nebo 12 instrukci. Vektorove to udelame na 4 instrukce AVX2.
					 //U AVX512 je to tak 16-24 vs 5 instrukci. U SSE by to bylo nastejno.
					 //vysvetleni nezahrnuje rezii kontroly Ci2 != CCi, ktera je u vektorove verze jen 1 instrukce navic

					 //Nejprve vypocteme masku rozdilu podle porovnani soucasnych hodnot
#ifdef __AVX2__
					 * reinterpret_cast<__m256i*>(&m_1[0]) = _mm256_cmpgt_epi64(*reinterpret_cast<__m256i*>(&VahM_offset[0]), *reinterpret_cast<__m256i*>(&GVahM_update[0]));
					 *reinterpret_cast<__m256i*>(&v_1[0]) = _mm256_cmpgt_epi64(*reinterpret_cast<__m256i*>(&Vah_offset[0]), *reinterpret_cast<__m256i*>(&GVah_update[0]));
					 *reinterpret_cast<__m256i*>(&p_1[0]) = _mm256_cmpgt_epi64(*reinterpret_cast<__m256i*>(&VahP_offset[0]), *reinterpret_cast<__m256i*>(&GVahP_update[0]));

					 *reinterpret_cast<__m256i*>(&m_c[0]) = _mm256_cmpgt_epi64(*reinterpret_cast<__m256i*>(&GVahM_update[0]), *reinterpret_cast<__m256i*>(&VahM_cnt[0]));
					 *reinterpret_cast<__m256i*>(&v_c[0]) = _mm256_cmpgt_epi64(*reinterpret_cast<__m256i*>(&GVah_update[0]), *reinterpret_cast<__m256i*>(&Vah_cnt[0]));
					 *reinterpret_cast<__m256i*>(&p_c[0]) = _mm256_cmpgt_epi64(*reinterpret_cast<__m256i*>(&GVahP_update[0]), *reinterpret_cast<__m256i*>(&VahP_cnt[0]));

#ifdef AVX512
					 * reinterpret_cast<__m256i*>(&m_1[4]) = _mm256_cmpgt_epi64(*reinterpret_cast<__m256i*>(&VahM_offset[4]), *reinterpret_cast<__m256i*>(&GVahM_update[4]));
					 *reinterpret_cast<__m256i*>(&v_1[4]) = _mm256_cmpgt_epi64(*reinterpret_cast<__m256i*>(&Vah_offset[4]), *reinterpret_cast<__m256i*>(&GVah_update[0]));
					 *reinterpret_cast<__m256i*>(&p_1[4]) = _mm256_cmpgt_epi64(*reinterpret_cast<__m256i*>(&VahP_offset[4]), *reinterpret_cast<__m256i*>(&GVahP_update[4]));


					 *reinterpret_cast<__m256i*>(&m_c[4]) = _mm256_cmpgt_epi64(*reinterpret_cast<__m256i*>(&GVahM_update[4]), *reinterpret_cast<__m256i*>(&VahM_cnt[4]));
					 *reinterpret_cast<__m256i*>(&v_c[4]) = _mm256_cmpgt_epi64(*reinterpret_cast<__m256i*>(&GVah_update[4]), *reinterpret_cast<__m256i*>(&Vah_cnt[4]));
					 *reinterpret_cast<__m256i*>(&p_c[4]) = _mm256_cmpgt_epi64(*reinterpret_cast<__m256i*>(&GVahP_update[4]), *reinterpret_cast<__m256i*>(&VahP_cnt[4]));
#endif
#else

					 for (size_t simd_idx = 0; simd_idx < simd_count; simd_idx++) {
						 m_1[simd_idx] = VahM_offset[simd_idx] > GVahM_update[simd_idx] ? static_cast<UI_64>(-1) : 0;
						 v_1[simd_idx] = Vah_offset[simd_idx] > GVah_update[simd_idx] ? static_cast<UI_64>(-1) : 0;
						 p_1[simd_idx] = VahP_offset[simd_idx] > GVahP_update[simd_idx] ? static_cast<UI_64>(-1) : 0;

						 m_c[simd_idx] = GVahM_update[simd_idx] > VahM_cnt[simd_idx] ? static_cast<UI_64>(-1) : 0;
						 v_c[simd_idx] = GVah_update[simd_idx] > Vah_cnt[simd_idx] ? static_cast<UI_64>(-1) : 0;
						 p_c[simd_idx] = GVahP_update[simd_idx] > VahP_cnt[simd_idx] ? static_cast<UI_64>(-1) : 0;
					 }
#endif

					 for (size_t simd_idx = 0; simd_idx < simd_count; simd_idx++) {
						 VahM_offset[simd_idx] -= m_1[simd_idx] & (VahM_offset[simd_idx] - GVahM_update[simd_idx]) & Ci2_NEQ_CCi[simd_idx];
						 Vah_offset[simd_idx] -= v_1[simd_idx] & (Vah_offset[simd_idx] - GVah_update[simd_idx]) & Ci2_NEQ_CCi[simd_idx];
						 VahP_offset[simd_idx] -= p_1[simd_idx] & (VahP_offset[simd_idx] - GVahP_update[simd_idx]) & Ci2_NEQ_CCi[simd_idx];

						 VahM_cnt[simd_idx] += m_c[simd_idx] & (GVahM_update[simd_idx] - VahM_cnt[simd_idx]) & Ci2_NEQ_CCi[simd_idx];
						 Vah_cnt[simd_idx] += v_c[simd_idx] & (GVah_update[simd_idx] - Vah_cnt[simd_idx]) & Ci2_NEQ_CCi[simd_idx];
						 VahP_cnt[simd_idx] += p_c[simd_idx] & (GVahP_update[simd_idx] - VahP_cnt[simd_idx]) & Ci2_NEQ_CCi[simd_idx];
					 }
				 }

				 //ted budeme inkrementovat pole vah a zaroven urcime
				 //shodnost Ci2 != CCi



				 if constexpr (!find_limits_no_write) {
#ifdef __INTEL_COMPILER	
#pragma novector
#elif _MSC_VER
#pragma loop(no_vector)
#endif
					 for (size_t simd_idx = 0; simd_idx < simd_count; simd_idx++) {
						 //protoze budeme testovat jenom ty, ktere jsou true
						  //uz nemusime znovu testovat na PocetCL
						 if (Ci2_NEQ_CCi[simd_idx] != 0) {
							 local_GVahM[GVahM_update[simd_idx]]++;
							 local_GVah[GVah_update[simd_idx]]++;
							 local_GVahP[GVahP_update[simd_idx]]++;
						 }
					 }
				 }


				 current_cycle_iteration++;
				 if (current_cycle_iteration >= least_common_cycle_count) {
					 //CPU se snadno nauci tuto podminku, ktera temer nikdy nebude splnena, az na par specialnich pripadu a nekolik malo zaverecnych kroku
					 //=> nasledujici instrukce se prakticky temer nebudou vykonavat a dtto branch missprediction 
					 
					 //pokud bychom nasledujici chteli vektorizovat, tak to nema smysl, protoze z predchozi implementace lze dovodit, ze nasledujici instrukce bychom stejne museli udelat
					 //respektive, navic mame jenom jeneden skok - else
					 //ale usetrime 4-5 vektorovych instrukci podle AVX2-AVX512

					 //tudiz navic, pro vyhodnoceni podminky tohoto bloku, mame 1xadd a 1xif, coz je stale mene instrukci nez puvodni vektorizovana verze
					 //=> hlavni prinos paralelizace je tak moznost rozdelit jednoho coset leadera na nekolik mensich
					 Ci2_NEQ_CCi_all = false;
					 for (size_t simd_idx = 0; simd_idx < simd_count; simd_idx++) {

						 if (current_cycle_iteration >= Cycle_Length[simd_idx]) {
							 Ci2_NEQ_CCi[simd_idx] = 0;
						 } else
							 Ci2_NEQ_CCi_all = true;
					 }
				 }


			 }			 

			 //najdeme globalni extremy
			 if constexpr (find_limits_no_write) {
				 for (size_t simd_idx = 0; simd_idx < simd_count; simd_idx++) {
					limits.VahM_1st = std::min(limits.VahM_1st, VahM_offset[simd_idx]);
					limits.VahM_cnt = std::max(limits.VahM_cnt, VahM_cnt[simd_idx]);

					limits.Vah_1st = std::min(limits.Vah_1st, Vah_offset[simd_idx]);
					limits.Vah_cnt = std::max(limits.Vah_cnt, Vah_cnt[simd_idx]);

					limits.VahP_1st = std::min(limits.VahP_1st, VahP_offset[simd_idx]);
					limits.VahP_cnt = std::max(limits.VahP_cnt, VahP_cnt[simd_idx]);
				 }
			 }

			 //Nyni je vhodny bod pro fail-safe, kdy se aktualni stav - leader ulozi na disk

			 completed_leaders++;

			 const auto ccl = completed_leaders.load();
			 const auto lcl = last_completed_leaders.load();

			 if ((ccl - lcl) > int_one_percent) {
				 double completed_percents = static_cast<double>(ccl) * inv_100_percent;
				 completed_percents = std::trunc(completed_percents * 1000.0);
				 completed_percents *= 0.1;
				 completed_percents = std::min(completed_percents, 100.0);


				 std::cout << completed_percents << "%... " << std::flush;

				 last_completed_leaders.store(ccl);
			 }

		 }

		  intermediates.push(std::move(local_intermediate));
	});

	std::cout << "100%" << std::endl;
}

CCalculator::CCalculator(const NMethod method, const UI_32 PocetCL, const UI_64 Gi, const UI_64_2 NMarkC)
	: mMethod(method), mPocetCL(PocetCL), mGi(Gi), mNMarkC{ NMarkC[0], NMarkC[1] } {
	//
}

tbb::concurrent_queue<std::unique_ptr<CGVah_Intermediate>> CCalculator::Allocate_Intermediate_Memory(const TVah_Limits &limits) const {

	tbb::concurrent_queue<std::unique_ptr<CGVah_Intermediate>> result;

	std::cout << "Alokuji pamet pro procesorova jadra..." << std::endl;

	size_t effective_concurrency = 0;
	for (; effective_concurrency < std::thread::hardware_concurrency(); effective_concurrency++) {
		std::unique_ptr<CGVah_Intermediate> mem = std::make_unique<CGVah_Intermediate>(limits);
		if (!mem->is_valid()) break;

		result.push(std::move(mem));
	}

	if (effective_concurrency >= std::thread::hardware_concurrency()) std::cout << "Uspesne alokovana pamet pro vsechny procesorova jadra." << std::endl;
	else {
		std::cout << "Podarilo se alokovat pamet jenom pro " << effective_concurrency << " procesorovych jader z celkovych " << std::thread::hardware_concurrency() << ".\n";
	}

	return result;
};


TVah_Limits CCalculator::Estimate_Range(const size_t desired_nvah,
										const TCoset_Leader* CosetLeader, const UI_8_256 BVah, const UI_32 N1) const {
	TVah_Limits result;

	auto fill_default_range = [&]() {
		result.VahM_1st = result.VahP_1st = result.Vah_1st = 0;
		result.VahM_cnt = result.VahP_cnt = result.Vah_cnt = desired_nvah+2;	//+2 is given by the algorithm
	};

	auto estimate_range = [&]() {
		tbb::concurrent_queue<std::unique_ptr<CGVah_Intermediate>> partial_limits;
		size_t effective_concurrency = 0;
		for (; effective_concurrency < std::thread::hardware_concurrency(); effective_concurrency++) {
			std::unique_ptr<CGVah_Intermediate> mem = std::make_unique<CGVah_Intermediate>(nullptr);
			partial_limits.push(std::move(mem));
		}

		tbb::global_control control{ tbb::global_control::max_allowed_parallelism, partial_limits.unsafe_size() };

		Calculate_Code_Single_Pass<true, true>(mPocetCL, mGi, mNMarkC, partial_limits, CosetLeader, BVah, N1);

		// a na zaver secteme vsechny mezivysledky
		std::unique_ptr<CGVah_Intermediate> default_intermediate, local_intermediate;
		partial_limits.try_pop(default_intermediate);	//prvni vysledek
		result = default_intermediate->found_limits;

		while (partial_limits.try_pop(local_intermediate))
			Join_Vah_Limits(result, local_intermediate->found_limits);

		//and convert ends to counts
		result.VahM_cnt -= result.VahM_1st - 1;
		result.Vah_cnt -= result.Vah_1st - 1;
		result.VahP_cnt -= result.VahP_1st - 1;
	};

	switch (mMethod) {
		case NMethod::two_passes:
			estimate_range();
			break;

		default:
			fill_default_range();
			break;

	}


	return result;
}

CGVah_Intermediate CCalculator::Calculate_Code(const TVah_Limits& limits,
	const TCoset_Leader* CosetLeader, const UI_8_256 BVah, const UI_32 N1) const {
	//alokujeme pracovni pamet
	auto intermediates = Allocate_Intermediate_Memory(limits);

	auto compute_only = [&] {
		CGVah_Intermediate result{ limits };
		if ((limits.VahM_1st == 0) && (limits.VahP_1st == 0) && (limits.Vah_1st == 0))
		Calculate_Code_Single_Pass<false, true>(mPocetCL, mGi, mNMarkC, intermediates, CosetLeader, BVah, N1);
		else Calculate_Code_Single_Pass<false, false>(mPocetCL, mGi, mNMarkC, intermediates, CosetLeader, BVah, N1);

		// a na zaver secteme vsechny mezivysledky
		std::unique_ptr<CGVah_Intermediate> local_intermediate;
		while (intermediates.try_pop(local_intermediate))
			result.add(*local_intermediate.get());

	

		return result;
	};

	auto compute_and_reduce = [&] {
		auto result = compute_only();

		result.GVahM->reduce();
		result.GVah->reduce();
		result.GVahP->reduce();		

		return result;
	};


#ifdef DDegrade_To_Serial
	const size_t effective_parallelism = 1;
#else
	const size_t effective_parallelism = intermediates.unsafe_size();	//we do no parallel op on it right now, so it is actually a safe size
#endif
	std::cout << "Omezuji planovac TBB na " << effective_parallelism << " jader z celkovych " << std::thread::hardware_concurrency() << "." << std::endl;

	//nastavime max. pocet vlaken, ktere nam ma TBB planovat 
	tbb::global_control control{ tbb::global_control::max_allowed_parallelism, effective_parallelism };
	
	CGVah_Intermediate result{ limits };
	if ((limits.VahM_1st == 0) && (limits.VahP_1st == 0) && (limits.Vah_1st == 0))
		Calculate_Code_Single_Pass<false, true>(mPocetCL, mGi, mNMarkC, intermediates, CosetLeader, BVah, N1);
	else Calculate_Code_Single_Pass<false, false>(mPocetCL, mGi, mNMarkC, intermediates, CosetLeader, BVah, N1);

	// a na zaver secteme vsechny mezivysledky
	std::unique_ptr<CGVah_Intermediate> local_intermediate;
	while (intermediates.try_pop(local_intermediate))
		result.add(*local_intermediate.get());

	//auto result = mMethod == NMethod::single_pass ? compute_and_reduce() : compute_only();	

	result.GVahM->reduce();
	result.GVah->reduce();
	result.GVahP->reduce();

	return result;
}
