結果
| 問題 | No.2464 To DAG |
| ユーザー |
 hitonanode
|
| 提出日時 | 2023-09-08 22:06:32 |
| 言語 | C++23 (gcc 12.3.0 + boost 1.83.0) |
| 結果 |
AC
|
| 実行時間 | 296 ms / 2,000 ms |
| コード長 | 59,343 bytes |
| コンパイル時間 | 2,164 ms |
| コンパイル使用メモリ | 195,060 KB |
| 実行使用メモリ | 52,544 KB |
| 最終ジャッジ日時 | 2024-06-26 15:09:53 |
| 合計ジャッジ時間 | 16,413 ms |
|
ジャッジサーバーID （参考情報） |
judge3 / judge2 |
（要ログイン）
テストケース
テストケース表示| 入力 | 結果 | 実行時間 実行使用メモリ |
|---|---|---|
| testcase_00 | AC | 2 ms
5,248 KB |
| testcase_01 | AC | 2 ms
5,248 KB |
| testcase_02 | AC | 296 ms
52,544 KB |
| testcase_03 | AC | 272 ms
47,204 KB |
| testcase_04 | AC | 251 ms
42,652 KB |
| testcase_05 | AC | 248 ms
41,784 KB |
| testcase_06 | AC | 256 ms
41,360 KB |
| testcase_07 | AC | 164 ms
29,604 KB |
| testcase_08 | AC | 191 ms
30,432 KB |
| testcase_09 | AC | 173 ms
26,860 KB |
| testcase_10 | AC | 172 ms
26,012 KB |
| testcase_11 | AC | 112 ms
19,664 KB |
| testcase_12 | AC | 91 ms
16,784 KB |
| testcase_13 | AC | 114 ms
19,024 KB |
| testcase_14 | AC | 171 ms
30,768 KB |
| testcase_15 | AC | 248 ms
39,812 KB |
| testcase_16 | AC | 105 ms
23,720 KB |
| testcase_17 | AC | 15 ms
6,528 KB |
| testcase_18 | AC | 59 ms
16,740 KB |
| testcase_19 | AC | 2 ms
5,376 KB |
| testcase_20 | AC | 2 ms
5,376 KB |
| testcase_21 | AC | 2 ms
5,376 KB |
| testcase_22 | AC | 86 ms
17,688 KB |
| testcase_23 | AC | 83 ms
17,520 KB |
| testcase_24 | AC | 86 ms
19,584 KB |
| testcase_25 | AC | 86 ms
19,596 KB |
| testcase_26 | AC | 40 ms
11,704 KB |
| testcase_27 | AC | 295 ms
48,576 KB |
| testcase_28 | AC | 283 ms
48,444 KB |
| testcase_29 | AC | 285 ms
47,068 KB |
| testcase_30 | AC | 180 ms
29,748 KB |
| testcase_31 | AC | 160 ms
26,472 KB |
| testcase_32 | AC | 146 ms
25,360 KB |
| testcase_33 | AC | 253 ms
42,608 KB |
| testcase_34 | AC | 112 ms
23,440 KB |
| testcase_35 | AC | 71 ms
18,144 KB |
| testcase_36 | AC | 73 ms
17,996 KB |
| testcase_37 | AC | 64 ms
17,380 KB |
| testcase_38 | AC | 66 ms
17,384 KB |
| testcase_39 | AC | 19 ms
22,016 KB |
| testcase_40 | AC | 20 ms
22,016 KB |
ソースコード
#include <algorithm>
#include <array>
#include <bitset>
#include <cassert>
#include <chrono>
#include <cmath>
#include <complex>
#include <deque>
#include <forward_list>
#include <fstream>
#include <functional>
#include <iomanip>
#include <ios>
#include <iostream>
#include <limits>
#include <list>
#include <map>
#include <memory>
#include <numeric>
#include <optional>
#include <queue>
#include <random>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <tuple>
#include <type_traits>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
using namespace std;
using lint = long long;
using pint = pair<int, int>;
using plint = pair<lint, lint>;
struct fast_ios { fast_ios(){ cin.tie(nullptr), ios::sync_with_stdio(false), cout << fixed << setprecision(20); }; } fast_ios_;
#define ALL(x) (x).begin(), (x).end()
#define FOR(i, begin, end) for(int i=(begin),i##_end_=(end);i<i##_end_;i++)
#define IFOR(i, begin, end) for(int i=(end)-1,i##_begin_=(begin);i>=i##_begin_;i--)
#define REP(i, n) FOR(i,0,n)
#define IREP(i, n) IFOR(i,0,n)
template <typename T> bool chmax(T &m, const T q) { return m < q ? (m = q, true) : false; }
template <typename T> bool chmin(T &m, const T q) { return m > q ? (m = q, true) : false; }
const std::vector<std::pair<int, int>> grid_dxs{{1, 0}, {-1, 0}, {0, 1}, {0, -1}};
int floor_lg(long long x) { return x <= 0 ? -1 : 63 - __builtin_clzll(x); }
template <class T1, class T2> T1 floor_div(T1 num, T2 den) { return (num > 0 ? num / den : -((-num + den - 1) / den)); }
template <class T1, class T2> std::pair<T1, T2> operator+(const std::pair<T1, T2> &l, const std::pair<T1, T2> &r) { return std::make_pair(l.first + r.first, l.second + r.second); }
template <class T1, class T2> std::pair<T1, T2> operator-(const std::pair<T1, T2> &l, const std::pair<T1, T2> &r) { return std::make_pair(l.first - r.first, l.second - r.second); }
template <class T> std::vector<T> sort_unique(std::vector<T> vec) { sort(vec.begin(), vec.end()), vec.erase(unique(vec.begin(), vec.end()), vec.end()); return vec; }
template <class T> int arglb(const std::vector<T> &v, const T &x) { return std::distance(v.begin(), std::lower_bound(v.begin(), v.end(), x)); }
template <class T> int argub(const std::vector<T> &v, const T &x) { return std::distance(v.begin(), std::upper_bound(v.begin(), v.end(), x)); }
template <class IStream, class T> IStream &operator>>(IStream &is, std::vector<T> &vec) { for (auto &v : vec) is >> v; return is; }
template <class OStream, class T> OStream &operator<<(OStream &os, const std::vector<T> &vec);
template <class OStream, class T, size_t sz> OStream &operator<<(OStream &os, const std::array<T, sz> &arr);
template <class OStream, class T, class TH> OStream &operator<<(OStream &os, const std::unordered_set<T, TH> &vec);
template <class OStream, class T, class U> OStream &operator<<(OStream &os, const pair<T, U> &pa);
template <class OStream, class T> OStream &operator<<(OStream &os, const std::deque<T> &vec);
template <class OStream, class T> OStream &operator<<(OStream &os, const std::set<T> &vec);
template <class OStream, class T> OStream &operator<<(OStream &os, const std::multiset<T> &vec);
template <class OStream, class T> OStream &operator<<(OStream &os, const std::unordered_multiset<T> &vec);
template <class OStream, class T, class U> OStream &operator<<(OStream &os, const std::pair<T, U> &pa);
template <class OStream, class TK, class TV> OStream &operator<<(OStream &os, const std::map<TK, TV> &mp);
template <class OStream, class TK, class TV, class TH> OStream &operator<<(OStream &os, const std::unordered_map<TK, TV, TH> &mp);
template <class OStream, class... T> OStream &operator<<(OStream &os, const std::tuple<T...> &tpl);
template <class OStream, class T> OStream &operator<<(OStream &os, const std::vector<T> &vec) { os << '['; for (auto v : vec) os << v << ','; os << ']'; return os; }
template <class OStream, class T, size_t sz> OStream &operator<<(OStream &os, const std::array<T, sz> &arr) { os << '['; for (auto v : arr) os << v << ','; os << ']'; return os; }
template <class... T> std::istream &operator>>(std::istream &is, std::tuple<T...> &tpl) { std::apply([&is](auto &&... args) { ((is >> args), ...);}, tpl); return is; }
template <class OStream, class... T> OStream &operator<<(OStream &os, const std::tuple<T...> &tpl) { os << '('; std::apply([&os](auto &&... args) { ((os << args << ','), ...);}, tpl); return os << ')'; }
template <class OStream, class T, class TH> OStream &operator<<(OStream &os, const std::unordered_set<T, TH> &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; }
template <class OStream, class T> OStream &operator<<(OStream &os, const std::deque<T> &vec) { os << "deq["; for (auto v : vec) os << v << ','; os << ']'; return os; }
template <class OStream, class T> OStream &operator<<(OStream &os, const std::set<T> &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; }
template <class OStream, class T> OStream &operator<<(OStream &os, const std::multiset<T> &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; }
template <class OStream, class T> OStream &operator<<(OStream &os, const std::unordered_multiset<T> &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; }
template <class OStream, class T, class U> OStream &operator<<(OStream &os, const std::pair<T, U> &pa) { return os << '(' << pa.first << ',' << pa.second << ')'; }
template <class OStream, class TK, class TV> OStream &operator<<(OStream &os, const std::map<TK, TV> &mp) { os << '{'; for (auto v : mp) os << v.first << "=>" << v.second << ','; os << '}'; return os; }
template <class OStream, class TK, class TV, class TH> OStream &operator<<(OStream &os, const std::unordered_map<TK, TV, TH> &mp) { os << '{'; for (auto v : mp) os << v.first << "=>" << v.second << ','; os << '}'; return os; }
#ifdef HITONANODE_LOCAL
const string COLOR_RESET = "\033[0m", BRIGHT_GREEN = "\033[1;32m", BRIGHT_RED = "\033[1;31m", BRIGHT_CYAN = "\033[1;36m", NORMAL_CROSSED = "\033[0;9;37m", RED_BACKGROUND = "\033[1;41m", NORMAL_FAINT = "\033[0;2m";
#define dbg(x) std::cerr << BRIGHT_CYAN << #x << COLOR_RESET << " = " << (x) << NORMAL_FAINT << " (L" << __LINE__ << ") " << __FILE__ << COLOR_RESET << std::endl
#define dbgif(cond, x) ((cond) ? std::cerr << BRIGHT_CYAN << #x << COLOR_RESET << " = " << (x) << NORMAL_FAINT << " (L" << __LINE__ << ") " << __FILE__ << COLOR_RESET << std::endl : std::cerr)
#else
#define dbg(x) ((void)0)
#define dbgif(cond, x) ((void)0)
#endif
#include <algorithm>
#include <cassert>
#include <chrono>
#include <cmath>
#include <limits>
#include <random>
#include <utility>
#include <vector>
// This program is the modificatiosn of the
// [lemon::NetworkSimplex](http://lemon.cs.elte.hu/pub/doc/latest-svn/a00404.html)
//
/* -*- mode: C++; indent-tabs-mode: nil; -*-
*
* This file is a part of LEMON, a generic C++ optimization library.
*
* Copyright (C) 2003-2013
* Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
* (Egervary Research Group on Combinatorial Optimization, EGRES).
*
* Permission to use, modify and distribute this software is granted
* provided that this copyright notice appears in all copies. For
* precise terms see the accompanying LICENSE file.
*
* This software is provided "AS IS" with no warranty of any kind,
* express or implied, and with no claim as to its suitability for any
* purpose.
*
*/
template <class Digraph, typename V = int, typename C = V> class NetworkSimplex {
public:
using Node = int;
using Arc = int;
static const int INVALID = -1;
typedef V Value; /// The type of the flow amounts, capacity bounds and supply values
typedef C Cost; /// The type of the arc costs
public:
enum ProblemType { INFEASIBLE, OPTIMAL, UNBOUNDED };
/// \brief Constants for selecting the type of the supply constraints.
///
/// Enum type containing constants for selecting the supply type,
/// i.e. the direction of the inequalities in the supply/demand
/// constraints of the \ref min_cost_flow "minimum cost flow problem".
///
/// The default supply type is \c GEQ, the \c LEQ type can be
/// selected using \ref supplyType().
/// The equality form is a special case of both supply types.
enum SupplyType {
/// This option means that there are <em>"greater or equal"</em>
/// supply/demand constraints in the definition of the problem.
GEQ,
/// This option means that there are <em>"less or equal"</em>
/// supply/demand constraints in the definition of the problem.
LEQ
};
/// \brief Constants for selecting the pivot rule.
///
/// Enum type containing constants for selecting the pivot rule for
/// the \ref run() function.
///
/// \ref NetworkSimplex provides five different implementations for
/// the pivot strategy that significantly affects the running time
/// of the algorithm.
/// According to experimental tests conducted on various problem
/// instances, \ref BLOCK_SEARCH "Block Search" and
/// \ref ALTERING_LIST "Altering Candidate List" rules turned out
/// to be the most efficient.
/// Since \ref BLOCK_SEARCH "Block Search" is a simpler strategy that
/// seemed to be slightly more robust, it is used by default.
/// However, another pivot rule can easily be selected using the
/// \ref run() function with the proper parameter.
enum PivotRule {
/// The \e First \e Eligible pivot rule.
/// The next eligible arc is selected in a wraparound fashion
/// in every iteration.
FIRST_ELIGIBLE,
/// The \e Best \e Eligible pivot rule.
/// The best eligible arc is selected in every iteration.
BEST_ELIGIBLE,
/// The \e Block \e Search pivot rule.
/// A specified number of arcs are examined in every iteration
/// in a wraparound fashion and the best eligible arc is selected
/// from this block.
BLOCK_SEARCH,
/// The \e Candidate \e List pivot rule.
/// In a major iteration a candidate list is built from eligible arcs
/// in a wraparound fashion and in the following minor iterations
/// the best eligible arc is selected from this list.
CANDIDATE_LIST,
/// The \e Altering \e Candidate \e List pivot rule.
/// It is a modified version of the Candidate List method.
/// It keeps only a few of the best eligible arcs from the former
/// candidate list and extends this list in every iteration.
ALTERING_LIST
};
private:
using IntVector = std::vector<int>;
using ValueVector = std::vector<Value>;
using CostVector = std::vector<Cost>;
using CharVector = std::vector<signed char>;
enum ArcState { STATE_UPPER = -1, STATE_TREE = 0, STATE_LOWER = 1 };
enum ArcDirection { DIR_DOWN = -1, DIR_UP = 1 };
private:
// Data related to the underlying digraph
const Digraph &_graph;
int _node_num;
int _arc_num;
int _all_arc_num;
int _search_arc_num;
// Parameters of the problem
bool _has_lower;
SupplyType _stype;
Value _sum_supply;
// Data structures for storing the digraph
IntVector _source;
IntVector _target;
// Node and arc data
ValueVector _lower;
ValueVector _upper;
ValueVector _cap;
CostVector _cost;
ValueVector _supply;
ValueVector _flow;
CostVector _pi;
// Data for storing the spanning tree structure
IntVector _parent;
IntVector _pred;
IntVector _thread;
IntVector _rev_thread;
IntVector _succ_num;
IntVector _last_succ;
CharVector _pred_dir;
CharVector _state;
IntVector _dirty_revs;
int _root;
// Temporary data used in the current pivot iteration
int in_arc, join, u_in, v_in, u_out, v_out;
Value delta;
const Value MAX;
public:
/// \brief Constant for infinite upper bounds (capacities).
///
/// Constant for infinite upper bounds (capacities).
/// It is \c std::numeric_limits<Value>::infinity() if available,
/// \c std::numeric_limits<Value>::max() otherwise.
const Value INF;
private:
// Implementation of the First Eligible pivot rule
class FirstEligiblePivotRule {
private:
// References to the NetworkSimplex class
const IntVector &_source;
const IntVector &_target;
const CostVector &_cost;
const CharVector &_state;
const CostVector &_pi;
int &_in_arc;
int _search_arc_num;
// Pivot rule data
int _next_arc;
public:
// Constructor
FirstEligiblePivotRule(NetworkSimplex &ns)
: _source(ns._source), _target(ns._target), _cost(ns._cost), _state(ns._state),
_pi(ns._pi), _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), _next_arc(0) {}
// Find next entering arc
bool findEnteringArc() {
Cost c;
for (int e = _next_arc; e != _search_arc_num; ++e) {
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
if (c < 0) {
_in_arc = e;
_next_arc = e + 1;
return true;
}
}
for (int e = 0; e != _next_arc; ++e) {
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
if (c < 0) {
_in_arc = e;
_next_arc = e + 1;
return true;
}
}
return false;
}
}; // class FirstEligiblePivotRule
// Implementation of the Best Eligible pivot rule
class BestEligiblePivotRule {
private:
// References to the NetworkSimplex class
const IntVector &_source;
const IntVector &_target;
const CostVector &_cost;
const CharVector &_state;
const CostVector &_pi;
int &_in_arc;
int _search_arc_num;
public:
// Constructor
BestEligiblePivotRule(NetworkSimplex &ns)
: _source(ns._source), _target(ns._target), _cost(ns._cost), _state(ns._state),
_pi(ns._pi), _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num) {}
// Find next entering arc
bool findEnteringArc() {
Cost c, min = 0;
for (int e = 0; e != _search_arc_num; ++e) {
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
if (c < min) {
min = c;
_in_arc = e;
}
}
return min < 0;
}
}; // class BestEligiblePivotRule
// Implementation of the Block Search pivot rule
class BlockSearchPivotRule {
private:
// References to the NetworkSimplex class
const IntVector &_source;
const IntVector &_target;
const CostVector &_cost;
const CharVector &_state;
const CostVector &_pi;
int &_in_arc;
int _search_arc_num;
// Pivot rule data
int _block_size;
int _next_arc;
public:
// Constructor
BlockSearchPivotRule(NetworkSimplex &ns)
: _source(ns._source), _target(ns._target), _cost(ns._cost), _state(ns._state),
_pi(ns._pi), _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), _next_arc(0) {
// The main parameters of the pivot rule
const double BLOCK_SIZE_FACTOR = 1.0;
const int MIN_BLOCK_SIZE = 10;
_block_size = std::max(
int(BLOCK_SIZE_FACTOR * std::sqrt(double(_search_arc_num))), MIN_BLOCK_SIZE);
}
// Find next entering arc
bool findEnteringArc() {
Cost c, min = 0;
int cnt = _block_size;
int e;
for (e = _next_arc; e != _search_arc_num; ++e) {
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
if (c < min) {
min = c;
_in_arc = e;
}
if (--cnt == 0) {
if (min < 0) goto search_end;
cnt = _block_size;
}
}
for (e = 0; e != _next_arc; ++e) {
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
if (c < min) {
min = c;
_in_arc = e;
}
if (--cnt == 0) {
if (min < 0) goto search_end;
cnt = _block_size;
}
}
if (min >= 0) return false;
search_end:
_next_arc = e;
return true;
}
}; // class BlockSearchPivotRule
// Implementation of the Candidate List pivot rule
class CandidateListPivotRule {
private:
// References to the NetworkSimplex class
const IntVector &_source;
const IntVector &_target;
const CostVector &_cost;
const CharVector &_state;
const CostVector &_pi;
int &_in_arc;
int _search_arc_num;
// Pivot rule data
IntVector _candidates;
int _list_length, _minor_limit;
int _curr_length, _minor_count;
int _next_arc;
public:
/// Constructor
CandidateListPivotRule(NetworkSimplex &ns)
: _source(ns._source), _target(ns._target), _cost(ns._cost), _state(ns._state),
_pi(ns._pi), _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), _next_arc(0) {
// The main parameters of the pivot rule
const double LIST_LENGTH_FACTOR = 0.25;
const int MIN_LIST_LENGTH = 10;
const double MINOR_LIMIT_FACTOR = 0.1;
const int MIN_MINOR_LIMIT = 3;
_list_length = std::max(
int(LIST_LENGTH_FACTOR * std::sqrt(double(_search_arc_num))), MIN_LIST_LENGTH);
_minor_limit = std::max(int(MINOR_LIMIT_FACTOR * _list_length), MIN_MINOR_LIMIT);
_curr_length = _minor_count = 0;
_candidates.resize(_list_length);
}
/// Find next entering arc
bool findEnteringArc() {
Cost min, c;
int e;
if (_curr_length > 0 && _minor_count < _minor_limit) {
// Minor iteration: select the best eligible arc from the
// current candidate list
++_minor_count;
min = 0;
for (int i = 0; i < _curr_length; ++i) {
e = _candidates[i];
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
if (c < min) {
min = c;
_in_arc = e;
} else if (c >= 0) {
_candidates[i--] = _candidates[--_curr_length];
}
}
if (min < 0) return true;
}
// Major iteration: build a new candidate list
min = 0;
_curr_length = 0;
for (e = _next_arc; e != _search_arc_num; ++e) {
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
if (c < 0) {
_candidates[_curr_length++] = e;
if (c < min) {
min = c;
_in_arc = e;
}
if (_curr_length == _list_length) goto search_end;
}
}
for (e = 0; e != _next_arc; ++e) {
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
if (c < 0) {
_candidates[_curr_length++] = e;
if (c < min) {
min = c;
_in_arc = e;
}
if (_curr_length == _list_length) goto search_end;
}
}
if (_curr_length == 0) return false;
search_end:
_minor_count = 1;
_next_arc = e;
return true;
}
}; // class CandidateListPivotRule
// Implementation of the Altering Candidate List pivot rule
class AlteringListPivotRule {
private:
// References to the NetworkSimplex class
const IntVector &_source;
const IntVector &_target;
const CostVector &_cost;
const CharVector &_state;
const CostVector &_pi;
int &_in_arc;
int _search_arc_num;
// Pivot rule data
int _block_size, _head_length, _curr_length;
int _next_arc;
IntVector _candidates;
CostVector _cand_cost;
// Functor class to compare arcs during sort of the candidate list
class SortFunc {
private:
const CostVector &_map;
public:
SortFunc(const CostVector &map) : _map(map) {}
bool operator()(int left, int right) { return _map[left] < _map[right]; }
};
SortFunc _sort_func;
public:
// Constructor
AlteringListPivotRule(NetworkSimplex &ns)
: _source(ns._source), _target(ns._target), _cost(ns._cost), _state(ns._state),
_pi(ns._pi), _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), _next_arc(0),
_cand_cost(ns._search_arc_num), _sort_func(_cand_cost) {
// The main parameters of the pivot rule
const double BLOCK_SIZE_FACTOR = 1.0;
const int MIN_BLOCK_SIZE = 10;
const double HEAD_LENGTH_FACTOR = 0.01;
const int MIN_HEAD_LENGTH = 3;
_block_size = std::max(
int(BLOCK_SIZE_FACTOR * std::sqrt(double(_search_arc_num))), MIN_BLOCK_SIZE);
_head_length = std::max(int(HEAD_LENGTH_FACTOR * _block_size), MIN_HEAD_LENGTH);
_candidates.resize(_head_length + _block_size);
_curr_length = 0;
}
// Find next entering arc
bool findEnteringArc() {
// Check the current candidate list
int e;
Cost c;
for (int i = 0; i != _curr_length; ++i) {
e = _candidates[i];
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
if (c < 0) {
_cand_cost[e] = c;
} else {
_candidates[i--] = _candidates[--_curr_length];
}
}
// Extend the list
int cnt = _block_size;
int limit = _head_length;
for (e = _next_arc; e != _search_arc_num; ++e) {
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
if (c < 0) {
_cand_cost[e] = c;
_candidates[_curr_length++] = e;
}
if (--cnt == 0) {
if (_curr_length > limit) goto search_end;
limit = 0;
cnt = _block_size;
}
}
for (e = 0; e != _next_arc; ++e) {
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
if (c < 0) {
_cand_cost[e] = c;
_candidates[_curr_length++] = e;
}
if (--cnt == 0) {
if (_curr_length > limit) goto search_end;
limit = 0;
cnt = _block_size;
}
}
if (_curr_length == 0) return false;
search_end:
// Perform partial sort operation on the candidate list
int new_length = std::min(_head_length + 1, _curr_length);
std::partial_sort(_candidates.begin(), _candidates.begin() + new_length,
_candidates.begin() + _curr_length, _sort_func);
// Select the entering arc and remove it from the list
_in_arc = _candidates[0];
_next_arc = e;
_candidates[0] = _candidates[new_length - 1];
_curr_length = new_length - 1;
return true;
}
}; // class AlteringListPivotRule
public:
NetworkSimplex(const Digraph &graph)
: _graph(graph), MAX(std::numeric_limits<Value>::max()),
INF(std::numeric_limits<Value>::has_infinity ? std::numeric_limits<Value>::infinity()
: MAX) {
// Check the number types
static_assert(-Value(1) < 0, "Value must be signed");
static_assert(-Cost(1) < 0, "Cost must be signed");
static_assert(std::numeric_limits<Value>::max() > 0, "max() must be greater than 0");
// Reset data structures
reset();
}
template <typename LowerMap> NetworkSimplex &lowerMap(const LowerMap &map) {
_has_lower = true;
for (Arc a = 0; a < _arc_num; a++) _lower[a] = map[a];
return *this;
}
template <typename UpperMap> NetworkSimplex &upperMap(const UpperMap &map) {
for (Arc a = 0; a < _arc_num; a++) _upper[a] = map[a];
return *this;
}
// Set costs of arcs (default value: 1)
template <typename CostMap> NetworkSimplex &costMap(const CostMap &map) {
for (Arc a = 0; a < _arc_num; a++) _cost[a] = map[a];
return *this;
}
template <typename SupplyMap> NetworkSimplex &supplyMap(const SupplyMap &map) {
for (Node n = 0; n < _node_num; n++) _supply[n] = map[n];
return *this;
}
NetworkSimplex &stSupply(const Node &s, const Node &t, Value k) { // set s-t flow
for (int i = 0; i != _node_num; ++i) _supply[i] = 0;
_supply[s] = k, _supply[t] = -k;
return *this;
}
/// \brief Set the type of the supply constraints.
///
/// This function sets the type of the supply/demand constraints.
/// If it is not used before calling \ref run(), the \ref GEQ supply
/// type will be used.
NetworkSimplex &supplyType(SupplyType supply_type) {
_stype = supply_type;
return *this;
}
/// @}
/// This function can be called more than once. All the given parameters
/// are kept for the next call, unless \ref resetParams() or \ref reset()
/// is used, thus only the modified parameters have to be set again.
/// If the underlying digraph was also modified after the construction
/// of the class (or the last \ref reset() call), then the \ref reset()
/// function must be called.
ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
if (!init()) return INFEASIBLE;
return start(pivot_rule);
}
/// \brief Reset all the parameters that have been given before.
///
/// This function resets all the paramaters that have been given
/// before using functions \ref lowerMap(), \ref upperMap(),
/// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
///
/// It is useful for multiple \ref run() calls. Basically, all the given
/// parameters are kept for the next \ref run() call, unless
/// \ref resetParams() or \ref reset() is used.
/// If the underlying digraph was also modified after the construction
/// of the class or the last \ref reset() call, then the \ref reset()
/// function must be used, otherwise \ref resetParams() is sufficient.
///
/// For example,
/// \code
/// NetworkSimplex<ListDigraph> ns(graph);
///
/// // First run
/// ns.lowerMap(lower).upperMap(upper).costMap(cost)
/// .supplyMap(sup).run();
///
/// // Run again with modified cost map (resetParams() is not called,
/// // so only the cost map have to be set again)
/// cost[e] += 100;
/// ns.costMap(cost).run();
///
/// // Run again from scratch using resetParams()
/// // (the lower bounds will be set to zero on all arcs)
/// ns.resetParams();
/// ns.upperMap(capacity).costMap(cost)
/// .supplyMap(sup).run();
/// \endcode
///
/// \return <tt>(*this)</tt>
///
/// \see reset(), run()
NetworkSimplex &resetParams() {
for (int i = 0; i != _node_num; ++i) { _supply[i] = 0; }
for (int i = 0; i != _arc_num; ++i) {
_lower[i] = 0;
_upper[i] = INF;
_cost[i] = 1;
}
_has_lower = false;
_stype = GEQ;
return *this;
}
/// \brief Reset the internal data structures and all the parameters
/// that have been given before.
///
/// This function resets the internal data structures and all the
/// paramaters that have been given before using functions \ref lowerMap(),
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
/// \ref supplyType().
///
/// It is useful for multiple \ref run() calls. Basically, all the given
/// parameters are kept for the next \ref run() call, unless
/// \ref resetParams() or \ref reset() is used.
/// If the underlying digraph was also modified after the construction
/// of the class or the last \ref reset() call, then the \ref reset()
/// function must be used, otherwise \ref resetParams() is sufficient.
///
/// See \ref resetParams() for examples.
///
/// \return <tt>(*this)</tt>
///
/// \see resetParams(), run()
NetworkSimplex &reset() {
// Resize vectors
_node_num = _graph.countNodes();
_arc_num = _graph.countArcs();
int all_node_num = _node_num + 1;
int max_arc_num = _arc_num + 2 * _node_num;
_source.resize(max_arc_num);
_target.resize(max_arc_num);
_lower.resize(_arc_num);
_upper.resize(_arc_num);
_cap.resize(max_arc_num);
_cost.resize(max_arc_num);
_supply.resize(all_node_num);
_flow.resize(max_arc_num);
_pi.resize(all_node_num);
_parent.resize(all_node_num);
_pred.resize(all_node_num);
_pred_dir.resize(all_node_num);
_thread.resize(all_node_num);
_rev_thread.resize(all_node_num);
_succ_num.resize(all_node_num);
_last_succ.resize(all_node_num);
_state.resize(max_arc_num);
for (int a = 0; a < _arc_num; ++a) {
_source[a] = _graph.source(a);
_target[a] = _graph.target(a);
}
// Reset parameters
resetParams();
return *this;
}
/// @}
template <typename Number = Cost> Number totalCost() const {
Number c = 0;
for (Arc a = 0; a < _arc_num; a++) c += Number(_flow[a]) * Number(_cost[a]);
return c;
}
Value flow(const Arc &a) const { return _flow[a]; }
template <typename FlowMap> void flowMap(FlowMap &map) const {
for (Arc a = 0; a < _arc_num; a++) { map.set(a, _flow[a]); }
}
ValueVector flowMap() const { return _flow; }
Cost potential(const Node &n) const { return _pi[n]; }
template <typename PotentialMap> void potentialMap(PotentialMap &map) const {
for (int n = 0; n < _graph.V; n++) { map.set(n, _pi[n]); }
}
CostVector potentialMap() const { return _pi; }
private:
// Initialize internal data structures
bool init() {
if (_node_num == 0) return false;
// Check the sum of supply values
_sum_supply = 0;
for (int i = 0; i != _node_num; ++i) { _sum_supply += _supply[i]; }
if (!((_stype == GEQ && _sum_supply <= 0) || (_stype == LEQ && _sum_supply >= 0)))
return false;
// Check lower and upper bounds
// LEMON_DEBUG(checkBoundMaps(), "Upper bounds must be greater or equal to the lower bounds");
// Remove non-zero lower bounds
if (_has_lower) {
for (int i = 0; i != _arc_num; ++i) {
Value c = _lower[i];
if (c >= 0) {
_cap[i] = _upper[i] < MAX ? _upper[i] - c : INF;
} else {
_cap[i] = _upper[i] < MAX + c ? _upper[i] - c : INF;
}
_supply[_source[i]] -= c;
_supply[_target[i]] += c;
}
} else {
for (int i = 0; i != _arc_num; ++i) { _cap[i] = _upper[i]; }
}
// Initialize artifical cost
Cost ART_COST;
if (std::numeric_limits<Cost>::is_exact) {
ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;
} else {
ART_COST = 0;
for (int i = 0; i != _arc_num; ++i) {
if (_cost[i] > ART_COST) ART_COST = _cost[i];
}
ART_COST = (ART_COST + 1) * _node_num;
}
// Initialize arc maps
for (int i = 0; i != _arc_num; ++i) {
_flow[i] = 0;
_state[i] = STATE_LOWER;
}
// Set data for the artificial root node
_root = _node_num;
_parent[_root] = -1;
_pred[_root] = -1;
_thread[_root] = 0;
_rev_thread[0] = _root;
_succ_num[_root] = _node_num + 1;
_last_succ[_root] = _root - 1;
_supply[_root] = -_sum_supply;
_pi[_root] = 0;
// Add artificial arcs and initialize the spanning tree data structure
if (_sum_supply == 0) {
// EQ supply constraints
_search_arc_num = _arc_num;
_all_arc_num = _arc_num + _node_num;
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
_parent[u] = _root;
_pred[u] = e;
_thread[u] = u + 1;
_rev_thread[u + 1] = u;
_succ_num[u] = 1;
_last_succ[u] = u;
_cap[e] = INF;
_state[e] = STATE_TREE;
if (_supply[u] >= 0) {
_pred_dir[u] = DIR_UP;
_pi[u] = 0;
_source[e] = u;
_target[e] = _root;
_flow[e] = _supply[u];
_cost[e] = 0;
} else {
_pred_dir[u] = DIR_DOWN;
_pi[u] = ART_COST;
_source[e] = _root;
_target[e] = u;
_flow[e] = -_supply[u];
_cost[e] = ART_COST;
}
}
} else if (_sum_supply > 0) {
// LEQ supply constraints
_search_arc_num = _arc_num + _node_num;
int f = _arc_num + _node_num;
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
_parent[u] = _root;
_thread[u] = u + 1;
_rev_thread[u + 1] = u;
_succ_num[u] = 1;
_last_succ[u] = u;
if (_supply[u] >= 0) {
_pred_dir[u] = DIR_UP;
_pi[u] = 0;
_pred[u] = e;
_source[e] = u;
_target[e] = _root;
_cap[e] = INF;
_flow[e] = _supply[u];
_cost[e] = 0;
_state[e] = STATE_TREE;
} else {
_pred_dir[u] = DIR_DOWN;
_pi[u] = ART_COST;
_pred[u] = f;
_source[f] = _root;
_target[f] = u;
_cap[f] = INF;
_flow[f] = -_supply[u];
_cost[f] = ART_COST;
_state[f] = STATE_TREE;
_source[e] = u;
_target[e] = _root;
_cap[e] = INF;
_flow[e] = 0;
_cost[e] = 0;
_state[e] = STATE_LOWER;
++f;
}
}
_all_arc_num = f;
} else {
// GEQ supply constraints
_search_arc_num = _arc_num + _node_num;
int f = _arc_num + _node_num;
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
_parent[u] = _root;
_thread[u] = u + 1;
_rev_thread[u + 1] = u;
_succ_num[u] = 1;
_last_succ[u] = u;
if (_supply[u] <= 0) {
_pred_dir[u] = DIR_DOWN;
_pi[u] = 0;
_pred[u] = e;
_source[e] = _root;
_target[e] = u;
_cap[e] = INF;
_flow[e] = -_supply[u];
_cost[e] = 0;
_state[e] = STATE_TREE;
} else {
_pred_dir[u] = DIR_UP;
_pi[u] = -ART_COST;
_pred[u] = f;
_source[f] = u;
_target[f] = _root;
_cap[f] = INF;
_flow[f] = _supply[u];
_state[f] = STATE_TREE;
_cost[f] = ART_COST;
_source[e] = _root;
_target[e] = u;
_cap[e] = INF;
_flow[e] = 0;
_cost[e] = 0;
_state[e] = STATE_LOWER;
++f;
}
}
_all_arc_num = f;
}
return true;
}
// Check if the upper bound is greater than or equal to the lower bound
// on each arc.
bool checkBoundMaps() {
for (int j = 0; j != _arc_num; ++j) {
if (_upper[j] < _lower[j]) return false;
}
return true;
}
// Find the join node
void findJoinNode() {
int u = _source[in_arc];
int v = _target[in_arc];
while (u != v) {
if (_succ_num[u] < _succ_num[v]) {
u = _parent[u];
} else {
v = _parent[v];
}
}
join = u;
}
// Find the leaving arc of the cycle and returns true if the
// leaving arc is not the same as the entering arc
bool findLeavingArc() {
// Initialize first and second nodes according to the direction
// of the cycle
int first, second;
if (_state[in_arc] == STATE_LOWER) {
first = _source[in_arc];
second = _target[in_arc];
} else {
first = _target[in_arc];
second = _source[in_arc];
}
delta = _cap[in_arc];
int result = 0;
Value c, d;
int e;
// Search the cycle form the first node to the join node
for (int u = first; u != join; u = _parent[u]) {
e = _pred[u];
d = _flow[e];
if (_pred_dir[u] == DIR_DOWN) {
c = _cap[e];
d = c >= MAX ? INF : c - d;
}
if (d < delta) {
delta = d;
u_out = u;
result = 1;
}
}
// Search the cycle form the second node to the join node
for (int u = second; u != join; u = _parent[u]) {
e = _pred[u];
d = _flow[e];
if (_pred_dir[u] == DIR_UP) {
c = _cap[e];
d = c >= MAX ? INF : c - d;
}
if (d <= delta) {
delta = d;
u_out = u;
result = 2;
}
}
if (result == 1) {
u_in = first;
v_in = second;
} else {
u_in = second;
v_in = first;
}
return result != 0;
}
// Change _flow and _state vectors
void changeFlow(bool change) {
// Augment along the cycle
if (delta > 0) {
Value val = _state[in_arc] * delta;
_flow[in_arc] += val;
for (int u = _source[in_arc]; u != join; u = _parent[u]) {
_flow[_pred[u]] -= _pred_dir[u] * val;
}
for (int u = _target[in_arc]; u != join; u = _parent[u]) {
_flow[_pred[u]] += _pred_dir[u] * val;
}
}
// Update the state of the entering and leaving arcs
if (change) {
_state[in_arc] = STATE_TREE;
_state[_pred[u_out]] = (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
} else {
_state[in_arc] = -_state[in_arc];
}
}
// Update the tree structure
void updateTreeStructure() {
int old_rev_thread = _rev_thread[u_out];
int old_succ_num = _succ_num[u_out];
int old_last_succ = _last_succ[u_out];
v_out = _parent[u_out];
// Check if u_in and u_out coincide
if (u_in == u_out) {
// Update _parent, _pred, _pred_dir
_parent[u_in] = v_in;
_pred[u_in] = in_arc;
_pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
// Update _thread and _rev_thread
if (_thread[v_in] != u_out) {
int after = _thread[old_last_succ];
_thread[old_rev_thread] = after;
_rev_thread[after] = old_rev_thread;
after = _thread[v_in];
_thread[v_in] = u_out;
_rev_thread[u_out] = v_in;
_thread[old_last_succ] = after;
_rev_thread[after] = old_last_succ;
}
} else {
// Handle the case when old_rev_thread equals to v_in
// (it also means that join and v_out coincide)
int thread_continue = old_rev_thread == v_in ? _thread[old_last_succ] : _thread[v_in];
// Update _thread and _parent along the stem nodes (i.e. the nodes
// between u_in and u_out, whose parent have to be changed)
int stem = u_in; // the current stem node
int par_stem = v_in; // the new parent of stem
int next_stem; // the next stem node
int last = _last_succ[u_in]; // the last successor of stem
int before, after = _thread[last];
_thread[v_in] = u_in;
_dirty_revs.clear();
_dirty_revs.push_back(v_in);
while (stem != u_out) {
// Insert the next stem node into the thread list
next_stem = _parent[stem];
_thread[last] = next_stem;
_dirty_revs.push_back(last);
// Remove the subtree of stem from the thread list
before = _rev_thread[stem];
_thread[before] = after;
_rev_thread[after] = before;
// Change the parent node and shift stem nodes
_parent[stem] = par_stem;
par_stem = stem;
stem = next_stem;
// Update last and after
last = _last_succ[stem] == _last_succ[par_stem] ? _rev_thread[par_stem]
: _last_succ[stem];
after = _thread[last];
}
_parent[u_out] = par_stem;
_thread[last] = thread_continue;
_rev_thread[thread_continue] = last;
_last_succ[u_out] = last;
// Remove the subtree of u_out from the thread list except for
// the case when old_rev_thread equals to v_in
if (old_rev_thread != v_in) {
_thread[old_rev_thread] = after;
_rev_thread[after] = old_rev_thread;
}
// Update _rev_thread using the new _thread values
for (int i = 0; i != int(_dirty_revs.size()); ++i) {
int u = _dirty_revs[i];
_rev_thread[_thread[u]] = u;
}
// Update _pred, _pred_dir, _last_succ and _succ_num for the
// stem nodes from u_out to u_in
int tmp_sc = 0, tmp_ls = _last_succ[u_out];
for (int u = u_out, p = _parent[u]; u != u_in; u = p, p = _parent[u]) {
_pred[u] = _pred[p];
_pred_dir[u] = -_pred_dir[p];
tmp_sc += _succ_num[u] - _succ_num[p];
_succ_num[u] = tmp_sc;
_last_succ[p] = tmp_ls;
}
_pred[u_in] = in_arc;
_pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
_succ_num[u_in] = old_succ_num;
}
// Update _last_succ from v_in towards the root
int up_limit_out = _last_succ[join] == v_in ? join : -1;
int last_succ_out = _last_succ[u_out];
for (int u = v_in; u != -1 && _last_succ[u] == v_in; u = _parent[u]) {
_last_succ[u] = last_succ_out;
}
// Update _last_succ from v_out towards the root
if (join != old_rev_thread && v_in != old_rev_thread) {
for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ; u = _parent[u]) {
_last_succ[u] = old_rev_thread;
}
} else if (last_succ_out != old_last_succ) {
for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ; u = _parent[u]) {
_last_succ[u] = last_succ_out;
}
}
// Update _succ_num from v_in to join
for (int u = v_in; u != join; u = _parent[u]) { _succ_num[u] += old_succ_num; }
// Update _succ_num from v_out to join
for (int u = v_out; u != join; u = _parent[u]) { _succ_num[u] -= old_succ_num; }
}
// Update potentials in the subtree that has been moved
void updatePotential() {
Cost sigma = _pi[v_in] - _pi[u_in] - _pred_dir[u_in] * _cost[in_arc];
int end = _thread[_last_succ[u_in]];
for (int u = u_in; u != end; u = _thread[u]) { _pi[u] += sigma; }
}
// Heuristic initial pivots
bool initialPivots() {
Value curr, total = 0;
std::vector<Node> supply_nodes, demand_nodes;
for (int u = 0; u < _node_num; ++u) {
curr = _supply[u];
if (curr > 0) {
total += curr;
supply_nodes.push_back(u);
} else if (curr < 0) {
demand_nodes.push_back(u);
}
}
if (_sum_supply > 0) total -= _sum_supply;
if (total <= 0) return true;
IntVector arc_vector;
if (_sum_supply >= 0) {
if (supply_nodes.size() == 1 && demand_nodes.size() == 1) {
// Perform a reverse graph search from the sink to the source
std::vector<char> reached(_node_num, false);
Node s = supply_nodes[0], t = demand_nodes[0];
std::vector<Node> stack;
reached[t] = true;
stack.push_back(t);
while (!stack.empty()) {
Node u, v = stack.back();
stack.pop_back();
if (v == s) break;
// for (InArcIt a(_graph, v); a != INVALID; ++a) {
for (auto a : _graph.in_eids[v]) {
if (reached[u = _graph.source(a)]) continue;
int j = a;
if (_cap[j] >= total) {
arc_vector.push_back(j);
reached[u] = true;
stack.push_back(u);
}
}
}
} else {
// Find the min. cost incoming arc for each demand node
for (int i = 0; i != int(demand_nodes.size()); ++i) {
Node v = demand_nodes[i];
Cost c, min_cost = std::numeric_limits<Cost>::max();
Arc min_arc = INVALID;
for (auto a : _graph.in_eids[v]) {
// for (InArcIt a(_graph, v); a != INVALID; ++a) {
c = _cost[a];
if (c < min_cost) {
min_cost = c;
min_arc = a;
}
}
if (min_arc != INVALID) { arc_vector.push_back(min_arc); }
}
}
} else {
// Find the min. cost outgoing arc for each supply node
for (Node u : supply_nodes) {
Cost c, min_cost = std::numeric_limits<Cost>::max();
Arc min_arc = INVALID;
for (auto a : _graph.out_eids[u]) {
c = _cost[a];
if (c < min_cost) {
min_cost = c;
min_arc = a;
}
}
if (min_arc != INVALID) { arc_vector.push_back(min_arc); }
}
}
// Perform heuristic initial pivots
for (int i = 0; i != int(arc_vector.size()); ++i) {
in_arc = arc_vector[i];
if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] - _pi[_target[in_arc]]) >= 0)
continue;
findJoinNode();
bool change = findLeavingArc();
if (delta >= MAX) return false;
changeFlow(change);
if (change) {
updateTreeStructure();
updatePotential();
}
}
return true;
}
// Execute the algorithm
ProblemType start(PivotRule pivot_rule) {
// Select the pivot rule implementation
switch (pivot_rule) {
case FIRST_ELIGIBLE: return start<FirstEligiblePivotRule>();
case BEST_ELIGIBLE: return start<BestEligiblePivotRule>();
case BLOCK_SEARCH: return start<BlockSearchPivotRule>();
case CANDIDATE_LIST: return start<CandidateListPivotRule>();
case ALTERING_LIST: return start<AlteringListPivotRule>();
}
return INFEASIBLE; // avoid warning
}
template <typename PivotRuleImpl> ProblemType start() {
PivotRuleImpl pivot(*this);
// Perform heuristic initial pivots
if (!initialPivots()) return UNBOUNDED;
// Execute the Network Simplex algorithm
while (pivot.findEnteringArc()) {
findJoinNode();
bool change = findLeavingArc();
if (delta >= MAX) return UNBOUNDED;
changeFlow(change);
if (change) {
updateTreeStructure();
updatePotential();
}
}
// Check feasibility
for (int e = _search_arc_num; e != _all_arc_num; ++e) {
if (_flow[e] != 0) return INFEASIBLE;
}
// Transform the solution and the supply map to the original form
if (_has_lower) {
for (int i = 0; i != _arc_num; ++i) {
Value c = _lower[i];
if (c != 0) {
_flow[i] += c;
_supply[_source[i]] += c;
_supply[_target[i]] -= c;
}
}
}
// Shift potentials to meet the requirements of the GEQ/LEQ type
// optimality conditions
if (_sum_supply == 0) {
if (_stype == GEQ) {
Cost max_pot = -std::numeric_limits<Cost>::max();
for (int i = 0; i != _node_num; ++i) {
if (_pi[i] > max_pot) max_pot = _pi[i];
}
if (max_pot > 0) {
for (int i = 0; i != _node_num; ++i) _pi[i] -= max_pot;
}
} else {
Cost min_pot = std::numeric_limits<Cost>::max();
for (int i = 0; i != _node_num; ++i) {
if (_pi[i] < min_pot) min_pot = _pi[i];
}
if (min_pot < 0) {
for (int i = 0; i != _node_num; ++i) _pi[i] -= min_pot;
}
}
}
return OPTIMAL;
}
}; // class NetworkSimplex
template <typename Capacity = long long, typename Weight = long long> struct mcf_graph_ns {
struct Digraph {
const int V;
int E;
std::vector<std::vector<int>> in_eids, out_eids;
std::vector<std::pair<int, int>> arcs;
Digraph(int V = 0) : V(V), E(0), in_eids(V), out_eids(V){};
int add_edge(int s, int t) {
assert(0 <= s and s < V);
assert(0 <= t and t < V);
in_eids[t].push_back(E), out_eids[s].push_back(E), arcs.emplace_back(s, t), E++;
return E - 1;
}
int countNodes() const noexcept { return V; }
int countArcs() const noexcept { return E; }
int source(int arcid) const { return arcs[arcid].first; }
int target(int arcid) const { return arcs[arcid].second; }
};
struct edge {
int eid;
int from, to;
Capacity lo, hi;
Weight weight;
};
int n;
std::vector<Capacity> bs;
bool infeasible;
std::vector<edge> Edges;
mcf_graph_ns(int V = 0) : n(V), bs(V), infeasible(false) {}
int add_edge(int from, int to, Capacity lower, Capacity upper, Weight weight) {
assert(from >= 0 and from < n);
assert(to >= 0 and to < n);
int idnow = Edges.size();
Edges.push_back({idnow, from, to, lower, upper, weight});
return idnow;
}
void add_supply(int v, Capacity b) {
assert(v >= 0 and v < n);
bs[v] += b;
}
std::vector<Capacity> flow;
std::vector<Capacity> potential;
template <typename RetVal = __int128> [[nodiscard]] RetVal solve() {
std::mt19937 rng(std::chrono::steady_clock::now().time_since_epoch().count());
std::vector<int> vid(n), eid(Edges.size());
std::iota(vid.begin(), vid.end(), 0);
std::shuffle(vid.begin(), vid.end(), rng);
std::iota(eid.begin(), eid.end(), 0);
std::shuffle(eid.begin(), eid.end(), rng);
flow.clear();
potential.clear();
Digraph graph(n + 1);
std::vector<Capacity> supplies(graph.countNodes());
std::vector<Capacity> lowers(Edges.size());
std::vector<Capacity> uppers(Edges.size());
std::vector<Weight> weights(Edges.size());
for (int i = 0; i < n; i++) supplies[vid[i]] = bs[i];
for (auto i : eid) {
const auto &e = Edges[i];
int arc = graph.add_edge(vid[e.from], vid[e.to]);
lowers[arc] = e.lo;
uppers[arc] = e.hi;
weights[arc] = e.weight;
}
NetworkSimplex<Digraph, Capacity, Weight> ns(graph);
auto status = ns.supplyMap(supplies).costMap(weights).lowerMap(lowers).upperMap(uppers).run(
decltype(ns)::BLOCK_SEARCH);
if (status == decltype(ns)::INFEASIBLE) {
return infeasible = true, 0;
} else {
flow.resize(Edges.size());
potential.resize(n);
for (int i = 0; i < int(Edges.size()); i++) flow[eid[i]] = ns.flow(i);
for (int i = 0; i < n; i++) potential[i] = ns.potential(vid[i]);
return ns.template totalCost<RetVal>();
}
}
};
// This implementation is based on https://ei1333.github.io/library/graph/others/eulerian-trail.cpp
template <bool directed> struct EulerianTrail {
using Vertex = int;
using Edge = int;
int M;
std::vector<int> deg;
std::vector<std::pair<Vertex, Vertex>> edges;
std::vector<std::vector<std::pair<Vertex, Edge>>> to;
std::vector<int> used_vertex, used_edge;
EulerianTrail(int V = 0) : M(0), deg(V, 0), to(V), used_vertex(V) {}
void add_edge(Vertex a, Vertex b) {
to[a].emplace_back(b, M);
if (directed) {
deg[a]++;
deg[b]--;
} else {
to[b].emplace_back(a, M);
deg[a]++;
deg[b]++;
}
edges.emplace_back(a, b);
M++;
used_edge.push_back(0);
}
std::vector<std::vector<Edge>> solve() {
// if (directed) {
// if (std::count_if(deg.begin(), deg.end(), [](int x) { return x != 0; })) return {};
// } else {
// if (std::count_if(deg.begin(), deg.end(), [](int x) { return (x & 1) != 0; }))
// return {};
// }
std::vector<std::vector<Edge>> ret;
for (int i = 0; i < int(to.size()); i++) {
if (to[i].empty() or used_vertex[i]) continue;
ret.emplace_back(go(i));
}
return ret;
}
std::vector<Edge> go(Vertex s) {
std::vector<std::pair<Vertex, Edge>> st;
std::vector<Edge> ord;
st.emplace_back(s, -1);
while (!st.empty()) {
Vertex now = st.back().first;
used_vertex[now] = true;
if (to[now].empty()) {
ord.emplace_back(st.back().second);
st.pop_back();
} else {
auto e = to[now].back();
to[now].pop_back();
if (used_edge[e.second]) continue;
used_edge[e.second] = true;
st.emplace_back(e);
}
}
ord.pop_back();
std::reverse(ord.begin(), ord.end());
return ord;
}
};
int main() {
int N, M;
cin >> N >> M;
vector<pint> edges(M);
vector<int> d(N);
EulerianTrail<true> et(N);
for (auto &[s, t] : edges) {
cin >> s >> t, --s, --t;
et.add_edge(s, t);
d.at(s)++;
d.at(t)--;
}
dbg(d);
vector<vector<pint>> to(N);
REP(e, M) to.at(edges.at(e).first).emplace_back(edges.at(e).second, e);
vector<int> is_erased_edge(M);
vector<int> isin(N);
vector<int> stv;
vector<int> ste;
auto push_v = [&](int v) {
assert(!isin.at(v));
isin.at(v) = 1;
stv.push_back(v);
};
auto pop_ev = [&] {
isin.at(stv.back()) = 0;
stv.pop_back();
ste.pop_back();
};
REP(i, N) {
push_v(i);
while (true) {
int now = stv.back();
if (to.at(now).size()) {
auto [nxt, nxteid] = to.at(now).back();
to.at(now).pop_back();
if (isin.at(nxt)) {
is_erased_edge.at(nxteid) = 1;
while (stv.back() != nxt) {
is_erased_edge.at(ste.back()) = 1;
pop_ev();
}
} else {
push_v(nxt);
ste.push_back(nxteid);
now = nxt;
}
} else if (stv.size() > 1) {
pop_ev();
} else {
break;
}
}
assert(ste.empty());
assert(stv.size() == 1);
ste.clear();
while (stv.size()) {
isin.at(stv.back()) = 0;
stv.pop_back();
}
dbg(isin);
dbg(stv);
dbg(ste);
}
vector<pint> ans;
REP(e, M) if (!is_erased_edge.at(e)) ans.push_back(edges.at(e));
cout << N << ' ' << ans.size() << '\n';
for (auto [s, t] : ans) cout << s + 1 << ' ' << t + 1 << '\n';
}