//
// Copyright (c) 2003 Kresimir Fresl
// 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)
//

#ifndef BOOST_NUMERIC_BINDINGS_DETAIL_COMPLEX_UTILS_HPP
#define BOOST_NUMERIC_BINDINGS_DETAIL_COMPLEX_UTILS_HPP

#include <iterator>
#include <boost/numeric/bindings/begin.hpp>
#include <boost/numeric/bindings/end.hpp>
#include <boost/numeric/bindings/is_complex.hpp>
#include <boost/numeric/bindings/remove_imaginary.hpp>
#include <boost/numeric/bindings/value_type.hpp>
#include <boost/numeric/bindings/vector_view.hpp>
#include <boost/utility/enable_if.hpp>

namespace boost {
namespace numeric {
namespace bindings {

namespace detail {

#ifdef BOOST_NUMERIC_BINDINGS_BY_THE_BOOK
template <typename It>
void inshuffle(It it, std::size_t n) {
  if (n==0) return;
  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 = it[n+i];
    it[n+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) {
      if (k1 < kmin) {
        kmin = k1;
        inshuffle(it+n, i+1);
      }
      else
        inshuffle(it+n+1, i);
    }
  }
  return inshuffle(it+n+(n%2), n/2);
}
#else
template <typename It>
void inshuffle(It it, std::size_t n) {
  while (n > 0) {
    std::size_t kmin = 1;
    while (kmin <= n)
      kmin *= 2;
    {
      std::size_t kk = kmin/2;
      It itn = it + n;
      for (std::size_t i = 0, s = (n+1)/2; i < s; ++i) {
        std::size_t k = (2*i+1)*kk;
        while (k > n) {
          k /= 2;
          kk /= 2;
        }
        // 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;
      }
    }
    // the optimized computation of k fails for n=2,
    // so skip the 'normalization' loop when possible
    if (n > 3) {
      std::size_t kk = kmin/4;
      for (std::size_t i = 1; 4*i < n+3; ++i) {
        std::size_t k = (2*i+1)*kk;
        if (k > n) {
          kk /= 2;
          if (k < kmin) {
            kmin = k;
            // if kmin is updated, do an in-shuffle
            inshuffle(it+n, i);
          }
          else
            // otherwise do an out-shuffle
            inshuffle(it+n+1, i-1);
        }
      }
    }
    // implement the tail recursion as an iteration
    it += n+(n%2);
    n /= 2;
  }
}
#endif

// Reorders a real array followed by an imaginary array to a true complex array
// where real and imaginary part of each number directly follow each other.
template <typename VectorW>
typename boost::enable_if< is_complex< typename bindings::value_type< VectorW >::type >, void >::type
interlace (VectorW& w) {
  typedef typename bindings::value_type< VectorW >::type value_type;
  typedef typename bindings::remove_imaginary< value_type >::type real_type;
  value_type* pw = bindings::begin_value(w);
  std::ptrdiff_t n = bindings::end_value(w) - pw;
  if (n < 2) return;
  inshuffle(reinterpret_cast<real_type*> (pw)+1, n-1);
}

} // namespace detail

namespace result_of {

template< typename VectorW >
struct real_part_view {
    typedef typename bindings::result_of::vector_view< typename
      bindings::remove_imaginary< typename
      bindings::value_type< VectorW >::type
      >::type >::type type;
};

template< typename VectorW >
struct imag_part_view {
    typedef typename bindings::result_of::vector_view< typename
      bindings::remove_imaginary< typename
      bindings::value_type< VectorW >::type
      >::type >::type type;
};

} // namespace result_of

namespace detail {

// Creates a real vector_view to the first half of the complex array,
// which is intended to be filled by the real part
template <typename VectorW>
typename boost::enable_if< is_complex< typename bindings::value_type< VectorW >::type >,
        typename result_of::real_part_view< VectorW >::type const >::type
real_part_view (VectorW& w) {
  typedef typename bindings::value_type< VectorW >::type value_type;
  typedef typename bindings::remove_imaginary< value_type >::type real_type;
  value_type* pw = bindings::begin_value(w);
  std::ptrdiff_t n = bindings::end_value(w) - pw;
  return bindings::vector_view(reinterpret_cast<real_type*> (pw), n);
}

// Creates a real vector_view to the second half of the complex array,
// which is intended to be filled by the imaginary part
template <typename VectorW>
typename boost::enable_if< is_complex< typename bindings::value_type< VectorW >::type >,
        typename result_of::imag_part_view< VectorW >::type const >::type
imag_part_view (VectorW& w) {
  typedef typename bindings::value_type< VectorW >::type value_type;
  typedef typename bindings::remove_imaginary< value_type >::type real_type;
  value_type* pw = bindings::begin_value(w);
  std::ptrdiff_t n = bindings::end_value(w) - pw;
  return bindings::vector_view(reinterpret_cast<real_type*> (pw)+n, n);
}

} // namespace detail

} // namespace bindings
} // namespace numeric
} // namespace boost

#endif