namespace dk {
	namespace implem {
		///----------------------------------------------------------------------
		///Equivalent to (parA + parB) % parDiv
		///----------------------------------------------------------------------
		inline CoordinateScalarType sum_mod (CoordinateScalarType parA, CoordinateScalarType parB, CoordinateScalarType parDiv) {
			//const auto sum = m_pos[index] + parValue;
			//const auto sz = m_size[index] + 1;
			//const auto remainder = sum - (sum / sz) * sz;

			//m_pos[index] = remainder;
			//parValue /= sz;

			const auto rounded_sum = sum_div(parA, parB, parDiv) * parDiv;
			return std::max(parA, parB) - rounded_sum + std::min(parA, parB);
		}

		///----------------------------------------------------------------------
		///Equivalent to (parA + parB) / parDiv
		///----------------------------------------------------------------------
		inline CoordinateScalarType sum_div (CoordinateScalarType parA, CoordinateScalarType parB, CoordinateScalarType parDiv) {
			DK_ASSERT(parA >= 0 and parB >= 0 and parDiv > 0);
			DK_ASSERT(std::numeric_limits<CoordinateScalarType>::max() - (parA % parDiv) >= (parB % parDiv));
			const auto x = ((parA % parDiv) + (parB % parDiv)) / parDiv;
			DK_ASSERT(std::numeric_limits<CoordinateScalarType>::max() - parA / parDiv >= parB / parDiv);
			DK_ASSERT(std::numeric_limits<CoordinateScalarType>::max() - parA / parDiv - parB / parDiv >= x);
			return parA / parDiv + parB / parDiv + x;
		}

		inline CoordinateScalarType sub_mod (CoordinateScalarType parA, CoordinateScalarType parB, CoordinateScalarType parDiv) {
			typedef std::make_unsigned<CoordinateScalarType>::type unsigned_coord;
			const unsigned_coord rounded_sum = static_cast<unsigned_coord>(sub_div(parA, parB, parDiv) * parDiv);
			const unsigned_coord ret = static_cast<unsigned_coord>(parA) - rounded_sum - static_cast<unsigned_coord>(parB);
			return static_cast<CoordinateScalarType>(ret);
		}

		inline CoordinateScalarType sub_div (CoordinateScalarType parA, CoordinateScalarType parB, CoordinateScalarType parDiv) {
			const auto x = ((parA % parDiv) + (parB % parDiv)) / parDiv;
			return parA / parDiv - parB / parDiv + x;
		}
	} //namespace implem

	template <uint32_t D>
	inline CoordinateDistType to_index (const Vector<D>& parPos, const Vector<D>& parUpper) {
		CoordinateDistType index = 0;
		for (CoordinateDistType d = 0; d < D; ++d) {
			auto pos = static_cast<CoordinateDistType>(parPos[D - 1 - d]);
			for (CoordinateDistType p = 0; p < D - 1 - d; ++p) {
				pos *= static_cast<CoordinateDistType>(parUpper[p] + 1);
			}

			index += pos;
		}
		return index;
	}

#if defined(NDEBUG)
	template <>
	inline CoordinateDistType to_index<2> (const Vector<2>& parPos, const Vector<2>& parUpper) {
		return parPos.y() * (parUpper.x() + 1) + parPos.x();
	}
#endif

	template <uint32_t D>
	inline CoordinateDistType to_index (const TileCoords<D>& parTC) {
		return to_index<D>(parTC.position(), parTC.upper());
	}

	template <uint32_t D>
	TileCoords<D>::TileCoords (const coords& parSize) :
		m_pos(CoordinateScalarType(0)),
		m_size(parSize)
	{
		DK_ASSERT(m_pos <= m_size);
		DK_ASSERT(m_size < coords(std::numeric_limits<CoordinateScalarType>::max()));
	}

	template <uint32_t D>
	TileCoords<D>::TileCoords (const coords& parValue, const coords& parSize) :
		m_pos(parValue),
		m_size(parSize)
	{
		DK_ASSERT(m_pos <= m_size);
		DK_ASSERT(m_size < coords(std::numeric_limits<CoordinateScalarType>::max()));
	}

	template <uint32_t D>
	TileCoords<D>& TileCoords<D>::operator++ () {
		const coords lower(0);
		for (std::size_t d = 0; d < static_cast<std::size_t>(D - 1); ++d) {
			++m_pos[d];
			if (m_pos[d] > m_size[d])
				m_pos[d] = lower[d];
			else
				return *this;
		}
		++m_pos[static_cast<std::size_t>(D - 1)];
		return *this;
	}

	template <uint32_t D>
	TileCoords<D> TileCoords<D>::operator++ (int) {
		TileCoords<D> retval(*this);
		++(*this);
		return retval;
	}

	template <uint32_t D>
	TileCoords<D>& TileCoords<D>::operator-- () {
		const coords lower(0);
		for (std::size_t d = 0; d < static_cast<std::size_t>(D); ++d) {
			if (m_pos[d] > lower[d]) {
				--m_pos[d];
				return;
			}
			else {
				m_pos[d] = m_size[d];
			}
		}
		++m_pos[static_cast<std::size_t>(D - 1)];
		return *this;
	}

	template <uint32_t D>
	TileCoords<D> TileCoords<D>::operator-- (int) {
		TileCoords<D> retval(*this);
		--(*this);
		return retval;
	}

	template <uint32_t D>
	const TileCoords<D>& TileCoords<D>::operator-= (CoordinateScalarType parValue) {
		std::size_t index = 0;
		if (parValue > 0) {
			while (parValue) {
				const auto t = parValue % (m_size[index] + 1);
				DK_ASSERT(t >= 0);
				const CoordinateScalarType r = (t > m_pos[index] ? 1 : 0);
				m_pos[index] = (m_size[index] + 1) * r + m_pos[index] - t;
				parValue /= (m_size[index] + 1);
				parValue += r;

				++index;
			}
		}
		else if (parValue < 0) {
			while (parValue) {
				const auto new_pos = implem::sub_mod(m_pos[index], parValue, m_size[index] + 1);
				parValue = -implem::sub_div(m_pos[index], parValue, m_size[index] + 1);
				m_pos[index] = new_pos;

				++index;
			}
		}

		return *this;
	}

	template <uint32_t D>
	const TileCoords<D>& TileCoords<D>::operator+= (CoordinateScalarType parValue) {
		std::size_t index = 0;
		if (parValue > 0) {
			while (parValue) {
				//naïve implementation
				//parValue += m_pos[index];
				//m_pos[index] = parValue % (m_size[index] + 1);
				//parValue /= (m_size[index] + 1);

				//overflow-aware implementation
				const auto new_pos = implem::sum_mod(m_pos[index], parValue, m_size[index] + 1);
				parValue = implem::sum_div(m_pos[index], parValue, m_size[index] + 1);
				m_pos[index] = new_pos;

				++index;
			}
		}
		else if (parValue < 0) {
			while (parValue) {
				//const auto t = (-parValue) % (m_size[index] + 1);
				const auto t = -(parValue % (m_size[index] + 1));
				DK_ASSERT(t >= 0);
				const CoordinateScalarType r = (t > m_pos[index] ? 1 : 0);
				m_pos[index] = (m_size[index] + 1) * r + m_pos[index] - t;
				parValue /= (m_size[index] + 1);
				parValue -= r;

				++index;
			}
		}

		return *this;
	}

	template <uint32_t D>
	bool TileCoords<D>::operator== (const TileCoords& parOther) const {
		return m_pos == parOther.m_pos and m_size == parOther.m_size;
	}

	template <uint32_t D>
	bool TileCoords<D>::operator!= (const TileCoords& parOther) const {
		return not(*this == parOther);
	}

	template <uint32_t D>
	const CoordinateScalarType& TileCoords<D>::operator[] (int parIndex) const {
		DK_ASSERT(parIndex >= 0);
		DK_ASSERT(static_cast<uint32_t>(parIndex) < D);
		return m_pos[parIndex];
	}
	template <uint32_t D>
	CoordinateScalarType& TileCoords<D>::operator[] (int parIndex) {
		DK_ASSERT(parIndex >= 0);
		DK_ASSERT(static_cast<uint32_t>(parIndex) < D);
		return m_pos[parIndex];
	}

	template <uint32_t D>
	auto TileCoords<D>::position() const -> const coords& {
		return m_pos;
	}

	template <uint32_t D>
	auto TileCoords<D>::upper() const -> const coords& {
		return m_size;
	}
} //namespace dk