結果
| 問題 | No.1775 Love Triangle 2 |
| ユーザー |
 Pachicobue
|
| 提出日時 | 2021-12-05 05:24:30 |
| 言語 | C++17 (gcc 12.3.0 + boost 1.83.0) |
| 結果 |
WA
|
| 実行時間 | - |
| コード長 | 61,639 bytes |
| コンパイル時間 | 3,645 ms |
| コンパイル使用メモリ | 244,744 KB |
| 実行使用メモリ | 5,524 KB |
| 最終ジャッジ日時 | 2023-09-21 13:17:48 |
| 合計ジャッジ時間 | 7,840 ms |
|
ジャッジサーバーID （参考情報） |
judge12 / judge13 |
（要ログイン）
テストケース
テストケース表示| 入力 | 結果 | 実行時間 実行使用メモリ |
|---|---|---|
| testcase_00 | AC | 2 ms
4,380 KB |
| testcase_01 | AC | 1 ms
4,380 KB |
| testcase_02 | AC | 1 ms
4,376 KB |
| testcase_03 | AC | 2 ms
4,388 KB |
| testcase_04 | WA | - |
| testcase_05 | WA | - |
| testcase_06 | WA | - |
| testcase_07 | WA | - |
| testcase_08 | WA | - |
| testcase_09 | AC | 3 ms
4,380 KB |
| testcase_10 | WA | - |
| testcase_11 | WA | - |
| testcase_12 | WA | - |
| testcase_13 | WA | - |
| testcase_14 | WA | - |
| testcase_15 | AC | 3 ms
4,380 KB |
| testcase_16 | AC | 3 ms
4,380 KB |
| testcase_17 | AC | 4 ms
4,376 KB |
| testcase_18 | AC | 5 ms
4,920 KB |
| testcase_19 | AC | 3 ms
4,376 KB |
| testcase_20 | WA | - |
| testcase_21 | WA | - |
| testcase_22 | WA | - |
| testcase_23 | AC | 2 ms
4,380 KB |
| testcase_24 | WA | - |
| testcase_25 | AC | 3 ms
4,376 KB |
| testcase_26 | AC | 3 ms
4,380 KB |
| testcase_27 | AC | 3 ms
4,380 KB |
| testcase_28 | AC | 4 ms
4,380 KB |
| testcase_29 | AC | 6 ms
5,172 KB |
| testcase_30 | AC | 5 ms
4,964 KB |
| testcase_31 | WA | - |
| testcase_32 | WA | - |
| testcase_33 | WA | - |
| testcase_34 | WA | - |
| testcase_35 | WA | - |
| testcase_36 | WA | - |
| testcase_37 | WA | - |
| testcase_38 | AC | 3 ms
4,376 KB |
| testcase_39 | AC | 6 ms
5,228 KB |
| testcase_40 | WA | - |
| testcase_41 | WA | - |
| testcase_42 | WA | - |
| testcase_43 | WA | - |
| testcase_44 | WA | - |
| testcase_45 | WA | - |
| testcase_46 | WA | - |
| testcase_47 | WA | - |
| testcase_48 | WA | - |
| testcase_49 | AC | 6 ms
5,216 KB |
| testcase_50 | AC | 5 ms
5,008 KB |
| testcase_51 | WA | - |
| testcase_52 | WA | - |
| testcase_53 | WA | - |
| testcase_54 | WA | - |
| testcase_55 | WA | - |
| testcase_56 | WA | - |
| testcase_57 | WA | - |
| testcase_58 | WA | - |
| testcase_59 | AC | 5 ms
5,008 KB |
| testcase_60 | WA | - |
| testcase_61 | WA | - |
| testcase_62 | WA | - |
| testcase_63 | WA | - |
| testcase_64 | WA | - |
| testcase_65 | AC | 3 ms
4,376 KB |
| testcase_66 | AC | 3 ms
4,380 KB |
| testcase_67 | WA | - |
| testcase_68 | WA | - |
| testcase_69 | WA | - |
| testcase_70 | WA | - |
| testcase_71 | WA | - |
| testcase_72 | WA | - |
| testcase_73 | WA | - |
| testcase_74 | WA | - |
| testcase_75 | WA | - |
| testcase_76 | WA | - |
| testcase_77 | WA | - |
| testcase_78 | WA | - |
| testcase_79 | WA | - |
| testcase_80 | WA | - |
| testcase_81 | WA | - |
| testcase_82 | WA | - |
| testcase_83 | WA | - |
| testcase_84 | WA | - |
| testcase_85 | WA | - |
| testcase_86 | WA | - |
| testcase_87 | WA | - |
| testcase_88 | WA | - |
| testcase_89 | WA | - |
ソースコード
#include <bits/stdc++.h>
#pragma region Header
using i32 = int;
using u32 = unsigned int;
using i64 = long long;
using u64 = unsigned long long;
using i128 = __int128_t;
using u128 = __uint128_t;
using f64 = double;
using f80 = long double;
using f128 = __float128;
constexpr i32 operator"" _i32(u64 v)
{
return v;
}
constexpr i32 operator"" _u32(u64 v)
{
return v;
}
constexpr i64 operator"" _i64(u64 v)
{
return v;
}
constexpr u64 operator"" _u64(u64 v)
{
return v;
}
constexpr f64 operator"" _f64(f80 v)
{
return v;
}
constexpr f80 operator"" _f80(f80 v)
{
return v;
}
using Istream = std::istream;
using Ostream = std::ostream;
using Str = std::string;
template<typename T>
using Lt = std::less<T>;
template<typename T>
using Gt = std::greater<T>;
template<typename T>
using IList = std::initializer_list<T>;
template<int n>
using BSet = std::bitset<n>;
template<typename T1, typename T2>
using Pair = std::pair<T1, T2>;
template<typename... Ts>
using Tup = std::tuple<Ts...>;
template<typename T, int N>
using Arr = std::array<T, N>;
template<typename... Ts>
using Deq = std::deque<Ts...>;
template<typename... Ts>
using Set = std::set<Ts...>;
template<typename... Ts>
using MSet = std::multiset<Ts...>;
template<typename... Ts>
using USet = std::unordered_set<Ts...>;
template<typename... Ts>
using UMSet = std::unordered_multiset<Ts...>;
template<typename... Ts>
using Map = std::map<Ts...>;
template<typename... Ts>
using MMap = std::multimap<Ts...>;
template<typename... Ts>
using UMap = std::unordered_map<Ts...>;
template<typename... Ts>
using UMMap = std::unordered_multimap<Ts...>;
template<typename... Ts>
using Vec = std::vector<Ts...>;
template<typename... Ts>
using Stack = std::stack<Ts...>;
template<typename... Ts>
using Queue = std::queue<Ts...>;
template<typename T>
using MaxHeap = std::priority_queue<T>;
template<typename T>
using MinHeap = std::priority_queue<T, Vec<T>, Gt<T>>;
using NSec = std::chrono::nanoseconds;
using USec = std::chrono::microseconds;
using MSec = std::chrono::milliseconds;
using Sec = std::chrono::seconds;
template<typename T>
constexpr T LIMMIN = std::numeric_limits<T>::min();
template<typename T>
constexpr T LIMMAX = std::numeric_limits<T>::max();
template<typename T>
constexpr T INF = (LIMMAX<T> - 1) / 2;
template<typename T>
constexpr T PI = T{3.141592653589793238462643383279502884};
template<typename T = u64>
constexpr T TEN(const int n)
{
return n == 0 ? T{1} : TEN<T>(n - 1) * T{10};
}
Ostream& operator<<(Ostream& os, i128 v)
{
bool minus = false;
if (v < 0) { minus = true, v = -v; }
Str ans;
if (v == 0) { ans = "0"; }
while (v) {
ans.push_back('0' + v % 10), v /= 10;
}
std::reverse(ans.begin(), ans.end());
return os << (minus ? "-" : "") << ans;
}
Ostream& operator<<(Ostream& os, u128 v)
{
Str ans;
if (v == 0) { ans = "0"; }
while (v) {
ans.push_back('0' + v % 10), v /= 10;
}
std::reverse(ans.begin(), ans.end());
return os << ans;
}
template<typename T>
bool chmin(T& a, const T& b)
{
if (a > b) {
a = b;
return true;
} else {
return false;
}
}
template<typename T>
bool chmax(T& a, const T& b)
{
if (a < b) {
a = b;
return true;
} else {
return false;
}
}
template<typename T>
constexpr T floorDiv(T x, T y)
{
if (y < T{}) { x = -x, y = -y; }
return x >= T{} ? x / y : (x - y + 1) / y;
}
template<typename T>
constexpr T ceilDiv(T x, T y)
{
if (y < T{}) { x = -x, y = -y; }
return x >= T{} ? (x + y - 1) / y : x / y;
}
template<typename T, typename I>
constexpr T modPower(T v, I n, T mod)
{
T ans = 1 % mod;
for (; n > 0; n >>= 1, (v *= v) %= mod) {
if (n % 2 == 1) { (ans *= v) %= mod; }
}
return ans;
}
template<typename T, typename I>
constexpr T power(T v, I n)
{
T ans = 1;
for (; n > 0; n >>= 1, v *= v) {
if (n % 2 == 1) { ans *= v; }
}
return ans;
}
template<typename T, typename I>
constexpr T power(T v, I n, const T& e)
{
T ans = e;
for (; n > 0; n >>= 1, v *= v) {
if (n % 2 == 1) { ans *= v; }
}
return ans;
}
template<typename T>
Vec<T> operator+=(Vec<T>& vs1, const Vec<T>& vs2)
{
vs1.insert(vs1.end(), vs2.begin(), vs2.end());
return vs1;
}
template<typename T>
Vec<T> operator+(const Vec<T>& vs1, const Vec<T>& vs2)
{
auto vs = vs1;
vs += vs2;
return vs;
}
template<typename Vs, typename V>
void fillAll(Vs& arr, const V& v)
{
if constexpr (std::is_convertible<V, Vs>::value) {
arr = v;
} else {
for (auto& subarr : arr) {
fillAll(subarr, v);
}
}
}
template<typename Vs>
void sortAll(Vs& vs)
{
std::sort(std::begin(vs), std::end(vs));
}
template<typename Vs, typename C>
void sortAll(Vs& vs, C comp)
{
std::sort(std::begin(vs), std::end(vs), comp);
}
template<typename Vs>
void reverseAll(Vs& vs)
{
std::reverse(std::begin(vs), std::end(vs));
}
template<typename V, typename Vs>
V sumAll(const Vs& vs)
{
if constexpr (std::is_convertible<Vs, V>::value) {
return static_cast<V>(vs);
} else {
V ans = 0;
for (const auto& v : vs) {
ans += sumAll<V>(v);
}
return ans;
}
}
template<typename Vs>
int minInd(const Vs& vs)
{
return std::min_element(std::begin(vs), std::end(vs)) - std::begin(vs);
}
template<typename Vs>
int maxInd(const Vs& vs)
{
return std::max_element(std::begin(vs), std::end(vs)) - std::begin(vs);
}
template<typename Vs, typename V>
int lbInd(const Vs& vs, const V& v)
{
return std::lower_bound(std::begin(vs), std::end(vs), v) - std::begin(vs);
}
template<typename Vs, typename V>
int ubInd(const Vs& vs, const V& v)
{
return std::upper_bound(std::begin(vs), std::end(vs), v) - std::begin(vs);
}
template<typename T, typename F>
Vec<T> genVec(int n, F gen)
{
Vec<T> ans;
std::generate_n(std::back_insert_iterator(ans), n, gen);
return ans;
}
Vec<int> iotaVec(int n, int offset = 0)
{
Vec<int> ans(n);
std::iota(ans.begin(), ans.end(), offset);
return ans;
}
constexpr int popcount(const u64 v)
{
return v ? __builtin_popcountll(v) : 0;
}
constexpr int log2p1(const u64 v)
{
return v ? 64 - __builtin_clzll(v) : 0;
}
constexpr int lsbp1(const u64 v)
{
return __builtin_ffsll(v);
}
constexpr int clog(const u64 v)
{
return v ? log2p1(v - 1) : 0;
}
constexpr u64 ceil2(const u64 v)
{
const int l = clog(v);
return (l == 64) ? 0_u64 : (1_u64 << l);
}
constexpr u64 floor2(const u64 v)
{
return v ? (1_u64 << (log2p1(v) - 1)) : 0_u64;
}
constexpr bool ispow2(const u64 v)
{
return (v > 0) and ((v & (v - 1)) == 0);
}
constexpr bool btest(const u64 mask, const int ind)
{
return (mask >> ind) & 1_u64;
}
template<typename F>
struct Fix : F
{
Fix(F&& f) : F{std::forward<F>(f)} {}
template<typename... Args>
auto operator()(Args&&... args) const
{
return F::operator()(*this, std::forward<Args>(args)...);
}
};
class irange
{
private:
struct itr
{
itr(i64 start = 0, i64 step = 1) : m_cnt{start}, m_step{step} {}
bool operator!=(const itr& it) const
{
return m_cnt != it.m_cnt;
}
int operator*()
{
return m_cnt;
}
itr& operator++()
{
m_cnt += m_step;
return *this;
}
i64 m_cnt, m_step;
};
i64 m_start, m_end, m_step;
public:
irange(i64 start, i64 end, i64 step = 1)
{
assert(step != 0);
const i64 d = std::abs(step);
const i64 l = (step > 0 ? start : end);
const i64 r = (step > 0 ? end : start);
int n = (r - l) / d + ((r - l) % d ? 1 : 0);
if (l >= r) { n = 0; }
m_start = start;
m_end = start + step * n;
m_step = step;
}
itr begin() const
{
return itr{m_start, m_step};
}
itr end() const
{
return itr{m_end, m_step};
}
};
irange rep(int end)
{
return irange(0, end, 1);
}
irange per(int rend)
{
return irange(rend - 1, -1, -1);
}
#pragma COMMENT("[REFS] Xoshiro: https://prng.di.unimi.it")
namespace xoshiro_impl {
u64 x;
u64 next()
{
uint64_t z = (x += 0x9e3779b97f4a7c15);
z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
return z ^ (z >> 31);
}
}
class Xoshiro32
{
public:
using result_type = u32;
using T = result_type;
Xoshiro32(T seed = 0)
{
xoshiro_impl::x = seed;
s[0] = xoshiro_impl::next();
s[1] = xoshiro_impl::next();
s[2] = xoshiro_impl::next();
s[3] = xoshiro_impl::next();
}
static constexpr T min()
{
return LIMMIN<T>;
}
static constexpr T max()
{
return LIMMAX<T>;
}
T operator()()
{
return next();
}
private:
static constexpr T rotl(const T x, int k)
{
return (x << k) | (x >> (32 - k));
}
T next()
{
const T ans = rotl(s[1] * 5, 7) * 9;
const T t = s[1] << 9;
s[2] ^= s[0];
s[3] ^= s[1];
s[1] ^= s[2];
s[0] ^= s[3];
s[2] ^= t;
s[3] = rotl(s[3], 11);
return ans;
}
T s[4];
};
class Xoshiro64
{
public:
using result_type = u64;
using T = result_type;
Xoshiro64(T seed = 0)
{
xoshiro_impl::x = seed;
s[0] = xoshiro_impl::next();
s[1] = xoshiro_impl::next();
s[2] = xoshiro_impl::next();
s[3] = xoshiro_impl::next();
}
static constexpr T min()
{
return LIMMIN<T>;
}
static constexpr T max()
{
return LIMMAX<T>;
}
T operator()()
{
return next();
}
private:
static constexpr T rotl(const T x, int k)
{
return (x << k) | (x >> (64 - k));
}
T next()
{
const T ans = rotl(s[1] * 5, 7) * 9;
const T t = s[1] << 17;
s[2] ^= s[0];
s[3] ^= s[1];
s[1] ^= s[2];
s[0] ^= s[3];
s[2] ^= t;
s[3] = rotl(s[3], 45);
return ans;
}
T s[4];
};
template<typename Rng>
class RNG
{
public:
using result_type = typename Rng::result_type;
using T = result_type;
static constexpr T min()
{
return Rng::min();
}
static constexpr T max()
{
return Rng::max();
}
RNG() : RNG(std::random_device{}()) {}
RNG(T seed) : m_rng(seed) {}
T operator()()
{
return m_rng();
}
template<typename T>
T val(T min, T max)
{
return std::uniform_int_distribution<T>(min, max)(m_rng);
}
template<typename T>
Pair<T, T> pair(T min, T max)
{
return std::minmax({val<T>(min, max), val<T>(min, max)});
}
template<typename T>
Vec<T> vec(int n, T min, T max)
{
return genVec<T>(n, [&]() { return val<T>(min, max); });
}
template<typename T>
Vec<Vec<T>> vvec(int n, int m, T min, T max)
{
return genVec<Vec<T>>(n, [&]() { return vec(m, min, max); });
}
private:
Rng m_rng;
};
RNG<std::mt19937> rng;
RNG<std::mt19937_64> rng64;
RNG<Xoshiro32> rng_xo;
RNG<Xoshiro64> rng_xo64;
class Scanner
{
public:
Scanner(Istream& is = std::cin) : m_is{is}
{
m_is.tie(nullptr)->sync_with_stdio(false);
}
template<typename T>
T val()
{
T v;
return m_is >> v, v;
}
template<typename T>
T val(T offset)
{
return val<T>() - offset;
}
template<typename T>
Vec<T> vec(int n)
{
return genVec<T>(n, [&]() { return val<T>(); });
}
template<typename T>
Vec<T> vec(int n, T offset)
{
return genVec<T>(n, [&]() { return val<T>(offset); });
}
template<typename T>
Vec<Vec<T>> vvec(int n, int m)
{
return genVec<Vec<T>>(n, [&]() { return vec<T>(m); });
}
template<typename T>
Vec<Vec<T>> vvec(int n, int m, const T offset)
{
return genVec<Vec<T>>(n, [&]() { return vec<T>(m, offset); });
}
template<typename... Args>
auto tup()
{
return Tup<Args...>{val<Args>()...};
}
template<typename... Args>
auto tup(const Args&... offsets)
{
return Tup<Args...>{val<Args>(offsets)...};
}
private:
Istream& m_is;
};
Scanner in;
class Printer
{
public:
Printer(Ostream& os = std::cout) : m_os{os}
{
m_os << std::fixed << std::setprecision(15);
}
template<typename... Args>
int operator()(const Args&... args)
{
dump(args...);
return 0;
}
template<typename... Args>
int ln(const Args&... args)
{
dump(args...), m_os << '\n';
return 0;
}
template<typename... Args>
int el(const Args&... args)
{
dump(args...), m_os << std::endl;
return 0;
}
private:
template<typename T>
void dump(const T& v)
{
m_os << v;
}
template<typename T>
void dump(const Vec<T>& vs)
{
for (const int i : rep(vs.size())) {
m_os << (i ? " " : ""), dump(vs[i]);
}
}
template<typename T>
void dump(const Vec<Vec<T>>& vss)
{
for (const int i : rep(vss.size())) {
m_os << (i ? "\n" : ""), dump(vss[i]);
}
}
template<typename T, typename... Ts>
int dump(const T& v, const Ts&... args)
{
dump(v), m_os << ' ', dump(args...);
return 0;
}
Ostream& m_os;
};
Printer out;
template<typename T, int n, int i = 0>
auto ndVec(int const (&szs)[n], const T x = T{})
{
if constexpr (i == n) {
return x;
} else {
return std::vector(szs[i], ndVec<T, n, i + 1>(szs, x));
}
}
template<typename T, typename F>
T binSearch(T ng, T ok, F check)
{
while (std::abs(ok - ng) > 1) {
const T mid = (ok + ng) / 2;
(check(mid) ? ok : ng) = mid;
}
return ok;
}
template<u32 mod_, u32 root_, u32 max2p_>
class modint
{
template<typename U = u32&>
static U modRef()
{
static u32 s_mod = 0;
return s_mod;
}
template<typename U = u32&>
static U rootRef()
{
static u32 s_root = 0;
return s_root;
}
template<typename U = u32&>
static U max2pRef()
{
static u32 s_max2p = 0;
return s_max2p;
}
public:
template<typename U = const u32>
static constexpr std::enable_if_t<mod_ != 0, U> mod()
{
return mod_;
}
template<typename U = const u32>
static std::enable_if_t<mod_ == 0, U> mod()
{
return modRef();
}
template<typename U = const u32>
static constexpr std::enable_if_t<mod_ != 0, U> root()
{
return root_;
}
template<typename U = const u32>
static std::enable_if_t<mod_ == 0, U> root()
{
return rootRef();
}
template<typename U = const u32>
static constexpr std::enable_if_t<mod_ != 0, U> max2p()
{
return max2p_;
}
template<typename U = const u32>
static std::enable_if_t<mod_ == 0, U> max2p()
{
return max2pRef();
}
template<typename U = u32>
static void setMod(std::enable_if_t<mod_ == 0, U> m)
{
modRef() = m;
}
template<typename U = u32>
static void setRoot(std::enable_if_t<mod_ == 0, U> r)
{
rootRef() = r;
}
template<typename U = u32>
static void setMax2p(std::enable_if_t<mod_ == 0, U> m)
{
max2pRef() = m;
}
constexpr modint() : m_val{0} {}
constexpr modint(i64 v) : m_val{normll(v)} {}
constexpr void setRaw(u32 v)
{
m_val = v;
}
constexpr modint operator-() const
{
return modint{0} - (*this);
}
constexpr modint& operator+=(const modint& m)
{
m_val = norm(m_val + m.val());
return *this;
}
constexpr modint& operator-=(const modint& m)
{
m_val = norm(m_val + mod() - m.val());
return *this;
}
constexpr modint& operator*=(const modint& m)
{
m_val = normll((i64)m_val * (i64)m.val() % (i64)mod());
return *this;
}
constexpr modint& operator/=(const modint& m)
{
return *this *= m.inv();
}
constexpr modint operator+(const modint& m) const
{
auto v = *this;
return v += m;
}
constexpr modint operator-(const modint& m) const
{
auto v = *this;
return v -= m;
}
constexpr modint operator*(const modint& m) const
{
auto v = *this;
return v *= m;
}
constexpr modint operator/(const modint& m) const
{
auto v = *this;
return v /= m;
}
constexpr bool operator==(const modint& m) const
{
return m_val == m.val();
}
constexpr bool operator!=(const modint& m) const
{
return not(*this == m);
}
friend Istream& operator>>(Istream& is, modint& m)
{
i64 v;
return is >> v, m = v, is;
}
friend Ostream& operator<<(Ostream& os, const modint& m)
{
return os << m.val();
}
constexpr u32 val() const
{
return m_val;
}
template<typename I>
constexpr modint pow(I n) const
{
return power(*this, n);
}
constexpr modint inv() const
{
return pow(mod() - 2);
}
static modint sinv(u32 n)
{
static Vec<modint> is{1, 1};
for (u32 i = (u32)is.size(); i <= n; i++) {
is.push_back(-is[mod() % i] * (mod() / i));
}
return is[n];
}
static modint fact(u32 n)
{
static Vec<modint> fs{1, 1};
for (u32 i = (u32)fs.size(); i <= n; i++) {
fs.push_back(fs.back() * i);
}
return fs[n];
}
static modint ifact(u32 n)
{
static Vec<modint> ifs{1, 1};
for (u32 i = (u32)ifs.size(); i <= n; i++) {
ifs.push_back(ifs.back() * sinv(i));
}
return ifs[n];
}
static modint comb(int n, int k)
{
return k > n or k < 0 ? modint{0} : fact(n) * ifact(n - k) * ifact(k);
}
private:
static constexpr u32 norm(u32 x)
{
return x < mod() ? x : x - mod();
}
static constexpr u32 normll(i64 x)
{
return norm(u32(x % (i64)mod() + (i64)mod()));
}
u32 m_val;
};
using modint_1000000007 = modint<1000000007, 5, 1>;
using modint_998244353 = modint<998244353, 3, 23>;
template<int id>
using modint_dynamic = modint<0, 0, id>;
template<typename T = int>
class Graph
{
struct Edge
{
Edge() = default;
Edge(int i, int t, T c) : id{i}, to{t}, cost{c} {}
int id;
int to;
T cost;
operator int() const
{
return to;
}
};
public:
Graph(int n) : m_v{n}, m_edges(n) {}
void addEdge(int u, int v, bool bi = false)
{
assert(0 <= u and u < m_v);
assert(0 <= v and v < m_v);
m_edges[u].emplace_back(m_e, v, 1);
if (bi) { m_edges[v].emplace_back(m_e, u, 1); }
m_e++;
}
void addEdge(int u, int v, const T& c, bool bi = false)
{
assert(0 <= u and u < m_v);
assert(0 <= v and v < m_v);
m_edges[u].emplace_back(m_e, v, c);
if (bi) { m_edges[v].emplace_back(m_e, u, c); }
m_e++;
}
const Vec<Edge>& operator[](const int u) const
{
assert(0 <= u and u < m_v);
return m_edges[u];
}
Vec<Edge>& operator[](const int u)
{
assert(0 <= u and u < m_v);
return m_edges[u];
}
int v() const
{
return m_v;
}
int e() const
{
return m_e;
}
friend Ostream& operator<<(Ostream& os, const Graph& g)
{
for (int u : rep(g.v())) {
for (const auto& [id, v, c] : g[u]) {
os << "[" << id << "]: ";
os << u << "->" << v << "(" << c << ")\n";
}
}
return os;
}
Vec<T> sizes(int root = 0) const
{
const int N = v();
assert(0 <= root and root < N);
Vec<T> ss(N, 1);
Fix([&](auto dfs, int u, int p) -> void {
for (const auto& [id, v, c] : m_edges[u]) {
static_cast<void>(id);
if (v == p) { continue; }
dfs(v, u);
ss[u] += ss[v];
}
})(root, -1);
return ss;
}
Vec<T> depths(int root = 0) const
{
const int N = v();
assert(0 <= root and root < N);
Vec<T> ds(N, 0);
Fix([&](auto dfs, int u, int p) -> void {
for (const auto& [id, v, c] : m_edges[u]) {
static_cast<void>(id);
if (v == p) { continue; }
ds[v] = ds[u] + c;
dfs(v, u);
}
})(root, -1);
return ds;
}
Vec<int> parents(int root = 0) const
{
const int N = v();
assert(0 <= root and root < N);
Vec<int> ps(N, -1);
Fix([&](auto dfs, int u, int p) -> void {
for (const auto& [id, v, c] : m_edges[u]) {
static_cast<void>(id);
if (v == p) { continue; }
ps[v] = u;
dfs(v, u);
}
})(root, -1);
return ps;
}
private:
int m_v;
int m_e = 0;
Vec<Vec<Edge>> m_edges;
};
#pragma endregion
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;
typedef C Cost;
public:
enum ProblemType { INFEASIBLE, OPTIMAL, UNBOUNDED };
enum SupplyType {
GEQ,
LEQ
};
enum PivotRule {
FIRST_ELIGIBLE,
BEST_ELIGIBLE,
BLOCK_SEARCH,
CANDIDATE_LIST,
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:
const Digraph& _graph;
int _node_num;
int _arc_num;
int _all_arc_num;
int _search_arc_num;
bool _has_lower;
SupplyType _stype;
Value _sum_supply;
IntVector _source;
IntVector _target;
ValueVector _lower;
ValueVector _upper;
ValueVector _cap;
CostVector _cost;
ValueVector _supply;
ValueVector _flow;
CostVector _pi;
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;
int in_arc, join, u_in, v_in, u_out, v_out;
Value delta;
const Value MAX;
public:
const Value INF;
private:
class FirstEligiblePivotRule
{
private:
const IntVector& _source;
const IntVector& _target;
const CostVector& _cost;
const CharVector& _state;
const CostVector& _pi;
int& _in_arc;
int _search_arc_num;
int _next_arc;
public:
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)
{}
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 BestEligiblePivotRule
{
private:
const IntVector& _source;
const IntVector& _target;
const CostVector& _cost;
const CharVector& _state;
const CostVector& _pi;
int& _in_arc;
int _search_arc_num;
public:
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)
{}
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 BlockSearchPivotRule
{
private:
const IntVector& _source;
const IntVector& _target;
const CostVector& _cost;
const CharVector& _state;
const CostVector& _pi;
int& _in_arc;
int _search_arc_num;
int _block_size;
int _next_arc;
public:
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)
{
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);
}
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 CandidateListPivotRule
{
private:
const IntVector& _source;
const IntVector& _target;
const CostVector& _cost;
const CharVector& _state;
const CostVector& _pi;
int& _in_arc;
int _search_arc_num;
IntVector _candidates;
int _list_length, _minor_limit;
int _curr_length, _minor_count;
int _next_arc;
public:
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)
{
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);
}
bool findEnteringArc()
{
Cost min, c;
int e;
if (_curr_length > 0 && _minor_count < _minor_limit) {
++_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;
}
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 AlteringListPivotRule
{
private:
const IntVector& _source;
const IntVector& _target;
const CostVector& _cost;
const CharVector& _state;
const CostVector& _pi;
int& _in_arc;
int _search_arc_num;
int _block_size, _head_length, _curr_length;
int _next_arc;
IntVector _candidates;
CostVector _cand_cost;
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:
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)
{
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;
}
bool findEnteringArc()
{
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];
}
}
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:
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);
_in_arc = _candidates[0];
_next_arc = e;
_candidates[0] = _candidates[new_length - 1];
_curr_length = new_length - 1;
return true;
}
};
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)
{
static_assert(std::numeric_limits<Value>::is_signed,
"Value must be signed");
static_assert(std::numeric_limits<Cost>::is_signed,
"Cost must be signed");
static_assert(std::numeric_limits<Value>::max() > 0,
"max() must be greater than 0");
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;
}
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)
{
for (int i = 0; i != _node_num; ++i)
_supply[i] = 0;
_supply[s] = k, _supply[t] = -k;
return *this;
}
NetworkSimplex& supplyType(SupplyType supply_type)
{
_stype = supply_type;
return *this;
}
ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH)
{
if (!init()) return INFEASIBLE;
return start(pivot_rule);
}
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;
}
NetworkSimplex& reset()
{
_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);
}
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:
bool init()
{
if (_node_num == 0) return false;
_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;
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];
}
}
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;
}
for (int i = 0; i != _arc_num; ++i) {
_flow[i] = 0;
_state[i] = STATE_LOWER;
}
_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;
if (_sum_supply == 0) {
_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) {
_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 {
_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;
}
bool checkBoundMaps()
{
for (int j = 0; j != _arc_num; ++j) {
if (_upper[j] < _lower[j]) return false;
}
return true;
}
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;
}
bool findLeavingArc()
{
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;
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;
}
}
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;
}
void changeFlow(bool change)
{
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;
}
}
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];
}
}
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];
if (u_in == u_out) {
_parent[u_in] = v_in;
_pred[u_in] = in_arc;
_pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
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 {
int thread_continue = old_rev_thread == v_in
? _thread[old_last_succ]
: _thread[v_in];
int stem = u_in;
int par_stem = v_in;
int next_stem;
int last = _last_succ[u_in];
int before, after = _thread[last];
_thread[v_in] = u_in;
_dirty_revs.clear();
_dirty_revs.push_back(v_in);
while (stem != u_out) {
next_stem = _parent[stem];
_thread[last] = next_stem;
_dirty_revs.push_back(last);
before = _rev_thread[stem];
_thread[before] = after;
_rev_thread[after] = before;
_parent[stem] = par_stem;
par_stem = stem;
stem = next_stem;
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;
if (old_rev_thread != v_in) {
_thread[old_rev_thread] = after;
_rev_thread[after] = old_rev_thread;
}
for (int i = 0; i != int(_dirty_revs.size()); ++i) {
int u = _dirty_revs[i];
_rev_thread[_thread[u]] = u;
}
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;
}
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;
}
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;
}
}
for (int u = v_in; u != join; u = _parent[u]) {
_succ_num[u] += old_succ_num;
}
for (int u = v_out; u != join; u = _parent[u]) {
_succ_num[u] -= old_succ_num;
}
}
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;
}
}
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) {
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 (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 {
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]) {
c = _cost[a];
if (c < min_cost) {
min_cost = c;
min_arc = a;
}
}
if (min_arc != INVALID) { arc_vector.push_back(min_arc); }
}
}
} else {
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); }
}
}
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;
}
ProblemType start(PivotRule pivot_rule)
{
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;
}
template<typename PivotRuleImpl>
ProblemType start()
{
PivotRuleImpl pivot(*this);
if (!initialPivots()) return UNBOUNDED;
while (pivot.findEnteringArc()) {
findJoinNode();
bool change = findLeavingArc();
if (delta >= MAX) return UNBOUNDED;
changeFlow(change);
if (change) {
updateTreeStructure();
updatePotential();
}
}
for (int e = _search_arc_num; e != _all_arc_num; ++e) {
if (_flow[e] != 0) return INFEASIBLE;
}
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;
}
}
}
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;
}
};
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;
friend Ostream& operator<<(Ostream& os, const edge& e)
{
return (os << e.from << "->" << e.to << ":" << e.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 set_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>();
}
}
};
int main()
{
const auto [N, M] = in.tup<int, int>();
const auto [X, Y, Z] = in.tup<int, int, int>(1, 1, 1);
auto g = ndVec({N, N}, 1);
for (int i : rep(N)) {
g[i][i] = 0;
}
for (int i : rep(M)) {
const auto [A, B] = in.tup<int, int>(1, 1);
g[A][B] = g[B][A] = 0;
}
auto vin = [&, X = X, Y = Y, Z = Z](int i) {
if (i == X) { return 2 * X; }
if (i == Y) { return 2 * Y; }
return 2 * i;
};
auto vout = [&, X = X, Y = Y, Z = Z, N = N](int i) {
if (i == X) { return 2 * X; }
if (i == Y) { return 2 * Y; }
return 2 * i + 1;
};
mcf_graph_ns<int, i64> flow(4 * N);
for (int i : rep(N)) {
if (g[X][i]) { flow.add_edge(vout(X), vin(i), 0, 1, 0); }
}
for (int i : rep(N)) {
if (g[i][Y]) {
flow.add_edge(vout(i), vin(Y), 0, 1, 0);
flow.add_edge(vout(i) + 2 * N, vin(Y) + 2 * N, 0, 1, 0);
}
}
for (int i : rep(N)) {
if (i == X or i == Y) { continue; }
if (i == Z) {
flow.add_edge(vin(i), vout(i) + 2 * N, 0, 1, 1);
} else {
flow.add_edge(vin(i), vout(i), 0, 1, 1);
flow.add_edge(vin(i) + 2 * N, vout(i) + 2 * N, 0, 1, 1);
}
}
for (int i : rep(N)) {
if (i == X or i == Y) { continue; }
for (int j : rep(N)) {
if (j == X or j == Y) { continue; }
if (g[i][j] == 0) { continue; }
flow.add_edge(vout(i), vin(j), 0, 1, 0);
flow.add_edge(vout(i) + 2 * N, vin(j) + 2 * N, 0, 1, 0);
}
}
flow.set_supply(vin(X), 2);
flow.set_supply(vin(Y), -1);
flow.set_supply(vin(Y) + 2 * N, -1);
i64 ans = flow.solve() + 2;
if (flow.infeasible) {
out.ln(-1);
} else {
void(0);
void(0);
void(0);
for (int i : rep(flow.flow.size())) {
if (flow.flow[i]) { void(0); }
}
out.ln(ans);
}
return 0;
}