Cod sursa(job #2717103)

Utilizator lflorin29Florin Laiu lflorin29 Data 6 martie 2021 13:39:02
Problema Ciclu hamiltonian de cost minim Scor 0
Compilator cpp-64 Status done
Runda Arhiva educationala Marime 6.06 kb
#include <bits/stdc++.h>
using namespace std;

const double EPS = 1e-6;

#define VI vector<int>
#define VVD vector<vector<double>>
#define VD vector<double>
#define DOUBLE double

// rezolva LP min c^T * x
// cu Ax = b
struct LPSolver
{
    int m, n;
    VI B, N;
    VVD D;


    LPSolver(const VVD &A, const VD &b, const VD &c):
        m(b.size()), n(c.size()), N(n + 1), B(m), D(m + 2, VD(n + 2))
    {
        for (int i = 0; i < m; i++)
            for (int j = 0; j < n; j++)
                D[i][j] = A[i][j];
        for (int i = 0; i < m; i++)
        {
            B[i] = n + i;
            D[i][n] = -1;
            D[i][n + 1] = b[i];
        }
        for (int j = 0; j < n; j++)
        {
            N[j] = j;
            D[m][j] = -c[j];
        }
        N[n] = -1;
        D[m + 1][n] = 1;
    }

    void Pivot(int r, int s)
    {
        DOUBLE inv = 1.0L / D[r][s];
        for (int i = 0; i < m + 2; i++)
            if (i != r && abs(D[i][s]) > EPS / 100)
                for (int j = 0; j < n + 2; j++)
                    if (j != s)
                        D[i][j] -= D[r][j] * D[i][s] * inv;
        for (int j = 0; j < n + 2; j++)
            if (j != s)
                D[r][j] *= inv;
        for (int i = 0; i < m + 2; i++)
            if (i != r)
                D[i][s] *= -inv;
        D[r][s] = inv;
        swap(B[r], N[s]);
    }

    bool Simplex(int phase)
    {
        int x = phase == 1 ? m + 1 : m;
        while (true)
        {
            int s = -1;
            for (int j = 0; j <= n; j++)
            {
                if (phase == 2 && N[j] == -1)
                    continue;
                if (s == -1 || D[x][j] < D[x][s] || D[x][j] == D[x][s] && N[j] < N[s])
                    s = j;
            }
            if (D[x][s] > -EPS)
                return true;
            int r = -1;
            for (int i = 0; i < m; i++)
            {
                if (D[i][s] < EPS)
                    continue;
                if (r == -1 || D[i][n + 1] / D[i][s] < D[r][n + 1] / D[r][s] ||
                        (D[i][n + 1] / D[i][s]) == (D[r][n + 1] / D[r][s]) && B[i] < B[r])
                    r = i;
            }
            if (r == -1)
                return false;
            Pivot(r, s);
        }
    }

    DOUBLE Solve(VD &x)
    {
        int r = 0;
        for (int i = 1; i < m; i++)
            if (D[i][n + 1] < D[r][n + 1])
                r = i;
        if (D[r][n + 1] < -EPS)
        {
            Pivot(r, n);
            if (!Simplex(1) || D[m + 1][n + 1] < -EPS)
                return -numeric_limits<DOUBLE>::infinity();
            for (int i = 0; i < m; i++)
                if (B[i] == -1)
                {
                    int s = -1;
                    for (int j = 0; j <= n; j++)
                        if (s == -1 || D[i][j] < D[i][s] || D[i][j] == D[i][s] && N[j] < N[s])
                            s = j;
                    Pivot(i, s);
                }
        }
        if (!Simplex(2))
            return numeric_limits<DOUBLE>::infinity();
        x = VD(n);
        for (int i = 0; i < m; i++)
            if (B[i] < n)
                x[B[i]] = D[i][n + 1];
        return D[m][n + 1];
    }
};

int main()
{
    ifstream cin("hamilton.in");
    ofstream cout("hamilton.out");

    int n, m;
    cin >> n >> m;
    vector<vector<int>>c(n, vector<int>(n,-1000));
    for(int i(0); i < m; ++i)
    {
        int x, y, z;
        cin>>x>>y>>z;
        c[x][y] = -z;
    }

    auto idx = [&] (int i, int j)
    {
        return i * n + j;
    };


    VD C(n * n);
    for(int i = 0; i < n; ++i)
        for(int j = 0; j < n; ++j)
            C[idx(i, j)] = c[i][j];

    VVD A(3 * n, VD(n * n));
    VD B(3 * n);
    int k = 0;
    for(int i = 0; i < n; ++i)
    {
        B[k] = -1;
        // adaugam constragerile pt i
        for(int j = 0; j < n; ++j)
        {
            if(j != i)
                A[k][idx(i, j)] = -1;
        }


        k++;
        B[k] = -1;
        for(int j = 0; j < n; ++j)
        {
            if(j != i)
                A[k][idx(j, i)] = -1;
        }
        k++;
        A[k][idx(i, i)] = 1;
        B[k] = 0;
        k++;
    }

    VD x;
    LPSolver S(A, B, C);


    auto test = [&] (VD &x)
    {
        vector<int>go(n);
        for(int i = 0; i < n * n; ++i)
        {
            if(x[i] == 1)
            {
                int a = i / n, b = i % n;
                go[a] = b;
            }
        }
        vector<pair<int, int> > sol;
        int i = 0;
        vector<bool>vis(n);
        bool ok = 1;
        for(int i = 0; i < n; ++i)
        {
            if(vis[i])
                continue;
            int been = 0;
            int sz = 0;
            while(!vis[i])
            {
                been |= 1 << i;
                vis[i] = 1;
                sol.push_back(make_pair(i, go[i]));
                i = go[i];
                sz++;
            }
            if(sz != n)
            {
                ok = 0;
                // adaugam o variabila slack, nu afecteaza ec anterioare
                for(int i = 0; i < A.size(); ++i)
                    A[i].push_back(0);
                C.push_back(0); // variabilele slack nu influenteaza obiectivul
                A.push_back(VD(A.back().size())); // adaug o ecuatie noua
                B.push_back(sz - 1); //
                // trebe inca 1 constragere.
                for(int n1 = 0; n1 < n; ++n1)
                    if(been & (1 << n1))
                       for(int n2 = 0; n2 < n; n2++)
                          if(been & (1 << n2))
                             A[k][idx(n1, n2)] = 1;
                A[k].back() = 1; // pt ecuatia noua,variabila slack are coef 1 !
                k++;
            }
        }
        return make_pair(ok, sol);
    };


    long long ans = -1;
    while(true)
    {
        S.Solve(x);
        if(test(x).first == 1) break;
        S = LPSolver(A, B, C);
    }
    cout << -S.Solve(x) << '\n';

}