Cod sursa(job #3031120)

Utilizator lucametehauDart Monkey lucametehau Data 18 martie 2023 19:23:38
Problema Flux maxim Scor 0
Compilator cpp-64 Status done
Runda Arhiva educationala Marime 20.25 kb
#include <iostream>
#include <cmath>
#include <vector>
#include <map>
#include <unordered_map>
#include <set>
#include <cstring>
#include <chrono>
#include <cassert>
#include <bitset>
#include <stack>
#include <queue>
#include <iomanip>
#include <random>
#include <string>
#include <complex>
#include <algorithm>
#include <numeric>

#ifdef _MSC_VER
#  include <intrin.h>
#  define __builtin_popcount __popcnt
#  define __builtin_popcountll __popcnt64
#endif

#pragma warning(disable : 4996)

#pragma GCC optimize("Ofast")
#define x first
#define y second
#define ld long double
#define ll long long
#define ull unsigned long long
#define us unsigned short
#define lsb(x) ((x) & (-(x)))
#define pii pair <int, int>
#define pll pair <ll, ll>

using namespace std;

//const int MOD = (int)1e9 + 7;

mt19937_64 gen(chrono::duration_cast<chrono::nanoseconds>(chrono::steady_clock::now().time_since_epoch()).count());
uniform_int_distribution <ull> rng;
const ld PI = acos(-1);

// the following namespace includes many useful things
// for solving linear recurrences
namespace recurrences {
    //const int MOD = (int)1e9 + 7;
    // binary exponentiation
    template <int MOD>
    int lgput(int n, int p) {
        int ans = 1, x = n;

        while (p) {
            if (p & 1)
                ans = 1LL * ans * x % MOD;
            x = 1LL * x * x % MOD;
            p >>= 1;
        }

        return ans;
    }

    // modular integer class
    template <int MOD>
    struct Int {
        int x;

        Int() {
            x = 0;
        }

        Int(int _x) {
            if (_x < 0)
                _x += MOD;
            if (_x >= MOD)
                _x -= MOD;
            x = _x;
        }

        friend ostream& operator << (ostream& os, const Int& X) {
            os << (false ? X.x - MOD : X.x);
            return os;
        }

        Int operator + (const Int& other) const {
            int val = x + other.x;

            return (val >= MOD ? val - MOD : val);
        }

        Int operator += (const Int& other) {
            return *this = *this + other;
        }

        Int operator - (const Int& other) const {
            int val = x - other.x;

            return (val < 0 ? val + MOD : val);
        }
        Int operator -= (const Int& other) {
            return *this = *this - other;
        }

        Int operator * (const Int& other) const {
            return 1LL * x * other.x % MOD;
        }

        Int operator *= (const Int& other) {
            return *this = *this * other;
        }

        Int operator / (const Int& other) const {
            return 1LL * x * other.inv() % MOD;
        }

        Int operator /= (const Int& other) {
            return *this = *this / other;
        }

        bool operator == (const Int& other) const {
            return x == other.x;
        }

        bool operator != (const Int& other) const {
            return x != other.x;
        }

        Int pow(int p) const {
            return lgput<MOD>(x, p);
        }

        int inv() const {
            return lgput<MOD>(x, MOD - 2);
        }

        void nice_print() const {
            for (int d = 2; d <= 1000; d++) {
                //cout << d << " " << (*this * d).x << "\n";
                if ((*this * d).x <= 1000) {
                    cout << (*this * d).x << "/" << d << "\n";
                    break;
                }
            }
        }
    };

    const bool slow_mult = false;

    template <int MOD>
    int get_primitive_root() {
        vector<int> fact;
        int phi = MOD - 1, n = phi;
        for (int i = 2; i * i <= n; ++i)
            if (n % i == 0) {
                fact.push_back(i);
                while (n % i == 0)
                    n /= i;
            }
        if (n > 1)
            fact.push_back(n);

        for (int res = 2; res <= MOD; ++res) {
            bool ok = true;
            for (size_t i = 0; i < fact.size() && ok; ++i)
                ok &= lgput<MOD>(res, phi / fact[i]) != 1;
            if (ok) {
                return res;
            }
        }
        return -1;
    }

    int get_smallest_power(int n) {
        int p = 1;
        while (p < n)
            p <<= 1;
        return p;
    }

    // modular 
    template <int MOD>
    struct Poly {
        vector <Int<MOD>> p;

        Poly() {
            p.clear();
        }

        Poly(vector <Int<MOD>> values) {
            p = values;
        }

        Poly(int val) {
            p = { val };
        }

        Int<MOD>& operator [] (int index) {
            assert(index < (int)p.size());
            return p[index];
        }

        void setDegree(int deg) {
            p.resize(deg + 1);
        }

        int deg() const {
            return p.size() - 1;
        }

        friend ostream& operator << (ostream& os, const Poly& P) {
            for (auto& i : P.p)
                os << i << " ";
            return os;
        }

        bool operator == (const Poly& other) const {
            if (deg() != other.deg())
                return 0;

            for (int i = 0; i <= deg(); i++) {
                if (p[i] != other.p[i])
                    return 0;
            }

            return 1;
        }

        Poly operator + (const Poly& other) const {
            Poly sum(p);

            sum.setDegree(max(deg(), other.deg()));

            for (int i = 0; i <= other.deg(); i++)
                sum[i] += other.p[i];

            return sum;
        }

        // add
        Poly operator += (const Poly& other) {
            return *this = *this + other;
        }

        Poly operator - (const Poly& other) const {
            Poly dif(p);

            dif.setDegree(max(deg(), other.deg()));

            for (int i = 0; i <= other.deg(); i++)
                dif[i] -= other.p[i];

            return dif;
        }

        // substract
        Poly operator -= (const Poly& other) {
            return *this = *this - other;
        }

        // scalar multiplication
        Poly operator * (const Int<MOD>& other) const {
            Poly mult(*this);

            for (auto& i : mult.p)
                i *= other;

            return mult;
        }
        // scalar multiplication
        Poly operator *= (const Int<MOD>& other) {
            return *this = *this * other;
        }

        // scalar division
        Poly operator / (const Int<MOD>& other) const {
            Poly mult(*this);
            Int<MOD> val = other.inv();

            for (auto& i : mult.p)
                i *= val;

            return mult;
        }

        // scalar division
        Poly operator /= (const Int<MOD>& other) {
            return *this = *this / other;
        }

        Poly fft(bool invert) {
            Poly Ans(p);

            static int primitive_root = get_primitive_root<MOD>();
            Int<MOD> root(primitive_root);

            int ind = 0, n = deg();
            for (int i = 1; i < n; i++) {
                int b;
                for (b = n / 2; ind & b; b >>= 1)
                    ind ^= b;
                ind ^= b;

                if (i < ind)
                    swap(Ans[i], Ans[ind]);
            }

            for (int l = 2, p = 0; l <= n; l <<= 1, p++) {
                assert((MOD - 1) % ((1 << (p + 1))) == 0);
                Int<MOD> val = root.pow((MOD - 1) >> (p + 1));
                Int<MOD> bw(!invert ? val : val.inv());

                for (int i = 0; i < n; i += l) {
                    Int<MOD> w(1);

                    for (int j = 0; j < l / 2; j++) {
                        int i1 = i + j, i2 = i + j + l / 2;
                        Int<MOD> val1(Ans[i1]), val2(Ans[i2] * w);
                        Ans[i1] = val1 + val2;
                        Ans[i2] = val1 - val2;
                        w *= bw;
                    }
                }
            }

            if (invert) {
                Int<MOD> inv = Int<MOD>(n).inv();

                for (int i = 0; i < n; i++)
                    Ans[i] *= inv;
            }

            return Ans;
        }

        // multiplies 2 polynomials
        Poly operator * (const Poly& other) const {
            //cout << get_primitive_root<MOD>() << " " << MOD << "\n";
            if (!get_primitive_root<MOD>()) { // for small sizes, use naive multiplication
                Poly mult;
                mult.setDegree(deg() + other.deg());

                for (int i = 0; i <= deg(); i++) {
                    for (int j = 0; j <= other.deg(); j++)
                        mult[i + j] += p[i] * other.p[j];
                }

                return mult;
            }
            Poly A(p), B(other.p);

            int sz = max(get_smallest_power(A.deg() + 1), get_smallest_power(B.deg() + 1)) * 2;

            A.setDegree(sz), B.setDegree(sz);

            A = A.fft(0);
            B = B.fft(0);

            for (int i = 0; i < sz; i++)
                A[i] *= B[i];

            A = A.fft(1);

            A.setDegree(deg() + other.deg());

            return A;
        }

        // p mod q
        Poly operator % (const Poly& other) const {
            int d = deg() - other.deg();

            if (d < 0)
                return *this;

            if (false) {
                Poly R2(p);
                for (int i = deg(); i >= other.deg(); i--) {
                    R2 -= (other * R2[i] / other.p[other.deg()]).shift(i - other.deg());
                }

                R2.setDegree(other.deg() - 1);

                //return R;
            }

            Poly A, B = other;

            for (int i = 0; i <= d; i++) {
                A.p.push_back(p[deg() - i]);
            }

            for (int i = 0; i <= B.deg() / 2; i++)
                swap(B[i], B[B.deg() - i]);

            Poly C = A * B.inverse(d);

            C.setDegree(d);

            for (int i = 0; i <= d / 2; i++)
                swap(C[i], C[d - i]);

            Poly R = *this - other * C;

            R.setDegree(other.deg() - 1);

            return R;
        }

        // *x^n
        Poly shift(int index) {
            Poly q = p;
            q.setDegree(deg() + index);
            for (int i = deg(); i >= 0; i--)
                q[i + index] = q[i];
            for (int i = index - 1; i >= 0; i--)
                q[i] = Int<MOD>(0);
            return q;
        }

        // derivate of P
        Poly derivative() {
            Poly D;

            D.setDegree(deg() - 1);

            for (int i = 1; i <= deg(); i++)
                D[i - 1] = Int<MOD>(i) * p[i];

            return D;
        }

        // integral of P
        Poly integral() {
            Poly I;

            I.setDegree(deg() + 1);

            for (int i = 0; i <= deg(); i++)
                I[i + 1] = p[i] / Int<MOD>(i + 1);

            return I;
        }

        // P^-1 mod x^n
        Poly inverse(int n) {
            Poly Inv(p[0].inv()), two(2);

            int power = 1;

            while ((power / 2) <= n) {
                Poly A;
                for (int i = 0; i <= power; i++)
                    A.p.push_back((i <= deg() ? p[i] : 0));

                Inv = Inv * (two - A * Inv);
                Inv.setDegree(power);

                power <<= 1;
            }

            Inv.setDegree(n);

            return Inv;
        }

        // log(P) mod x^n
        Poly log(int n) {
            Poly Log(derivative());

            Log = Log * this->inverse(n);

            return Log.integral();
        }

        // e^P mod x^n
        Poly exp(int n) {
            Poly Exp(1);

            int power = 1;

            while ((power / 2) <= n) {
                Poly A(p);
                A.setDegree(power);
                Exp += Exp * A - Exp * Exp.log(power);
                Exp.setDegree(power);
                power <<= 1;
            }

            Exp.setDegree(n);
            return Exp;
        }

        // p^power mod mod, where mod is a polynomial
        Poly pow(uint64_t power, Poly mod) {
            Poly Pow(1), X(p);

            while (power) {
                if (power & 1)
                    Pow = Pow * X % mod;
                X = X * X % mod;
                power >>= 1;
            }

            return Pow;
        }

        Poly pow(uint64_t power) {
            Poly Pow(1), X(p);

            while (power) {
                if (power & 1)
                    Pow = Pow * X;
                X = X * X;
                power >>= 1;
                //cout << Pow << "\n";
            }

            return Pow;
        }
    };

    // berlekamp-massey algorithm
    template <int MOD>
    Poly<MOD> berlekamp_massey(vector <Int<MOD>> values) {
        Poly<MOD> P(values), B, C, Pol(-1);
        int n = values.size(), lst = -1;
        Int<MOD> lstError;

        for (int i = 0; i < n; i++) {
            Int<MOD> error = P[i];

            for (int j = 1; j <= C.deg() + 1; j++)
                error -= C[j - 1] * P[i - j];

            if (error == Int<MOD>(0))
                continue;

            if (lst == -1) {
                C = C.shift(i + 1);
                lst = i;
                lstError = P[i];
                continue;
            }

            Poly<MOD> D = (Pol * error / lstError);
            Poly<MOD> t = C;

            // instead of shifting D with i - lst - 1 positions
            // do the substraction on the last D.deg() coeficients
            if (i - lst - 1 + D.deg() > C.deg())
                C.setDegree(i - lst - 1 + D.deg());

            for (int j = 0; j <= D.deg(); j++)
                C[j + i - lst - 1] -= D[j];

            if (i - t.deg() > lst - B.deg()) {
                B = t;
                Pol = B;
                Pol = Pol.shift(1);
                Pol[0] = Int<MOD>(-1);
                lst = i;
                lstError = error;
            }

            //C -= D;
        }

        return C;
    }

    // find kth term based on terms of recurrence
    // assuming first term has index 1
    template <int MOD>
    Int<MOD> kth_term(vector <Int<MOD>> v, uint64_t k) {
        vector <Int<MOD>> x = { Int<MOD>(0), Int<MOD>(1) };
        Poly<MOD> CP = berlekamp_massey<MOD>(v);
        Poly<MOD> X(x);

        // characteristic polynomial is of form
        // x^n - sigma(i = 1..n, c[i] * x^(n-i))
        // that's why we need to reverse the recurrence
        // from berlekamp-massey

        CP *= Int<MOD>(-1);
        CP = CP.shift(1);
        for (int i = 0; i <= CP.deg() / 2; i++)
            swap(CP[i], CP[CP.deg() - i]);
        CP[CP.deg()] = 1;

        X = X.pow(k - 1, CP);

        Int<MOD> term(0);

        for (int i = 0; i <= X.deg(); i++)
            term += X[i] * v[i];

        return term;
    }

    template <int MOD>
    void berlekamp_massey_test() {
        int n;
        vector<Int<MOD>> v;

        cin >> n;
        for (int i = 0; i < n; i++) {
            int x;
            cin >> x;
            v.push_back(x);
        }

        Poly<MOD> rec = berlekamp_massey<MOD>(v);

        // sanity check

        for (int k = rec.deg() + 2; k <= 30; k++) {
            Int<MOD> val(0);

            for (int i = k - rec.deg() - 1; i < k; i++)
                val += v[i - 1] * rec[k - 1 - i];

            assert(val == kth_term(v, k));
            if (k > n)
                v.push_back(val);
        }

        uint64_t k;
        cin >> k;
        cout << kth_term<MOD>(v, k) << "\n";
    }

    template <int MOD>
    void berlekamp_massey_speed_test() {
        int n;
        vector<Int<MOD>> v;

        n = 1000;
        for (int i = 0; i < n; i++) {
            int x;
            //cin >> x;
            x = rng(gen);
            v.push_back(x);
        }

        ld t1 = clock();

        for (int k = 1; k <= 200; k++)
            kth_term<MOD>(v, k);

        ld t2 = clock();

        cout << (t2 - t1) / CLOCKS_PER_SEC << "s\n";
    }

    template <int MOD>
    void inverse_test() {
        int n;
        cin >> n;

        vector <Int<MOD>> v;
        for (int i = 1; i <= n; i++) {
            int x;
            cin >> x;

            v.push_back(Int<MOD>(x));
        }

        Poly<MOD> P(v), Inv, Exp, Log;
        int deg = 10;

        Inv = P.inverse(deg);

        cout << "Inverse: " << Inv << "\n";

        Exp = P.exp(deg);

        cout << "Exponential: " << Exp << "\n";

        Log = P.log(deg);

        cout << "Logarithm: " << Log << '\n';
    }
};

using namespace recurrences;

class Edge {
public:
    int to, cap, init_cap, ind;
};

class MaxFlow {
public:
    vector <vector <Edge>> g;
    vector <int> lvl, ind, path;
    vector <vector <int>> paths;

    void init(int n) {
        g.resize(n + 1);
        lvl.resize(n + 1);
        ind.resize(n + 1);
    }

    void add_edge(int x, int y, int cap = 1) {
        g[x].push_back({ y, cap, cap, (int)g[y].size() });
        g[y].push_back({ x, 0, 0, (int)g[x].size() - 1 });
    }

    bool bfs(int src, int dst) {
        queue <int> q;
        fill(lvl.begin(), lvl.end(), 0);
        lvl[src] = 1;
        q.push(src);

        while (!q.empty()) {
            int node = q.front();
            q.pop();

            for (auto& [son, cap, _, __] : g[node]) {
                if (!lvl[son] && cap > 0) {
                    lvl[son] = lvl[node] + 1;
                    q.push(son);
                }
            }
        }

        //cout << lvl[dst] << "\n";

        return lvl[dst] > 0;
    }

    int dfs(int node, int dst, int flow = (int)1e9) {
        //cout << node << " " << ind[node] << " " << flow << " heya\n";
        if (node == dst) {
            return flow;
        }


        for (; ind[node] < (int)g[node].size(); ind[node]++) {
            auto& [son, cap, _, idx] = g[node][ind[node]];
            if (cap <= 0 || lvl[son] != lvl[node] + 1)
                continue;

            int curr_flow = dfs(son, dst, min(flow, cap));

            if (curr_flow) {
                cap -= curr_flow;
                g[son][idx].cap += curr_flow;
                return curr_flow;
            }

            lvl[son] = 0;
        }

        return 0;
    }

    int max_flow(int src, int dst) {
        int max_flow = 0, flow;
        while (bfs(src, dst)) {
            //cout << src << " " << dst << "\n";
            fill(ind.begin(), ind.end(), 0);
            while (flow = dfs(src, dst)) {
                max_flow += flow;
            }
        }
        return max_flow;
    }

    void dfs_fill(int node, int dst) {
        path.push_back(node);
        if (node == dst) {
            return;
        }
        for (auto &[son, cap, init_cap, idx] : g[node]) {
            if (cap < init_cap) {
                cap++;
                g[son][idx].cap--;
                dfs_fill(son, dst);
                return;
            }
        }
    }

    void reconstruct(int src, int dst, int max_flow) {
        while (true) {
            int* a = new int[10000005];
            delete[] a;
        }
        while (max_flow) {
            path.clear();
            dfs_fill(src, dst);
            cout << path.size() << "\n";
            for (auto& node : path)
                cout << node << " ";
            cout << "\n";
            max_flow--;
        }
    }
} network;

void solve(int test = 0) {
    int n, m;
    cin >> n >> m;
    network.init(n);

    for (; m; m--) {
        int x, y, z;
        cin >> x >> y >> z;
        network.add_edge(x, y, z);
    }

    int flow = network.max_flow(1, n);
    cout << flow << "\n";
    network.reconstruct(1, n, flow);
}

int main() {
    ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);
    srand(time(0));

#ifdef LOCAL
    freopen("test.in", "r", stdin);
    freopen("test.out", "w", stdout);
#else
    freopen("maxflow.in", "r", stdin);
    freopen("maxflow.out", "w", stdout);
#endif

    int T = 1;

    //cin >> T;

    for (int t = 1; t <= T; t++) {
        solve(t);
    }

    return 0;
}