qbpp.hpp

Go to the documentation of this file.

 
#ifndef QUBOPP_HPP
#define QUBOPP_HPP
 
#define VERSION "2025-01-15"
 
#include <tbb/blocked_range.h>
#include <tbb/parallel_for.h>
#include <tbb/parallel_reduce.h>
#include <tbb/parallel_sort.h>
 
#include <algorithm>
#include <boost/container/small_vector.hpp>
#include <boost/container/static_vector.hpp>
#include <boost/endian/conversion.hpp>
#include <boost/multi_array.hpp>
#include <boost/multiprecision/cpp_int.hpp>
#include <chrono>
#include <cmath>
#include <initializer_list>
#include <iomanip>
#include <iostream>
#include <limits>
#include <list>
#include <memory>
#include <mutex>
#include <optional>
#include <sstream>
#include <string>
#include <thread>
#include <type_traits>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <variant>
#include <vector>
 
 
 
namespace qbpp {
 
class Inf;
template <typename T>
class Vector;
class Var;
class Term;
class Expr;
class VarIntCore;
class VarInt;
class Sol;
class Model;
class QuadModel;
class SolHolder;
 
namespace impl {
class VarSet;
class BitVector;
class IndexVarMapper;
struct var_hash;
struct vars_hash;
}  
 
 
 
template <typename... Args>
std::string file_line(const char *file, int line, Args... args) {
  std::ostringstream oss;
  oss << file << "(" << line << ") ";
  (oss << ... << args);
  return oss.str();
}
 
#define THROW_MESSAGE(...) file_line(__FILE__, __LINE__, __VA_ARGS__)
 
#define TERM_CAPACITY 4
 
 
#ifndef MAXDEG
#define MAXDEG 0
#endif
 
#if MAXDEG == 0
using Vars = boost::container::small_vector<Var, 2>;
#else
 
using Vars = boost::container::static_vector<Var, MAXDEG>;
#endif
 
using Terms = std::vector<Term>;
 
using int128_t = boost::multiprecision::int128_t;
using uint128_t = boost::multiprecision::uint128_t;
using int256_t = boost::multiprecision::int256_t;
using uint256_t = boost::multiprecision::uint256_t;
using int512_t = boost::multiprecision::int512_t;
using uint512_t = boost::multiprecision::uint512_t;
using int1024_t = boost::multiprecision::int1024_t;
using uint1024_t = boost::multiprecision::uint1024_t;
using cpp_int = boost::multiprecision::cpp_int;
 
#ifndef COEFF_TYPE
#define COEFF_TYPE int32_t
#endif
 
#ifndef ENERGY_TYPE
#define ENERGY_TYPE int64_t
#endif
 
using vindex_t = uint32_t;
 
const vindex_t vindex_limit = std::numeric_limits<vindex_t>::max();
 
using var_val_t = int8_t;
 
using coeff_t = COEFF_TYPE;
 
using energy_t = ENERGY_TYPE;
 
 
 
using MapList = std::list<std::pair<std::variant<Var, VarInt>, Expr>>;
 
using MapDict = std::unordered_map<Var, Expr, impl::var_hash>;
 
using VarValMap = std::unordered_map<Var, var_val_t, impl::var_hash>;
 
using VarIndexMap = std::unordered_map<Var, vindex_t, impl::var_hash>;
 
using VarExprMap = std::unordered_map<Var, Expr, impl::var_hash>;
 
using VarsCoeffMap = std::unordered_map<Vars, coeff_t, impl::vars_hash>;
 
namespace impl {
 
struct var_hash {
  size_t operator()(const Var var) const;
};
 
struct vars_hash {
  size_t operator()(const Vars &vars) const;
};
 
}  
 
 
 
 
Var var(const std::string &var_str);
Var var();
 
std::string str(Var var);
std::string str(const char *param);
std::string str(const Vars &vars);
std::string str(const Term &term);
std::string str(const Expr &expr, const std::string &prefix);
std::string str(const Model &model);
std::string str(const QuadModel &quad_model);
std::string str(const Sol &sol);
std::string str(const MapList &map_list);
std::string str_short(const Expr &expr);
 
std::ostream &operator<<(std::ostream &os, Var var);
std::ostream &operator<<(std::ostream &os, const Term &term);
std::ostream &operator<<(std::ostream &os, const Expr &expr);
std::ostream &operator<<(std::ostream &os, const Model &model);
std::ostream &operator<<(std::ostream &os, const QuadModel &quad_model);
std::ostream &operator<<(std::ostream &os, const Sol &sol);
 
void sort_vars_in_place(Vars &vars);
Vars sort_vars(const Vars &vars);
Vars sort_vars_as_binary(const Vars &vars);
Vars sort_vars_as_spin(const Vars &vars);
Expr simplify(const Expr &expr, Vars (*sort_vars_func)(const Vars &));
Expr simplify_as_binary(const Expr &expr);
Expr simplify_as_spin(const Expr &expr);
bool is_simplified(const Expr &expr);
bool is_binary(const Expr &expr);
 
template <typename T>
Vector<Expr> sqr(const Vector<T> &arg);
template <typename T>
auto sqr(const Vector<Vector<T>> &arg) -> Vector<decltype(sqr(arg[0]))>;
 
energy_t eval(const Expr &expr, const Sol &sol);
energy_t eval(const Expr &expr, const MapList &map_list);
Expr reduce(const Expr &expr);
 
Expr binary_to_spin(const Expr &expr);
Expr spin_to_binary(const Expr &expr);
 
template <typename T>
Vector<Expr> binary_to_spin(const Vector<T> &arg);
template <typename T>
auto binary_to_spin(const Vector<Vector<T>> &arg)
    -> Vector<decltype(binary_to_spin(arg[0]))>;
 
template <typename T>
Vector<Expr> spin_to_binary(const Vector<T> &arg);
template <typename T>
auto spin_to_binary(const Vector<Vector<T>> &arg)
    -> Vector<decltype(spin_to_binary(arg[0]))>;
 
Vector<Expr> get_row(const Vector<Vector<Expr>> &vec, vindex_t index);
Vector<Expr> get_col(const Vector<Vector<Expr>> &vec, size_t index);
 
double get_time();
vindex_t all_var_count();
 
template <typename T>
Expr sum(const Vector<T> &items);
 
template <typename T>
Vector<Vector<Expr>> transpose(const Vector<Vector<T>> &vec);
 
Expr operator-(const Expr &expr);
Expr operator-(Expr &&expr);
 
energy_t toInt(const Expr &expr);
VarValMap list_to_var_val(const MapList &map_list);
 
 
inline std::string str(int8_t val) {
  return std::to_string(static_cast<int32_t>(val));
}
 
template <typename T>
auto str(const T &val) -> decltype(std::to_string(val)) {
  return std::to_string(val);
}
 
template <typename T>
auto str(T val) -> decltype(val.str()) {
  return val.str();
}
 
template <typename T>
T abs(T val) {
  return val < 0 ? -val : val;
}
 
inline energy_t gcd(energy_t a, energy_t b) {
  a = abs(a);
  b = abs(b);
  if (a == 0) return b;
  if (b == 0) return a;
  while (b != 0) {
    energy_t r = a % b;
    a = b;
    b = r;
  }
  return a;
}
 
 
class Inf {
  bool positive_ = true;
 
 public:
  Inf() = default;
 
  bool is_positive() const { return positive_; }
 
  bool is_negative() const { return !positive_; }
 
  Inf operator+() const { return *this; }
 
  Inf operator-() const {
    Inf inf;
    inf.positive_ = !positive_;
    return inf;
  }
};
 
const Inf inf;
 
template <typename T>
class Vector {
  std::vector<T> data_;
 
  template <typename U, typename Op>
  Vector<T> &vector_operation(const Vector<U> &rhs, Op operation) {
    if (size() != rhs.size()) {
      throw std::out_of_range(
          THROW_MESSAGE("Vector size mismatch: ", size(), " != ", rhs.size()));
    }
    tbb::parallel_for(tbb::blocked_range<size_t>(0, size()),
                      [&](const tbb::blocked_range<size_t> &range) {
                        for (size_t i = range.begin(); i < range.end(); ++i) {
                          operation((*this)[i], rhs[i]);
                          (*this)[i] = static_cast<T>((*this)[i]);
                        }
                      });
    return *this;
  }
 
  template <typename U, typename Op>
  Vector<T> &vector_operation(Vector<U> &&rhs, Op operation) {
    if (size() != rhs.size()) {
      throw std::out_of_range(
          THROW_MESSAGE("Vector size mismatch: ", size(), " != ", rhs.size()));
    }
 
    tbb::parallel_for(tbb::blocked_range<size_t>(0, size()),
                      [&](const tbb::blocked_range<size_t> &range) {
                        for (size_t i = range.begin(); i < range.end(); ++i) {
                          operation((*this)[i], std::move(rhs[i]));
                          (*this)[i] = static_cast<T>((*this)[i]);
                        }
                      });
 
    return *this;
  }
 
  template <typename Op>
  Vector<T> &vector_operation(const Expr &rhs, Op operation) {
    tbb::parallel_for(tbb::blocked_range<size_t>(0, size()),
                      [&](const tbb::blocked_range<size_t> &range) {
                        for (size_t i = range.begin(); i < range.end(); ++i) {
                          operation((*this)[i], rhs);
                        }
                      });
 
    return *this;
  }
 
 public:
  Vector() = default;
 
  Vector(const Vector<T> &) = default;
 
  Vector(Vector<T> &&) = default;
 
  Vector(std::initializer_list<T> init) : data_(init) {}
 
  Vector &operator=(const Vector<T> &) = default;
 
  Vector &operator=(Vector<T> &&) = default;
 
  explicit Vector(size_t size) : data_(size) {}
 
  template <typename U, typename = typename std::enable_if<
                            std::is_integral<U>::value>::type>
  Vector(U size, const T &value) : data_(static_cast<size_t>(size), value) {}
 
  template <typename U, typename = typename std::enable_if<
                            !std::is_integral<U>::value>::type>
  Vector(U begin, U end) : data_(begin, end) {}
 
  template <typename U>
  Vector(const Vector<U> &rhs) {
    data_.resize(rhs.size());
    tbb::parallel_for(size_t(0), rhs.size(),
                      [&](size_t i) { data_[i] = rhs[i]; });
  }
 
  const std::vector<T> get_data() const { return data_; }
 
  std::vector<T> get_data() { return data_; }
 
  void resize(size_t size) { data_.resize(size); }
 
  void push_back(const T &t) { data_.push_back(t); }
 
  void emplace_back(T &&t) { data_.emplace_back(std::forward<T>(t)); }
 
  void reserve(size_t size) { data_.reserve(size); }
 
  T &operator[](size_t i) { return data_[i]; }
 
  const T &operator[](size_t i) const { return data_[i]; }
 
  size_t size() const { return data_.size(); }
 
  auto begin() const { return data_.begin(); }
 
  auto end() const { return data_.end(); }
 
  auto begin() { return data_.begin(); }
 
  auto end() { return data_.end(); }
 
  typename std::vector<T>::iterator erase(
      typename std::vector<T>::iterator pos) {
    return data_.erase(pos);
  }
 
  typename std::vector<T>::iterator erase(
      typename std::vector<T>::iterator first,
      typename std::vector<T>::iterator last) {
    return data_.erase(first, last);
  }
 
 
  Vector<T> &operator=(const Expr &rhs) {
    tbb::parallel_for(size_t(0), size(), [&](size_t i) { (*this)[i] = rhs; });
    return *this;
  }
 
  Vector<T> &operator=(energy_t rhs) {
    tbb::parallel_for(size_t(0), size(), [&](size_t i) { (*this)[i] = rhs; });
    return *this;
  }
 
  template <typename U>
  Vector<T> &operator+=(const Vector<U> &rhs) {
    return vector_operation(rhs, [](T &larg, const U &rarg) { larg += rarg; });
  }
 
  template <typename U>
  Vector<T> &operator+=(Vector<U> &&rhs) {
    return vector_operation(std::move(rhs),
                            [](T &lhs_elem, U rhs_elem) mutable {
                              lhs_elem += std::move(rhs_elem);
                            });
  }
 
  template <typename U>
  Vector<T> &operator-=(const Vector<U> &rhs) {
    return vector_operation(rhs, [](T &larg, const U &rarg) { larg -= rarg; });
  }
 
  template <typename U>
  Vector<T> &operator-=(Vector<U> &&rhs) {
    return vector_operation(std::move(rhs),
                            [](T &lhs_elem, U rhs_elem) mutable {
                              lhs_elem -= std::move(rhs_elem);
                            });
  }
 
  template <typename U>
  Vector<T> &operator*=(const Vector<U> &rhs) {
    return vector_operation(rhs, [](T &larg, const U &rarg) { larg *= rarg; });
  }
 
  template <typename U>
  Vector<T> &operator*=(Vector<U> &&rhs) {
    return vector_operation(std::move(rhs),
                            [](T &lhs_elem, U rhs_elem) mutable {
                              lhs_elem *= std::move(rhs_elem);
                            });
  }
 
  Vector<T> &operator+=(const Expr &expr) {
    return vector_operation(expr,
                            [](T &elem, const Expr &rhs) { elem += rhs; });
  }
 
  Vector<T> &operator-=(const Expr &expr) {
    return vector_operation(expr,
                            [](T &elem, const Expr &rhs) { elem -= rhs; });
  }
 
  Vector<T> &operator*=(const Expr &expr) {
    return vector_operation(expr,
                            [](T &elem, const Expr &rhs) { elem *= rhs; });
  }
 
  Vector<T> &operator/=(energy_t val) {
    tbb::parallel_for(tbb::blocked_range<size_t>(0, size()),
                      [&](const tbb::blocked_range<size_t> &range) {
                        for (size_t i = range.begin(); i < range.end(); ++i) {
                          (*this)[i] /= val;
                        }
                      });
    return *this;
  }
 
  Vector<T> &sqr() {
    *this *= *this;
    return *this;
  }
 
  Vector<T> &simplify(Vars (*sort_vars_func)(const Vars &) = sort_vars);
 
  Vector<T> &simplify_as_binary() {
    *this = simplify(sort_vars_as_binary);
    return *this;
  }
 
  Vector<T> &simplify_as_spin() {
    *this = simplify(sort_vars_as_spin);
    return *this;
  }
 
  Vector<T> &binary_to_spin() {
    *this = qbpp::binary_to_spin(*this);
    return *this;
  }
 
  Vector<T> &spin_to_binary() {
    *this = qbpp::spin_to_binary(*this);
    return *this;
  }
 
  Vector<T> &replace(const MapList &map_list);
 
  Vector<T> &reduce();
 
  Vector<T> &transpose() {
    *this = transpose(*this);
    return *this;
  }
};
 
class Var {
  vindex_t index_;
 
  Var() = delete;
 
 public:
  Var(const Var &var) = default;
 
  Var &operator=(const Var &var) = default;
 
  explicit Var(vindex_t index) : index_(index) {}
  vindex_t get_index() const { return index_; }
 
  bool operator==(Var var) const { return index_ == var.index_; }
 
  bool operator!=(Var var) const { return index_ != var.index_; }
  bool operator<(Var var) const { return index_ < var.index_; }
 
  bool operator>(Var var) const { return index_ > var.index_; }
 
  bool operator<=(Var var) const { return index_ <= var.index_; }
 
  size_t size() const { return 1; }
};  
 
namespace impl {
 
 
class VarSet {
  std::vector<std::string> var_str_;
 
  vindex_t unnamed_vars{0};
 
  VarSet() {}
 
  VarSet(const impl::VarSet &) = delete;
 
  VarSet &operator=(const impl::VarSet &) = delete;
 
 public:
  static VarSet &get_instance() {
    static std::mutex mutex_;
    static VarSet singleton_var_set;
    std::lock_guard<std::mutex> lock(mutex_);
    return singleton_var_set;
  }
 
  static vindex_t all_var_count() {
    return static_cast<vindex_t>(get_instance().var_str_.size());
  }
 
  static Var var(const std::string &var_str) {
    static std::mutex mutex_;
    std::lock_guard<std::mutex> lock(mutex_);
    vindex_t new_index = get_instance().all_var_count();
    get_instance().var_str_.push_back(var_str);
    return Var(new_index);
  }
 
  static std::string new_unnamed_var_str() {
    static std::mutex mutex_;
    std::lock_guard<std::mutex> lock(mutex_);
    return "{" + std::to_string(get_instance().unnamed_vars++) + "}";
  }
 
  static std::string str(Var var) {
    return get_instance().var_str_[var.get_index()];
  }
};  
 
class IndexVarMapper {
  const std::vector<Var> index_var_;
  const std::vector<vindex_t> var_index_;
 
  std::vector<Var> compute_index_var(const Expr &expr);
 
  std::vector<vindex_t> compute_var_index(const std::vector<Var> &index_var);
 
  IndexVarMapper() = delete;
 
 public:
  explicit IndexVarMapper(const Expr &expr)
      : index_var_(compute_index_var(expr)),
        var_index_(compute_var_index(index_var_)) {}
 
  const std::vector<Var> &get_index_var() const { return index_var_; }
 
  const std::vector<vindex_t> &get_var_index() const { return var_index_; }
 
  vindex_t var_count() const {
    return static_cast<vindex_t>(index_var_.size());
  }
 
  vindex_t get_index(Var var) const {
    return static_cast<vindex_t>(var_index_[var.get_index()]);
  }
 
  Var get_var(vindex_t index) const { return get_index_var()[index]; }
 
  bool has(Var var) const {
    return var_index_[var.get_index()] != vindex_limit;
  }
};
}  
 
class Term {
  coeff_t coeff_{1};
 
  Vars vars_;
 
 public:
 
  Term() = default;
 
  Term(const Term &) = default;
 
  Term(Term &&) noexcept = default;
 
  explicit Term(coeff_t coeff) : coeff_(coeff) {}
 
  explicit Term(Var var, coeff_t coeff = 1) : coeff_(coeff), vars_({var}) {}
 
 
  explicit Term(const Vars &vars, coeff_t coeff = 1)
      : coeff_(coeff), vars_(vars) {}
 
 
  explicit Term(Vars &&vars, coeff_t coeff = 1) noexcept
      : coeff_(coeff), vars_(std::move(vars)) {}
 
 
  Term &operator=(const Term &) = default;
 
  Term &operator=(Term &&) noexcept = default;
 
  coeff_t get_coeff() const { return coeff_; }
 
  void set_coeff(coeff_t coeff) { coeff_ = coeff; }
 
  Vars &get_vars() { return vars_; }
 
  const Vars &get_vars() const { return vars_; }
 
  vindex_t var_count() const { return static_cast<vindex_t>(vars_.size()); }
 
  bool operator==(const Term &term) const {
    return (coeff_ == term.coeff_) && (vars_ == term.vars_);
  }
 
  bool operator!=(const Term &term) const { return !(*this == term); }
 
  bool operator<(const Term &term) const {
    if (var_count() < term.var_count())
      return true;
    else if (var_count() > term.var_count())
      return false;
    return vars_ < term.vars_;
  }
 
  bool operator<=(const Term &term) const {
    return *this < term || *this == term;
  }
 
  Term &operator*=(coeff_t val) {
    coeff_ *= val;
    return *this;
  }
 
 
  Term &operator*=(const Term &term) {
    coeff_ *= term.coeff_;
    vars_.reserve(var_count() + term.var_count());
    vars_.insert(vars_.end(), term.vars_.begin(), term.vars_.end());
    return *this;
  }
 
  Term &operator/=(energy_t val) {
    if constexpr (std::is_integral<coeff_t>::value) {
      if ((coeff_ / val) * val != coeff_) {
        throw std::runtime_error(
            THROW_MESSAGE("Indivisible division occurred."));
      }
    }
    coeff_ /= static_cast<coeff_t>(val);
    return *this;
  }
 
  Term operator*(energy_t val) const {
    Term result = *this;
    result.coeff_ *= static_cast<coeff_t>(val);
    return result;
  }
 
  Term operator*(const Term &term) const {
    Term result = *this;
    result.coeff_ *= term.coeff_;
    result.vars_.insert(result.vars_.end(), term.vars_.begin(),
                        term.vars_.end());
    return result;
  }
 
  Term &operator-() && {
    coeff_ = -coeff_;
    return *this;
  }
 
  Term operator-() const & {
    Term result = *this;
    result.coeff_ = -result.coeff_;
    return result;
  }
};
 
class Expr {
 
  energy_t constant_{0};
 
  Terms terms_;
 
  Expr(bool b) = delete;
 
 
 public:
 
  Expr() { terms_.reserve(TERM_CAPACITY); };
 
  Expr(const Expr &expr, size_t term_capacity = TERM_CAPACITY) {
    constant_ = expr.constant_;
    terms_.reserve(term_capacity);
    terms_ = expr.terms_;
  };
 
  Expr(Expr &&expr) noexcept = default;
 
  template <typename T,
            typename std::enable_if<std::is_convertible<T, energy_t>::value,
                                    int>::type = 0>
  Expr(T value, size_t term_capacity = TERM_CAPACITY) {
    constant_ = static_cast<energy_t>(value);
    terms_.reserve(term_capacity);
  }
 
  Expr(Var var) : Expr() { terms_.emplace_back(Term(var)); }
 
  Expr(const Term &term, size_t term_capacity = TERM_CAPACITY) {
    terms_.reserve(term_capacity);
    terms_.push_back(term);
  }
 
  Expr(Term &&term, size_t term_capacity = TERM_CAPACITY) noexcept {
    terms_.reserve(term_capacity);
    terms_.push_back(std::move(term));
  }
 
  Expr(energy_t constant, Terms &&terms) noexcept
      : constant_(constant), terms_(std::move(terms)) {}
 
  Expr(energy_t constant, const Terms &terms) noexcept
      : constant_(constant), terms_(terms) {}
 
  Expr(const VarInt &var_int);
 
  Expr(VarInt &&var_int);
 
 
  Expr &operator=(const Expr &expr) {
    if (this != &expr) {
      constant_ = expr.constant_;
      terms_ = expr.terms_;
    }
    return *this;
  }
 
  Expr &operator=(Expr &&expr) noexcept {
    if (this != &expr) {
      constant_ = expr.constant_;
      terms_ = std::move(expr.terms_);
    }
    return *this;
  }
 
 
  template <typename T,
            typename std::enable_if<std::is_convertible<T, energy_t>::value,
                                    int>::type = 0>
  Expr &operator*=(T val) {
    if (val == 0) {
      constant_ = 0;
      terms_.clear();
    } else {
      constant_ *= val;
      if (terms_.size() < 1000) {
        for (auto &t : terms_) {
          t *= static_cast<coeff_t>(val);
        }
      } else {
        tbb::parallel_for(size_t(0), terms_.size(), [&](size_t i) {
          terms_[i] *= static_cast<coeff_t>(val);
        });
      }
    }
    return *this;
  }
 
  Expr &operator*=(const Term &term) {
    if (term.var_count() == 0) {
      return *this *= term.get_coeff();
    }
    if (terms_.size() < 1000) {
      for (auto &t : terms_) {
        t *= term;
      }
    } else {
      tbb::parallel_for(size_t(0), terms_.size(),
                        [&](size_t i) { terms_[i] *= term; });
    }
    if (constant_ != 0) {
      *this += term * constant_;
      constant_ = 0;
    }
    return *this;
  }
 
 
 
  Expr &operator*=(const Expr &expr) {
    if (expr.term_count() == 0) {
      if (term_count() < 1000) {
        return *this *= expr.constant_;
      }
      constant_ *= expr.constant_;
      tbb::parallel_for(size_t(0), terms_.size(), [&](size_t i) {
        terms_[i] *= static_cast<coeff_t>(expr.constant_);
      });
      return *this;
    }
 
    if (expr.term_count() == 1 && expr.constant_ == 0) {
      return *this *= expr.terms_[0];
    }
 
    const Expr &expr_small =
        this->term_count() <= expr.term_count() ? *this : expr;
    const Expr &expr_large =
        this->term_count() <= expr.term_count() ? expr : *this;
 
    Expr result(expr_small.constant_ * expr_large.constant_);
 
    size_t total_terms = expr_large.term_count() * expr_small.term_count();
 
    if (total_terms < 1000) {
      result.terms_.reserve(total_terms);
      for (const auto &term1 : expr_large.terms_) {
        for (const auto &term2 : expr_small.terms_) {
          result.terms_.push_back(term1 * term2);
        }
      }
      if (expr_large.constant_ != 0) {
        result.terms_.reserve(total_terms += expr_small.term_count());
        for (const auto &term : expr_small.terms_) {
          result.terms_.push_back(term * expr_large.constant_);
        }
      }
 
      if (expr_small.constant_ != 0) {
        result.terms_.reserve(total_terms += expr_large.term_count());
        for (const auto &term : expr_large.terms_) {
          result.terms_.push_back(term * expr_small.constant_);
        }
      }
      *this = std::move(result);
      return *this;
    }
 
    result.terms_.resize(total_terms);
    tbb::parallel_for(size_t(0), expr_large.term_count(), [&](size_t i) {
      for (size_t j = 0; j < expr_small.term_count(); ++j) {
        result.terms_[i * expr_small.term_count() + j] =
            expr_large.terms_[i] * expr_small.terms_[j];
      }
    });
 
    if (expr_large.constant_ != 0) {
      result.terms_.resize(total_terms + expr_small.term_count());
      tbb::parallel_for(size_t(0), expr_small.term_count(), [&](size_t i) {
        result.terms_[total_terms + i] =
            expr_small.terms_[i] * expr_large.constant_;
      });
      total_terms += expr_small.term_count();
    }
 
    if (expr_small.constant_ != 0) {
      result.terms_.resize(total_terms + expr_large.term_count());
      tbb::parallel_for(size_t(0), expr_large.term_count(), [&](size_t i) {
        result.terms_[total_terms + i] =
            expr_large.terms_[i] * expr_small.constant_;
      });
    }
 
    *this = std::move(result);
    return *this;
  }
 
  Expr &operator+=(const Expr &expr) {
    constant_ += expr.constant_;
    auto expr2 = expr;
    terms_.insert(terms_.end(), expr2.terms_.begin(), expr2.terms_.end());
    return *this;
  }
 
  Expr &operator-=(const Expr &expr) {
    auto expr2 = -expr;
    *this += std::move(expr2);
    return *this;
  }
 
  Expr &operator/=(energy_t val) {
    if (val == 0) {
      throw std::runtime_error(THROW_MESSAGE("Division by zero."));
    }
    tbb::parallel_for(size_t(0), terms_.size(),
                      [&](size_t i) { terms_[i] /= val; });
    if ((constant_ / val) * val != constant_) {
      throw std::runtime_error(THROW_MESSAGE("Indivisible division occurred."));
    }
    constant_ /= val;
    return *this;
  }
 
 
 
  Expr &sqr() {
    *this *= *this;
    return *this;
  }
 
 
  Expr &simplify(Vars (*sort_vars)(const Vars &));
 
  Expr &simplify_as_binary();
 
  Expr &simplify_as_spin();
 
 
 
 
  energy_t get_constant() const { return constant_; }
 
  void set_constant(energy_t constant) { constant_ = constant; }
 
  const Terms &get_terms() const { return terms_; }
 
  Terms &get_terms() { return terms_; }
 
  size_t term_count() const { return terms_.size(); }
 
  size_t size() const { return term_count(); }
 
 
 
  Expr &replace(const MapList &map_list);
 
  Expr operator()(const MapList &map_list) const {
    auto result = *this;
    return result.replace(map_list);
  }
 
 
 
 
  Expr &reduce() { return *this = qbpp::reduce(*this); }
 
 
 
  Expr &binary_to_spin() { return *this = qbpp::binary_to_spin(*this); }
 
  Expr &spin_to_binary() { return *this = qbpp::spin_to_binary(*this); }
 
  energy_t pos_sum() const {
    Expr temp = qbpp::simplify_as_binary(*this);
    energy_t sum = temp.constant_;
    for (const auto &term : temp.terms_) {
      if (term.get_coeff() > 0) {
        sum += term.get_coeff();
      }
    }
    return sum;
  }
 
  energy_t neg_sum() const {
    Expr temp = qbpp::simplify_as_binary(*this);
    energy_t sum = temp.constant_;
    for (const auto &term : temp.terms_) {
      if (term.get_coeff() < 0) {
        sum += term.get_coeff();
      }
    }
    return sum;
  }
};
 
class VarIntCore {
 public:
  const std::string var_str_;
 
  explicit VarIntCore(const std::string &var_str) : var_str_(var_str) {}
 
  VarIntCore() : var_str_(impl::VarSet::new_unnamed_var_str()) {}
};
 
class VarInt {
 
  const std::string var_str_;
 
  const energy_t min_val_;
  const energy_t max_val_;
  const std::shared_ptr<std::vector<coeff_t>> coeffs_ptr_;
  const Vector<Var> vars_;
  const Expr expr_;
 
  const Vector<Var> new_vars() const {
    Vector<Var> vars;
    if (coeffs_ptr_->size() == 1) {
      vars.push_back(var(var_str_));
    } else {
      for (size_t i = 0; i < coeffs_ptr_->size(); i++) {
        vars.push_back(var(var_str_ + "[" + std::to_string(i) + "]"));
      }
    }
    return vars;
  }
 
  const Expr new_expr() const {
    Expr result = min_val_;
    for (size_t i = 0; i < coeffs_ptr_->size(); i++) {
      result += Term(vars_[i], (*coeffs_ptr_)[i]);
    }
    return result;
  }
 
 public:
  VarInt(const std::string &var_str, energy_t min_val, energy_t max_val,
         std::shared_ptr<std::vector<coeff_t>> coeffs_ptr)
      : var_str_(var_str),
        min_val_(min_val),
        max_val_(max_val),
        coeffs_ptr_(coeffs_ptr),
        vars_(new_vars()),
        expr_(new_expr()) {}
 
  VarInt(const VarInt &) = default;
 
  std::string get_name() const { return var_str_; }
 
  vindex_t var_count() const { return static_cast<vindex_t>(vars_.size()); }
 
  energy_t min_val() const { return min_val_; }
 
  energy_t max_val() const { return max_val_; }
 
  coeff_t get_coeff(vindex_t i) const { return (*coeffs_ptr_)[i]; }
 
  Var get_var(vindex_t i) const { return vars_[i]; }
 
  Var operator[](vindex_t i) const { return vars_[i]; }
 
  const Expr &get_expr() const { return expr_; }
 
  VarValMap get_val_map(energy_t val) const {
    if (val < min_val_ || val > max_val_) {
      throw std::out_of_range(
          THROW_MESSAGE("Value (", str(val), " out of range."));
    }
    val -= min_val_;
    VarValMap result;
    for (size_t i = var_count(); i >= 1; i--) {
      if (val >= (*coeffs_ptr_)[static_cast<size_t>(i) - 1]) {
        result[vars_[i - 1]] = 1;
        val -= (*coeffs_ptr_)[i - 1];
      } else {
        result[vars_[i - 1]] = 0;
      }
    }
    return result;
  }
};
 
class Model {
 protected:
  const std::shared_ptr<const Expr> expr_ptr_;
 
  const std::shared_ptr<const impl::IndexVarMapper> index_var_ptr_;
 
  static const Expr &check_expr(const Expr &expr) {
    if (is_simplified(expr)) {
      return expr;
    } else {
      throw std::runtime_error(THROW_MESSAGE(
          "The expression must be simplified to be a QUBO expression."));
    }
  }
 
  static Expr check_expr(Expr &&expr) {
    if (is_simplified(expr)) {
      return std::move(expr);
    } else {
      throw std::runtime_error(THROW_MESSAGE(
          "The expression must be simplified to be a QUBO expression."));
    }
  }
 
  Model() = delete;
 
  Model &operator=(const Model &) = delete;
 
  Model &operator=(Model &&) = delete;
 
 public:
  Model(const Model &model)
      : expr_ptr_(model.expr_ptr_), index_var_ptr_(model.index_var_ptr_) {}
 
 
  Model(const Expr &expr)
      : expr_ptr_(std::make_shared<const Expr>(expr)),
        index_var_ptr_(
            std::make_shared<const impl::IndexVarMapper>(*expr_ptr_)) {}
 
  Model(Model &&) = default;
 
  Model(Expr &&expr) noexcept
      : expr_ptr_(std::make_shared<const Expr>(std::move(expr))),
        index_var_ptr_(
            std::make_shared<const impl::IndexVarMapper>(*expr_ptr_)) {}
 
  virtual ~Model() = default;
 
  vindex_t var_count() const { return index_var_ptr_->var_count(); }
 
  Var get_var(vindex_t index) const { return index_var_ptr_->get_var(index); }
 
  vindex_t get_index(Var var) const { return index_var_ptr_->get_index(var); }
 
  bool has(Var var) const { return index_var_ptr_->has(var); }
 
  virtual size_t term_count() const { return expr_ptr_->term_count(); }
 
  const Expr &get_expr() const { return *expr_ptr_; }
 
  operator const Expr &() const { return *expr_ptr_; }
 
  energy_t get_constant() const { return expr_ptr_->get_constant(); }
 
  const Terms &get_terms() const { return expr_ptr_->get_terms(); }
 
  std::vector<Var> get_index_var() const {
    return index_var_ptr_->get_index_var();
  }
 
};
 
class QuadModel : public Model {
 private:
 
  struct Impl {
    std::vector<coeff_t> linear_;
 
    std::vector<vindex_t> degree_;
 
    std::vector<std::vector<std::pair<vindex_t, coeff_t>>> quadratic_;
 
    coeff_t min_coeff_;
 
    coeff_t max_coeff_;
 
    explicit Impl(const qbpp::Model &model);
  };
 
  const std::shared_ptr<const Impl> pimpl_;
 
  const std::vector<coeff_t> &linear_;
 
  const std::vector<vindex_t> &degree_;
 
  const std::vector<std::vector<std::pair<vindex_t, coeff_t>>> &quadratic_;
 
  QuadModel &operator=(QuadModel &&) = delete;
 
  QuadModel &operator=(const QuadModel &) = delete;
 
 public:
 
  QuadModel(const QuadModel &) = default;
 
  QuadModel(QuadModel &&) noexcept = default;
 
  QuadModel(const Model &model)
      : Model(model),
        pimpl_(std::make_shared<Impl>(model)),
        linear_(pimpl_->linear_),
        degree_(pimpl_->degree_),
        quadratic_(pimpl_->quadratic_) {}
 
  QuadModel(Model &&model)
      : Model(std::move(model)),
        pimpl_(std::make_shared<Impl>(*this)),
        linear_(pimpl_->linear_),
        degree_(pimpl_->degree_),
        quadratic_(pimpl_->quadratic_) {}
 
  QuadModel(const Expr &expr)
      : Model(expr),
        pimpl_(std::make_shared<Impl>(Impl(*this))),
        linear_(pimpl_->linear_),
        degree_(pimpl_->degree_),
        quadratic_(pimpl_->quadratic_) {}
 
  QuadModel(Expr &&expr)
      : Model(std::move(expr)),
        pimpl_(std::make_shared<Impl>(Impl(*this))),
        linear_(pimpl_->linear_),
        degree_(pimpl_->degree_),
        quadratic_(pimpl_->quadratic_) {}
 
  const std::vector<coeff_t> &get_linear() const { return linear_; }
 
  const std::vector<vindex_t> &get_degree() const { return degree_; }
 
  const std::vector<std::vector<std::pair<vindex_t, coeff_t>>> &get_quadratic()
      const {
    return quadratic_;
  }
 
  coeff_t get_linear(vindex_t index) const { return linear_[index]; }
 
  vindex_t get_degree(vindex_t index) const { return degree_[index]; }
 
  std::pair<vindex_t, coeff_t> get_quadratic(vindex_t i, vindex_t j) const {
    return quadratic_[i][j];
  }
 
  size_t term_count(vindex_t deg) const {
    if (deg == 1) {
      return linear_.size();
    } else if (deg == 2) {
      return std::accumulate(degree_.begin(), degree_.end(), 0UL);
    } else {
      return 0;
    }
  }
 
  size_t term_count() const override { return term_count(1) + term_count(2); }
 
  coeff_t get_min_coeff() const { return pimpl_->min_coeff_; }
 
  coeff_t get_max_coeff() const { return pimpl_->max_coeff_; }
 
};
 
namespace impl {
class BitVector {
  const vindex_t bit_count_;
 
  const vindex_t size64_;
 
  std::unique_ptr<uint64_t[]> bits_;
 
 public:
  explicit BitVector(vindex_t bit_count)
      : bit_count_(bit_count),
        size64_((bit_count + 63) / 64),
        bits_(std::make_unique<uint64_t[]>(size64_)) {}
 
  explicit BitVector(const BitVector &bit_vector)
      : bit_count_(bit_vector.bit_count_),
        size64_(bit_vector.size64_),
        bits_(std::make_unique<uint64_t[]>(size64_)) {
    std::copy(bit_vector.bits_.get(), bit_vector.bits_.get() + size64_,
              bits_.get());
  }
 
  BitVector(BitVector &&bit_vector)
      : bit_count_(bit_vector.bit_count_),
        size64_(bit_vector.size64_),
        bits_(std::move(bit_vector.bits_)) {};
 
  BitVector &operator=(const BitVector &bit_vector) {
    if (this == &bit_vector) {
      return *this;
    }
    if (size64_ != bit_vector.size64_) {
      throw std::runtime_error(
          THROW_MESSAGE("BitVector: Different size of bit vector."));
    }
    std::copy(bit_vector.bits_.get(), bit_vector.bits_.get() + size64_,
              bits_.get());
    return *this;
  }
 
  BitVector &operator=(BitVector &bit_vector) {
    if (this == &bit_vector) {
      return *this;
    }
    if (size64_ != bit_vector.size64_) {
      throw std::runtime_error(
          THROW_MESSAGE("BitVector: Different size of bit vector."));
    }
    bits_ = std::move(bit_vector.bits_);
    return *this;
  }
 
  BitVector &operator=(BitVector &&bit_vector) {
    if (this == &bit_vector) {
      return *this;
    }
    if (bit_count_ != bit_vector.bit_count_) {
      throw std::runtime_error(
          THROW_MESSAGE("BitVector: Different size of bit vector."));
    }
    bits_ = std::move(bit_vector.bits_);
    return *this;
  };
 
  bool operator==(const BitVector &bit_vector) const {
    if (bit_count_ != bit_vector.bit_count_) {
      throw std::runtime_error(
          THROW_MESSAGE("BitVector: Different size of bit vector."));
    }
    return std::equal(bits_.get(), bits_.get() + size64_,
                      bit_vector.bits_.get());
  }
 
  bool operator!=(const BitVector &bit_vector) const {
    return !this->operator==(bit_vector);
  }
 
  bool operator<(const BitVector &bit_vector) const {
    auto reverse_bits = [](uint64_t n) -> uint64_t {
#if defined(__has_builtin)
#if __has_builtin(__builtin_bitreverse64)
      return __builtin_bitreverse64(n);
#endif
#endif
      n = ((n >> 1) & 0x5555555555555555) | ((n & 0x5555555555555555) << 1);
      n = ((n >> 2) & 0x3333333333333333) | ((n & 0x3333333333333333) << 2);
      n = ((n >> 4) & 0x0f0f0f0f0f0f0f0f) | ((n & 0x0f0f0f0f0f0f0f0f) << 4);
      n = ((n >> 8) & 0x00ff00ff00ff00ff) | ((n & 0x00ff00ff00ff00ff) << 8);
      n = ((n >> 16) & 0x0000ffff0000ffff) | ((n & 0x0000ffff0000ffff) << 16);
      n = ((n >> 32) & 0x00000000ffffffff) | ((n & 0x00000000ffffffff) << 32);
      return n;
    };
    for (vindex_t i = 0; i < size64_; i++) {
      if (reverse_bits(bits_[i]) < reverse_bits(bit_vector.bits_[i]))
        return true;
      if (reverse_bits(bits_[i]) > reverse_bits(bit_vector.bits_[i]))
        return false;
    }
    return false;
  }
 
  vindex_t size() const { return bit_count_; }
 
  vindex_t size64() const { return size64_; }
 
  bool get(vindex_t index) const {
    if (index >= bit_count_) {
      throw std::out_of_range(THROW_MESSAGE("Index out of range."));
    }
    return bits_[index / 64] & (1ULL << (index % 64));
  }
 
  void set(vindex_t index, bool value) {
    if (index >= bit_count_) {
      throw std::out_of_range(THROW_MESSAGE("Index out of range."));
    }
    if (value)
      bits_[index / 64] |= 1ULL << (index % 64);
    else
      bits_[index / 64] &= ~(1ULL << (index % 64));
  }
 
  void flip(vindex_t index) {
    if (index >= bit_count_) {
      throw std::out_of_range(THROW_MESSAGE("Index out of range."));
    }
    bits_[index / 64] ^= 1ULL << (index % 64);
  }
 
  uint64_t get64(vindex_t index64) const {
    if (index64 >= size64_) {
      throw std::out_of_range(THROW_MESSAGE("Index out of range."));
    }
    return bits_[index64];
  }
 
  uint64_t set64(vindex_t index64, uint64_t value) {
    if (index64 >= size64_) {
      throw std::out_of_range(THROW_MESSAGE("Index out of range."));
    }
    if (index64 < size64_ - 1 || index64 * 64 == bit_count_)
      return bits_[index64] = value;
    else
      return bits_[index64] = value & ((1ULL << (bit_count_ % 64)) - 1);
  }
 
  void clear() { std::fill(bits_.get(), bits_.get() + size64_, 0); }
 
  vindex_t popcount() const {
    vindex_t count = 0;
    for (vindex_t i = 0; i < size64_; i++) {
      count += static_cast<vindex_t>(__builtin_popcountll(bits_[i]));
    }
    return count;
  }
 
  std::string str() const {
    std::string result;
    for (vindex_t i = 0; i < bit_count_; i++) {
      result += get(i) ? "1" : "0";
    }
    return result;
  }
};
}  
 
class Sol {
 protected:
  const QuadModel quad_model_;
 
  impl::BitVector bit_vector_;
 
  mutable std::optional<energy_t> energy_ = std::nullopt;
 
  Sol() = delete;
 
  energy_t comp_energy() const;
 
  MapList get_map_list() const;
 
 public:
 
  Sol(const Sol &sol) = default;
 
  Sol(Sol &&sol) = default;
 
  Sol(const QuadModel &quad_model)
      : quad_model_(quad_model),
        bit_vector_(quad_model.var_count()),
        energy_(quad_model.get_constant()) {}
 
  virtual ~Sol() = default;
 
  Sol &operator=(const Sol &sol) {
    if (this == &sol) {
      return *this;
    }
    bit_vector_ = sol.bit_vector_;
    energy_ = sol.energy_;
    return *this;
  }
 
  Sol &operator=(Sol &&sol) {
    if (this == &sol) {
      return *this;
    }
    bit_vector_ = std::move(sol.bit_vector_);
    energy_ = std::move(sol.energy_);
    return *this;
  }
 
  bool operator==(const Sol &sol) const {
    if (get_energy() != sol.get_energy()) return false;
    return bit_vector_ == sol.bit_vector_;
  }
  bool operator<(const Sol &sol) const {
    if (get_energy() < sol.get_energy()) return true;
    if (get_energy() == sol.get_energy())
      return bit_vector_ < sol.bit_vector_;
    else
      return false;
  }
 
 
  var_val_t get(vindex_t index) const { return bit_vector_.get(index); }
 
  var_val_t get(Var var) const { return get(quad_model_.get_index(var)); }
 
  energy_t get(const Expr &expr) const { return eval(expr, *this); }
 
  bool has(Var var) const { return quad_model_.has(var); }
 
  void clear() {
    bit_vector_.clear();
    energy_ = quad_model_.get_constant();
  }
 
  template <typename T>
  auto get(const Vector<T> &vec) const {
    Vector<decltype(get(std::declval<T>()))> result;
    result.reserve(vec.size());
    for (const auto &elem : vec) {
      result.push_back(get(elem));
    }
    return result;
  }
 
  void set(vindex_t index, bool value) {
    energy_.reset();
    bit_vector_.set(index, value);
  }
 
  void set(Var var, bool value) {
    energy_.reset();
    set(quad_model_.get_index(var), value);
  }
 
  virtual void flip(vindex_t index) {
    energy_.reset();
    bit_vector_.flip(index);
  }
 
  void flip(Var var) {
    energy_.reset();
    bit_vector_.flip(quad_model_.get_index(var));
  }
 
  vindex_t popcount() const { return bit_vector_.popcount(); }
 
  energy_t get_energy() const {
    if (!energy_) energy_ = comp_energy();
    return *energy_;
  }
 
  void set_energy(energy_t energy) { energy_ = energy; }
 
  const impl::BitVector &get_bit_vector() const { return bit_vector_; }
 
 
 
  vindex_t var_count() const { return quad_model_.var_count(); }
 
  energy_t get_constant() const { return quad_model_.get_constant(); }
 
  Var get_var(vindex_t index) const { return quad_model_.get_var(index); }
 
  vindex_t get_index(Var var) const { return quad_model_.get_index(var); }
 
  operator MapList() const { return get_map_list(); }
 
  Sol set(const Sol &sol) {
    for (auto var : sol.quad_model_.get_index_var()) {
      set(var, sol.get(var));
    }
    return *this;
  }
 
  Sol set(const MapList &map_list) {
    VarValMap var_val_map = list_to_var_val(map_list);
    for (const auto &[var, val] : var_val_map) {
      set(var, val);
    }
    return *this;
  }
};
 
class SolHolder {
 protected:
  Sol sol_;
  std::optional<energy_t> bound_ = std::nullopt;
  double startx_ = get_time();
  double tts_{0.0};
  std::string solver_ = "N/A";
  mutable std::mutex mtx_;
 
 public:
  explicit SolHolder(const QuadModel &quad_model) : sol_(quad_model) {}
 
  SolHolder(const SolHolder &other) = delete;
 
 
  virtual ~SolHolder() = default;
 
 
  virtual std::optional<double> set_if_better(const qbpp::Sol &new_sol,
                                              const std::string &solver = "") {
    if (new_sol.get_energy() >= sol_.get_energy()) {
      return std::nullopt;
    }
    std::lock_guard<std::mutex> lock(mtx_);
    if (new_sol.get_energy() < sol_.get_energy()) {
      tts_ = get_time() - startx_;
      sol_ = new_sol;
      solver_ = solver;
      return std::make_optional(tts_);
    } else {
      return std::nullopt;
    }
  }
 
  virtual std::optional<qbpp::Sol> get_if_better(energy_t my_energy) {
    if (sol_.get_energy() >= my_energy) return std::nullopt;
    std::lock_guard<std::mutex> lock(mtx_);
    if (sol_.get_energy() < my_energy) {
      return sol_;
    } else {
      return std::nullopt;
    }
  }
 
  std::optional<qbpp::Sol> get_if_better(const Sol &my_sol) {
    return get_if_better(my_sol.get_energy());
  }
 
  std::tuple<Sol, double, std::string> get_sol_tts_solver() const {
    std::lock_guard<std::mutex> lock(mtx_);
    return std::make_tuple(sol_, tts_, solver_);
  }
 
  vindex_t var_count() const { return sol_.var_count(); }
 
  energy_t get_energy() const { return sol_.get_energy(); }
 
  const Sol &get_sol() const { return sol_; }
 
  Sol copy_sol() const {
    std::lock_guard<std::mutex> lock(mtx_);
    return sol_;
  }
 
  void clear() {
    std::lock_guard<std::mutex> lock(mtx_);
    sol_.clear();
    tts_ = 0;
    solver_ = "N/A";
  }
 
  double get_tts() const { return tts_; }
 
  std::string get_solver() const { return solver_; }
 
  void set_bound(energy_t bound) { bound_ = bound; }
 
  std::optional<energy_t> get_bound() const { return bound_; }
 
  operator Sol() const { return sol_; }
 
};
 
 
inline std::string str_impl(const Vars &vars,
                            std::function<std::string(Var)> str) {
  std::string result;
  for (const auto &item : vars) {
    if (!result.empty()) {
      result += "*";
    }
    result += str(item);
  }
  return result;
}
 
inline std::string str_impl(const Term &term,
                            std::function<std::string(Var)> str) {
  std::string result;
  if (term.var_count() == 0) return qbpp::str(term.get_coeff());
  if (term.get_coeff() == -1)
    result += "-";
  else if (term.get_coeff() != 1)
    result += qbpp::str(term.get_coeff()) + "*";
  result += str_impl(term.get_vars(), str);
  return result;
}
 
inline std::string str_impl(const Expr &expr,
                            std::function<std::string(Var)> str) {
  bool first = true;  
  std::string result;
  if (expr.get_constant() != 0) {
    result += qbpp::str(expr.get_constant());
    first = false;
  } else if (expr.term_count() == 0)
    return "0";
  for (const auto &term : expr.get_terms()) {
    if (first) {
      result += str_impl(term, str);
      first = false;
    } else {
      if (term.get_coeff() < 0)
        result += " " + str_impl(term, str);
      else
        result += " +" + str_impl(term, str);
    }
  }
  return result;
}
 
 
 
inline std::string str(Var var) { return impl::VarSet::str(var); }
 
inline std::string str(const Vars &vars) {
  return str_impl(vars, static_cast<std::string (*)(Var)>(str));
}
 
inline std::string str(const impl::BitVector &bit_vector) {
  std::string result;
  for (vindex_t i = 0; i < bit_vector.size(); i++) {
    result += bit_vector.get(i) ? "1" : "0";
  }
  return result;
}
 
 
inline std::string str(const char *param) {
  std::string param_str(param);
  if (param_str == "compilation_date") {
    return __DATE__;
  } else if (param_str == "author") {
    return "Koji Nakano";
  } else if (param_str == "version") {
    return VERSION;
  } else if (param_str == "license") {
    return "QUBO++: Only for non-commercial use, evaluation, and research "
           "purposes with no warranties. Redistribution is strictly "
           "prohibited.";
  } else if (param_str == "all_vars") {
    std::string result;
    for (vindex_t i = 0; i < impl::VarSet::get_instance().all_var_count();
         i++) {
      result += "{" + std::to_string(i) + "," + str(Var(i)) + "}";
    }
    return result;
  } else
    throw std::runtime_error(
        THROW_MESSAGE("Unknown parameter (", param, ") for str()."));
}
 
inline std::string str(const Term &term) {
  return str_impl(term, static_cast<std::string (*)(Var)>(str));
}
 
inline std::string str(const Expr &expr) {
  return str_impl(expr, static_cast<std::string (*)(Var)>(str));
}
 
inline std::string str(const Expr &expr, const std::string &prefix,
                       const std::string &separator [[maybe_unused]] = ",") {
  if (prefix == "") {
    return "{" + str_impl(expr, static_cast<std::string (*)(Var)>(str)) + "}";
  } else {
    return "{" + prefix + "," +
           str_impl(expr, static_cast<std::string (*)(Var)>(str)) + "}";
  }
}
 
inline std::string str(const Model &model) {
  auto str = [&model](Var var) -> std::string {
    return "<" + std::to_string(model.get_index(var)) + ">";
  };
  return str_impl(model.get_expr(), str);
}
 
inline std::string str(const QuadModel &quad_model) {
  auto str = [&quad_model](Var var) -> std::string {
    return "<" + qbpp::str(quad_model.get_index(var)) + ">";
  };
  std::string result = qbpp::str(quad_model.get_constant());
  for (vindex_t i = 0; i < quad_model.var_count(); i++) {
    result += " +" + str_impl(quad_model.get_var(i), str) + "*(" +
              qbpp::str(quad_model.get_linear(i));
    for (vindex_t j = 0; j < quad_model.get_degree(i); j++) {
      auto [index, coeff] = quad_model.get_quadratic(i, j);
      if (coeff == 1)
        result += " +";
      else if (coeff == -1)
        result += " -";
      else if (coeff > 1)
        result += " +" + qbpp::str(coeff) + "*";
      else
        result += " " + qbpp::str(coeff) + "*";
      result += str_impl(quad_model.get_var(index), str) + "/2";
    }
    result += ")";
  }
  return result;
}
 
inline std::string str(const MapList &map_list) {
  std::ostringstream oss;
  oss << "{";
  bool first = true;
 
  for (const auto &[key, val] : map_list) {
    if (!first) {
      oss << ",";
    }
    first = false;
 
    auto str_val = str(val);
    std::visit(
        [&](const auto &arg) {
          using T = std::decay_t<decltype(arg)>;
          oss << "{";
          if constexpr (std::is_same_v<T, VarInt>) {
            oss << arg.get_name();
          } else {
            oss << str(arg);
          }
          oss << "," << str_val << "}";
        },
        key);
  }
 
  oss << "}";
  return oss.str();
}
 
inline std::string str(const Sol &sol) {
  return qbpp::str(sol.get_energy()) + ":" + str(static_cast<MapList>(sol));
}
 
template <typename T>
std::string str(const Vector<T> &vec, const std::string &prefix = "",
                const std::string &separator = ",") {
  std::string result;
  for (size_t i = 0; i < vec.size(); i++) {
    result += str(vec[i], prefix + "[" + qbpp::str(i) + "]", separator);
    if (i < vec.size() - 1) result += separator;
  }
  return result;
}
 
inline std::ostream &operator<<(std::ostream &os, Var var) {
  return os << str(var);
}
 
inline std::ostream &operator<<(std::ostream &os, const VarInt &var_int) {
  return os << qbpp::Expr(var_int);
}
 
inline std::ostream &operator<<(std::ostream &os, const Term &term) {
  return os << str(term);
}
 
inline std::ostream &operator<<(std::ostream &os, const Expr &expr) {
  return os << str(expr);
}
 
inline std::ostream &operator<<(std::ostream &os, const Model &model) {
  return os << str(model);
}
 
inline std::ostream &operator<<(std::ostream &os, const QuadModel &quad_model) {
  return os << str(quad_model);
}
 
inline std::ostream &operator<<(std::ostream &os, const MapList &map_list) {
  return os << str(map_list);
}
 
inline std::ostream &operator<<(std::ostream &os, const Sol &sol) {
  return os << str(sol);
}
 
template <typename T>
std::ostream &operator<<(std::ostream &os, const Vector<T> &vec) {
  os << str(vec);
  return os;
}
 
template <typename T>
std::ostream &operator<<(std::ostream &os,
                         const std::pair<std::string, Vector<T>> &vec) {
  os << str(vec.second, vec.first);
  return os;
}
 
inline std::string str_short(const Expr &expr) {
  if (expr.term_count() <= 5) return str(expr);
  std::string result;
  if (expr.get_constant() != 0)
    result += qbpp::str(expr.get_constant()) + " + ";
  auto it = expr.get_terms().begin();
  result += str(*it) + " + ";
  ++it;
  result += str(*it) + " + ... + ";
  it = expr.get_terms().end();
  --it;
  --it;
  result += str(*it) + " + ";
  ++it;
  result += str(*it);
  result += " (" + qbpp::str(expr.term_count()) + " terms)";
  return result;
}
 
 
namespace impl {
inline size_t var_hash::operator()(const Var var) const {
  std::hash<vindex_t> hasher;
  return hasher(var.get_index());
}
 
inline size_t vars_hash::operator()(const Vars &vars) const {
  std::hash<vindex_t> hasher;
  size_t seed = 0;
  for (const auto &item : vars) {
    seed ^= hasher(item.get_index()) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
  }
  return seed;
}
}  
 
 
namespace impl {
template <typename T>
Expr equal(const T &lhs, const T &rhs) {
  return lhs - rhs == 0;
}
 
template <typename T>
Expr equal(const std::vector<T> &lhs, const std::vector<T> &rhs) {
  if (lhs.size() != rhs.size()) {
    throw std::runtime_error(THROW_MESSAGE("Vector dimension mismatch"));
  }
  Expr result = 0;
  for (size_t i = 0; i < lhs.size(); i++) {
    result += equal(lhs[i], rhs[i]);
  }
  return result;
}
 
}  
 
 
inline Expr &Expr::simplify(Vars (*sort_vars_func)(const Vars &) = sort_vars) {
  return *this = qbpp::simplify(*this, sort_vars_func);
}
 
inline Expr &Expr::simplify_as_binary() {
  return *this = qbpp::simplify_as_binary(*this);
}
 
inline Expr &Expr::simplify_as_spin() {
  return *this = qbpp::simplify_as_spin(*this);
}
 
inline Expr::Expr(const VarInt &var_int)
    : constant_(var_int.get_expr().get_constant()),
      terms_(var_int.get_expr().get_terms()) {
  if (var_int.max_val() == 0) {
    throw std::runtime_error(
        THROW_MESSAGE("Wrong VarInt object. Max value is not specified."));
  }
}
 
inline Expr::Expr(VarInt &&var_int)
    : constant_(var_int.get_expr().get_constant()),
      terms_(std::move(var_int.get_expr().get_terms())) {
  if (var_int.max_val() == 0) {
    throw std::runtime_error(
        THROW_MESSAGE("Wrong VarInt object. Max value is not specified."));
  }
}
 
 
inline VarValMap list_to_var_val(const MapList &map_list) {
  VarValMap var_int_map;
  for (const auto &[key, val] : map_list) {
    if (val.term_count() >= 1) {
      throw std::runtime_error(THROW_MESSAGE("Value must be an integer."));
    }
    auto constant = val.get_constant();
    std::visit(
        [&](auto &&arg) {
          using T = std::decay_t<decltype(arg)>;
          if constexpr (std::is_same_v<T, VarInt>) {
            for (const auto &[var, bin_val] :
                 arg.get_val_map(static_cast<coeff_t>(constant))) {
              var_int_map[var] = bin_val;
            }
          } else {
            var_int_map[arg] = static_cast<var_val_t>(constant);
          }
        },
        key);
  }
  return var_int_map;
}
 
inline MapDict list_to_dict(const MapList &map_list) {
  MapDict map_dict;
  for (const auto &[key, val] : map_list) {
    auto val_local = val;
    std::visit(
        [&](auto &&arg) {
          using T = std::decay_t<decltype(arg)>;
          if constexpr (std::is_same_v<T, VarInt>) {
            if (val_local.term_count() >= 1) {
              throw std::runtime_error(
                  THROW_MESSAGE("VarInt must have a constant."));
            }
            for (const auto &[var, bin_val] : arg.get_val_map(
                     static_cast<coeff_t>(val_local.get_constant()))) {
              map_dict[var] = bin_val;
            }
          } else {
            map_dict[arg] = val_local;
          }
        },
        key);
  }
  return map_dict;
}
 
inline Expr replace(const Term &term, const MapDict &map_dict) {
  Expr result = term.get_coeff();
  for (auto item : term.get_vars()) {
    auto it = map_dict.find(item);
    if (it == map_dict.end()) {
      result *= item;
    } else {
      result *= it->second;
    }
  }
  return result;
}
 
inline Expr &Expr::replace(const MapList &map_list) {
  MapDict map_dict = list_to_dict(map_list);
  Vector<Expr> exprs;
  exprs.resize(terms_.size());
  tbb::parallel_for(size_t(0), terms_.size(), [&](size_t i) {
    exprs[i] = qbpp::replace(terms_[i], map_dict);
  });
  Expr result = qbpp::sum(exprs);
  result += constant_;
  *this = std::move(result);
  return *this;
}
 
 
inline energy_t eval(const Term &term, const VarValMap &var_val_map) {
  energy_t result = term.get_coeff();
  for (const auto &item : term.get_vars()) {
    try {
      result *= var_val_map.at(item);
    } catch (const std::out_of_range &e) {
      throw std::runtime_error(THROW_MESSAGE("Variable ", qbpp::str(item),
                                             " is not found in the map."));
    }
  }
  return result;
}
 
inline energy_t eval(const Term &term, const Sol &sol) {
  energy_t result = term.get_coeff();
  for (const auto &item : term.get_vars()) {
    result *= sol.get(item);
  }
  return result;
}
 
 
inline Expr replace(const Expr &expr, const MapList &map_list) {
  auto result = expr;
  return result.replace(map_list);
}
 
inline energy_t eval_var_val_map(const Expr &expr,
                                 const VarValMap &var_val_map) {
  energy_t result = expr.get_constant();
 
  result += tbb::parallel_reduce(
      tbb::blocked_range<size_t>(0, expr.get_terms().size()),
      static_cast<energy_t>(0),
      [&](const tbb::blocked_range<size_t> &range, energy_t local_result) {
        for (size_t i = range.begin(); i < range.end(); ++i) {
          local_result += eval(expr.get_terms()[i], var_val_map);
        }
        return local_result;
      },
      [](energy_t lhs, energy_t rhs) { return lhs + rhs; });
 
  return result;
}
 
inline energy_t eval(const Expr &expr, const MapList &map_list) {
  VarValMap var_val_map = list_to_var_val(map_list);
  return eval_var_val_map(expr, var_val_map);
}
 
inline energy_t eval(const Expr &expr, const Sol &sol) {
  energy_t result = expr.get_constant();
  result += tbb::parallel_reduce(
      tbb::blocked_range<size_t>(0, expr.get_terms().size()),
      static_cast<energy_t>(0),
      [&](const tbb::blocked_range<size_t> &range, energy_t local_result) {
        for (size_t i = range.begin(); i < range.end(); ++i) {
          local_result += eval(expr.get_terms()[i], sol);
        }
        return local_result;
      },
      [](energy_t lhs, energy_t rhs) { return lhs + rhs; });
 
  return result;
}
 
 
 
inline std::vector<Var> impl::IndexVarMapper::compute_index_var(
    const Expr &expr) {
  std::vector<std::pair<int8_t, vindex_t>> var_set(all_var_count(), {1, 0});
  if (expr.term_count() == 0) return {};
 
  tbb::parallel_for(size_t(0), expr.get_terms().size(), [&](size_t i) {
    const auto &term = expr.get_terms()[i];
    for (const auto &var : term.get_vars()) {
      var_set[var.get_index()] = {0, var.get_index()};
    }
  });
 
  tbb::parallel_sort(var_set.begin(), var_set.end());
 
  vindex_t non_zero_count = all_var_count();
  for (vindex_t i = 1; i < all_var_count(); ++i) {
    if (var_set[i - 1].first == 0 && var_set[i].first == 1) {
      non_zero_count = i;
      break;
    }
  }
 
  std::vector<Var> index_var(non_zero_count, Var(vindex_limit));
  for (vindex_t i = 0; i < non_zero_count; ++i) {
    index_var[i] = Var(var_set[i].second);
  }
 
  return index_var;
}
 
inline std::vector<vindex_t> impl::IndexVarMapper::compute_var_index(
    const std::vector<Var> &index_var) {
  std::vector<vindex_t> var_index(all_var_count(), vindex_limit);
 
  if (index_var.empty()) return var_index;
 
  tbb::parallel_for(
      static_cast<vindex_t>(0), static_cast<vindex_t>(index_var.size()),
      [&](vindex_t i) { var_index[index_var[i].get_index()] = i; });
 
  return var_index;
}
 
 
 
inline Var var(const std::string &var_str) {
  auto name = var_str;
  return impl::VarSet::var(name);
}
 
inline Var var() { return var(impl::VarSet::new_unnamed_var_str()); }
 
template <typename T, typename... Args,
          typename = typename std::enable_if<std::is_integral<T>::value>::type>
auto var(const std::string &var_str, T size, Args... args) {
  if constexpr (sizeof...(Args) == 0) {
    Vector<Var> vars;
    if (size == 1) {
      vars.emplace_back(var(var_str));
    } else {
      for (vindex_t i = 0; i < static_cast<size_t>(size); ++i) {
        vars.emplace_back(var(var_str + "[" + qbpp::str(i) + "]"));
      }
    }
    return vars;
  } else {
    Vector<decltype(var(var_str, args...))> vars;
    vars.reserve(static_cast<size_t>(size));
    for (size_t i = 0; i < static_cast<size_t>(size); ++i) {
      vars.emplace_back(var(var_str + "[" + qbpp::str(i) + "]",
                            static_cast<size_t>(args)...));
    }
    return vars;
  }
}
 
template <typename T, typename... Args,
          typename = typename std::enable_if<std::is_integral<T>::value>::type>
auto var(T size, Args... args) {
  return var(impl::VarSet::new_unnamed_var_str(), size, args...);
}
 
inline vindex_t all_var_count() {
  return impl::VarSet::get_instance().all_var_count();
}
 
 
 
 
 
inline VarIntCore var_int(const std::string &var_str) {
  return VarIntCore(var_str);
}
 
inline VarIntCore var_int() {
  return var_int(impl::VarSet::new_unnamed_var_str());
}
 
 
template <typename T, typename... Args>
auto var_int(const std::string &var_str, T size, Args... args) {
  if constexpr (sizeof...(args) == 0) {
    Vector<VarIntCore> vars;
    vars.reserve(static_cast<size_t>(size));
    for (size_t i = 0; i < static_cast<size_t>(size); ++i) {
      vars.emplace_back(VarIntCore(var_str + "[" + qbpp::str(i) + "]"));
    }
    return vars;
  } else {
    Vector<decltype(var_int(var_str, args...))> vars;
    vars.reserve(static_cast<size_t>(size));
    for (size_t i = 0; i < static_cast<size_t>(size); ++i) {
      vars.emplace_back(var_int(var_str + "[" + qbpp::str(i) + "]",
                                static_cast<size_t>(args)...));
    }
    return vars;
  }
}
 
 
 
 
inline std::vector<coeff_t> comp_coeffs(energy_t min_val, energy_t max_val,
                                        energy_t base_coeff = 1) {
  std::vector<coeff_t> coeffs;
  coeff_t val = static_cast<coeff_t>(max_val - min_val);
  if (val <= 0) {
    throw std::invalid_argument(THROW_MESSAGE(
        "max_val(", max_val, ") must be greater than min_val(", min_val, ")"));
  }
  coeff_t current_coeff = static_cast<coeff_t>(base_coeff);
  while (val > 0) {
    if (current_coeff < val) {
      coeffs.push_back(current_coeff);
    } else {
      coeffs.push_back(val);
    }
    val -= current_coeff;
    current_coeff *= 2;
  }
  return coeffs;
}
 
inline VarInt new_var_int(const std::string &var_str, energy_t min_val,
                          energy_t max_val, energy_t base_coeff = 1) {
  std::shared_ptr<std::vector<coeff_t>> coeffs_ptr =
      std::make_shared<std::vector<coeff_t>>(
          comp_coeffs(min_val, max_val, base_coeff));
  return VarInt(var_str, min_val, max_val, coeffs_ptr);
}
 
inline VarInt new_var_int(const VarIntCore &var_int_core, energy_t min_val,
                          energy_t max_val, energy_t base_coeff = 1) {
  return new_var_int(var_int_core.var_str_, min_val, max_val, base_coeff);
}
 
inline auto operator<=(energy_t lhs, VarIntCore &&rhs) {
  return std::make_pair(lhs, std::move(rhs));
}
 
template <typename T>
auto operator<=(energy_t lhs, Vector<T> &&rhs) {
  return std::make_pair(lhs, std::forward<Vector<T>>(rhs));
}
 
template <typename T>
auto operator<=(const std::pair<energy_t, T> lhs, energy_t rhs) {
  return new_var_int(lhs.second.var_str_, lhs.first, rhs);
}
 
inline std::pair<energy_t, Vector<VarIntCore>> operator<=(
    energy_t lhs, const Vector<VarIntCore> &rhs) {
  return std::make_pair(lhs, rhs);
}
 
inline auto operator<=(const std::pair<energy_t, Vector<VarIntCore>> lhs,
                       energy_t rhs) {
  Vector<VarInt> result;
  result.reserve(lhs.second.size());
  for (const auto &item : lhs.second) {
    result.push_back(new_var_int(item, lhs.first, rhs));
  }
  return result;
}
 
void operator<=(energy_t , Var ) = delete;
 
void operator<=(Var , energy_t ) = delete;
 
 
 
inline Expr expr() { return Expr(); }
 
inline Vector<Expr> expr(vindex_t size) { return Vector<Expr>(size); }
 
template <typename T, typename... Args>
auto expr(T size, Args... args) {
  if constexpr (sizeof...(Args) == 0) {
    return Vector<Expr>(static_cast<size_t>(size));
  } else {
    Vector<decltype(expr(args...))> result(static_cast<size_t>(size));
    for (size_t i = 0; i < static_cast<size_t>(size); ++i)
      result[i] = expr(static_cast<size_t>(args)...);
    return result;
  }
}
 
template <typename T>
Expr toExpr(const T &arg) {
  return Expr(arg);
}
 
inline Expr toExpr(const Expr &arg) { return arg; }
 
inline Vector<Expr> toExpr(const Vector<Expr> &arg) { return arg; }
 
template <typename T>
auto toExpr(const Vector<T> &arg) {
  Vector<decltype(toExpr(arg[0]))> result;
  result.reserve(arg.size());
  for (const auto &item : arg) {
    result.push_back(toExpr(item));
  }
  return result;
}
 
template <typename T>
auto toExpr(const std::vector<T> &arg) {
  Vector<decltype(toExpr(arg[0]))> result;
  result.reserve(arg.size());
  for (const auto &item : arg) {
    result.push_back(toExpr(item));
  }
  return result;
}
 
template <typename T>
auto toExpr(const Vector<T> &arg) ->
    typename std::enable_if<std::is_same<T, decltype(toExpr(arg[0]))>::value,
                            Vector<T>>::type {
  return arg;
}
 
inline Vector<Expr> toExpr(const std::initializer_list<Expr> &list) {
  return Vector<Expr>(list);
}
 
template <typename T>
auto toExpr(const std::initializer_list<T> &list) {
  Vector<decltype(toExpr(std::declval<T>()))> result;
  result.reserve(list.size());
  for (const auto &item : list) {
    result.push_back(toExpr(item));
  }
  return result;
}
 
inline Vector<Vector<Expr>> toExpr(
    const std::initializer_list<std::initializer_list<Expr>> &list) {
  Vector<Vector<Expr>> result;
  result.reserve(list.size());
  for (const auto &item : list) {
    result.push_back(toExpr(item));
  }
  return result;
}
 
inline Vector<Vector<Vector<Expr>>> toExpr(
    const std::initializer_list<
        std::initializer_list<std::initializer_list<Expr>>> &list) {
  Vector<Vector<Vector<Expr>>> result;
  result.reserve(list.size());
  for (const auto &item : list) {
    result.push_back(toExpr(item));
  }
  return result;
}
 
inline Vector<Vector<Vector<Vector<Expr>>>> toExpr(
    const std::initializer_list<std::initializer_list<
        std::initializer_list<std::initializer_list<Expr>>>> &list) {
  Vector<Vector<Vector<Vector<Expr>>>> result;
  result.reserve(list.size());
  for (const auto &item : list) {
    result.push_back(toExpr(item));
  }
  return result;
}
 
 
 
 
inline Expr &&operator*(Expr &&lhs, Expr &&rhs) {
  if (lhs.term_count() >= rhs.term_count()) {
    lhs *= rhs;
    return std::move(lhs);
  } else {
    rhs *= lhs;
    return std::move(rhs);
  }
}
 
inline Expr &&operator*(Expr &&lhs, const Expr &rhs) {
  lhs *= rhs;
  return std::move(lhs);
}
 
inline Expr &&operator*(const Expr &lhs, Expr &&rhs) {
  rhs *= lhs;
  return std::move(rhs);
}
 
inline Expr operator*(const Expr &lhs, const Expr &rhs) {
  auto result = lhs;
  result *= rhs;
  return result;
}
 
inline Expr &&operator+(Expr &&lhs, Expr &&rhs) {
  lhs += rhs;
  return std::move(lhs);
}
 
inline Expr &&operator+(Expr &&lhs, const Expr &rhs) {
  lhs += rhs;
  return std::move(lhs);
}
 
inline Expr &&operator+(const Expr &lhs, Expr &&rhs) {
  rhs += lhs;
  return std::move(rhs);
}
 
inline Expr operator+(const Expr &lhs, const Expr &rhs) {
  auto result = lhs;
  return result += rhs;
}
 
inline Expr &&operator-(Expr &&lhs, Expr &&rhs) {
  lhs -= rhs;
  return std::move(lhs);
}
 
inline Expr &&operator-(Expr &&lhs, const Expr &rhs) {
  lhs -= rhs;
  return std::move(lhs);
}
 
inline Expr &&operator-(const Expr &lhs, Expr &&rhs) {
  rhs *= -1;
  rhs += lhs;
  return std::move(rhs);
}
 
inline Expr operator-(const Expr &lhs, const Expr &rhs) {
  auto result = lhs;
  return result -= rhs;
}
 
inline Expr &&operator/(Expr &&expr, energy_t val) {
  expr /= val;
  return std::move(expr);
}
 
inline Expr operator/(const Expr &expr, energy_t val) {
  auto temp = expr;
  return temp /= val;
}
 
 
inline Expr operator+(Expr &&expr) { return std::move(expr); }
 
inline const Expr &operator+(const Expr &expr) { return expr; }
 
inline Expr operator-(Expr &&expr) {
  expr.set_constant(-expr.get_constant());
  for (auto &term : expr.get_terms()) {
    term = -std::move(term);  
  }
  return std::move(expr);  
}
 
inline Expr operator-(const Expr &expr) {
  Expr result{-expr.get_constant()};
  for (const auto &term : expr.get_terms()) {
    result += -term;
  }
  return result;
}
 
 
 
 
inline Expr sqr(const Expr &expr) { return expr * expr; }
 
inline energy_t gcd(const Expr &expr) {
  const auto &terms = expr.get_terms();
  energy_t result = tbb::parallel_reduce(
        tbb::blocked_range<size_t>(0, terms.size()), 
        expr.get_constant(),                       
        [&](const tbb::blocked_range<size_t> &range, energy_t local_result) {
        for (size_t i = range.begin(); i < range.end(); ++i) {
            const auto &term = terms[i];
            local_result = gcd(local_result, term.get_coeff());
        }
        return local_result;
    },                               
        [](energy_t a, energy_t b) {                
            return gcd(a, b);
        }
    );
    return abs(result);
}
 
 
 
inline Expr operator==(const Expr &expr, energy_t val) {
  return sqr(expr - val);
}
 
inline Expr comparison(const Expr &expr, energy_t minimum, energy_t maximum) {
  energy_t val = maximum - minimum;
  Expr result;
  if (val < 0) {
    throw std::runtime_error(THROW_MESSAGE("RHS (", maximum,
                                           ") must be greater than LHS (",
                                           minimum, ") in operator <=."));
  } else if (val == 0) {
    result = (expr - minimum) == 0;
  } else if (val == 1) {
    result = (expr - minimum) * (expr - maximum);
  } else {
    VarInt slack_var = new_var_int(impl::VarSet::new_unnamed_var_str(), minimum,
                                   maximum - 1, 2);
    result = (expr - slack_var) * (expr - (slack_var + 1));
  }
  return result;
}
 
inline std::pair<energy_t, Expr> operator<=(energy_t min_val,
                                            const Expr &expr) {
  return std::make_pair(min_val, expr);
}
 
inline Expr operator<=(const std::pair<energy_t, Expr> &pair,
                       energy_t max_val) {
  energy_t min_val = pair.first;
  const Expr &result = pair.second;
  return comparison(result, min_val, max_val);
}
 
inline Expr operator<=(const std::pair<energy_t, Expr> &pair, Inf val) {
  energy_t min_val = pair.first;
  const Expr &result = pair.second;
  if (val.is_negative()) {
    throw std::runtime_error(
        THROW_MESSAGE("RHS of operator <= must not be -inf."));
  }
  return comparison(result, min_val, result.pos_sum());
}
 
inline std::pair<energy_t, Expr> operator<=(Inf val, const Expr &expr) {
  if (val.is_positive()) {
    throw std::runtime_error(
        THROW_MESSAGE("LHS of operator <= must not be inf."));
  }
  return std::make_pair(expr.neg_sum(), expr);
}
 
 
 
inline void sort_vars_in_place(Vars &vars) {
  if (vars.size() <= 1) return;
  if (vars.size() == 2)
    if (vars[0] > vars[1]) std::swap(vars[0], vars[1]);
  std::sort(vars.begin(), vars.end());
}
 
inline Vars sort_vars(const Vars &vars) {
  auto result = vars;
  sort_vars_in_place(result);
  return result;
}
 
inline Vars sort_vars_as_binary(const Vars &vars) {
  auto result = vars;
  sort_vars_in_place(result);
  if (result.size() >= 2)
    result.erase(std::unique(result.begin(), result.end()), result.end());
  return result;
}
 
inline Vars sort_vars_as_spin(const Vars &vars) {
  auto result = vars;
  sort_vars_in_place(result);
  if (result.size() >= 2) {
    for (auto it = result.begin();
         it != result.end() && it + 1 != result.end();) {
      if (*it == *(it + 1)) {
        it = result.erase(it, it + 2);
      } else {
        ++it;
      }
    }
  }
  return result;
}
 
 
inline Expr simplify_seq(const Expr &expr,
                         Vars (*sort_vars_func)(const Vars &) = sort_vars) {
  VarsCoeffMap vars_val_map;
  Expr result = expr.get_constant();
  for (const auto &term : expr.get_terms()) {
    auto new_vars = (*sort_vars_func)(term.get_vars());
    if (new_vars.size() == 0) {
      result += term.get_coeff();
    } else {
      vars_val_map[new_vars] += term.get_coeff();
    }
  }
  for (const auto &[vars, val] : vars_val_map) {
    if (val != 0) {
      result += Term(vars, val);
    }
  }
  std::sort(result.get_terms().begin(), result.get_terms().end());
  return result;
}
 
inline Expr simplify(const Expr &expr,
                     Vars (*sort_vars_func)(const Vars &) = sort_vars) {
  if (expr.size() <= 1000) return simplify_seq(expr, sort_vars_func);
 
  Terms simplified_terms;
  simplified_terms.resize(expr.get_terms().size());
 
  tbb::parallel_for(size_t(0), expr.get_terms().size(), [&](size_t i) {
    simplified_terms[i] = Term(sort_vars_func(expr.get_terms()[i].get_vars()),
                               expr.get_terms()[i].get_coeff());
  });
 
  tbb::parallel_sort(simplified_terms.begin(), simplified_terms.end());
 
  const size_t size = simplified_terms.size();
  constexpr size_t min_chunk_size = 32;
  size_t max_chunk_count = std::thread::hardware_concurrency();
  size_t chunk_count =
      std::min((size + min_chunk_size - 1) / min_chunk_size, max_chunk_count);
  size_t chunk_size = (size + chunk_count - 1) / chunk_count;
 
  std::vector<Terms> chunk(chunk_count);
  std::vector<energy_t> constant(chunk_count, 0);
  std::vector<size_t> indices(chunk_count);
 
  tbb::parallel_for(
      tbb::blocked_range<size_t>(0, chunk_count),
      [&](const tbb::blocked_range<size_t> &range) {
        for (size_t i = range.begin(); i < range.end(); ++i) {
          size_t start = i * chunk_size;
          size_t end = std::min(start + chunk_size, size);
 
          for (size_t j = start; j < end; ++j) {
            if (simplified_terms[j].get_coeff() == 0) {
              continue;
            }
            if (simplified_terms[j].var_count() == 0) {
              constant[i] += simplified_terms[j].get_coeff();
              continue;
            }
            if (chunk[i].empty()) {
              chunk[i].push_back(std::move(simplified_terms[j]));
              continue;
            }
            if (chunk[i].back().get_vars() == simplified_terms[j].get_vars()) {
              chunk[i].back().set_coeff(chunk[i].back().get_coeff() +
                                        simplified_terms[j].get_coeff());
              if (chunk[i].back().get_coeff() == 0) {
                chunk[i].pop_back();
              }
              continue;
            }
            chunk[i].emplace_back(std::move(simplified_terms[j]));
          }
        }
      });
 
  for (size_t i = 1; i < chunk_count; ++i) {
    if (!chunk[i - 1].empty()) {
      if (chunk[i].empty()) {
        chunk[i].push_back(std::move(chunk[i - 1].back()));
        chunk[i - 1].pop_back();
        continue;
      }
      if (chunk[i - 1].back().get_vars() == chunk[i].front().get_vars()) {
        chunk[i].front().set_coeff(chunk[i].front().get_coeff() +
                                   chunk[i - 1].back().get_coeff());
        chunk[i - 1].pop_back();
        if (chunk[i].front().get_coeff() == 0) {
          chunk[i].erase(chunk[i].begin());
        }
      }
    }
  }
 
  std::vector<size_t> prefix(chunk_count + 1, 0);
  energy_t result_constant = expr.get_constant();
  for (size_t i = 0; i < chunk_count; ++i) {
    prefix[i + 1] = prefix[i] + chunk[i].size();
    result_constant += constant[i];
  }
 
  std::vector<Term> result_terms;
  result_terms.resize(prefix.back());
  tbb::parallel_for(size_t(0), size_t(chunk_count), [&](size_t i) {
    std::copy(chunk[i].begin(), chunk[i].end(),
              result_terms.begin() + static_cast<std::ptrdiff_t>(prefix[i]));
  });
 
  return Expr(result_constant, std::move(result_terms));
}
 
template <typename T>
Vector<T> simplify(const Vector<T> &vec,
                   Vars (*sort_vars_func)(const Vars &) = sort_vars) {
  Vector<T> result(vec.size());
  tbb::parallel_for(size_t(0), vec.size(), [&](size_t i) {
    result[i] = simplify(vec[i], sort_vars_func);
  });
  return result;
}
 
inline Expr simplify_as_binary(const Expr &expr) {
  return simplify(expr, sort_vars_as_binary);
}
 
template <typename T>
auto simplify_as_binary(const Vector<T> &arg) {
  return simplify(arg, sort_vars_as_binary);
}
 
inline Expr simplify_as_spin(const Expr &expr) {
  return simplify(expr, sort_vars_as_spin);
}
 
template <typename T>
auto simplify_as_spin(const Vector<T> &arg) {
  return simplify(arg, sort_vars_as_spin);
}
 
inline bool is_simplified(const Expr &expr) {
  for (auto it = expr.get_terms().begin(); it != expr.get_terms().end(); ++it) {
    if (it->get_coeff() == 0) {
      return false;
    }
    for (auto it_var = it->get_vars().begin(); it_var != it->get_vars().end();
         ++it_var) {
      if (it_var != it->get_vars().begin() && *it_var < *(std::prev(it_var))) {
        return false;
      }
    }
    if (it != expr.get_terms().begin() && *it < *(std::prev(it))) {
      return false;
    }
  }
  return true;
}
 
inline bool is_binary(const Expr &expr) {
  for (const auto &term : expr.get_terms()) {
    if (term.var_count() > 2) {
      return false;
    }
    if (term.var_count() == 2 && term.get_vars()[0] == term.get_vars()[1]) {
      return false;
    }
  }
  return true;
}
 
 
 
inline energy_t Sol::comp_energy() const {
  const auto &terms = quad_model_.get_expr().get_terms();
  energy_t constant = quad_model_.get_constant();
 
  energy_t term_sum = tbb::parallel_reduce(
      tbb::blocked_range<size_t>(0, terms.size()), static_cast<energy_t>(0),
      [&](const tbb::blocked_range<size_t> &range, energy_t acc) {
        for (size_t i = range.begin(); i < range.end(); ++i) {
          const auto &term = terms[i];
          energy_t term_energy = term.get_coeff();
          for (auto item : term.get_vars()) {
            if (get(item) == 0) {
              term_energy = 0;
              break;
            }
          }
          acc += term_energy;
        }
        return acc;
      },
      std::plus<energy_t>());
 
  return term_sum + constant;
}
 
inline MapList Sol::get_map_list() const {
  MapList result;
  for (Var var : quad_model_.get_index_var()) {
    result.push_back({var, static_cast<var_val_t>(get(var))});
  }
  return result;
}
 
 
inline Expr reduce_sum(const Term &term) {
  coeff_t num_var = static_cast<coeff_t>(term.var_count());
  if (num_var <= 2) {  
    return term;
  }
  Expr sum_vars =
      std::accumulate(term.get_vars().begin(), term.get_vars().end(), Expr(0));
  if (term.get_coeff() < 0) {  
    return term.get_coeff() * qbpp::var() * (sum_vars - (num_var - 1));
  }
  VarInt aux_var =
      new_var_int(impl::VarSet::new_unnamed_var_str(), 0, num_var - 2, 2);
  Expr result =
      term.get_coeff() * (sum_vars - aux_var) * (sum_vars - (aux_var + 1));
  return simplify_as_binary(result) / 2;
}
 
inline Expr reduce_cascade(const Term &term) {
  coeff_t var_count = static_cast<coeff_t>(term.var_count());
  if (var_count <= 2) {  
    return term;
  }
  Term temp = term;
  auto aux_var = qbpp::var();
  auto tail1 = temp.get_vars().back();
  temp.get_vars().pop_back();
  auto tail0 = temp.get_vars().back();
  temp.get_vars().pop_back();
  temp.get_vars().push_back(aux_var);
  return reduce_cascade(temp) +
         abs(term.get_coeff()) *
             (tail0 * tail1 - 2 * (tail0 + tail1) * aux_var + 3 * aux_var);
}
 
 
inline Expr reduce(const Expr &expr) {
  Expr result = expr.get_constant();
  for (const auto &term : expr.get_terms()) {
    if (term.get_coeff() < 0 || term.var_count() <= 4) {
      result += reduce_sum(term);
    } else
      result += reduce_cascade(term);
  }
  return result;
}
 
 
 
 
inline Expr binary_to_spin(const Expr &expr) {
  Expr result = 4 * expr.get_constant();
  for (const auto &term : expr.get_terms()) {
    if (term.var_count() == 0) {
      result += term.get_coeff();
    } else if (term.var_count() == 1) {
      result += 2 * term.get_coeff() * (term.get_vars()[0] + 1);
    } else if (term.var_count() == 2) {
      result += term.get_coeff() * (term.get_vars()[0] + 1) *
                (term.get_vars()[1] + 1);
    } else {
      throw std::runtime_error(
          THROW_MESSAGE("The term must be linear or quadratic, but degree is ",
                        term.var_count(), "."));
    }
  }
  return result.simplify_as_spin();
}
 
inline Expr spin_to_binary(const Expr &expr) {
  Expr result = expr.get_constant();
  for (const auto &term : expr.get_terms()) {
    if (term.var_count() == 0) {
      result += term.get_coeff();
    } else if (term.var_count() == 1) {
      result += term.get_coeff() * (2 * term.get_vars()[0] - 1);
    } else if (term.var_count() == 2) {
      result += term.get_coeff() * (2 * term.get_vars()[0] - 1) *
                (2 * term.get_vars()[1] - 1);
    } else {
      throw std::runtime_error(
          THROW_MESSAGE("The term must be linear or quadratic, but degree is ",
                        term.var_count(), "."));
    }
  }
  return result.simplify_as_binary();
}
 
 
 
 
 
template <typename T, typename U>
Vector<Expr> operator+(const Vector<T> &lhs, const Vector<U> &rhs) {
  Vector<Expr> result = lhs;
  result += rhs;
  return result;
}
 
template <typename T, typename U>
auto operator+(const Vector<Vector<T>> &lhs, const Vector<Vector<U>> &rhs) {
  Vector<decltype(lhs[0] + rhs[0])> result = lhs;
  result += rhs;
  return result;
}
 
template <typename T>
Vector<Expr> operator+(const Vector<T> &lhs, const Expr &rhs) {
  Vector<Expr> result = lhs;
  result += rhs;
  return result;
}
 
template <typename T>
Vector<Expr> operator+(Vector<T> &&lhs, const Expr &rhs) {
  lhs += rhs;
  return std::move(lhs);
}
 
template <typename T>
auto operator+(const Vector<Vector<T>> &lhs, const Expr &rhs) {
  Vector<decltype(lhs[0] + rhs)> result = lhs;
  result += rhs;
  return result;
}
 
template <typename T>
auto operator+(Vector<Vector<T>> &&lhs, const Expr &rhs) {
  lhs += rhs;
  return std::move(lhs);
}
 
template <typename T>
auto operator+(const Expr &lhs, const Vector<T> &rhs) {
  return rhs + lhs;
}
 
template <typename T>
auto operator+(const Expr &lhs, Vector<T> &&rhs) {
  rhs += lhs;
  return std::move(rhs);
}
 
 
 
template <typename T, typename U>
Vector<Expr> operator*(const Vector<T> &lhs, const Vector<U> &rhs) {
  Vector<Expr> result = lhs;
  result *= rhs;
  return result;
}
 
template <typename T, typename U>
auto operator*(const Vector<Vector<T>> &lhs, const Vector<Vector<U>> &rhs) {
  Vector<decltype(lhs[0] * rhs[0])> result = lhs;
  result *= rhs;
  return result;
}
 
template <typename T>
Vector<Expr> operator*(const Vector<T> &lhs, const Expr &rhs) {
  Vector<Expr> result = lhs;
  result *= rhs;
  return result;
}
 
template <typename T>
Vector<Expr> operator*(Vector<T> &&lhs, const Expr &rhs) {
  lhs *= rhs;
  return std::move(lhs);
}
 
template <typename T>
auto operator*(const Vector<Vector<T>> &lhs, const Expr &rhs) {
  Vector<decltype(lhs[0] * rhs)> result = lhs;
  result *= rhs;
  return result;
}
 
template <typename T>
auto operator*(Vector<Vector<T>> &&lhs, const Expr &rhs) {
  lhs *= rhs;
  return lhs;
}
 
template <typename T>
auto operator*(const Expr &lhs, const Vector<T> &rhs) {
  return rhs * lhs;
}
 
template <typename T>
auto operator*(const Expr &lhs, Vector<T> &&rhs) {
  return std::move(rhs) * lhs;
}
 
 
 
template <typename T, typename U>
Vector<Expr> operator-(const Vector<T> &lhs, const Vector<U> &rhs) {
  Vector<Expr> result = lhs;
  result -= rhs;
  return result;
}
 
template <typename T, typename U>
auto operator-(const Vector<Vector<T>> &lhs, const Vector<Vector<U>> &rhs) {
  Vector<decltype(lhs[0] - rhs[0])> result = lhs;
  result -= rhs;
  return result;
}
 
template <typename T>
Vector<Expr> operator-(const Vector<T> &lhs, const Expr &rhs) {
  Vector<Expr> result = lhs;
  result -= rhs;
  return result;
}
 
template <typename T>
Vector<Expr> operator-(Vector<T> &&lhs, const Expr &rhs) {
  lhs -= rhs;
  return std::move(lhs);
}
 
template <typename T>
auto operator-(const Vector<Vector<T>> &lhs, const Expr &rhs) {
  Vector<decltype(lhs[0] - rhs)> result = lhs;
  result -= rhs;
  return result;
}
 
template <typename T>
auto operator-(Vector<Vector<T>> &&lhs, const Expr &rhs) {
  lhs -= rhs;
  return std::move(lhs);
}
 
template <typename T>
auto operator-(const Expr &lhs, const Vector<T> &rhs) {
  return -(rhs - lhs);
}
 
 
 
template <typename T>
auto operator/(const Vector<T> &lhs, energy_t rhs) {
  Vector<decltype(lhs[0] / rhs)> result = lhs;
  result /= rhs;
  return result;
}
 
 
 
template <typename T>
auto operator+(const Vector<T> &lhs) {
  return qbpp::toExpr(lhs);
}
 
template <typename T>
auto operator-(const Vector<T> &lhs) {
  return lhs * (-1);
}
 
 
 
template <typename T>
auto operator==(const Vector<T> &lhs, energy_t rhs) {
  Vector<decltype(lhs[0] == rhs)> result = qbpp::sqr(lhs - rhs);
  return result;
}
 
template <typename T>
auto operator<=(energy_t lhs, const Vector<T> &rhs) {
  return std::make_pair(lhs, rhs);
}
 
template <typename T>
auto operator<=(Inf lhs, const Vector<T> &rhs) {
  return std::make_pair(lhs, rhs);
}
 
template <typename T>
auto operator<=(const std::pair<energy_t, Vector<T>> &lhs, energy_t rhs) {
  Vector<Expr> result;
  result.resize(lhs.second.size());
  tbb::parallel_for(size_t(0), lhs.second.size(), [&](size_t i) {
    result[i] = (lhs.first <= lhs.second[i] <= rhs);
  });
  return result;
}
 
template <typename T>
auto operator<=(const std::pair<energy_t, Vector<T>> &lhs, Inf rhs) {
  Vector<Expr> result;
  result.resize(lhs.second.size());
  tbb::parallel_for(size_t(0), lhs.second.size(), [&](size_t i) {
    result[i] = (lhs.first <= lhs.second[i] <= rhs);
  });
  return result;
}
 
template <typename T>
auto operator<=(const std::pair<Inf, Vector<T>> &lhs, energy_t rhs) {
  Vector<Expr> result;
  result.resize(lhs.second.size());
  tbb::parallel_for(size_t(0), lhs.second.size(), [&](size_t i) {
    result[i] = (lhs.first <= lhs.second[i] <= rhs);
  });
  return result;
}
 
template <typename T>
auto operator<=(energy_t lhs, const Vector<Vector<T>> &rhs) {
  return std::make_pair(lhs, rhs);
}
 
template <typename T>
auto operator<=(const std::pair<energy_t, Vector<Vector<T>>> &lhs,
                energy_t rhs) {
  Vector<decltype(lhs.first <= lhs.second[0] <= rhs)> result;
  result.resize(lhs.second.size());
  tbb::parallel_for(size_t(0), lhs.second.size(), [&](size_t i) {
    result[i] = (lhs.first <= lhs.second[i] <= rhs);
  });
  return result;
}
 
template <typename T>
auto operator<=(const std::pair<energy_t, Vector<Vector<T>>> &lhs, Inf rhs) {
  Vector<Vector<T>> result;
  result.resize(lhs.second.size());
  tbb::parallel_for(size_t(0), lhs.second.size(), [&](size_t i) {
    result[i] = (lhs.first <= lhs.second[i] <= rhs);
  });
  return result;
}
 
template <typename T>
auto operator<=(const std::pair<Inf, Vector<Vector<T>>> &lhs, energy_t rhs) {
  Vector<decltype(lhs.second[0] <= rhs)> result;
  result.resize(lhs.second.size());
  tbb::parallel_for(size_t(0), lhs.second.size(), [&](size_t i) {
    result[i] = (lhs.first <= lhs.second[i] <= rhs);
  });
  return result;
}
 
 
 
template <typename T>
Expr sum(const T &arg [[maybe_unused]]) {
  throw std::runtime_error(
      "qbpp::sum cannot have a scalar argment. It must be a vector.");
  return toExpr(0);
}
 
template <typename T>
Expr sum(const Vector<Vector<T>> &arg [[maybe_unused]]) {
  throw std::runtime_error(
      "qbpp::sum cannot have a 2-d or higher array as an argument. It "
      "must be a 1-d vector. For 2-d or higher arrays, use qbpp::total_sum "
      "or "
      "qbpp::vector_sum.");
  return toExpr(0);
}
 
template <typename T>
Expr sum(const Vector<T> &items) {
  std::vector<Expr> result;
  result.reserve(items.size());
  for (const auto &item : items) {
    result.push_back(toExpr(item));
  }
  std::vector<size_t> indices(result.size() + 1);
  energy_t constant = 0;
  for (size_t i = 0; i < result.size(); ++i) {
    indices[i + 1] = indices[i] + result[i].term_count();
    constant += result[i].get_constant();
  }
  Terms terms;
  terms.resize(indices.back());
  tbb::parallel_for(size_t(0), result.size(), [&](size_t i) {
    std::copy(result[i].get_terms().begin(), result[i].get_terms().end(),
              terms.begin() + static_cast<std::ptrdiff_t>(indices[i]));
  });
  return Expr(constant, std::move(terms));
}
 
template <typename T>
Expr total_sum_impl(const T &item) {
  return toExpr(item);
}
 
template <typename T>
Expr total_sum_impl(const Vector<T> &items) {
  std::vector<Expr> result;
  result.resize(items.size());
  tbb::parallel_for(size_t(0), items.size(),
                    [&](size_t i) { result[i] = total_sum_impl(items[i]); });
  std::vector<size_t> indices(result.size() + 1);
  energy_t constant = 0;
 
  for (size_t i = 0; i < result.size(); ++i) {
    indices[i + 1] = indices[i] + result[i].term_count();
    constant += result[i].get_constant();
  }
  Terms terms;
  terms.resize(indices.back());
  tbb::parallel_for(
      tbb::blocked_range<size_t>(0, result.size()),
      [&](const tbb::blocked_range<size_t> &range) {
        for (size_t i = range.begin(); i < range.end(); ++i) {
          std::copy(result[i].get_terms().begin(), result[i].get_terms().end(),
                    terms.begin() + static_cast<std::ptrdiff_t>(indices[i]));
        }
      });
  return Expr(constant, std::move(terms));
}
 
template <typename T>
qbpp::Expr total_sum(const T &arg [[maybe_unused]]) {
  throw std::runtime_error(
      "qbpp::total_sum cannot have a scalar argment. It must be a vector.");
  return toExpr(0);
}
 
template <typename T>
qbpp::Expr total_sum(const Vector<T> &arg [[maybe_unused]]) {
  throw std::runtime_error(
      "qbpp::total_sum cannot have a 1-d vector as an argument. For "
      "1-d vector, use qbpp::sum.");
  return toExpr(0);
}
 
template <typename T>
Expr total_sum(const Vector<Vector<T>> &items) {
  return total_sum_impl(items);
}
 
template <typename T>
Expr vector_sum(const T &items [[maybe_unused]]) {
  throw std::runtime_error(
      "qbpp::vector_sum function does not support scalar values.");
  return toExpr(0);
}
 
template <typename T>
Expr vector_sum(const Vector<T> &items [[maybe_unused]]) {
  throw std::runtime_error(
      "qbpp::vector_sum function does not support 1-d vectors. Use "
      "qbpp::sum instead.");
  return toExpr(0);
}
 
template <typename T>
auto vector_sum_impl(const T &items) {
  return toExpr(items);
}
 
template <typename T>
auto vector_sum_impl(const Vector<T> &items) {
  return sum(items);
}
 
template <typename T>
auto vector_sum_impl(const Vector<Vector<T>> &items) {
  using SubResultType = decltype(vector_sum_impl(std::declval<Vector<T>>()));
  Vector<SubResultType> result(items.size());
  tbb::parallel_for(size_t(0), items.size(),
                    [&](size_t i) { result[i] = vector_sum_impl(items[i]); });
  return result;
}
 
template <typename T>
auto vector_sum(const Vector<Vector<T>> &items) {
  return vector_sum_impl(items);
}
 
template <typename T>
Expr product(const T &arg [[maybe_unused]]) {
  throw std::runtime_error(
      "qbpp::product cannot have a scalar argment. It must be a vector.");
  return toExpr(0);
}
 
template <typename T>
Expr product(const Vector<Vector<T>> &arg [[maybe_unused]]) {
  throw std::runtime_error(
      "qbpp::product cannot have a 2-d or higher array as an argument. It "
      "must be a 1-d vector. For 2-d or higher arrays, use "
      "qbpp::total_product "
      "or "
      "qbpp::vector_product.");
  return toExpr(0);
}
 
template <typename T>
Expr product(const Vector<T> &items) {
  return tbb::parallel_reduce(
      tbb::blocked_range<size_t>(0, items.size()), Expr{1},
      [&](const tbb::blocked_range<size_t> &range, Expr acc) {
        for (size_t i = range.begin(); i < range.end(); ++i) {
          acc *= items[i];
        }
        return acc;
      },
      [](const Expr &a, const Expr &b) { return a * b; });
}
 
template <typename T>
Expr total_product_impl(const T &item) {
  return toExpr(item);
}
 
template <typename T>
Expr total_product_impl(const Vector<T> &items) {
  Expr result = tbb::parallel_reduce(
      tbb::blocked_range<size_t>(0, items.size()), Expr{1},
      [&](const tbb::blocked_range<size_t> &range, Expr acc) {
        for (size_t i = range.begin(); i < range.end(); ++i) {
          acc *= total_product_impl(items[i]);
        }
        return acc;
      },
      [](const Expr &a, const Expr &b) { return a * b; });
  return result;
}
 
template <typename T>
Expr total_product(const Vector<Vector<T>> &items) {
  return total_product_impl(items);
}
 
template <typename T>
Expr total_product(const T &arg [[maybe_unused]]) {
  throw std::runtime_error(
      "qbpp::total_product cannot have a scalar argment. It must be a "
      "vector.");
  return toExpr(0);
}
 
template <typename T>
Expr total_product(const Vector<T> &arg [[maybe_unused]]) {
  throw std::runtime_error(
      "qbpp::total_product cannot have a 1-d vector as an argument. For "
      "1-d vector, use qbpp::product.");
  return toExpr(0);
}
 
template <typename T>
Expr vector_product(const T &items [[maybe_unused]]) {
  throw std::runtime_error(
      "qbpp::vector_product function does not support scalar values.");
  return toExpr(0);
}
 
template <typename T>
Expr vector_product(const Vector<T> &items [[maybe_unused]]) {
  throw std::runtime_error(
      "qbpp::vector_product function does not support 1-d vectors. Use "
      "qbpp::product instead.");
  return toExpr(0);
}
 
template <typename T>
auto vector_product_impl(const T &items) {
  return toExpr(items);
}
 
template <typename T>
auto vector_product_impl(const Vector<T> &items) {
  return product(items);
}
 
template <typename T>
auto vector_product_impl(const Vector<Vector<T>> &items) {
  using SubResultType = decltype(vector_product(std::declval<Vector<T>>()));
  Vector<SubResultType> result(items.size());
  tbb::parallel_for(tbb::blocked_range<size_t>(0, items.size()),
                    [&](const tbb::blocked_range<size_t> &range) {
                      for (size_t i = range.begin(); i < range.end(); ++i) {
                        result[i] = vector_product_impl(items[i]);
                      }
                    });
 
  return result;
}
 
template <typename T>
auto vector_product(const Vector<Vector<T>> &items) {
  return vector_product_impl(items);
}
 
 
 
inline energy_t toInt(const Expr &expr) {
  if (expr.term_count() == 0) {
    return static_cast<energy_t>(expr.get_constant());
  } else
    throw std::runtime_error(
        THROW_MESSAGE("The expression is not a constant."));
}
 
template <typename T>
auto toInt(const Vector<T> &arg) {
  using ResultType = decltype(toInt(arg[0]));
  Vector<ResultType> result(arg.size());
 
  tbb::parallel_for(tbb::blocked_range<size_t>(0, arg.size()),
                    [&](const tbb::blocked_range<size_t> &range) {
                      for (size_t i = range.begin(); i < range.end(); ++i) {
                        result[i] = toInt(arg[i]);
                      }
                    });
 
  return result;
}
 
template <typename T>
Vector<energy_t> eval_var_val_map(const Vector<T> &arg,
                                  const VarValMap &var_val_map) {
  Vector<energy_t> result(arg.size());
 
  tbb::parallel_for(tbb::blocked_range<size_t>(0, arg.size()),
                    [&](const tbb::blocked_range<size_t> &range) {
                      for (size_t i = range.begin(); i < range.end(); ++i) {
                        result[i] = eval_var_val_map(arg[i], var_val_map);
                      }
                    });
 
  return result;
}
 
template <typename T>
Vector<energy_t> eval(const Vector<T> &arg, const MapList &map_list) {
  VarValMap var_val_map = list_to_var_val(map_list);
  return eval_var_val_map(arg, var_val_map);
}
 
template <typename T>
auto eval_var_val_map(const Vector<Vector<T>> &arg,
                      const VarValMap &var_val_map) {
  using ResultType = decltype(eval_var_val_map(arg[0], var_val_map));
  Vector<ResultType> result(arg.size());
 
  tbb::parallel_for(tbb::blocked_range<size_t>(0, arg.size()),
                    [&](const tbb::blocked_range<size_t> &range) {
                      for (size_t i = range.begin(); i < range.end(); ++i) {
                        result[i] = eval_var_val_map(arg[i], var_val_map);
                      }
                    });
 
  return result;
}
 
template <typename T>
auto eval(const Vector<Vector<T>> &arg, const MapList &map_list) {
  VarValMap var_val_map = list_to_var_val(map_list);
  return eval_var_val_map(arg, var_val_map);
}
 
template <typename T>
Vector<Expr> replace(const Vector<T> &arg, const MapList &map_list) {
  Vector<Expr> result(arg.size());
 
  tbb::parallel_for(tbb::blocked_range<size_t>(0, arg.size()),
                    [&](const tbb::blocked_range<size_t> &range) {
                      for (size_t i = range.begin(); i < range.end(); ++i) {
                        result[i] = replace(arg[i], map_list);
                      }
                    });
 
  return result;
}
 
template <typename T>
auto replace(const Vector<Vector<T>> &arg, const MapList &map_list) {
  using ResultType = decltype(replace(arg[0], map_list));
  Vector<ResultType> result(arg.size());
 
  tbb::parallel_for(tbb::blocked_range<size_t>(0, arg.size()),
                    [&](const tbb::blocked_range<size_t> &range) {
                      for (size_t i = range.begin(); i < range.end(); ++i) {
                        result[i] = replace(arg[i], map_list);
                      }
                    });
 
  return result;
}
 
template <typename T, typename F>
auto element_wise(const Vector<T> &arg, F func) {
  Vector<decltype(func(arg[0]))> result(arg.size());
 
  tbb::parallel_for(tbb::blocked_range<size_t>(0, arg.size()),
                    [&](const tbb::blocked_range<size_t> &range) {
                      for (size_t i = range.begin(); i < range.end(); ++i) {
                        result[i] = func(arg[i]);
                      }
                    });
 
  return result;
}
 
template <typename T>
Vector<Expr> sqr(const Vector<T> &arg) {
  return element_wise(arg, [](const T &item) { return sqr(item); });
}
 
template <typename T>
auto sqr(const Vector<Vector<T>> &arg) -> Vector<decltype(sqr(arg[0]))> {
  return element_wise(arg, [](const Vector<T> &item) { return sqr(item); });
}
 
template <typename T>
Vector<Expr> reduce(const Vector<T> &arg) {
  return element_wise(arg, [](const T &item) { return reduce(item); });
}
 
template <typename T>
auto reduce(const Vector<Vector<T>> &arg) {
  return element_wise(arg, [](const Vector<T> &item) { return reduce(item); });
}
 
template <typename T>
Vector<Expr> binary_to_spin(const Vector<T> &arg) {
  return element_wise(arg, [](const T &item) { return binary_to_spin(item); });
}
 
template <typename T>
auto binary_to_spin(const Vector<Vector<T>> &arg)
    -> Vector<decltype(binary_to_spin(arg[0]))> {
  return element_wise(
      arg, [](const Vector<T> &item) { return binary_to_spin(item); });
}
 
template <typename T>
Vector<Expr> spin_to_binary(const Vector<T> &arg) {
  return element_wise(arg, [](const T &item) { return spin_to_binary(item); });
}
 
template <typename T>
auto spin_to_binary(const Vector<Vector<T>> &arg)
    -> Vector<decltype(spin_to_binary(arg[0]))> {
  return element_wise(
      arg, [](const Vector<T> &item) { return spin_to_binary(item); });
}
 
inline Vector<Expr> get_row(const Vector<Vector<Expr>> &vec, vindex_t index) {
  Vector<Expr> result = vec[index];
  return result;
}
 
inline Vector<Expr> get_col(const Vector<Vector<Expr>> &vec, size_t index) {
  Vector<Expr> result(vec.size());
 
  tbb::parallel_for(size_t(0), vec.size(),
                    [&](size_t i) { result[i] = vec[i][index]; });
 
  return result;
}
 
template <typename T>
energy_t gcd(const Vector<T> &vec) {
  return tbb::parallel_reduce(
      tbb::blocked_range<size_t>(0, vec.size()), static_cast<energy_t>(0),
      [&](const tbb::blocked_range<size_t> &range, energy_t acc) {
        for (size_t i = range.begin(); i < range.end(); ++i) {
          acc = gcd(acc, gcd(vec[i]));
        }
        return acc;
      },
      [](energy_t a, energy_t b) { return gcd(a, b); });
}
 
template <typename T>
Vector<Vector<Expr>> transpose(const Vector<Vector<T>> &vec) {
  Vector<Vector<Expr>> result(vec[0].size());
 
  tbb::parallel_for(size_t(0), vec[0].size(),
                    [&](size_t j) { result[j] = get_col(vec, j); });
 
  return result;
}
 
inline int32_t onehot_to_int(const Vector<var_val_t> &vec) {
  int32_t total = 0;
  int32_t result = -1;
  for (vindex_t i = 0; i < vec.size(); i++) {
    if (vec[i] == 1) {
      result = static_cast<int32_t>(i);
    }
    total += vec[i];
  }
  if (total != 0) {
    return result;
  } else {
    return -1;
  }
}
 
template <typename T>
auto onehot_to_int(const Vector<Vector<T>> &vec) {
  using ResultType = decltype(onehot_to_int(vec[0]));
  Vector<ResultType> result(vec.size());
 
  tbb::parallel_for(size_t(0), vec.size(),
                    [&](size_t i) { result[i] = onehot_to_int(vec[i]); });
 
  return result;
}
 
 
 
template <typename T>
Vector<T> &Vector<T>::simplify(Vars (*sort_vars_func)(const Vars &)) {
  *this = qbpp::simplify(*this, sort_vars_func);
  return *this;
}
 
template <typename T>
Vector<T> &Vector<T>::replace(const MapList &map_list) {
  *this = qbpp::replace(*this, map_list);
  return *this;
}
 
template <typename T>
Vector<T> &Vector<T>::reduce() {
  *this = qbpp::reduce(*this);
  return *this;
}
 
 
inline std::tuple<bool, size_t, size_t, coeff_t, coeff_t>
check_if_simplified_as_binary(const Expr &expr) {
  if (expr.term_count() == 0) {
    return {false, 0, 0, 0, 0};
  }
  const auto &terms = expr.get_terms();
  bool error = false;
  size_t linear_count = 0;
  coeff_t min_coeff = terms[0].get_coeff();
  coeff_t max_coeff = terms[0].get_coeff();
  std::mutex mtx;
 
  tbb::parallel_for(
      tbb::blocked_range<size_t>(0, terms.size()),
      [&](const tbb::blocked_range<size_t> &range) {
        coeff_t local_min = min_coeff;
        coeff_t local_max = max_coeff;
        for (size_t i = range.begin(); i < range.end(); ++i) {
          const auto &term = terms[i];
 
          if (local_min > term.get_coeff()) local_min = term.get_coeff();
          if (local_max < term.get_coeff()) local_max = term.get_coeff();
 
          if (term.var_count() == 0 || term.var_count() >= 3) {
            error = true;
            return;
          }
          if (term.var_count() == 2 &&
              term.get_vars()[1] <= term.get_vars()[0]) {
            error = true;
            return;
          }
 
          if (i > 0) {
            if (term <= terms[i - 1]) {
              error = true;
              return;
            }
            if (terms[i - 1].var_count() == 1 && term.var_count() == 2) {
              linear_count = i;
            }
          }
        }
 
        std::lock_guard<std::mutex> lock(mtx);
        if (local_min < min_coeff) min_coeff = local_min;
        if (local_max > max_coeff) max_coeff = local_max;
      });
 
  if (terms.back().var_count() == 1) {
    linear_count = terms.size();
  }
  return {error, linear_count, terms.size() - linear_count, min_coeff,
          max_coeff};
}
 
inline QuadModel::Impl::Impl(const qbpp::Model &model)
    : linear_(model.var_count(), 0),
      degree_(model.var_count(), 0),
      quadratic_(model.var_count()) {
  const auto var_count = model.var_count();
  const Expr &expr = model.get_expr();
  const Terms &terms = expr.get_terms();
  auto [error, linear_count, quadratic_count, min_coeff, max_coeff] =
      check_if_simplified_as_binary(expr);
  if (error) {
    throw std::runtime_error(
        THROW_MESSAGE("The expression must be simplified as binary or spin."));
  }
  min_coeff_ = min_coeff;
  max_coeff_ = max_coeff;
 
  if (model.get_expr().get_terms().size() == 0) {
    return;
  }
 
  if (linear_count > 0) {
    tbb::parallel_for(size_t(0), linear_count, [&](size_t i) {
      linear_[model.get_index(terms[i].get_vars()[0])] = terms[i].get_coeff();
    });
  }
 
  if (quadratic_count == 0) {
    return;
  }
 
  std::vector<Term> quadratic(quadratic_count * 2);
 
  auto linear_count_local = linear_count;  
  tbb::parallel_for(size_t(0), quadratic_count, [&](size_t i) {
    quadratic[2 * i] = terms[i + linear_count_local];
    quadratic[2 * i + 1].set_coeff(terms[i + linear_count_local].get_coeff());
    quadratic[2 * i + 1].get_vars().push_back(
        terms[i + linear_count_local].get_vars()[1]);
    quadratic[2 * i + 1].get_vars().push_back(
        terms[i + linear_count_local].get_vars()[0]);
  });
 
  tbb::parallel_sort(quadratic.begin(), quadratic.end());
 
  std::vector<size_t> quadratic_index(var_count + 1, 0);
  tbb::parallel_for(
      tbb::blocked_range<size_t>(1, quadratic_count * 2),
      [&](const tbb::blocked_range<size_t> &range) {
        for (size_t i = range.begin(); i < range.end(); ++i) {
          if (quadratic[i - 1].get_vars()[0] != quadratic[i].get_vars()[0]) {
            size_t start = model.get_index(quadratic[i - 1].get_vars()[0]) + 1;
            size_t end = model.get_index(quadratic[i].get_vars()[0]) + 1;
 
            for (size_t j = start; j < end; ++j) {
              quadratic_index[j] = i;
            }
          }
        }
      });
 
  quadratic_index[model.get_index(
                      quadratic[quadratic_count * 2 - 1].get_vars()[0]) +
                  1] = quadratic_count * 2;
 
  tbb::parallel_for(static_cast<vindex_t>(0), var_count, [&](size_t i) {
    degree_[i] =
        static_cast<vindex_t>(quadratic_index[i + 1] - quadratic_index[i]);
  });
 
  tbb::parallel_for(tbb::blocked_range<size_t>(0, var_count),
                    [&](const tbb::blocked_range<size_t> &range) {
                      for (size_t i = range.begin(); i < range.end(); ++i) {
                        quadratic_[i].reserve(degree_[i]);
                        for (size_t j = quadratic_index[i];
                             j < quadratic_index[i + 1]; ++j) {
                          quadratic_[i].push_back(
                              {model.get_index(quadratic[j].get_vars()[1]),
                               quadratic[j].get_coeff()});
                        }
                      }
                    });
}
 
 
inline double get_time() {
  static const auto start_ = std::chrono::high_resolution_clock::now();
  auto elapsed = std::chrono::high_resolution_clock::now() - start_;
  return std::chrono::duration<double>(elapsed).count();
}
 
}  
#endif