Cod sursa(job #2920682)

Utilizator iulianarsenoiuArsenoiu Iulian iulianarsenoiu Data 25 august 2022 12:57:03
Problema Componente tare conexe Scor 0
Compilator cpp-64 Status done
Runda Arhiva educationala Marime 5.44 kb
#include <bits/stdc++.h>

using namespace std;

ifstream f("network.in");
ofstream g("network.out");

const double limit = 1e-9;

int n,m;

vector<pair<int,pair<double,double>>> G[5005],B[5005];
vector<int> GC[5005],c[5005];

int componenta[5005];
int comp;

bool sel[5005];

vector<int> l,ltot;

double dp[5005][2];

int grad[5005];

double a[205][205];
int poz[205];

stack<int> st;

vector<pair<int,int>> e;

void dfs_sort(int nod)
{
    sel[nod] = true;
    for(auto it : GC[nod])
    {
        if(sel[it])
        {
            continue;
        }
        dfs_sort(it);
    }
    l.push_back(nod);
}

void dfs_sort_tot(int nod)
{
    sel[nod] = true;
    for(auto it : G[nod])
    {
        if(sel[it.first])
        {
            continue;
        }
        dfs_sort_tot(it.first);
    }
    ltot.push_back(nod);
}

void ctc(int nod)
{
    c[comp].push_back(nod);
    componenta[nod] = comp;
    sel[nod] = true;
    for(auto it : B[nod])
    {
        if(sel[it.first])
        {
            continue;
        }
        ctc(it.first);
    }
}

void top_sort()
{
    for(int i=1;i<=n;i++)
    {
        if(!sel[i])
        {
            dfs_sort_tot(i);
        }
    }
    reverse(ltot.begin(),ltot.end());
}


void gauss(int cp, int bad)
{
    int nr = c[cp].size();
    for(int i=0; i<2*nr; i++)
    {
        if(i==bad || i==bad + nr)
        {
            continue;
        }
        poz[i] = 0;
        for(int j=0; j<=2*nr; j++)
        {
            if(j==bad || j==bad + nr)
            {
                continue;
            }
            if(abs(a[i][j]) > limit)
            {
                poz[i] = j;
                break;
            }
        }
        for(int l=0; l<2*nr; l++)
        {
            if(i==l)
            {
                continue;
            }
            if(l==bad || l==bad + nr)
            {
                continue;
            }
            if(abs(a[l][poz[i]]) < limit)
            {
                continue;
            }
            double fact = 1.0 * a[l][poz[i]] / a[i][poz[i]];
            for(int j=0; j<=2*nr; j++)
            {
                if(j==bad || j==bad+nr)
                {
                    continue;
                }
                a[l][j] -= fact * a[i][j];
            }
        }
    }
    for(int i=0; i<2*nr; i++)
    {
        if(i==bad || i==bad + nr)
        {
            continue;
        }
        if(poz[i]<nr)
        {
            int nod = c[cp][poz[i]];
            dp[nod][0] = 1.0 * a[i][2 * nr] / a[i][poz[i]];
        }
        else
        {
            int nod = c[cp][poz[i] - nr];
            dp[nod][1] = 1.0 * a[i][2 * nr] / a[i][poz[i]];
        }
    }
}

void solve(int state)
{
    for(auto cp : l)
    {
        int nr = c[cp].size();
        int bad = - nr - 1;
        for(int i=0; i<c[cp].size(); i++)
        {
            int nod = c[cp][i];
            if(nod==1)
            {
                dp[nod][state] = 1;
                dp[nod][state ^ 1] = 0;
                bad = i;
                continue;
            }
            a[i][i] = -1;
            a[i + nr][i + nr] = -1;
            for(auto it : B[nod])
            {
                double p0 = it.second.first;
                double p1 = it.second.second;
                if(componenta[it.first]==cp)
                {
                    int poz = 0;
                    for(int j=0; j<c[cp].size(); j++)
                    {
                        if(c[cp][i]==it.first)
                        {
                            poz = j;
                            break;
                        }
                    }
                    a[i][poz] += (1.0 / grad[it.first]) * p0;
                    a[i][poz + nr] += (1.0 / grad[it.first]) * (1 - p1);
                    a[i + nr][poz] += (1.0 / grad[it.first]) * (1 - p0);
                    a[i + nr][poz + nr] += (1.0 / grad[it.first]) * p1;
                }
                else
                {
                    a[i][2 * nr] -= (dp[it.first][0] * p0 + dp[it.first][1] * (1 - p1)) * (1.0 / grad[it.first]);
                    a[i + nr][2 * nr] -= (dp[it.first][1] * p1 + dp[it.first][0] * (1 - p0)) * (1.0 / grad[it.first]);
                }
            }
        }
        gauss(cp,bad);
        for(int i=0; i<=2 * nr; i++)
        {
            for(int j=0; j<=2 * nr; j++)
            {
                a[i][j] = 0;
            }
        }
    }
    g<<fixed<<setprecision(7);
    g<<dp[n][state]<<'\n';
}

int main()
{
    f>>n>>m;
    for(int i=1; i<=m; i++)
    {
        int x,y;
        double p0,p1;
        f>>x>>y>>p0>>p1;
        ++x;
        ++y;
        G[x].push_back({y,{p0,p1}});
        B[y].push_back({x,{p0,p1}});
        e.push_back({x,y});
        ++grad[x];
    }
    ctc(1);
    memset(sel,0,sizeof(sel));
    top_sort();
    memset(sel,0,sizeof(sel));
    for(auto it : ltot)
    {
        if(sel[it])
        {
            continue;
        }
        ++comp;
        ctc(it);
    }
    for(auto it : e)
    {
        int x = it.first;
        int y = it.second;
        GC[componenta[x]].push_back(componenta[y]);
    }
    memset(sel,0,sizeof(sel));
    for(int i=1; i<=comp; i++)
    {
        if(!sel[i])
        {
            dfs_sort(i);
        }
    }
    reverse(l.begin(),l.end());
    solve(0);
    solve(1);
    return 0;
}