/*
								+--------------------------------+
								|                                |
								|*** Elementary combinatorics ***|
								|                                |
								|  Copyright © -tHE SWINe- 2009  |
								|                                |
								|           Combinat.h           |
								|                                |
								+--------------------------------+
*/

#pragma once
#ifndef __COMBINAT_INCLUDED
#define __COMBINAT_INCLUDED

/**
 *	@file Combinat.h
 *	@brief elementary combinatorics methods
 *	@date 2009
 *	@author -tHE SWINe-
 */

#include "MinMax.h"
#include "Integer.h"
#include <algorithm>
#include <vector>
#include <set>
#ifdef _OPENMP
#include <omp.h>
#endif // _OPENMP

/**
 *	@def __PERMUTATION_LSB_FIRST
 *	@brief if defined, CPermutation and CSeekablePermutation store the least significant
 *		item under index 0 (otherwise there's the most significant item)
 *
 *	If defined, CPermutation and CSeekablePermutation store items in such a manner
 *		that if the first item (index 0) is interpreted as the least significant and
 *		the last item as the most significant, the generated values (when calling Next())
 *		form monotonically increasing sequence.
 *
 *	@note This has no effect on CMixedRadixInteger or CCombination, which are always "LSB first".
 */
#define __PERMUTATION_LSB_FIRST

/**
 *	@brief a simple combination enumeration class
 *
 *	A simplie class that can be used to enumerate combination (of a set of n contiguous
 *	integers, starting with 0).
 *
 *	The documentation refers to "the set", which is the set of items being permuted.
 *	The set contains "items", which are represented by (zero-based) "indices".
 */
class CCombination {
protected:
	size_t m_n_size; /**< @brief set size */
	std::vector<size_t> m_combo; /**< @brief the combination */

public:
	/**
	 *	@brief default constructor; generates the first combination
	 *
	 *	@param[in] n_set_size is number of items in a set to choose from
	 *	@param[in] n_combination_size is number of items to choose (number of outputs)
	 *
	 *	@note This function can throw the std::bad_alloc exception.
	 */
	CCombination(size_t n_set_size, size_t n_combination_size); // throw(std::bad_alloc)

	/**
	 *	@brief gets size of the set the items are chosen from
	 *	@return Returns size of the set the items are chosen from (n in "n choose k").
	 */
	inline size_t n_Set_Size() const
	{
		return m_n_size;
	}

	/**
	 *	@brief gets the number of items that are to be chosen
	 *	@return Returns the number of items that are to be chosen
	 *		(k in "n choose k", number of outputs).
	 */
	inline size_t n_Combination_Size() const
	{
		return m_combo.size();
	}

	/**
	 *	@brief calculates the number of different combinations
	 *	@return Returns the number of different combinations, or 0 on overflow
	 *		(the value is larger than UINT64_MAX).
	 */
	uint64_t n_Combination_Num() const;

	/**
	 *	@brief gets the current combination
	 *	@return Returns the list of (zero-based) indices of the chosen items.
	 */
	inline const size_t *p_Indices() const
	{
		return &m_combo.front();
	}

	/**
	 *	@brief gets a single item from the current combination
	 *	@param[in] n_index is zero-based index of the item in combination (0 to n_Combination_Size() - 1)
	 *	@return Returns (zero-based) indicex of the chosen item.
	 */
	inline size_t operator [](size_t n_index) const
	{
		_ASSERTE(n_index < m_combo.size());
		return m_combo[n_index];
	}

	/**
	 *	@brief rewinds to the first combination
	 */
	void First();

	/**
	 *	@brief generates the next combination in lexicographic order
	 *	@return Returns true on success, false on failure (the combination
	 *		remains unchanged in such case).
	 *	@note The combination is automatically rewound to First() after reaching
	 *		the last one. As a consequence, this function returns false only once
	 *		in a while.
	 */
	bool Next();
};

/**
 *	@brief mixed radix positional integer representation with conversion to native values
 *	@note This can be used to enumerate combinations in a simple and clean manner.
 */
class CMixedRadixInteger {
protected:
	std::vector<size_t> m_comb; /**< @brief list of the individual digits */
	std::vector<size_t> m_base; /**< @brief list of radix values per each position */

public:
	/**
	 *	@brief default constructor; initializes the value to zero
	 *
	 *	@param[in] p_radix_list is list of n_radix_num radixes, the first being applied
	 *		to the first (the least significant) digit
	 *	@param[in] n_radix_num is number of digits
	 *
	 *	@note This function can throw the std::bad_alloc exception.
	 */
	CMixedRadixInteger(const size_t *p_radix_list, size_t n_radix_num); // throws(std::bad_alloc)

	/**
	 *	@brief calculates the maximal representable value
	 *	@return Returns the maximal representable value, or 0 on overflow
	 *		(the value is larger than UINT64_MAX).
	 */
	inline uint64_t n_Max_Value() const
	{
		uint64_t n_max_value;
		bool b_result = Get_Max_Value(n_max_value);
		//_ASSERTE(b_result); // 0 can be max value too, in extreme case. better assert this.
		return (b_result)? n_max_value : 0;
	}

	/**
	 *	@brief calculates the maximal representable value
	 *	@param[out] r_n_max_value is the maximal representable value
	 *	@return Returns true if the value fit in 64-bit unsigned integer,
	 *		otherwise returns false.
	 */
	bool Get_Max_Value(uint64_t &r_n_max_value) const;

	/**
	 *	@brief gets the number of digits
	 *	@return Returns the number of digits. 
	 */
	inline const size_t n_Digit_Num() const
	{
		return m_comb.size();
	}

	/**
	 *	@brief gets a specified number digit
	 *	@param[in] n_digit_index is zero-based index of the digit to retrieve
	 *	@return Returns number digit specified by the index.
	 */
	inline const size_t operator [](size_t n_digit_index) const
	{
		return m_comb[n_digit_index];
	}

	/**
	 *	@brief gets number digits
	 *	@return Returns list of n_Digit_Num() number digits,
	 *		the first-one being the least significant.
	 */
	inline const size_t *p_Digits() const
	{
		return (m_comb.empty())? 0 : &m_comb.front();
	}

	/**
	 *	@brief gets number radixes
	 *	@return Returns list of n_Digit_Num() number of radixes,
	 *		the first-one corresponding to the least significant digit.
	 */
	inline const size_t *p_Radixes() const
	{
		return (m_base.empty())? 0 : &m_base.front();
	}

	/**
	 *	@brief calculates the value as a native integer
	 *	@return Returns the value as a native integer, or 0 on overflow
	 *		(the value is larger than UINT64_MAX).
	 */
	uint64_t n_Value() const;

	/**
	 *	@brief sets the value
	 *	@param[in] n_value is the value to set
	 *	@return Returns true on success, false on failure
	 *		(the value was greater than n_Max_Value()).
	 */
	bool SetValue(uint64_t n_value);

	/**
	 *	@brief increments the number
	 *	@return Returns true on success, false on failure
	 *		(overflow, the value is now 0).
	 */
	bool Increment();

	/**
	 *	@brief decrements the number
	 *	@return Returns true on success, false on failure
	 *		(underflow, the value is now n_Max_Value()).
	 */
	bool Decrement();

	/**
	 *	@copydoc CMixedRadixInteger::Increment()
	 */
	inline bool operator ++()
	{
		return Increment();
	}

	/**
	 *	@copydoc CMixedRadixInteger::Decrement()
	 */
	inline bool operator --()
	{
		return Decrement();
	}
};

/**
 *	@brief a simple permutation enumeration class
 *
 *	A simplie class that can be used to enumerate permutations (of a set of n contiguous
 *	integers, starting with 0).
 *
 *	The documentation refers to "the set", which is the set of items being permuted.
 *	The set contains "items", which are represented by (zero-based) "indices".
 *
 *	@note This uses the Narayana Pandita algorithm, which is efficient in case
 *		the permutation size is not too much smaller than the set size. Otherwise
 *		finding the next permutation takes (n_set_size - n_permutation_size)! steps,
 *		and it is probably better to use CSeekablePermutation which will be
 *		more efficient in such cases.
 *	@note The set size is limited to 22 items as 22! is more than UINT64_MAX,
 *		CSeekablePermutation can be used to enumerate permutations of sets of
 *		more than 22 items.
 *	@note The order and values of generated permutations match the ones of CSeekablePermutation
 *		(given the same input paramters).
 *	@note This cannot seek to a particular permutation, but CSeekablePermutation
 *		can seek and it is in turn possible to copy item indices to CPermutation
 *		and iterate to the next permutations. The following code snippet illustrates
 *		such use:
 *
 *	@code
 *	int n = 5, k = 4;
 *	// n choose k (selecting subset of 4 different items from set of 5 items)
 *
 *	CPermutation perm(n, k);
 *	// this gives the first permutation (items 0, 1, 2, 3)
 *
 *	CSeekablePermutation perm_seek(n, k);
 *	perm_seek.Seek(3); // skip to the fourth permutation (the index is zero-based)
 *	perm.Set_Permutation(perm_seek.p_Indices()); // copy the permutation to perm
 *	// seek to the fourth permutation using CSeekablePermutation
 *
 *	do {
 *		printf("%d %d %d %d\n", perm[0], perm[1], perm[2], perm[3]);
 *	} while(perm.Next());
 *	// can efficiently use CPermutation::Next() to enumerate the following permutations
 *	@endcode
 */
class CPermutation {
protected:
	size_t m_n_size; /**< @brief permutation size */
	size_t m_n_off; /**< @brief set size - permutation size */
	std::vector<size_t> m_perm; /**< @brief the current permutation (contains set size elements) */
	uint64_t m_n_skip_ratio; /**< @brief (set size - permutation size)! */

public:
	/**
	 *	@brief default constructor; generates the first permutation
	 *
	 *	@param[in] n_set_size is number of items in a set to choose from
	 *	@param[in] n_permutation_size is number of items to choose (number of outputs)
	 *
	 *	@note This function can throw the std::bad_alloc exception.
	 *	@note This is only good for set sizes below 22. Use CSeekablePermutation for larger sets.
	 */
	CPermutation(size_t n_set_size, size_t n_permutation_size); // throws(std::bad_alloc)

	/**
	 *	@brief gets size of the set the items are chosen from
	 *	@return Returns size of the set the items are chosen from (n in "n choose k").
	 */
	inline size_t n_Set_Size() const
	{
		return m_perm.size();
	}

	/**
	 *	@brief gets the number of items that are to be chosen
	 *	@return Returns the number of items that are to be chosen
	 *		(k in "n choose k", number of outputs).
	 */
	inline size_t n_Permutation_Size() const
	{
		return m_n_size;
	}

	/**
	 *	@brief sets a particular permutation as current
	 *	@param[in] p_index_list is the list of n_Permutation_Size() item indices
	 *		(a permutation; the indices must be some of integers from 0 to n_Set_Size() - 1
	 *		and none of them shall appear more than once)
	 *	@return Returns true on success, false on failure (the sequence is either malformed or
	 *		it contains indices of items that do not belong to the set the items are chosen from).
	 *	@note Calling Next() will generate the next permutation (in lexicographic order)
	 *		to the one being set here (in case such permutation does exist).
	 */
	bool Set_Permutation(const size_t *p_index_list);

	/**
	 *	@brief gets the current permutation
	 *	@return Returns the list of (zero-based) indices of the chosen items.
	 *	@note The ordering of the returned indices can be controlled
	 *		using __PERMUTATION_LSB_FIRST.
	 */
	inline const size_t *p_Indices() const
	{
#ifdef __PERMUTATION_LSB_FIRST
		return (m_perm.empty())? 0 : &m_perm[m_n_off];
#else // __PERMUTATION_LSB_FIRST
		return (m_perm.empty())? 0 : &m_perm.front();
#endif // __PERMUTATION_LSB_FIRST
	}

	/**
	 *	@brief gets a single item from the current permutation
	 *	@param[in] n_index is zero-based index of the item in permutation (0 to n_Permutation_Size() - 1)
	 *	@return Returns (zero-based) indicex of the chosen item.
	 */
	inline size_t operator [](size_t n_index) const
	{
		_ASSERTE(n_index < m_n_size);
#ifdef __PERMUTATION_LSB_FIRST
		n_index += m_n_off;
#endif // __PERMUTATION_LSB_FIRST
		return m_perm[n_index];
	}

	/**
	 *	@brief calculates the number of different permutations
	 *	@return Returns the number of different permutations, or 0 on overflow
	 *		(the value is larger than UINT64_MAX).
	 */
	inline uint64_t n_Permutation_Num() const
	{
		uint64_t n_permutation_num;
		bool b_result = Get_Permutation_Num(n_permutation_num);
		return (b_result)? n_permutation_num : 0;
	}

	/**
	 *	@brief calculates the number of different permutations
	 *	@param[out] r_n_permutation_num is the number of different permutations
	 *	@return Returns true if the value fit in 64-bit unsigned integer,
	 *		otherwise returns false.
	 */
	bool Get_Permutation_Num(uint64_t &r_n_permutation_num) const;

	/**
	 *	@brief rewinds to the first permutation
	 */
	void First();

	/**
	 *	@brief generates the next permutation in lexicographic order
	 *
	 *	@return Returns true on success, false on failure (the permutation
	 *		remains unchanged in such case).
	 *
	 *	@note The order and values of generated permutations match the ones of CSeekablePermutation
	 *		(given the same input paramters).
	 *	@note Unlike CSeekablePermutation, calling Next() on the last permutation
	 *		fails and the permutation is not automatically changed to the first permutaiton.
	 *		Therefore after the enumeration finished, calling Next() always fails and
	 *		has no effect, until First() is called.
	 */
	bool Next();
};

/**
 *	@brief a simple permutation enumeration class
 *
 *	A simple class that can be used to enumerate permutations (of a set of n contiguous
 *	integers, starting with 0).
 *
 *	The documentation refers to "the set", which is the set of items being permuted.
 *	The set contains "items", which are represented by (zero-based) "indices".
 *
 *	@note This uses a factorial number system and Lehmer code to generate permutations.
 *		While this enables seeking to any particular permutation, the algorithm is slow.
 *		In some situations it may be better to use CPermutation (see it's class documentation).
 *		This implementation is particularly efficient in cases when the permutation size
 *		is significantly smaller than the set size.
 *	@note The order and values of generated permutations match the ones of CPermutation
 *		(given the same input paramters).
 */
class CSeekablePermutation {
private:
	size_t m_n_size; /**< @brief permutation size */
	size_t m_n_off; /**< @brief set size - permutation size */
	std::vector<size_t> m_perm; /**< @brief the permutation (size is always set size) */
	std::vector<size_t> m_count; /**< @brief Lehmer code */

public:
	/**
	 *	@copydoc CPermutation::CPermutation()
	 */
	CSeekablePermutation(size_t n_set_size, size_t n_permutation_size); // throws(std::bad_alloc)

	/**
	 *	@copydoc CPermutation::n_Set_Size()
	 */
	inline size_t n_Set_Size() const
	{
		return m_perm.size();
	}

	/**
	 *	@copydoc CPermutation::n_Permutation_Size()
	 */
	inline size_t n_Permutation_Size() const
	{
		return m_n_size;
	}

	/**
	 *	@copydoc CPermutation::p_Indices()
	 */
	inline const size_t *p_Indices() const
	{
#ifdef __PERMUTATION_LSB_FIRST
		_ASSERTE(m_perm.size() >= m_n_off);
		return (m_perm.size() == m_n_off)? 0 : &m_perm[m_n_off];
#else // __PERMUTATION_LSB_FIRST
		return (m_perm.empty())? 0 : &m_perm.front();
#endif // __PERMUTATION_LSB_FIRST
	}

	/**
	 *	@copydoc CPermutation::operator[]()
	 */
	inline size_t operator [](size_t n_index) const
	{
		_ASSERTE(n_index < m_n_size);
#ifdef __PERMUTATION_LSB_FIRST
		n_index += m_n_off;
#endif // __PERMUTATION_LSB_FIRST
		return m_perm[n_index];
	}

	/**
	 *	@copydoc CPermutation::n_Permutation_Num()
	 */
	inline uint64_t n_Permutation_Num() const
	{
		uint64_t n_permutation_num;
		bool b_result = Get_Permutation_Num(n_permutation_num);
		return (b_result)? n_permutation_num : 0;
	}

	/**
	 *	@copydoc CPermutation::Get_Permutation_Num()
	 */
	bool Get_Permutation_Num(uint64_t &r_n_permutation_num) const;

	/**
	 *	@brief seeks to a particular permutation
	 *	@param[in] n_permutation_index is zero-based index of the permutation to seek to
	 *		(counted from the first permutation, in increasing lexicographical order)
	 *	@return Returns true on success, false on failure (n_permutation_index is greater
	 *		than (or equal to) the number of permutations).
	 */
	bool Seek(uint64_t n_permutation_index);

	/**
	 *	@copydoc CPermutation::First()
	 */
	void First();

	/**
	 *	@brief generates the next permutation in lexicographic order
	 *	@return Returns true on success, false on failure (the permutation
	 *		returns to the first permutation in such case, the next call to Next()
	 *		succeeds and generates the second permutation).
	 *	@note The order and values of generated permutations match the ones of CSeekablePermutation
	 *		(given the same input paramters).
	 */
	bool Next();

protected:
	void FinishPermutation();
};

/**
 *	@brief a naive unordered partition enumerator
 *
 *	This enumerates ways of distributing N items into P bins,
 *	where the order of the items or bins is not significant.
 */
class CUnorderedPartition {
protected:
	size_t m_n_item_num; /**< @brief number of items to be partitioned */
	size_t m_n_bin_num; /**< @brief number of bins */
	std::vector<size_t> m_counters; /**< @brief internal counters, representing the partition */
	std::vector<size_t> m_bin_values; /**< @brief bin states */

public:
	/**
	 *	@brief default constructor
	 *
	 *	@param[in] n_item_num is number of items to be partitioned
	 *	@param[in] n_bin_num is number of bins to partinion items to
	 *
	 *	@note This function can throw the std::bad_alloc exception.
	 */
	inline CUnorderedPartition(size_t n_item_num, size_t n_bin_num) // throw(std::bad_alloc)
		:m_n_item_num(n_item_num), m_n_bin_num(n_bin_num),
		m_counters(std::max(size_t(1), n_bin_num) - 1, 0), m_bin_values(n_bin_num, 0)
	{
		m_bin_values.back() = n_item_num; // the last one is number of items
		//First(); // no need, already there
	}

	/**
	 *	@brief calculates number of partitions
	 *	@return Returns number of unordered partitions of n_Item_Num() items into n_Bin_Num() bins.
	 *	@note This uses recursive approach.
	 */
	inline uint64_t n_Partition_Num() const
	{
		return n_UP_Size(m_n_item_num, m_n_bin_num);
	}

	/**
	 *	@brief gets number of bins
	 *	@return Returns number of bins.
	 */
	inline size_t n_Bin_Num() const
	{
		return m_n_bin_num;
	}

	/**
	 *	@brief gets number of items
	 *	@return Returns number of items.
	 */
	inline size_t n_Item_Num() const
	{
		return m_n_item_num;
	}

	/**
	 *	@brief gets the internal state vector
	 *	@return Returns const reference to the vector, containing internal state.
	 */
	inline const std::vector<size_t> &r_Internal_State() const
	{
		return m_counters;
	}

	/**
	 *	@brief gets numbers of items in bins
	 *	@return Returns const reference to vector of bin values.
	 *	@note The values of bins are in nondecreasing order.
	 */
	inline const std::vector<size_t> &r_Bin_Values() const
	{
		return m_bin_values;
	}

	/**
	 *	@brief gets numbers of items in bins
	 *	@return Returns const pointer to the array of n_Bin_Num() bin values.
	 */
	inline const size_t *p_Bin_Values() const
	{
		return &m_bin_values.front();
	}

	/**
	 *	@brief gets a single bin state from the current partition
	 *	@param[in] n_index is zero-based index of the bin (0 to n_Bin_Num() - 1)
	 *	@return Returns number of items in the selected bin, for the current partition.
	 */
	inline size_t operator [](size_t n_index) const
	{
		return m_bin_values[n_index];
	}

	/**
	 *	@brief sets the state to the first partition
	 */
	void First();

	/**
	 *	@brief increments to the next partition
	 *	@return Returns true if there is a next combination, otherwise returns false.
	 *	@todo Specify what happens when the last permutation is reached. Is it reset? Is it kept? Is it damaged to something else?
	 */
	bool Next();

	/**
	 *	@brief enumerates unordered partitions in parallel
	 *
	 *	@tparam CListener is listener object (needs to implement function operator
	 *		that gets vector of bins, number of threads and the current thread id)
	 *
	 *	@param[in] n_item_num is number of items to be partitioned
	 *	@param[in] n_bin_num is number of bins to partinion items to
	 *	@param[in] parallel_listener is parallel listener object
	 *
	 *	@return Returns number of unordered partitions.
	 */
	template <class CListener>
	static uint64_t n_Fast_Enumerate(size_t n_item_num, size_t n_bin_num, CListener &parallel_listener)
	{
		uint64_t n_combo_num = 0;
		if(n_bin_num <= 2) {
			CUnorderedPartition mypart(n_item_num, n_bin_num);

			do {
				++ n_combo_num;
				parallel_listener(mypart.m_bin_values, 1, 0);
				// report
			} while(mypart._b_Next(1));
			// loop for all combinations in this listener

			return n_combo_num;
		}
		// very small problems can't be parallelized this way

#ifdef _OPENMP
		int n_cpu_num;
		#pragma omp parallel
		{
			if(!omp_get_thread_num())
				n_cpu_num = omp_get_num_threads();
		}
#else // _OPENMP
		int n_cpu_num = 1;
#endif // _OPENMP
		// get number of CPUs

		size_t n_max_first_counter = (n_item_num + n_bin_num - 1) / n_bin_num;
		int _n_max_first_counter = int(n_max_first_counter);

		#pragma	omp parallel
		{
#ifdef _OPENMP
			int n_tid = omp_get_thread_num();
#else // _OPENMP
			int n_tid = 0;
#endif // _OPENMP
			// each thread knows it's id

			CUnorderedPartition mypart(n_item_num, n_bin_num);
			// create a smaller partition problem

			#pragma	omp for schedule(dynamic, 1) reduction(+: n_combo_num)
			for(int i = 0; i <= _n_max_first_counter; ++ i) {
				size_t n_counter_0 = i;
				size_t n_remainder = n_item_num - n_counter_0;

				mypart.m_counters.front() = n_counter_0;
				for(size_t j = 0; j < n_bin_num - 2; ++ j)
					mypart.m_counters[j + 1] = 2 * mypart.m_counters[j] - ((j)? mypart.m_counters[j - 1] : 0);
				// initialize the counters in the lower

				{
					const size_t *p_i = &mypart.m_counters.front(); // the smaller problem
					const size_t n_block_size = n_item_num;
					const size_t n_bins = n_bin_num;
					const size_t n_last = n_bins - 2; // last one in the smaller problem
					// n_last > 0 ~ n_bins - 2 > 0 ~ n_bins > 2
					size_t n_prev = (n_bins > 2/*n_last > 0*/)? p_i[n_bins - 3/*n_last - 1*/] : 0; // can change to compile-time constant
					_ASSERTE(n_last + 1 > n_bins - 2); // condition below is always false, can skip the branch
					size_t n_next = /*(n_last + 1 <= n_bins - 2)? p_i[n_last + 1] :*/ n_block_size;
					// get prev and next counter

					if((p_i[n_last] - n_prev > n_next - p_i[n_last]))
						continue; // the first counter too high
				}
				// make sure that this will generate good sequence

				mypart.m_bin_values.front() = mypart.m_counters.front();
				for(size_t j = 1; j < n_bin_num - 1; ++ j)
					mypart.m_bin_values[j] = mypart.m_counters[j] - mypart.m_counters[j - 1];
				mypart.m_bin_values[n_bin_num - 1] = n_item_num - mypart.m_counters[n_bin_num - 2];
				// fixup the histogram

				do {
					++ n_combo_num;
					parallel_listener(mypart.m_bin_values, n_cpu_num, n_tid);
					// report
				} while(mypart._b_Next(1));
				// loop for all combinations in this listener
			}
		}

		return n_combo_num;
	}

protected:
	inline bool _b_Next(size_t n_min_depth = 0);
	static uint64_t n_UP2_Size(size_t n_item_num, size_t n_bin_num);
	static uint64_t n_UP_Size(size_t n_item_num, size_t n_bin_num);
};

#endif // !__COMBINAT_INCLUDED
