text: #pragma once
  #ifndef BITMANIP_HPP
  #define BITMANIP_HPP
  
  #include <cstddef>
  #include <cstdint>
  #include <type_traits>
  
  // ===== USER CONFIGURATION ============================================================================================
  
  // Uncomment the following #define if you prefer std::includes over the custom solutions provided in this section.
  // There are few functions that use <array> and <algorithm> (of which only std::min is used), so you can save on
  // compilation time by using the definitions in this section.
  // #define BITMANIP_PREFER_STD_INCLUDES
  
  // Uncomment the following #define if you want a bitmanip:: namespace.
  // #define BITMANIP_PREFER_NAMESPACE
  
  // ====== CONDITIONALLY CONFIGURED CODE ================================================================================
  #ifdef PREFER_STD_INCLUDES
  #include <algorithm>
  #include <array>
  #include <limits>
  
  #ifdef BITMANIP_PREFER_NAMESPACE
  namespace bitmanip {
  #endif
  
  template <typename T, size_t N>
  using arr = std::array<T, N>;
  
  using std::min;
  
  /**
   * @brief templated variable which contains the number of bits for any given integer type.
   * Example: bits_v<uint32_t> = 32
   */
  template <typename Int, std::enable_if_t<std::is_integral_v<Int>, int> = 0>
  constexpr size_t bits_v = std::numeric_limits<Int>::digits;
  #else
  
  #ifdef BITMANIP_PREFER_NAMESPACE
  namespace bitmanip {
  #endif
  /**
   * @brief Simple implementation of min to avoid including the entirety of <algorithm>
   * @param x the first value
   * @param y the second value
   * @param the minimum of the values
   */
  template <typename T>
  constexpr T min(T x, T y)
  {
      return x < y ? x : y;
  }
  
  template <typename T, size_t N>
  struct arr {
      T data[N];
  };
  
  /**
   * @brief templated variable which contains the number of bits for any given integer type.
   * Example: bits_v<uint32_t> = 32
   */
  template <typename Int, std::enable_if_t<std::is_integral_v<Int>, int> = 0>
  constexpr size_t bits_v = sizeof(Int) * 8;    
  #endif
  // ==== UTILITY ========================================================================================================
  
  /**
   * @brief Integer division but with ceiling (rounding up to the next integer) instead of rounding down as usual.
   * @param the dividend
   * @param the divisor
   * @param the quotient, rounded up to the nearest integer
   */
  template <typename Int, std::enable_if_t<(std::is_integral_v<Int>), int> = 0>
  constexpr Int div_ceil(Int x, Int y)
  {
      return 1 + ((x - 1) / y);
  }
  
  // ===== BIT MANIPULATION ==============================================================================================
  
  /**
   * @brief Returns whether an unsigned integer is a power of 2 or zero.
   * Note that this test is faster than having to test if val is a power of 2.
   * @param val the parameter to test
   * @return true if val is a power of 2 or if val is zero
   */
  template <typename Uint, std::enable_if_t<std::is_unsigned_v<Uint>, int> = 0>
  constexpr bool is_pow2_or_zero(Uint val)
  {
      return (val & (val - 1)) == 0;
  }
  
  /**
   * @brief Returns whether an unsigned integer is a power of 2.
   * @param val the parameter to test
   * @return true if val is a power of 2
   * @see is_pow2_or_zero
   */
  template <typename Uint, std::enable_if_t<std::is_unsigned_v<Uint>, int> = 0>
  constexpr bool is_pow2(Uint val)
  {
      return val != 0 && is_pow2_or_zero(val);
  }
  
  /**
   * @brief Naive implementation of log2 using repeated single-bit rightshifting.
   */
  template <typename Uint, std::enable_if_t<std::is_unsigned_v<Uint>, int> = 0>
  constexpr Uint log2_floor_naive(Uint val) noexcept
  {
      Uint result = 0;
      while (val >>= 1) {
          ++result;
      }
      return result;
  }
  
  /**
   * @brief Fast implementation of log2.
   * See https://graphics.stanford.edu/~seander/bithacks.html#IntegerLog.
   *
   * Unrolled 32-bit version:
   * unsigned shift = (v > 0xFFFF) << 4;
      v >>= shift;
      r |= shift;
  
      shift = (v > 0xFF  ) << 3;
      v >>= shift;
      r |= shift;
  
      shift = (v > 0xF   ) << 2;
      v >>= shift;
      r |= shift;
  
      shift = (v > 0x3   ) << 1;
      v >>= shift;
      r |= shift;
  
      shift = (v > 1) << 0;
      r >>= shift;
      r |= shift;
   */
  template <typename Uint, std::enable_if_t<std::is_unsigned_v<Uint>, int> = 0>
  constexpr Uint log2_floor_fast(Uint v) noexcept
  {
      constexpr size_t iterations = log2_floor_naive(bits_v<Uint>);
      unsigned result = 0;
  
      for (size_t i = iterations; i != 0; --i) {
          unsigned compBits = 1 << (i - 1);
          Uint compShift = Uint{1} << compBits;
          unsigned shift = unsigned{v >= compShift} << (i - 1);
          v >>= shift;
          result |= shift;
      }
  
      return result;
  }
  
  /**
   * @brief log2_floor implementation using De Bruijn multiplication.
   * See https://graphics.stanford.edu/~seander/bithacks.html#IntegerLogDeBruijn.
   * @param val the value
   */
  constexpr uint32_t log2_floor_debruijn(uint32_t val) noexcept
  {
      constexpr uint32_t MultiplyDeBruijnBitPosition[32] = {0, 9,  1,  10, 13, 21, 2,  29, 11, 14, 16, 18, 22, 25, 3, 30,
                                                            8, 12, 20, 28, 15, 17, 24, 7,  19, 27, 23, 6,  26, 5,  4, 31};
  
      // first round up to one less than a power of 2
      // this step is not necessary if val is a power of 2
      // Note: the link in @see says that this is rounding down, but it is actually rounding up
      val |= val >> 1;
      val |= val >> 2;
      val |= val >> 4;
      val |= val >> 8;
      val |= val >> 16;
  
      return MultiplyDeBruijnBitPosition[(val * uint32_t{0x07C4ACDD}) >> 27];
  }
  
  /**
   * @brief Computes the floored binary logarithm of a given integer.
   * Example: log2_floor(100) = 6
   *
   * This templated function will choose the best available method depending on the type of the integer.
   * It may fail if a negative signed integer is given.
   *
   * @param v the value
   * @return the floored binary logarithm
   */
  template <typename Int, std::enable_if_t<std::is_integral_v<Int>, int> = 0>
  constexpr Int log2_floor(Int v) noexcept
  {
      if constexpr (std::is_signed_v<Int>) {
          using Uint = std::make_unsigned_t<Int>;
          DEBUG_ASSERT_GE(v, 0);
          return log2_floor_fast<Uint>(static_cast<Uint>(v));
      }
      else {
          return log2_floor_fast<Int>(v);
      }
  }
  
  template <>
  constexpr uint32_t log2_floor<uint32_t>(uint32_t v) noexcept
  {
      return log2_floor_debruijn(v);
  }
  
  /**
   * @brief Computes the ceiled binary logarithm of a given integer.
   * Example: log2_ceil(100) = 7
   *
   * This templated function will choose the best available method depending on the type of the integer.
   * It may fail if a negative signed integer is given.
   *
   * @param v the value
   * @return the floored binary logarithm
   */
  template <typename Int, std::enable_if_t<std::is_integral_v<Int>, int> = 0>
  constexpr Int log2_ceil(Int val) noexcept
  {
      const Int result = log2_floor(val);
      const bool inputIsntPow2 = val != (static_cast<Int>(1) << result);
      return result + inputIsntPow2;
  }
  
  /**
   * @brief Rounds up an unsigned integer to the next power of 2.
   * Powers of two are not affected.
   * Examples: 100 -> 128, 1 -> 1, 3 -> 4, 3000 -> 4096
   * @param v the value to round up
   */
  template <typename Uint, std::enable_if_t<std::is_unsigned_v<Uint>, int> = 0>
  constexpr Uint ceil_pow2(Uint v)
  {
      DEBUG_ASSERT_NE(v, 0);
      constexpr size_t iterations = log2_floor_naive(bits_v<Uint>);
      // decrement is necessary so that powers of 2 don't get increased
      v--;
      for (size_t i = 0; i < iterations; ++i) {
          // after all iterations, all bits right to the msb will be filled with 1
          v |= v >> (1 << i);
      }
      // this step turns the number back to a power of two
      return v + 1;
  }
  
  /**
   * @brief Rounds down an unsigned integer to the next power of 2.
   * Powers of 2 are not affected.
   * Examples: 100 -> 64, 1 -> 1, 3 -> 2, 3000 -> 2048
   * @param v the value to round down
   */
  template <typename Uint, std::enable_if_t<std::is_unsigned_v<Uint>, int> = 0>
  constexpr Uint floor_pow2(Uint v)
  {
      DEBUG_ASSERT_NE(v, 0);
      constexpr size_t iterations = log2_floor_naive(bits_v<Uint>);
      // leftshift is necessary to floor instead of ceil
      v >>= 1;
      for (size_t i = 0; i < iterations; ++i) {
          // after all iterations, all bits right to the msb will be filled with 1
          v |= v >> (1 << i);
      }
      // this step turns the number back to a power of two
      return v + 1;
  }
  
  /**
   * @brief Interleaves an input number with <bits> zero-bits per input bit. Example: 0b11 -> 0b0101
   * @param input the input number
   * @param bits the amount of zero bits per input bit to interleave
   * @return the input number interleaved with input-bits
   */
  constexpr uint64_t ileave_zeros_naive(uint32_t input, size_t bits)
  {
      const auto lim = min<size_t>(32, div_ceil<size_t>(64, bits + 1));
  
      uint64_t result = 0;
      for (size_t i = 0, b_out = 0; i < lim; ++i) {
          result |= static_cast<uint64_t>(input & 1) << b_out;
          input >>= 1;
          b_out += bits + 1;
      }
  
      return result;
  }
  
  /**
   * @brief Removes each interleaved <bits> bits. Example: 0b010101 --rem 1--> 0b111
   * @param input the input number
   * @param bits input bits per output bit
   * @return the the output with removed bits
   */
  constexpr uint64_t rem_ileaved_bits_naive(uint64_t input, size_t bits) noexcept
  {
      // increment once to avoid modulo divisions by 0
      // this way our function is noexcept and safe for all inputs
      ++bits;
      uint64_t result = 0;
      for (size_t i = 0, b_out = 0; i < 64; ++i) {
          if (i % bits == 0) {
              result |= (input & 1) << b_out++;
          }
          input >>= 1;
      }
  
      return result;
  }
  
  /**
   * @brief Duplicates each input bit <bits> times. Example: 0b101 -> 0b110011
   * @param input the input number
   * @param out_bits_per_in_bits the output bits per input bits
   * @return the output number or 0 if the bits parameter was zero
   */
  constexpr uint64_t dupl_bits_naive(uint64_t input, size_t out_bits_per_in_bits)
  {
      if (out_bits_per_in_bits == 0) {
          return 0;
      }
      const auto lim = div_ceil<size_t>(64, out_bits_per_in_bits);
  
      uint64_t result = 0;
      for (size_t i = 0, b_out = 0; i < lim; ++i) {
          for (size_t j = 0; j < out_bits_per_in_bits; ++j, ++b_out) {
              result |= static_cast<uint64_t>((input >> i) & 1) << b_out;
          }
      }
  
      return result;
  }
  
  /**
   * @brief Interleaves <BITS> zero-bits inbetween each input bit.
   * @param input the input number
   * @tparam BITS the number of bits to be interleaved, must be > 0
   * @return the input interleaved with <BITS> bits
   */
  template <size_t BITS>
  constexpr uint64_t ileave_zeros(uint32_t input)
  {
      if constexpr (BITS == 0) {
          return input;
      }
      else {
          constexpr uint64_t MASKS[] = {
              dupl_bits_naive(ileave_zeros_naive(~uint32_t(0), BITS), 1),
              dupl_bits_naive(ileave_zeros_naive(~uint32_t(0), BITS), 2),
              dupl_bits_naive(ileave_zeros_naive(~uint32_t(0), BITS), 4),
              dupl_bits_naive(ileave_zeros_naive(~uint32_t(0), BITS), 8),
              dupl_bits_naive(ileave_zeros_naive(~uint32_t(0), BITS), 16),
          };
          // log2_floor(0) == 0 so this is always safe, even for 1 bit
          constexpr int start = 4 - (log2_floor(BITS >> 1));
  
          uint64_t n = input;
          for (int i = start; i != -1; --i) {
              size_t shift = BITS * (1 << i);
              n |= n << shift;
              n &= MASKS[i];
          }
  
          return n;
      }
  }
  
  /**
   * @brief Removes each interleaved <BITS> bits. Example: 0b010101 --rem 1--> 0b111
   * If BITS is zero, no bits are removed and the input is returned.
   * @param input the input number
   * @param bits input bits per output bit
   * @return the the output with removed bits
   */
  template <size_t BITS>
  constexpr uint64_t rem_ileaved_bits(uint64_t input) noexcept
  {
      if constexpr (BITS == 0) {
          return input;
      }
      else {
          constexpr uint64_t MASKS[] = {
              dupl_bits_naive(ileave_zeros_naive(~uint32_t(0), BITS), 1),
              dupl_bits_naive(ileave_zeros_naive(~uint32_t(0), BITS), 2),
              dupl_bits_naive(ileave_zeros_naive(~uint32_t(0), BITS), 4),
              dupl_bits_naive(ileave_zeros_naive(~uint32_t(0), BITS), 8),
              dupl_bits_naive(ileave_zeros_naive(~uint32_t(0), BITS), 16),
              dupl_bits_naive(ileave_zeros_naive(~uint32_t(0), BITS), 32),
          };
          // log2_floor(0) == 0 so this is always safe, even for 1 bit
          constexpr size_t iterations = 5 - (log2_floor(BITS >> 1));
  
          input &= MASKS[0];
  
          for (size_t i = 0; i < iterations; ++i) {
              size_t rshift = (1 << i) * BITS;
              input |= input >> rshift;
              input &= MASKS[i + 1];
          }
  
          return input;
      }
  }
  
  /**
   * @brief Interleaves 2 integers, where hi comprises the upper bits of each bit pair and lo the lower bits.
   * Examle: ileave2(0b111, 0b000) = 0b101010
   *
   * This is also referred to as a Morton Code in scientific literature.
   *
   * @param hi the high bits
   * @param lo the low bits
   * @return the interleaved bits
   */
  constexpr uint64_t ileave2(uint32_t hi, uint32_t lo)
  {  // TODO implement like ileav3
      constexpr uint64_t MASKS[] = {0x5555'5555'5555'5555,
                                    0x3333'3333'3333'3333,
                                    0x0F0F'0F0F'0F0F'0F0F,
                                    0x00FF'00FF'00FF'00FF,
                                    0x0000'FFFF'0000'FFFF};
  
      uint64_t result = 0;
      uint32_t *nums[] = {&hi, &lo};
  
      for (size_t i = 0; i < 2; ++i) {
          uint64_t n = *nums[i];
          for (int i = 4; i != -1; --i) {
              n |= n << (1 << i);
              n &= MASKS[i];
          }
          result |= n << (1 - i);
      }
  
      return result;
  }
  
  /**
   * @brief Interleaves 3 integers, where x comprises the uppermost bits of each bit triple and z the lowermost bits.
   *
   * This is also referred to as a Morton Code in scientific literature.
   *
   * @param x the highest bits
   * @param y the middle bits
   * @param z the lowest bits
   * @return the interleaved bits
   */
  constexpr uint64_t ileave3(uint32_t x, uint32_t y, uint32_t z)
  {
      return (ileave_zeros<2>(x) << 2) | (ileave_zeros<2>(y) << 1) | ileave_zeros<2>(z);
  }
  
  /**
   * @brief Deinterleaves 3 integers which are interleaved in a single number.
   * Visualization: abcdefghi -> (adg, beh, cfi)
   *
   * This is also referred to as a Morton Code in scientific literature.
   *
   * @param n the number
   * @return the interleaved bits
   */
  constexpr arr<uint32_t, 3> dileave3(uint64_t n)
  {
      uint32_t x = static_cast<uint32_t>(rem_ileaved_bits<2>(n >> 2));
      uint32_t y = static_cast<uint32_t>(rem_ileaved_bits<2>(n >> 1));
      uint32_t z = static_cast<uint32_t>(rem_ileaved_bits<2>(n >> 0));
      return {x, y, z};
  }
  
  #ifdef BITMANIP_PREFER_NAMESPACE
  }
  #endif
  
  #endif  // BITMANIP_HPP
  