/*
  Copyright (C) 2018 David Robillard <d@drobilla.net>
  Copyright (C) 2006-2007 Chris Hamilton <chamilton@cs.dal.ca>

  This program is free software: you can redistribute it and/or modify it under
  the terms of the GNU General Public License as published by the Free Software
  Foundation, either version 2 of the License, or (at your option) any later
  version.

  This program is distributed in the hope that it will be useful, but WITHOUT
  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
  details.

  You should have received a copy of the GNU General Public License along with
  this program.  If not, see <https://www.gnu.org/licenses/>.
*/

#ifndef CHILBERT_GRAYCODERANK_HPP
#define CHILBERT_GRAYCODERANK_HPP

#include <cassert>
#include <cstddef>
#include <numeric>

namespace chilbert {

// This is the bulk of the cost in calculating
// a Compact Hilbert Index.  It compresses a previously
// calculated index when provided the rotation
// at each level of precision.
template <class H, class HC>
inline void
compact_index(const size_t* const ms,
              const size_t* const ds,
              const size_t        n,
              const size_t        m,
              H&                  h,
              HC&                 hc)
{
	assert(h.size() >= n * m);
	assert(hc.size() >= std::accumulate(ms, ms + n, size_t(0)));

	reset_bits(hc);

	auto hm  = h.mask(0);
	auto hcm = hc.mask(0);

	// Run through the levels of precision
	for (size_t i = 0; i < m; i++) {
		// Run through the dimensions
		size_t j = ds[i];
		do {
			// This dimension contributes a bit?
			if (ms[j] > i) {
				if (h.test(hm)) {
					hc.set(hcm);
				}
				++hcm;
			}

			if (++j == n) {
				j = 0;
			}
			++hm;
		} while (j != ds[i]);
	}
}

template <class I>
inline void
gray_code_rank(const I& mask, const I& gi, const size_t n, I& r)
{
	assert(mask.size() == n);
	assert(gi.size() == n);
	assert(r.size() == n);

	r.reset();

	auto mi  = mask.begin();
	auto gii = gi.begin();
	auto ri  = r.begin();

	for (size_t i = 0; i < n; ++i, ++mi, ++gii) {
		if (*mi) {
			if (*gii) {
				ri.set();
			}
			++ri;
		}
	}
}

template <class I>
inline void
gray_code_rank_inv(const I&     mask,
                   const I&     ptrn,
                   const I&     r,
                   const size_t n,
                   const size_t b,
                   I&           g,
                   I&           gi)
{
	assert(mask.size() == n);
	assert(ptrn.size() == n);
	assert(r.size() == n);
	assert(g.size() == n);
	assert(gi.size() == n);

	g.reset();
	gi.reset();

	auto m  = mask.mask(n - 1);
	auto ri = r.begin(b - 1);

	typename I::Rack gi0 = 0;
	typename I::Rack gi1 = 0;
	typename I::Rack g0  = 0;

	for (size_t i = 0; i < n; ++i) {
		if (mask.test(m)) { // Unconstrained bit
			gi1 = gi0;
			gi0 = *ri;
			g0  = gi0 ^ gi1;
			if (gi0) {
				gi.set(m);
			}
			if (g0) {
				g.set(m);
			}
			--ri;
		} else { // Constrained bit
			g0  = (ptrn.test(m) > 0);
			gi1 = gi0;
			gi0 = g0 ^ gi1;
			if (gi0) {
				gi.set(m);
			}
			if (g0) {
				g.set(m);
			}
		}

		--m;
	}
}

template <class I>
inline void
extract_mask(const size_t* const ms,
             const size_t        n,
             const size_t        d,
             const size_t        i,
             I&                  mask,
             size_t&             b)
{
	assert(0 <= d && d < n);
	assert(mask.size() == n);

	mask.reset();
	b = 0;

	auto   mi = mask.begin();
	size_t j  = d; // #D j = (d==n-1) ? 0 : d+1;
	do {
		if (ms[j] > i) {
			mi.set();
			++b;
		}

		++mi;
		if (++j == n) {
			j = 0;
		}
	} while (j != d);
}

} // namespace chilbert

#endif