結果

問題 No.1775 Love Triangle 2
ユーザー Pachicobue
提出日時 2021-12-05 05:24:30
言語 C++17
(gcc 13.3.0 + boost 1.87.0)
結果
WA  
実行時間 -
コード長 61,639 bytes
コンパイル時間 3,339 ms
コンパイル使用メモリ 239,032 KB
最終ジャッジ日時 2025-01-26 03:59:42
ジャッジサーバーID
（参考情報） 		judge2 / judge3
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 24 WA * 66
権限があれば一括ダウンロードができます

ソースコード

diff #
プレゼンテーションモードにする

#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;
}
הההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההה
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
0