//
// Copyright (c) 2010 Thomas Klimpel
//
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
//

#include <iostream>
#include <vector>
#include <iterator>
#include <cassert>

#if 1

template <typename It>
void inshuffle(It it, std::size_t n) {
  while (n > 0) {
    It itn = it + n;
    std::size_t kmin = 1;
    while (2*kmin <= n)
      kmin *= 2;
    std::size_t kk = kmin;
    // compute k for the first iteration here, so that we can use "if" instead of "while"
    std::size_t k = kk;
    kk /= 2;
    for (std::size_t i = 0, s = (n+1)/2; i < s; ++i) {
      // compute k for the next iteration, let's call it k1
      std::size_t k1 = (2*(i+1)+1)*kk;
      assert (i == 0 || k1 == k + 2*kk);
      if (k1 > n) {
        k1 /= 2;
        kk /= 2;
      }
      assert(k1 <= n);
      // apply the cyclic permutation
      typename std::iterator_traits<It>::value_type tmp = itn[i];
      itn[i] = it[k-1];
      while (k % 2 == 0) {
        it[k-1] = it[(k/2)-1];
        k /= 2;
      }
      it[k-1] = tmp;
      // update k with k1
      k = k1;
    }
    // the optimized computation of k1 fails for n=2,
    // so skip the 'normalization' loop when possible
    if (n > 3) {
      kk = kmin;
      // compute k for the first iteration here, so that we can use "if" instead of "while"
      k = kk;
      kk /= 2;
      for (std::size_t i = 0, s = (n+2)/4; i < s; ++i) {
        // compute k for the next iteration, let's call it k1
        std::size_t k1 = (2*(i+1)+1)*kk;
        assert (i == 0 || k1 == k + 2*kk);
        if (k1 > n) {
          k1 /= 2;
          kk /= 2;
        }
        assert(k1 <= n);
        if (k > k1) {
          std::size_t outshuffle = 1;
          if (k1 < kmin) {
            kmin = k1;
            // if kmin is updated, do an in-shuffle, otherwise do an out-shuffle.
            outshuffle = 0;
          }
          inshuffle(itn+outshuffle,i+1-outshuffle);
        }
        // update k with k1
        k = k1;
      }
    }
    // implement the tail recursion as an iteration
    it = itn + (n % 2);
    n /= 2;
  }
}

#else

template <typename It>
void inshuffle(It it, std::size_t n) {
  if (n == 0)
    return;
  It itn = it + n;
  for (std::size_t i = 0; 2*i < n; ++i) {
    std::size_t k = 2*i + 1;
    while (2*k <= n)
      k *= 2;
    typename std::iterator_traits<It>::value_type tmp = itn[i];
    itn[i] = it[k-1];
    while (k % 2 == 0) {
      it[k-1] = it[(k/2)-1];
      k /= 2;
    }
    it[k-1] = tmp;
  }
  std::size_t kmin = 1;
  while (2*kmin <= n)
    kmin *= 2;
  for (std::size_t i = 0; (4*i+1) < n; ++i) {
    std::size_t k = 2*i + 1;
    while (2*k <= n)
      k *= 2;
    std::size_t k1 = 2*(i+1) + 1;
    while (2*k1 <= n)
      k1 *= 2;
    if (k > k1) {
      std::size_t outshuffle = 1;
      if (k1 < kmin) {
        kmin = k1;
        outshuffle = 0;
      }
      inshuffle(itn+outshuffle,i+1-outshuffle);
    }
  }
  return inshuffle(itn + (n % 2), n/2);
}

#endif

int main() {
  for (std::size_t n = 1; n <= 14; ++n) {
    std::cout << "Test interlace for n=" << n << std::endl;
    std::vector<std::size_t> data(2*n);
    for (std::size_t i = 0; i < n; ++i) {
      data[i] = 2*i+1;
      data[i+n] = 2*i;
    }
    inshuffle(&data[0],n);
    for (std::size_t i = 0; i < 2*n; ++i)
      std::cout << " " << data[i];
    std::cout << std::endl;
    for (std::size_t i = 1; i < 2*n; ++i)
      if (data[i-1]+1 != data[i]) {
        std::cout << "logic error" << std::endl;
        return 255;
      }
  }
  return 0;
}