text: /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
  /*                                                                       */
  /*    This file is part of the HiGHS linear optimization suite           */
  /*                                                                       */
  /*    Written and engineered 2008-2022 at the University of Edinburgh    */
  /*                                                                       */
  /*    Available as open-source under the MIT License                     */
  /*                                                                       */
  /*    Authors: Julian Hall, Ivet Galabova, Leona Gottwald and Michael    */
  /*    Feldmeier                                                          */
  /*                                                                       */
  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
  #ifndef HIGHS_UTIL_HASH_H_
  #define HIGHS_UTIL_HASH_H_
  
  #include <array>
  #include <cassert>
  #include <cmath>
  #include <cstddef>
  #include <cstdint>
  #include <cstring>
  #include <functional>
  #include <iterator>
  #include <memory>
  #include <type_traits>
  #include <utility>
  #include <vector>
  
  #include "util/HighsInt.h"
  
  #ifdef HIGHS_HAVE_BITSCAN_REVERSE
  #include <intrin.h>
  #pragma intrinsic(_BitScanReverse)
  #pragma intrinsic(_BitScanReverse64)
  #endif
  
  #if __GNUG__ && __GNUC__ < 5
  #define IS_TRIVIALLY_COPYABLE(T) __has_trivial_copy(T)
  #else
  #define IS_TRIVIALLY_COPYABLE(T) std::is_trivially_copyable<T>::value
  #endif
  
  template <typename T>
  struct HighsHashable : std::integral_constant<bool, IS_TRIVIALLY_COPYABLE(T)> {
  };
  
  template <typename U, typename V>
  struct HighsHashable<std::pair<U, V>>
      : public std::integral_constant<bool, HighsHashable<U>::value &&
                                                HighsHashable<V>::value> {};
  
  template <typename U, typename V>
  struct HighsHashable<std::tuple<U, V>> : public HighsHashable<std::pair<U, V>> {
  };
  template <typename U, typename V, typename W, typename... Args>
  struct HighsHashable<std::tuple<U, V, W, Args...>>
      : public std::integral_constant<
            bool, HighsHashable<U>::value &&
                      HighsHashable<std::tuple<V, W, Args...>>::value> {};
  
  struct HighsHashHelpers {
    using u8 = std::uint8_t;
    using i8 = std::int8_t;
  
    using u16 = std::uint16_t;
    using i16 = std::int16_t;
  
    using u32 = std::uint32_t;
    using i32 = std::int32_t;
  
    using u64 = std::uint64_t;
    using i64 = std::uint64_t;
  
    static constexpr u64 c[] = {
        u64{0xc8497d2a400d9551}, u64{0x80c8963be3e4c2f3}, u64{0x042d8680e260ae5b},
        u64{0x8a183895eeac1536}, u64{0xa94e9c75f80ad6de}, u64{0x7e92251dec62835e},
        u64{0x07294165cb671455}, u64{0x89b0f6212b0a4292}, u64{0x31900011b96bf554},
        u64{0xa44540f8eee2094f}, u64{0xce7ffd372e4c64fc}, u64{0x51c9d471bfe6a10f},
        u64{0x758c2a674483826f}, u64{0xf91a20abe63f8b02}, u64{0xc2a069024a1fcc6f},
        u64{0xd5bb18b70c5dbd59}, u64{0xd510adac6d1ae289}, u64{0x571d069b23050a79},
        u64{0x60873b8872933e06}, u64{0x780481cc19670350}, u64{0x7a48551760216885},
        u64{0xb5d68b918231e6ca}, u64{0xa7e5571699aa5274}, u64{0x7b6d309b2cfdcf01},
        u64{0x04e77c3d474daeff}, u64{0x4dbf099fd7247031}, u64{0x5d70dca901130beb},
        u64{0x9f8b5f0df4182499}, u64{0x293a74c9686092da}, u64{0xd09bdab6840f52b3},
        u64{0xc05d47f3ab302263}, u64{0x6b79e62b884b65d6}, u64{0xa581106fc980c34d},
        u64{0xf081b7145ea2293e}, u64{0xfb27243dd7c3f5ad}, u64{0x5211bf8860ea667f},
        u64{0x9455e65cb2385e7f}, u64{0x0dfaf6731b449b33}, u64{0x4ec98b3c6f5e68c7},
        u64{0x007bfd4a42ae936b}, u64{0x65c93061f8674518}, u64{0x640816f17127c5d1},
        u64{0x6dd4bab17b7c3a74}, u64{0x34d9268c256fa1ba}, u64{0x0b4d0c6b5b50d7f4},
        u64{0x30aa965bc9fadaff}, u64{0xc0ac1d0c2771404d}, u64{0xc5e64509abb76ef2},
        u64{0xd606b11990624a36}, u64{0x0d3f05d242ce2fb7}, u64{0x469a803cb276fe32},
        u64{0xa4a44d177a3e23f4}, u64{0xb9d9a120dcc1ca03}, u64{0x2e15af8165234a2e},
        u64{0x10609ba2720573d4}, u64{0xaa4191b60368d1d5}, u64{0x333dd2300bc57762},
        u64{0xdf6ec48f79fb402f}, u64{0x5ed20fcef1b734fa}, u64{0x4c94924ec8be21ee},
        u64{0x5abe6ad9d131e631}, u64{0xbe10136a522e602d}, u64{0x53671115c340e779},
        u64{0x9f392fe43e2144da}};
  
    /// mersenne prime 2^61 - 1
    static constexpr u64 M61() { return u64{0x1fffffffffffffff}; };
  
  #ifdef HIGHS_HAVE_BUILTIN_CLZ
    static int log2i(uint64_t n) { return 63 - __builtin_clzll(n); }
  
    static int log2i(unsigned int n) { return 31 - __builtin_clz(n); }
  
  #elif defined(HIGHS_HAVE_BITSCAN_REVERSE)
    static int log2i(uint64_t n) {
      unsigned long result;
      _BitScanReverse64(&result, n);
      return result;
    }
  
    static int log2i(unsigned int n) {
      unsigned long result;
      _BitScanReverse64(&result, (unsigned long)n);
      return result;
    }
  #else
    // integer log2 algorithm without floating point arithmetic. It uses an
    // unrolled loop and requires few instructions that can be well optimized.
    static int log2i(uint64_t n) {
      int x = 0;
  
      auto log2Iteration = [&](int p) {
        if (n >= uint64_t{1} << p) {
          x += p;
          n >>= p;
        }
      };
  
      log2Iteration(32);
      log2Iteration(16);
      log2Iteration(8);
      log2Iteration(4);
      log2Iteration(2);
      log2Iteration(1);
  
      return x;
    }
  
    static int log2i(uint32_t n) {
      int x = 0;
  
      auto log2Iteration = [&](int p) {
        if (n >= 1u << p) {
          x += p;
          n >>= p;
        }
      };
  
      log2Iteration(16);
      log2Iteration(8);
      log2Iteration(4);
      log2Iteration(2);
      log2Iteration(1);
  
      return x;
    }
  
  #endif
  
    /// compute a * b mod 2^61-1
    static u64 multiply_modM61(u64 a, u64 b) {
      u64 ahi = a >> 32;
      u64 bhi = b >> 32;
      u64 alo = a & 0xffffffffu;
      u64 blo = b & 0xffffffffu;
  
      // compute the different order terms with adicities 2^64, 2^32, 2^0
      u64 term_64 = ahi * bhi;
      u64 term_32 = ahi * blo + bhi * alo;
      u64 term_0 = alo * blo;
  
      // Partially reduce term_0 and term_32 modulo M61() individually to not deal
      // with a possible carry bit (thanks @https://github.com/WTFHCN for catching
      // the bug with this). We do not need to completely reduce by an additional
      // check for the range of the resulting term as this is done in the end in
      // any case and the reduced sizes do not cause troubles with the available
      // 64 bits.
      term_0 = (term_0 & M61()) + (term_0 >> 61);
      term_0 += ((term_32 >> 29) + (term_32 << 32)) & M61();
  
      // The lower 61 bits of term_0 are now the lower 61 bits of the result that
      // we need. Now extract the upper 61 of the result so that we can compute
      // the result of the multiplication modulo M61()
      u64 ab61 = (term_64 << 3) | (term_0 >> 61);
  
      // finally take the result modulo M61 which is computed by exploiting
      // that M61 is a mersenne prime, particularly, if a * b = q * 2^61 + r
      // then a * b = (q + r) (mod 2^61 - 1)
      u64 result = (term_0 & M61()) + ab61;
      if (result >= M61()) result -= M61();
      return result;
    }
  
    static u64 modexp_M61(u64 a, u64 e) {
      // the exponent need to be greater than zero
      assert(e > 0);
      u64 result = a;
  
      while (e != 1) {
        // square
        result = multiply_modM61(result, result);
  
        // multiply with a if exponent is odd
        if (e & 1) result = multiply_modM61(result, a);
  
        // shift to next bit
        e = e >> 1;
      }
  
      return result;
    }
  
    /// mersenne prime 2^31 - 1
    static constexpr u64 M31() { return u32{0x7fffffff}; };
  
    /// compute a * b mod 2^31-1
    static u32 multiply_modM31(u32 a, u32 b) {
      u64 result = u64(a) * u64(b);
      result = (result >> 31) + (result & M31());
      if (result >= M31()) result -= M31();
      return result;
    }
  
    static u32 modexp_M31(u32 a, u64 e) {
      // the exponent need to be greater than zero
      assert(e > 0);
      u32 result = a;
  
      while (e != 1) {
        // square
        result = multiply_modM31(result, result);
  
        // multiply with a if exponent is odd
        if (e & 1) result = multiply_modM31(result, a);
  
        // shift to next bit
        e = e >> 1;
      }
  
      return result;
    }
  
    template <HighsInt k>
    static u64 pair_hash(u32 a, u32 b) {
      return (a + c[2 * k]) * (b + c[2 * k + 1]);
    }
  
    static void sparse_combine(u64& hash, HighsInt index, u64 value) {
      // we take each value of the sparse hash as coefficient for a polynomial
      // of the finite field modulo the mersenne prime 2^61-1 where the monomial
      // for a sparse entry has the degree of its index. We evaluate the
      // polynomial at a random constant. This allows to compute the hashes of
      // sparse vectors independently of each others nonzero contribution and
      // therefore allows to use the order of best access patterns for cache
      // performance. E.g. we can compute a strong hash value for parallel row and
      // column detection and only need to loop over the nonzeros once in
      // arbitrary order. This comes at the expense of more expensive hash
      // calculations as it would be more efficient to evaluate the polynomial
      // with horners scheme, but allows for parallelization and arbitrary order.
      // Since we have 64 random constants available, we slightly improve
      // the scheme by using a lower degree polynomial with 64 variables
      // which we evaluate at the random vector of 64.
  
      // make sure input value is never zero and at most 61bits are used
      value = ((value << 1) & M61()) | 1;
  
      // make sure that the constant has at most 61 bits, as otherwise the modulo
      // algorithm for multiplication mod M61 might not work properly due to
      // overflow
      u64 a = c[index & 63] & M61();
      HighsInt degree = (index >> 6) + 1;
  
      hash += multiply_modM61(value, modexp_M61(a, degree));
      hash = (hash >> 61) + (hash & M61());
      if (hash >= M61()) hash -= M61();
      assert(hash < M61());
    }
  
    static void sparse_inverse_combine(u64& hash, HighsInt index, u64 value) {
      // same hash algorithm as sparse_combine(), but for updating a hash value to
      // the state before it was changed with a call to sparse_combine(). This is
      // easily possible as the hash value just uses finite field arithmetic. We
      // can simply add the additive inverse of the previous hash value. This is a
      // very useful routine for symmetry detection. During partition refinement
      // the hashes do not need to be recomputed but can be updated with this
      // procedure.
  
      // make sure input value is never zero and at most 61bits are used
      value = ((value << 1) & M61()) | 1;
  
      u64 a = c[index & 63] & M61();
      HighsInt degree = (index >> 6) + 1;
      // add the additive inverse (M61() - hashvalue) instead of the hash value
      // itself
      hash += M61() - multiply_modM61(value, modexp_M61(a, degree));
      hash = (hash >> 61) + (hash & M61());
      if (hash >= M61()) hash -= M61();
      assert(hash < M61());
    }
  
    /// overload that is not taking a value and saves one multiplication call
    /// useful for sparse hashing of bit vectors
    static void sparse_combine(u64& hash, HighsInt index) {
      u64 a = c[index & 63] & M61();
      HighsInt degree = (index >> 6) + 1;
  
      hash += modexp_M61(a, degree);
      hash = (hash >> 61) + (hash & M61());
      if (hash >= M61()) hash -= M61();
      assert(hash < M61());
    }
  
    /// overload that is not taking a value and saves one multiplication call
    /// useful for sparse hashing of bit vectors
    static void sparse_inverse_combine(u64& hash, HighsInt index) {
      // same hash algorithm as sparse_combine(), but for updating a hash value to
      // the state before it was changed with a call to sparse_combine(). This is
      // easily possible as the hash value just uses finite field arithmetic. We
      // can simply add the additive inverse of the previous hash value. This is a
      // very useful routine for symmetry detection. During partition refinement
      // the hashes do not need to be recomputed but can be updated with this
      // procedure.
  
      u64 a = c[index & 63] & M61();
      HighsInt degree = (index >> 6) + 1;
      // add the additive inverse (M61() - hashvalue) instead of the hash value
      // itself
      hash += M61() - modexp_M61(a, degree);
      hash = (hash >> 61) + (hash & M61());
      if (hash >= M61()) hash -= M61();
      assert(hash < M61());
    }
  
    static void sparse_combine32(u32& hash, HighsInt index, u64 value) {
      // we take each value of the sparse hash as coefficient for a polynomial
      // of the finite field modulo the mersenne prime 2^61-1 where the monomial
      // for a sparse entry has the degree of its index. We evaluate the
      // polynomial at a random constant. This allows to compute the hashes of
      // sparse vectors independently of each others nonzero contribution and
      // therefore allows to use the order of best access patterns for cache
      // performance. E.g. we can compute a strong hash value for parallel row and
      // column detection and only need to loop over the nonzeros once in
      // arbitrary order. This comes at the expense of more expensive hash
      // calculations as it would be more efficient to evaluate the polynomial
      // with horners scheme, but allows for parallelization and arbitrary order.
      // Since we have 16 random constants available, we slightly improve
      // the scheme by using a lower degree polynomial with 16 variables
      // which we evaluate at the random vector of 16.
  
      // make sure input value is never zero and at most 31bits are used
      value = (pair_hash<0>(value, value >> 32) >> 33) | 1;
  
      // make sure that the constant has at most 31 bits, as otherwise the modulo
      // algorithm for multiplication mod M31 might not work properly due to
      // overflow
      u32 a = c[index & 63] & M31();
      HighsInt degree = (index >> 6) + 1;
  
      hash += multiply_modM31(value, modexp_M31(a, degree));
      hash = (hash >> 31) + (hash & M31());
      if (hash >= M31()) hash -= M31();
      assert(hash < M31());
    }
  
    static void sparse_inverse_combine32(u32& hash, HighsInt index, u64 value) {
      // same hash algorithm as sparse_combine(), but for updating a hash value to
      // the state before it was changed with a call to sparse_combine(). This is
      // easily possible as the hash value just uses finite field arithmetic. We
      // can simply add the additive inverse of the previous hash value. This is a
      // very useful routine for symmetry detection. During partition refinement
      // the hashes do not need to be recomputed but can be updated with this
      // procedure.
  
      // make sure input value is never zero and at most 31bits are used
      value = (pair_hash<0>(value, value >> 32) >> 33) | 1;
  
      u32 a = c[index & 63] & M31();
      HighsInt degree = (index >> 6) + 1;
      // add the additive inverse (M31() - hashvalue) instead of the hash value
      // itself
      hash += M31() - multiply_modM31(value, modexp_M31(a, degree));
      hash = (hash >> 31) + (hash & M31());
      if (hash >= M31()) hash -= M31();
      assert(hash < M31());
    }
  
    static constexpr u64 fibonacci_muliplier() { return u64{0x9e3779b97f4a7c15}; }
  
    template <typename T,
              typename std::enable_if<HighsHashable<T>::value, int>::type = 0>
    static u64 vector_hash(const T* vals, size_t numvals) {
      std::array<u32, 2> pair{};
      u64 hash = 0;
      HighsInt k = 0;
  
      const char* dataptr = (const char*)vals;
      const char* dataend = (const char*)(vals + numvals);
  
      while (dataptr != dataend) {
        using std::size_t;
        size_t numBytes = std::min(size_t(dataend - dataptr), size_t{256});
        size_t numPairs = (numBytes + 7) / 8;
        size_t lastPairBytes = numBytes - (numPairs - 1) * 8;
        u64 chunkhash[] = {u64{0}, u64{0}};
  
  #define HIGHS_VECHASH_CASE_N(N, B)                         \
    std::memcpy(&pair[0], dataptr, B);                       \
    chunkhash[N & 1] += pair_hash<32 - N>(pair[0], pair[1]); \
    dataptr += B;
  
        switch (numPairs) {
          case 32:
            if (hash != 0) {
              // make sure hash is reduced mod M61() before multiplying with the
              // next random constant. For vectors at most 240 bytes we never
              // get here and only use the fast pair hashing scheme
              // for vectors with 240 bytes to 256 bytes we do have the one
              // additional check for hash != 0 above which will return false
              // and only for longer vectors we ever reduce modulo M61
              if (hash >= M61()) hash -= M61();
              hash = multiply_modM61(hash, c[(k++) & 63] & M61());
            }
            HIGHS_VECHASH_CASE_N(32, 8)
            // fall through
          case 31:
            HIGHS_VECHASH_CASE_N(31, 8)
            // fall through
          case 30:
            HIGHS_VECHASH_CASE_N(30, 8)
            // fall through
          case 29:
            HIGHS_VECHASH_CASE_N(29, 8)
            // fall through
          case 28:
            HIGHS_VECHASH_CASE_N(28, 8)
            // fall through
          case 27:
            HIGHS_VECHASH_CASE_N(27, 8)
            // fall through
          case 26:
            HIGHS_VECHASH_CASE_N(26, 8)
            // fall through
          case 25:
            HIGHS_VECHASH_CASE_N(25, 8)
            // fall through
          case 24:
            HIGHS_VECHASH_CASE_N(24, 8)
            // fall through
          case 23:
            HIGHS_VECHASH_CASE_N(23, 8)
            // fall through
          case 22:
            HIGHS_VECHASH_CASE_N(22, 8)
            // fall through
          case 21:
            HIGHS_VECHASH_CASE_N(21, 8)
            // fall through
          case 20:
            HIGHS_VECHASH_CASE_N(20, 8)
            // fall through
          case 19:
            HIGHS_VECHASH_CASE_N(19, 8)
            // fall through
          case 18:
            HIGHS_VECHASH_CASE_N(18, 8)
            // fall through
          case 17:
            HIGHS_VECHASH_CASE_N(17, 8)
            // fall through
          case 16:
            HIGHS_VECHASH_CASE_N(16, 8)
            // fall through
          case 15:
            HIGHS_VECHASH_CASE_N(15, 8)
            // fall through
          case 14:
            HIGHS_VECHASH_CASE_N(14, 8)
            // fall through
          case 13:
            HIGHS_VECHASH_CASE_N(13, 8)
            // fall through
          case 12:
            HIGHS_VECHASH_CASE_N(12, 8)
            // fall through
          case 11:
            HIGHS_VECHASH_CASE_N(11, 8)
            // fall through
          case 10:
            HIGHS_VECHASH_CASE_N(10, 8)
            // fall through
          case 9:
            HIGHS_VECHASH_CASE_N(9, 8)
            // fall through
          case 8:
            HIGHS_VECHASH_CASE_N(8, 8)
            // fall through
          case 7:
            HIGHS_VECHASH_CASE_N(7, 8)
            // fall through
          case 6:
            HIGHS_VECHASH_CASE_N(6, 8)
            // fall through
          case 5:
            HIGHS_VECHASH_CASE_N(5, 8)
            // fall through
          case 4:
            HIGHS_VECHASH_CASE_N(4, 8)
            // fall through
          case 3:
            HIGHS_VECHASH_CASE_N(3, 8)
            // fall through
          case 2:
            HIGHS_VECHASH_CASE_N(2, 8)
            // fall through
          case 1:
            HIGHS_VECHASH_CASE_N(1, lastPairBytes)
        }
  
        hash += (chunkhash[0] >> 3) ^ (chunkhash[1] >> 32);
      }
  
  #undef HIGHS_VECHASH_CASE_N
  
      return hash * fibonacci_muliplier();
    }
  
    template <typename T,
              typename std::enable_if<HighsHashable<T>::value &&
                                          (sizeof(T) <= 8) && (sizeof(T) >= 1),
                                      int>::type = 0>
    static u64 hash(const T& val) {
      std::array<u32, 2> bytes;
      if (sizeof(T) < 4) bytes[0] = 0;
      if (sizeof(T) < 8) bytes[1] = 0;
      std::memcpy(&bytes[0], &val, sizeof(T));
      return pair_hash<1>(bytes[0], bytes[1]) ^
             pair_hash<0>(bytes[0], bytes[1]) >> 32;
    }
  
    template <typename T,
              typename std::enable_if<HighsHashable<T>::value &&
                                          (sizeof(T) >= 9) && (sizeof(T) <= 16),
                                      int>::type = 0>
    static u64 hash(const T& val) {
      std::array<u32, 4> bytes;
      if (sizeof(T) < 12) bytes[2] = 0;
      if (sizeof(T) < 16) bytes[3] = 0;
      std::memcpy(&bytes[0], &val, sizeof(T));
      return (pair_hash<0>(bytes[0], bytes[1]) ^
              (pair_hash<1>(bytes[2], bytes[3]) >> 32)) *
             fibonacci_muliplier();
    }
  
    template <typename T,
              typename std::enable_if<HighsHashable<T>::value &&
                                          (sizeof(T) >= 17) && (sizeof(T) <= 24),
                                      int>::type = 0>
    static u64 hash(const T& val) {
      std::array<u32, 6> bytes;
      if (sizeof(T) < 20) bytes[4] = 0;
      if (sizeof(T) < 24) bytes[5] = 0;
      std::memcpy(&bytes[0], &val, sizeof(T));
      return (pair_hash<0>(bytes[0], bytes[1]) ^
              ((pair_hash<1>(bytes[2], bytes[3]) +
                pair_hash<2>(bytes[4], bytes[5])) >>
               32)) *
             fibonacci_muliplier();
    }
  
    template <typename T,
              typename std::enable_if<HighsHashable<T>::value &&
                                          (sizeof(T) >= 25) && (sizeof(T) <= 32),
                                      int>::type = 0>
    static u64 hash(const T& val) {
      std::array<u32, 8> bytes;
      if (sizeof(T) < 28) bytes[6] = 0;
      if (sizeof(T) < 32) bytes[7] = 0;
      std::memcpy(&bytes[0], &val, sizeof(T));
      return ((pair_hash<0>(bytes[0], bytes[1]) +
               pair_hash<1>(bytes[2], bytes[3])) ^
              ((pair_hash<2>(bytes[4], bytes[5]) +
                pair_hash<3>(bytes[6], bytes[7])) >>
               32)) *
             fibonacci_muliplier();
    }
  
    template <typename T,
              typename std::enable_if<HighsHashable<T>::value &&
                                          (sizeof(T) >= 33) && (sizeof(T) <= 40),
                                      int>::type = 0>
    static u64 hash(const T& val) {
      std::array<u32, 10> bytes;
      if (sizeof(T) < 36) bytes[8] = 0;
      if (sizeof(T) < 40) bytes[9] = 0;
      std::memcpy(&bytes[0], &val, sizeof(T));
      return ((pair_hash<0>(bytes[0], bytes[1]) +
               pair_hash<1>(bytes[2], bytes[3])) ^
              ((pair_hash<2>(bytes[4], bytes[5]) +
                pair_hash<3>(bytes[6], bytes[7]) +
                pair_hash<4>(bytes[8], bytes[9])) >>
               32)) *
             fibonacci_muliplier();
    }
  
    template <typename T,
              typename std::enable_if<HighsHashable<T>::value &&
                                          (sizeof(T) >= 41) && (sizeof(T) <= 48),
                                      int>::type = 0>
    static u64 hash(const T& val) {
      std::array<u32, 12> bytes;
      if (sizeof(T) < 44) bytes[10] = 0;
      if (sizeof(T) < 48) bytes[11] = 0;
      std::memcpy(&bytes[0], &val, sizeof(T));
      return ((pair_hash<0>(bytes[0], bytes[1]) +
               pair_hash<1>(bytes[2], bytes[3]) +
               pair_hash<2>(bytes[4], bytes[5])) ^
              ((pair_hash<3>(bytes[6], bytes[7]) +
                pair_hash<4>(bytes[8], bytes[9]) +
                pair_hash<5>(bytes[10], bytes[11])) >>
               32)) *
             fibonacci_muliplier();
    }
  
    template <typename T,
              typename std::enable_if<HighsHashable<T>::value &&
                                          (sizeof(T) >= 49) && (sizeof(T) <= 56),
                                      int>::type = 0>
    static u64 hash(const T& val) {
      std::array<u32, 14> bytes;
      if (sizeof(T) < 52) bytes[12] = 0;
      if (sizeof(T) < 56) bytes[13] = 0;
      std::memcpy(&bytes[0], &val, sizeof(T));
      return ((pair_hash<0>(bytes[0], bytes[1]) +
               pair_hash<1>(bytes[2], bytes[3]) +
               pair_hash<2>(bytes[4], bytes[5])) ^
              ((pair_hash<3>(bytes[6], bytes[7]) +
                pair_hash<4>(bytes[8], bytes[9]) +
                pair_hash<5>(bytes[10], bytes[11]) +
                pair_hash<6>(bytes[12], bytes[13])) >>
               32)) *
             fibonacci_muliplier();
    }
  
    template <typename T,
              typename std::enable_if<HighsHashable<T>::value &&
                                          (sizeof(T) >= 57) && (sizeof(T) <= 64),
                                      int>::type = 0>
    static u64 hash(const T& val) {
      std::array<u32, 16> bytes;
      if (sizeof(T) < 60) bytes[14] = 0;
      if (sizeof(T) < 64) bytes[15] = 0;
      std::memcpy(&bytes[0], &val, sizeof(T));
      return ((pair_hash<0>(bytes[0], bytes[1]) +
               pair_hash<1>(bytes[2], bytes[3]) +
               pair_hash<2>(bytes[4], bytes[5]) +
               pair_hash<3>(bytes[6], bytes[7])) ^
              ((pair_hash<4>(bytes[8], bytes[9]) +
                pair_hash<5>(bytes[10], bytes[11]) +
                pair_hash<6>(bytes[12], bytes[13]) +
                pair_hash<7>(bytes[14], bytes[15])) >>
               32)) *
             fibonacci_muliplier();
    }
  
    template <typename T,
              typename std::enable_if<HighsHashable<T>::value && (sizeof(T) > 64),
                                      int>::type = 0>
    static u64 hash(const T& val) {
      return vector_hash(&val, 1);
    }
  
    template <typename T,
              typename std::enable_if<HighsHashable<T>::value, int>::type = 0>
    static u64 hash(const std::vector<T>& val) {
      return vector_hash(val.data(), val.size());
    }
  
    template <typename T, typename std::enable_if<
                              std::is_same<decltype(*reinterpret_cast<T*>(0) ==
                                                    *reinterpret_cast<T*>(0)),
                                           bool>::value,
                              int>::type = 0>
    static bool equal(const T& a, const T& b) {
      return a == b;
    }
  
    template <typename T,
              typename std::enable_if<HighsHashable<T>::value, int>::type = 0>
    static bool equal(const std::vector<T>& a, const std::vector<T>& b) {
      if (a.size() != b.size()) return false;
      return std::memcmp(a.data(), b.data(), sizeof(T) * a.size()) == 0;
    }
  
    static constexpr double golden_ratio_reciprocal() {
      return 0.61803398874989484;
    }
  
    static u32 double_hash_code(double val) {
      // we multiply by some irrational number, so that the buckets in which we
      // put the real numbers do not break on a power of two pattern. E.g.
      // consider the use case for detecting parallel rows when we have two
      // parallel rows scaled to have their largest coefficient 1.0 and another
      // coefficient which is 0.5
      // +- epsilon. Clearly we want to detect those rows as parallel and give
      // them the same hash value for small enough epsilon. The exponent,
      // however will switch to -2 for the value just below 0.5 and the hashcodes
      // will differ. when multiplying with the reciprocal of the golden ratio the
      // exact 0.5 will yield 0.30901699437494742 and 0.5 - 1e-9 will yield
      // 0.3090169937569134 which has the same exponent and matches in the most
      // significant bits. Hence it yields the same hashcode. Obviously there will
      // now be different values which exhibit the same pattern as the 0.5 case,
      // but they do not have a small denominator like 1/2 in their rational
      // representation but are power of two multiples of the golden ratio and
      // therefore irrational, which we do not expect in non-artifical input data.
      int exponent;
      double hashbits = std::frexp(val * golden_ratio_reciprocal(), &exponent);
  
      // some extra casts to be more verbose about what is happening.
      // We want the exponent to use only 16bits so that the remaining 16 bits
      // are used for the most significant bits of the mantissa and the sign bit.
      // casting to unsigned 16bits first ensures that the value after the cast is
      // defined to be UINT16_MAX - |exponent| when the exponent is negative.
      // casting the exponent to a uint32_t directly would give wrong promotion
      // of negative exponents as UINT32_MAX - |exponent| and take up to many bits
      // or possibly loose information after the 16 bit shift. For the mantissa we
      // take the 15 most significant bits, even though we could squeeze out a few
      // more of the exponent. We don't need more bits as this would make the
      // buckets very small and might miss more values that are equal within
      // epsilon. Therefore the most significant 15 bits of the mantissa and the
      // sign is encoded in the 16 lower bits of the hashcode and the upper 16bits
      // encode the sign and value of the exponent.
      u32 hashvalue = (u32)(u16)(i16)exponent;
      hashvalue = (hashvalue << 16) | (u32)(u16)(i16)std::ldexp(hashbits, 15);
  
      return hashvalue;
    }
  };
  
  struct HighsHasher {
    template <typename T>
    size_t operator()(const T& x) const {
      return HighsHashHelpers::hash(x);
    }
  };
  
  struct HighsVectorHasher {
    template <typename T>
    size_t operator()(const std::vector<T>& vec) const {
      return HighsHashHelpers::vector_hash(vec.data(), vec.size());
    }
  };
  
  struct HighsVectorEqual {
    template <typename T>
    bool operator()(const std::vector<T>& vec1,
                    const std::vector<T>& vec2) const {
      if (vec1.size() != vec2.size()) return false;
      return std::equal(vec1.begin(), vec1.end(), vec2.begin());
    }
  };
  
  template <typename K, typename V>
  struct HighsHashTableEntry {
   private:
    K key_;
    V value_;
  
   public:
    template <typename K_>
    HighsHashTableEntry(K_&& k) : key_(k), value_() {}
    template <typename K_, typename V_>
    HighsHashTableEntry(K_&& k, V_&& v) : key_(k), value_(v) {}
  
    const K& key() const { return key_; }
    const V& value() const { return value_; }
    V& value() { return value_; }
  };
  
  template <typename T>
  struct HighsHashTableEntry<T, void> {
   private:
    T value_;
  
   public:
    template <typename... Args>
    HighsHashTableEntry(Args&&... args) : value_(std::forward<Args>(args)...) {}
  
    const T& key() const { return value_; }
    const T& value() const { return value_; }
  };
  
  template <typename K, typename V = void>
  class HighsHashTable {
    struct OpNewDeleter {
      void operator()(void* ptr) { ::operator delete(ptr); }
    };
  
   public:
    using u8 = std::uint8_t;
    using i8 = std::int8_t;
  
    using u16 = std::uint16_t;
    using i16 = std::int16_t;
  
    using u32 = std::uint32_t;
    using i32 = std::int32_t;
  
    using u64 = std::uint64_t;
    using i64 = std::int64_t;
  
    using Entry = HighsHashTableEntry<K, V>;
    using KeyType = K;
    using ValueType = typename std::remove_reference<decltype(
        reinterpret_cast<Entry*>(0)->value())>::type;
  
    std::unique_ptr<Entry, OpNewDeleter> entries;
    std::unique_ptr<u8[]> metadata;
    u64 tableSizeMask;
    u64 numHashShift;
    u64 numElements = 0;
  
    template <typename IterType>
    class HashTableIterator {
      u8* pos;
      u8* end;
      Entry* entryEnd;
  
     public:
      using difference_type = std::ptrdiff_t;
      using value_type = IterType;
      using pointer = IterType*;
      using reference = IterType&;
      using iterator_category = std::forward_iterator_tag;
      HashTableIterator(u8* pos, u8* end, Entry* entryEnd)
          : pos(pos), end(end), entryEnd(entryEnd) {}
      HashTableIterator() = default;
  
      HashTableIterator<IterType> operator++(int) {
        // postfix
        HashTableIterator<IterType> oldpos = *this;
        for (++pos; pos != end; ++pos)
          if ((*pos) & 0x80u) break;
  
        return oldpos;
      }
  
      HashTableIterator<IterType>& operator++() {
        // prefix
        for (++pos; pos != end; ++pos)
          if ((*pos) & 0x80u) break;
  
        return *this;
      }
  
      reference operator*() const { return *(entryEnd - (end - pos)); }
      pointer operator->() const { return (entryEnd - (end - pos)); }
      HashTableIterator<IterType> operator+(difference_type v) const {
        for (difference_type k = 0; k != v; ++k) ++(*this);
      }
  
      bool operator==(const HashTableIterator<IterType>& rhs) const {
        return pos == rhs.pos;
      }
      bool operator!=(const HashTableIterator<IterType>& rhs) const {
        return pos != rhs.pos;
      }
    };
  
    using const_iterator = HashTableIterator<const Entry>;
    using iterator = HashTableIterator<Entry>;
  
    HighsHashTable() { makeEmptyTable(128); }
    HighsHashTable(u64 minCapacity) {
      u64 initCapacity = u64{1} << (u64)std::ceil(
                             std::log2(std::max(128.0, 8 * minCapacity / 7.0)));
      makeEmptyTable(initCapacity);
    }
  
    iterator end() {
      u64 capacity = tableSizeMask + 1;
      return iterator{metadata.get() + capacity, metadata.get() + capacity,
                      entries.get() + capacity};
    };
  
    const_iterator end() const {
      u64 capacity = tableSizeMask + 1;
      return const_iterator{metadata.get() + capacity, metadata.get() + capacity,
                            entries.get() + capacity};
    };
  
    const_iterator begin() const {
      if (numElements == 0) return end();
      u64 capacity = tableSizeMask + 1;
      const_iterator iter{metadata.get(), metadata.get() + capacity,
                          entries.get() + capacity};
      if (!occupied(metadata[0])) ++iter;
  
      return iter;
    };
  
    iterator begin() {
      if (numElements == 0) return end();
      u64 capacity = tableSizeMask + 1;
      iterator iter{metadata.get(), metadata.get() + capacity,
                    entries.get() + capacity};
      if (!occupied(metadata[0])) ++iter;
  
      return iter;
    };
  
   private:
    u8 toMetadata(u64 hash) const { return (hash >> numHashShift) | 0x80u; }
  
    static constexpr u64 maxDistance() { return 127; }
  
    void makeEmptyTable(u64 capacity) {
      tableSizeMask = capacity - 1;
      numHashShift = 64 - HighsHashHelpers::log2i(capacity);
      assert(capacity == (u64{1} << (64 - numHashShift)));
      numElements = 0;
  
      metadata = decltype(metadata)(new u8[capacity]{});
      entries =
          decltype(entries)((Entry*)::operator new(sizeof(Entry) * capacity));
    }
  
    bool occupied(u8 meta) const { return meta & 0x80; }
  
    u64 distanceFromIdealSlot(u64 pos) const {
      // we store 7 bits of the hash in the metadata. Assuming a decent
      // hashfunction it is practically never happening that an item travels more
      // then 127 slots from its ideal position, therefore, we can compute the
      // distance from the ideal position just as it would normally be done
      // assuming there is at most one overflow. Consider using 3 bits which gives
      // values from 0 to 7. When an item is at a position with lower bits 7 and
      // is placed 3 positions after its ideal position, the lower bits of the
      // hash value will overflow and yield the value 2. With the assumption that
      // an item never cycles through one full cycle of the range 0 to 7, its
      // position would never be placed in a position with lower bits 7 other than
      // its ideal position. This allows us to compute the distance from its ideal
      // position by simply ignoring an overflow. In our case the correct answer
      // would be 3, but we get (2 - 7)=-5. This, however, is the correct result 3
      // when promoting to an unsigned value and looking at the lower 3 bits.
  
      return ((pos - metadata[pos])) & 0x7f;
    }
  
    void growTable() {
      decltype(entries) oldEntries = std::move(entries);
      decltype(metadata) oldMetadata = std::move(metadata);
      u64 oldCapactiy = tableSizeMask + 1;
  
      makeEmptyTable(2 * oldCapactiy);
  
      for (u64 i = 0; i != oldCapactiy; ++i)
        if (occupied(oldMetadata[i])) insert(std::move(oldEntries.get()[i]));
    }
  
    void shrinkTable() {
      decltype(entries) oldEntries = std::move(entries);
      decltype(metadata) oldMetadata = std::move(metadata);
      u64 oldCapactiy = tableSizeMask + 1;
  
      makeEmptyTable(oldCapactiy / 2);
  
      for (u64 i = 0; i != oldCapactiy; ++i)
        if (occupied(oldMetadata[i])) insert(std::move(oldEntries.get()[i]));
    }
  
    bool findPosition(const KeyType& key, u8& meta, u64& startPos, u64& maxPos,
                      u64& pos) const {
      u64 hash = HighsHashHelpers::hash(key);
      startPos = hash >> numHashShift;
      maxPos = (startPos + maxDistance()) & tableSizeMask;
      meta = toMetadata(hash);
  
      const Entry* entryArray = entries.get();
      pos = startPos;
      do {
        if (!occupied(metadata[pos])) return false;
        if (metadata[pos] == meta &&
            HighsHashHelpers::equal(key, entryArray[pos].key()))
          return true;
  
        u64 currentDistance = (pos - startPos) & tableSizeMask;
        if (currentDistance > distanceFromIdealSlot(pos)) return false;
  
        pos = (pos + 1) & tableSizeMask;
      } while (pos != maxPos);
  
      return false;
    }
  
   public:
    void clear() {
      if (numElements) {
        u64 capacity = tableSizeMask + 1;
        for (u64 i = 0; i < capacity; ++i)
          if (occupied(metadata[i])) entries.get()[i].~Entry();
        makeEmptyTable(128);
      }
    }
  
    const ValueType* find(const KeyType& key) const {
      u64 pos, startPos, maxPos;
      u8 meta;
      if (findPosition(key, meta, startPos, maxPos, pos))
        return &(entries.get()[pos].value());
  
      return nullptr;
    }
  
    ValueType* find(const KeyType& key) {
      u64 pos, startPos, maxPos;
      u8 meta;
      if (findPosition(key, meta, startPos, maxPos, pos))
        return &(entries.get()[pos].value());
  
      return nullptr;
    }
  
    ValueType& operator[](const KeyType& key) {
      Entry* entryArray = entries.get();
      u64 pos, startPos, maxPos;
      u8 meta;
      if (findPosition(key, meta, startPos, maxPos, pos))
        return entryArray[pos].value();
  
      if (numElements == ((tableSizeMask + 1) * 7) / 8 || pos == maxPos) {
        growTable();
        return (*this)[key];
      }
  
      using std::swap;
      ValueType& insertLocation = entryArray[pos].value();
      Entry entry(key, ValueType());
      ++numElements;
  
      do {
        if (!occupied(metadata[pos])) {
          metadata[pos] = meta;
          new (&entryArray[pos]) Entry{std::move(entry)};
          return insertLocation;
        }
  
        u64 currentDistance = (pos - startPos) & tableSizeMask;
        u64 distanceOfCurrentOccupant = distanceFromIdealSlot(pos);
        if (currentDistance > distanceOfCurrentOccupant) {
          // steal the position
          swap(entry, entryArray[pos]);
          swap(meta, metadata[pos]);
  
          startPos = (pos - distanceOfCurrentOccupant) & tableSizeMask;
          maxPos = (startPos + maxDistance()) & tableSizeMask;
        }
        pos = (pos + 1) & tableSizeMask;
      } while (pos != maxPos);
  
      growTable();
      insert(std::move(entry));
      return (*this)[key];
    }
  
    template <typename... Args>
    bool insert(Args&&... args) {
      Entry entry(std::forward<Args>(args)...);
  
      u64 pos, startPos, maxPos;
      u8 meta;
      if (findPosition(entry.key(), meta, startPos, maxPos, pos)) return false;
  
      if (numElements == ((tableSizeMask + 1) * 7) / 8 || pos == maxPos) {
        growTable();
        return insert(std::move(entry));
      }
  
      using std::swap;
      Entry* entryArray = entries.get();
      ++numElements;
  
      do {
        if (!occupied(metadata[pos])) {
          metadata[pos] = meta;
          new (&entryArray[pos]) Entry{std::move(entry)};
          return true;
        }
  
        u64 currentDistance = (pos - startPos) & tableSizeMask;
        u64 distanceOfCurrentOccupant = distanceFromIdealSlot(pos);
        if (currentDistance > distanceOfCurrentOccupant) {
          // steal the position
          swap(entry, entryArray[pos]);
          swap(meta, metadata[pos]);
  
          startPos = (pos - distanceOfCurrentOccupant) & tableSizeMask;
          maxPos = (startPos + maxDistance()) & tableSizeMask;
        }
        pos = (pos + 1) & tableSizeMask;
      } while (pos != maxPos);
  
      growTable();
      insert(std::move(entry));
      return true;
    }
  
    bool erase(const KeyType& key) {
      u64 pos, startPos, maxPos;
      u8 meta;
      if (!findPosition(key, meta, startPos, maxPos, pos)) return false;
      // delete element at position pos
      Entry* entryArray = entries.get();
      entryArray[pos].~Entry();
      metadata[pos] = 0;
  
      // retain at least a quarter of slots occupied, otherwise shrink the table
      // if its not at its minimum size already
      --numElements;
      u64 capacity = tableSizeMask + 1;
      if (capacity != 128 && numElements < capacity / 4) {
        shrinkTable();
        return true;
      }
  
      // shift elements after pos backwards
      while (true) {
        u64 shift = (pos + 1) & tableSizeMask;
        if (!occupied(metadata[shift])) return true;
  
        u64 dist = distanceFromIdealSlot(shift);
        if (dist == 0) return true;
  
        entryArray[pos] = std::move(entryArray[shift]);
        metadata[pos] = metadata[shift];
        metadata[shift] = 0;
        pos = shift;
      }
    }
  
    i64 size() const { return numElements; }
  
    HighsHashTable(HighsHashTable<K, V>&&) = default;
    HighsHashTable<K, V>& operator=(HighsHashTable<K, V>&&) = default;
  
    ~HighsHashTable() {
      if (metadata) {
        u64 capacity = tableSizeMask + 1;
        for (u64 i = 0; i < capacity; ++i) {
          if (occupied(metadata[i])) entries.get()[i].~Entry();
        }
      }
    }
  };
  
  #endif
  