std::vector<std::list<int>> P( graph.size() );
std::vector<int> sigma( graph.size() );
std::vector<int> d( graph.size() );
std::vector<int> delta( graph.size() );

for ( auto i = a; i < b; ++i ) {
    if ( graph[i].first == -1 ) continue; // omit nonexistent nodes
    int s = graph[i].first;
    std::stack<int> S;

    for ( auto j = 0; j < graph.size(); ++j ) {
        int w = graph[j].first;
        if ( w == -1 ) continue; // omit nonexistent nodes
        P[w] = std::list<int>();
        sigma[w] = 0;
        d[w] = -1;
        delta[w] = 0;
    }

    sigma[s] = 1;
    d[s] = 0;
    std::queue<int> Q;
    Q.push( s );

    while ( !Q.empty() ) {
        int v = Q.front();
        Q.pop();
        S.push( v );

        for ( auto w : graph[v].second ) {
            if ( d[w] < 0 ) {
                Q.push( w );
                d[w] = d[v] + 1;
            }
            if ( d[w] == d[v] + 1 ) {
                sigma[w] += sigma[v];
                P[w].push_back(v);
            }
        }
    }

    while ( !S.empty() ) {
        int w = S.top();
        S.pop();

        for ( auto v : P[w] ) {
            delta[v] += ( sigma[v] / sigma[w] ) * ( 1 + delta[w] );
        }

        if ( w != s ) {
            cout << "bsp [" << a <<", " << b << "): Increasing BC[" << w << "] by delta[" << w << "] = " << delta[w] << endl;
            BC[w].second += delta[w];
        }
    }
}