factorization.cpp

Go to the documentation of this file.

 
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/miller_rabin.hpp>
#include <boost/program_options.hpp>
#include <chrono>
#include <iostream>
#include <string>
#include <thread>
 
#include "qbpp.hpp"
#include "qbpp_abs2.hpp"
#include "qbpp_easy_solver.hpp"
#include "qbpp_grb.hpp"
#include "qbpp_misc.hpp"
#include "qbpp_multiplier.hpp"
 
namespace qbpp {
namespace factorization {
 
std::string toBinaryString(uint64_t num) {
  if (num == 0) return "0";
  std::string binary = "";
  while (num > 0) {
    binary = (num % 2 == 0 ? "0" : "1") + binary;
    num /= 2;
  }
  return binary;
}
 
uint64_t as_uint64(const qbpp::Sol &sol, const qbpp::Vector<qbpp::Var> &var,
                   bool fix_lsb) {
  uint64_t result = 0;
  for (int i = var.size() - 1; i >= 1; --i) {
    result = (result << 1) + sol.get(var[i]);
  }
  if (fix_lsb) {
    result = (result << 1) + 1;
  } else {
    result = (result << 1) + sol.get(var[0]);
  }
  return result;
}
 
uint64_t as_uint64(const qbpp::SolHolder &sol_holder,
                   const qbpp::Vector<qbpp::Var> &var, bool fix_lsb) {
  return as_uint64(sol_holder.copy_sol(), var, fix_lsb);
}
 
bool is_prime(uint32_t num) {
  return boost::multiprecision::miller_rabin_test(
      boost::multiprecision::uint128_t(num), 25);
}
 
uint32_t rand_prime(int size) {
  uint32_t num;
  do {
    num = qbpp::misc::RandomGenerator::gen(1 << (size - 2)) |
          (1 << (size - 1)) | 1;
  } while (is_prime(num) == false);
  return num;
}
 
 
class SolHolder : public qbpp::SolHolder {
  const bool factorization_mode = true;
  const bool fix_lsb = false;
  const qbpp::Vector<qbpp::Var> &var_x;
  const qbpp::Vector<qbpp::Var> &var_y;
  const qbpp::Vector<qbpp::Var> &var_z;
  std::optional<int64_t> bound = std::nullopt;
  mutable std::mutex mtx;
 
 public:
  SolHolder(const qbpp::QuadModel &quad_model, bool fix_lsb,
            const qbpp::Vector<qbpp::Var> &var_x,
            const qbpp::Vector<qbpp::Var> &var_y)
      : qbpp::SolHolder(quad_model),
        fix_lsb(fix_lsb),
        var_x(var_x),
        var_y(var_y),
        var_z(std::move(qbpp::Vector<qbpp::Var>{})) {}
 
  SolHolder(const qbpp::QuadModel &quad_model, bool fix_lsb,
            const qbpp::Vector<qbpp::Var> &var_z)
      : qbpp::SolHolder(quad_model),
        factorization_mode(false),
        fix_lsb(fix_lsb),
        var_x(std::move(qbpp::Vector<qbpp::Var>{})),
        var_y(std::move(qbpp::Vector<qbpp::Var>{})),
        var_z(var_z) {}
 
  SolHolder(const SolHolder &sol_holder) = default;
 
  void set_bound(int64_t bound) { this->bound = bound; }
 
  std::optional<int64_t> get_bound() const { return bound; }
 
  std::string get_bound_str() const {
    if (bound.has_value()) {
      return std::to_string(bound.value());
    } else {
      return "N/A";
    }
  }
 
  std::optional<double> set_if_better(const qbpp::Sol &new_sol,
                                      const std::string &solver) override {
    std::lock_guard<std::mutex> lock(mtx);
    std::optional<double> tts = qbpp::SolHolder::set_if_better(new_sol, solver);
    if (tts.has_value()) {
      std::cerr << std::left << std::setw(4) << get_solver() << " "
                << std::fixed << std::setprecision(3) << tts.value()
                << "s Energy = " << new_sol.get_energy()
                << " Bound = " << get_bound_str() << " : ";
      if (factorization_mode) {
        std::cerr << toBinaryString(as_uint64(new_sol, var_x, fix_lsb)) << " * "
                  << toBinaryString(as_uint64(new_sol, var_y, fix_lsb)) << " = "
                  << toBinaryString(as_uint64(new_sol, var_x, fix_lsb) *
                                    as_uint64(new_sol, var_y, fix_lsb))
                  << " (" << as_uint64(new_sol, var_x, fix_lsb) << " * "
                  << as_uint64(new_sol, var_y, fix_lsb) << " = "
                  << as_uint64(new_sol, var_x, fix_lsb) *
                         as_uint64(new_sol, var_y, fix_lsb)
                  << ")" << std::endl;
      } else {
        std::cout << toBinaryString(as_uint64(new_sol, var_z, fix_lsb)) << " ("
                  << as_uint64(new_sol, var_z, fix_lsb) << ")" << std::endl;
      }
    }
    return tts;
  }
};
 
class ABS2Callback : public qbpp_abs2::Callback {
  //  SolHolder &sol_holder;
  const std::shared_ptr<SolHolder> sol_holder_ptr_;
 
 public:
  ABS2Callback(const qbpp_abs2::QuadModel &quad_model,
               std::shared_ptr<SolHolder> sol_holder_ptr)
      : qbpp_abs2::Callback(quad_model), sol_holder_ptr_(sol_holder_ptr) {}
 
  void callback(const std::string &event) override {
    if (event == "init") {
      // Enable the callback for every new best solution.
      qbpp_abs2::Callback::set("new");
    } else if (event == "new") {
      auto sol = get_sol();
      sol_holder_ptr_->set_if_better(sol, "ABS2");
      auto hint = sol_holder_ptr_->get_if_better(sol);
      if (hint.has_value()) set_hint(hint.value());
    }
  }
};  // namespace factorization
 
class GRB_Callback : public qbpp_grb::Callback {
  //  SolHolder &sol_holder;
  std::shared_ptr<SolHolder> sol_holder_ptr_;
 
 public:
  GRB_Callback(qbpp_grb::QuadModel &quad_model,
               std::shared_ptr<SolHolder> sol_holder_ptr)
      : qbpp_grb::Callback(quad_model), sol_holder_ptr_(sol_holder_ptr) {}
 
  void callback() override {
    if (where == GRB_CB_MIPSOL) {
      qbpp_grb::Sol sol = get_sol();
      sol_holder_ptr_->set_if_better(sol, "GRB");
    }
    if (where == GRB_CB_MIP) {
      sol_holder_ptr_->set_bound(
          std::round(getDoubleInfoPublic(GRB_CB_MIP_OBJBND)));
      auto hint = sol_holder_ptr_->get_if_better(
          std::round(getDoubleInfoPublic(GRB_CB_MIP_OBJBST)));
      if (hint.has_value()) {
        for (qbpp::vindex_t i = 0; i < quad_model.var_count(); ++i) {
          setSolution(quad_model.get_grb_var(i), hint.value().get(i));
        }
        useSolution();
      }
    }
    if (where == GRB_CB_MIPNODE) {
      sol_holder_ptr_->set_bound(
          std::round(getDoubleInfoPublic(GRB_CB_MIPNODE_OBJBND)));
      auto hint = sol_holder_ptr_->get_if_better(
          std::round(getDoubleInfoPublic(GRB_CB_MIPNODE_OBJBST)));
      if (hint.has_value()) {
        for (qbpp::vindex_t i = 0; i < quad_model.var_count(); ++i) {
          setSolution(quad_model.get_grb_var(i), hint.value().get(i));
        }
        useSolution();
      }
    }
    abort_if_target_energy(sol_holder_ptr_->get_energy());
  }
};
 
}  // namespace factorization
}  // namespace qbpp
 
using namespace qbpp::factorization;
namespace po = boost::program_options;
 
int main(int argc, char *argv[]) {
  // clang-format off
  po::options_description desc(
      "Factorization/Multiplication z = x * y using QUBO solvers");
  desc.add_options()("help,h", "produce help message")
     ("multiplication,m", "solve multiplication problem (default: factorization)")
     ("x,x", po::value<uint32_t>(),"set value of x")
     ("y,y", po::value<uint32_t>(), "set value of y")
     ("random,r", po::value<uint32_t>(), "set bit size of random prime numbers for x and y (default: 32)")
     ("time,t", po::value<uint32_t>(), "set time limit in seconds")
     ("fix,f", "fix the LSBs of x, y, and z")
     ("seed,s", po::value<uint32_t>(), "set random seed")
     ("abs2,a", "use ABS2 GPU QUBO solver")
     ("gurobi,g", "use Gurobi optimizer")
     ("easy,e", "use QUBO++ Easy solver");
  // clang-format on
 
  po::variables_map vm;
  try {
    po::store(po::parse_command_line(argc, argv, desc), vm);
  } catch (const std::exception &e) {
    std::cout << "Wrong arguments. Please use -h/--help option to see the "
                 "usage.\n";
    return 1;
  }
  po::notify(vm);
 
  if (vm.count("help")) {
    std::cout << desc << std::endl;
    return 0;
  }
 
  // Set the random seed for deterministic behavior if seed is provided.
  if (vm.count("seed")) {
    qbpp::misc::RandomGenerator::set_seed(vm["seed"].as<uint32_t>());
  }
 
  bool fix_lsb = vm.count("fix");
  bool use_abs2 = vm.count("abs2");
  bool use_grb = vm.count("gurobi");
  bool use_easy = vm.count("easy");
  if (!use_abs2 && !use_grb && !use_easy) {
    use_abs2 = true;
    use_grb = true;
    use_easy = true;
  }
 
  // Set the time limit to 60 seconds by default.
  int time_limit;
  if (vm.count("time")) {
    time_limit = vm["time"].as<uint32_t>();
  } else {
    time_limit = 60;
  }
 
  // Set the problem type to factorization or multiplication.
  bool factorization_mode = true;
  if (vm.count("multiplication")) factorization_mode = false;
 
  // The input values X and Y are randomly generated by default.
  uint64_t X = 0, Y = 0;
  if (vm.count("x") && vm.count("y")) {
    X = vm["x"].as<uint32_t>();
    Y = vm["y"].as<uint32_t>();
    if (is_prime(X) == false) {
      std::cerr << "x must be a prime number." << std::endl;
      return 1;
    }
    if (is_prime(Y) == false) {
      std::cerr << "y must be a prime number." << std::endl;
      return 1;
    }
  } else {
    int size = 32;
    if (vm.count("random")) size = vm["random"].as<uint32_t>();
    do {
      X = rand_prime(size);
      Y = rand_prime(size);
    } while (X == Y);
  }
 
  // The product Z of X and Y is calculated.
  uint64_t Z = X * Y;
  // The number of bits of X, Y, and Z are calculated.
  uint32_t Xsize = 32 - __builtin_clz(X);
  uint32_t Ysize = 32 - __builtin_clz(Y);
  uint32_t Zsize = Xsize + Ysize;
 
  // Create variables x, y, and z.
  auto var_x = qbpp::var("x", Xsize);
  auto var_y = qbpp::var("y", Ysize);
  auto var_z = qbpp::var("z", Zsize);
 
  // Stores the input values of z
  //  qbpp::VarIntMap x_var_map, y_var_map, z_var_map;
  //  qbpp::VarIntMap expected_var_map;
  qbpp::MapList x_var_map, y_var_map, z_var_map;
  qbpp::MapList opt_sol;
 
  for (uint32_t i = 0; i < Xsize; ++i) {
    x_var_map.push_back({var_x[i], (X >> i) & 1});
  }
  for (uint32_t i = 0; i < Ysize; ++i) {
    y_var_map.push_back({var_y[i], (Y >> i) & 1});
  }
 
  for (uint32_t i = 0; i < Zsize; ++i) {
    z_var_map.push_back({var_z[i], (Z >> i) & 1});
  }
 
  opt_sol.insert(opt_sol.end(), x_var_map.begin(), x_var_map.end());
  opt_sol.insert(opt_sol.end(), y_var_map.begin(), y_var_map.end());
  opt_sol.insert(opt_sol.end(), z_var_map.begin(), z_var_map.end());
 
  // Returns QUBO quad_model for multiplier of var_z = var_x * var_y
  qbpp::Expr expr = multiplier(var_x, var_y, var_z, opt_sol);
 
  // Assignment by assign performed and then simplified as a binary
  // expression.
  if (factorization_mode) {
    expr.replace(z_var_map);
  } else {
    expr.replace(x_var_map);
    expr.replace(y_var_map);
  }
 
  if (fix_lsb) {
    expr.replace({{var_x[0], 1}, {var_y[0], 1}, {var_z[0], 1}});
  }
 
  expr.simplify_as_binary();
 
  // Create a QUBO quad_model from the expression.
  auto quad_model = qbpp::QuadModel(expr);
 
  // Create a holder to store the best solution by ABS2 and Gurobi.
  std::shared_ptr<SolHolder> sol_holder_ptr;
  if (factorization_mode) {
    sol_holder_ptr =
        std::make_shared<SolHolder>(quad_model, fix_lsb, var_x, var_y);
  } else {
    sol_holder_ptr = std::make_shared<SolHolder>(quad_model, fix_lsb, var_z);
  }
  std::cout << "QUBO : Variables = " << quad_model.var_count()
            << " Linear terms = " << quad_model.term_count(1)
            << " Quad terms = " << quad_model.term_count(2) << std::endl;
  std::cout << "Target : " << toBinaryString(X) << " * " << toBinaryString(Y)
            << " = " << toBinaryString(Z) << " (" << X << " * " << Y << " = "
            << Z << ")\n";
 
  std::thread easy_thread;
  std::thread abs2_thread;
  std::thread grb_thread;
 
  std::shared_ptr<qbpp::easy_solver::EasySolver> easy_solver_ptr;
 
  if (use_easy) {
    easy_solver_ptr = std::make_shared<qbpp::easy_solver::EasySolver>(
        quad_model, sol_holder_ptr);
    easy_solver_ptr->set_time_limit(time_limit);
    easy_solver_ptr->set_target_energy(0);
    easy_thread = std::thread([&easy_solver_ptr, &sol_holder_ptr]() {
      auto sol = easy_solver_ptr->search();
      sol_holder_ptr->set_if_better(sol, "EASY");
    });
  }
 
  std::shared_ptr<qbpp_abs2::Solver> abs2_solver_ptr;
  std::shared_ptr<qbpp_abs2::QuadModel> abs2_quad_model_ptr;
  std::shared_ptr<qbpp_abs2::Param> abs2_param_ptr;
  std::shared_ptr<ABS2Callback> abs2_callback_ptr;
 
  if (use_abs2) {
    // Initialize ABS2 solver
    abs2_solver_ptr = std::make_shared<qbpp_abs2::Solver>();
    // Create a quad_model for ABS2 solver from
    // quad_model
    abs2_quad_model_ptr = std::make_shared<qbpp_abs2::QuadModel>(quad_model);
    // Create a callback function for ABS2 solver.
    abs2_callback_ptr =
        std::make_shared<ABS2Callback>(*abs2_quad_model_ptr, sol_holder_ptr);
    // Create a parameters object to store solver
    // parameters.
    abs2_param_ptr = std::make_shared<qbpp_abs2::Param>();
    // Set a time limit
    abs2_param_ptr->set_time_limit(time_limit);
    // Set the target energy to 0.
    abs2_param_ptr->set_target_energy(0);
    // Set the callback function
    abs2_param_ptr->set(*abs2_callback_ptr);
    abs2_thread = std::thread([&abs2_solver_ptr, &abs2_quad_model_ptr,
                               &abs2_param_ptr, &sol_holder_ptr]() {
      auto sol = (*abs2_solver_ptr)(*abs2_quad_model_ptr, *abs2_param_ptr);
      sol_holder_ptr->set_if_better(sol, "ABS2");
    });
  }
 
  std::shared_ptr<qbpp_grb::QuadModel> grb_quad_model_ptr;
  std::shared_ptr<GRB_Callback> grb_callback_ptr;
  if (use_grb) {
    // Create a Gurobi quad_model from quad_model.
    grb_quad_model_ptr = std::make_shared<qbpp_grb::QuadModel>(quad_model);
    // Set the time limit for the Gurobi solver.
    grb_quad_model_ptr->set_time_limit(time_limit);
    // Create a callback function for Gurobi solver.
    grb_callback_ptr =
        std::make_shared<GRB_Callback>(*grb_quad_model_ptr, sol_holder_ptr);
    // Set the target energy to 0.
    grb_callback_ptr->set_target_energy(0);
    // Set the callback function
    grb_quad_model_ptr->set(*grb_callback_ptr);
    grb_thread = std::thread([&grb_quad_model_ptr, &sol_holder_ptr]() {
      auto sol = grb_quad_model_ptr->optimize();
      sol_holder_ptr->set_if_better(sol, "GRB");
    });
  }
 
  if (use_easy) easy_thread.join();
  if (use_abs2) abs2_thread.join();
  if (use_grb) grb_thread.join();
 
  std::cout << "Target : " << toBinaryString(X) << " * " << toBinaryString(Y)
            << " = " << toBinaryString(Z) << " (" << X << " * " << Y << " = "
            << Z << ")\n";
}