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) {
			const auto x = ((parA % parDiv) + (parB % parDiv)) / parDiv;
			return parA / parDiv + parB / parDiv + x;
		}
	} //namespace implem

	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;
		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;
		}
		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];
	}
} //namespace dk