digraph.h

Go to the documentation of this file.

// Copyright (C) 1997-2010 by Pawel Pilarczyk.
//
// This file is part of the Homology Library.  This library is free software;
// you can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation;
// either version 2 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this software; see the file "license.txt".  If not, write to the
// Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
// MA 02111-1307, USA.

// Started in January 2006. Last revision: May 29, 2010.


#ifndef _CHOMP_STRUCT_DIGRAPH_H_
#define _CHOMP_STRUCT_DIGRAPH_H_

#include <iostream>
#include <new>
#include <stack>
#include <queue>
#include <set>
#include <memory>
#include <vector>
#include <algorithm>

#include "chomp/struct/multitab.h"
#include "chomp/struct/hashsets.h"
#include "chomp/struct/flatmatr.h"
#include "chomp/struct/bitfield.h"
#include "chomp/struct/bitsets.h"
#include "chomp/struct/fibheap.h"
#include "chomp/system/textfile.h"
#include "chomp/system/timeused.h"

namespace chomp {
namespace homology {


// --------------------------------------------------
// ----------------- DUMMY ROUNDING -----------------
// --------------------------------------------------

template <class numType>
class dummyRounding
{
public:
        static inline numType add_down (const numType &x, const numType &y)
        { return x + y; }

        static inline numType add_up (const numType &x, const numType &y)
        { return x + y; }

        static inline numType sub_down (const numType &x, const numType &y)
        { return x - y; }

        static inline numType sub_up (const numType &x, const numType &y)
        { return x - y; }

        static inline numType mul_down (const numType &x, const numType &y)
        { return x * y; }

        static inline numType mul_up (const numType &x, const numType &y)
        { return x * y; }

        static inline numType div_down (const numType &x, const numType &y)
        { return x / y; }

        static inline numType div_down (const numType &x, int_t y)
        { return x / y; }

        static inline numType div_up (const numType &x, const numType &y)
        { return x / y; }

private:
}; /* class dummyRounding */


// --------------------------------------------------
// ------------------ DUMMY ARRAY -------------------
// --------------------------------------------------

class dummyArray
{
public:
        dummyArray () {n = 0;}

        int_t &operator [] (int_t) {return n;}

private:
        int_t n;

}; /* class dummyArray */


// --------------------------------------------------
// ----------------- DIRECTED GRAPH -----------------
// --------------------------------------------------

// #define DIGRAPH_DEBUG

template <class wType = int>
class diGraph
{
public:
        typedef wType weight_type;
        
        diGraph ();

        ~diGraph ();

        void swap (diGraph<wType> &g);

        void addVertex (void);

        void addEdge (int_t target);

        void addEdge (int_t target, const wType &weight);

        void setWeight (int_t vertex, int_t i, const wType &weight);

        void setWeight (int_t edge, const wType &weight);
        
        template <class Table>
        void setWeights (const Table &tab);

        void removeVertex (void);

        void removeVertex (int_t vertex, bool updateweights = false);

        int_t countVertices (void) const;

        int_t countEdges (void) const;

        int_t countEdges (int_t vertex) const;
        
        int_t getEdge (int_t vertex, int_t i) const;

        const wType &getWeight (int_t vertex, int_t i) const;

        const wType &getWeight (int_t edge) const;

        template <class Table>
        void getWeights (Table &tab) const;

        template <class Table>
        void writeEdges (Table &tab) const;

        template <class wType1>
        void transpose (diGraph<wType1> &result,
                bool copyweights = false) const;

        template <class Table, class wType1>
        void subgraph (diGraph<wType1> &result, const Table &tab,
                bool copyweights = false) const;

        template <class Table, class Color>
        void DFScolor (Table &tab, const Color &color,
                int_t vertex = 0) const;

        template <class Table, class Color>
        void DFScolorRecurrent (Table &tab, const Color &color,
                int_t vertex = 0) const;

        template <class Table, class Color>
        void DFScolorStack (Table &tab, const Color &color,
                int_t vertex = 0) const;

        template <class Table>
        void DFSfinishTime (Table &tab) const;

        template <class Table1, class Table2, class Table3>
        int_t DFSforest (const Table1 &ordered, Table2 &compVertices,
                Table3 &compEnds, bool nontrivial = false,
                diGraph<wType> *sccGraph = 0) const;

        int_t shortestPath (int_t source, int_t destination) const;

        int_t shortestLoop (int_t origin) const;

        template <class lenTable, class weightsType, class roundType>
        void Dijkstra (const roundType &rounding, int_t source,
                lenTable &len, weightsType &edgeWeights) const;

        template <class lenTable, class roundType>
        void Dijkstra (const roundType &rounding, int_t source,
                lenTable &len) const;

        template <class lenTable>
        void Dijkstra (int_t source, lenTable &len) const;

        template <class lenTable, class predTable, class roundType>
        bool BellmanFord (const roundType &rounding, int_t source,
                lenTable &len, wType *infinity, predTable pred) const;

        template <class lenTable, class predTable>
        bool BellmanFord (int_t source, lenTable &len, wType *infinity,
                predTable pred) const;

        template <class roundType>
        bool BellmanFord (const roundType &rounding, int_t source) const;

        bool BellmanFord (int_t source) const;

        wType Edmonds () const;

        wType EdmondsOld () const;

        template <class arrayType, class roundType>
        wType FloydWarshall (const roundType &rounding, arrayType &arr,
                bool setInfinity = true, bool ignoreNegLoop = false) const;

        template <class arrayType>
        wType FloydWarshall (arrayType &arr, bool setInfinity = true,
                bool ignoreNegLoop = false) const;

        template <class arrayType, class roundType>
        wType Johnson (const roundType &rounding, arrayType &arr,
                bool setInfinity = true, bool ignoreNegLoop = false) const;

        template <class arrayType>
        wType Johnson (arrayType &arr, bool setInfinity = true,
                bool ignoreNegLoop = false) const;

        template <class roundType>
        wType minPathWeight (const roundType &rounding,
                bool ignoreNegLoop = false, int sparseGraph = -1) const;

        wType minPathWeight (bool ignoreNegLoop = false,
                int sparseGraph = -1) const;

        wType minMeanCycleWeight (diGraph<wType> *transposed = 0) const;

        template <class roundType>
        wType minMeanCycleWeight (const roundType &rounding,
                diGraph<wType> *transposed) const;

        template <class arrayType, class roundType>
        wType minMeanPathWeight (const roundType &rounding,
                const arrayType &starting, int_t n) const;

        template <class arrayType>
        wType minMeanPathWeight (const arrayType &starting, int_t n) const;

        template <class wType1, class wType2>
        friend bool operator == (const diGraph<wType1> &g1,
                const diGraph<wType2> &g2);

        template <class outType>
        outType &show (outType &out, bool showWeights = false) const;

protected:
        int_t nVertices;

        multitable<int_t> edgeEnds;

        multitable<int_t> edges;

        multitable<wType> weights;

private:
        template <class Table>
        void DFSfinishTimeRecurrent (Table &tab, int_t vertex,
                int_t &counter) const;

        template <class Table>
        void DFSfinishTimeStack (Table &tab, int_t vertex,
                int_t &counter) const;

        template <class Table1, class Table2>
        bool DFSforestRecurrent (Table1 &tab, Table1 &ntab, int_t vertex,
                int_t treeNumber, int_t countTrees, Table2 &compVertices,
                int_t &curVertex, diGraph *sccGraph, int_t *sccEdgeAdded)
                const;

        template <class Table1, class Table2>
        bool DFSforestStack (Table1 &tab, Table1 &ntab, int_t vertex,
                int_t treeNumber, int_t countTrees, Table2 &compVertices,
                int_t &curVertex, diGraph *sccGraph, int_t *sccEdgeAdded)
                const;

        struct edgeTriple
        {
                int_t vert1;
                int_t vert2;
                
                wType weight;

                friend bool operator < (const edgeTriple &x,
                        const edgeTriple &y)
                {
                        return (x. weight < y. weight);
                }
        };

        class posWeight
        {
        public:
                posWeight ()
                {
                        value = 0;
                        return;
                }

                explicit posWeight (const wType &_value)
                {
                        if (_value < 0)
                                value = 0;
                        else
                                value = _value;
                        return;
                }

                void setInfinity ()
                {
                        value = -1;
                        return;
                }

                bool isInfinity () const
                {
                        return (value < 0);
                }

                const wType &getValue () const
                {
                        return value;
                }

                friend bool operator < (const posWeight &x,
                        const posWeight &y)
                {
                        if (y. isInfinity ())
                                return !(x. isInfinity ());
                        else if (x. isInfinity ())
                                return false;
                        return (x. value < y. value);
                }

                friend std::ostream &operator << (std::ostream &out,
                        const posWeight &x)
                {
                        if (x. isInfinity ())
                                out << "+oo";
                        else
                                out << x. getValue ();
                        return out;
                }

        private:
                wType value;
        };

}; /* class diGraph */

// --------------------------------------------------

template <class wType>
inline diGraph<wType>::diGraph (): nVertices (0),
        edgeEnds (1024), edges (4096), weights (4096)
{
        return;
} /* diGraph::diGraph */

template <class wType>
inline diGraph<wType>::~diGraph ()
{
        return;
} /* diGraph::~diGraph */

// --------------------------------------------------

template <class wType1, class wType2>
inline bool operator == (const diGraph<wType1> &g1,
        const diGraph<wType2> &g2)
{
        if (g1. nVertices != g2. nVertices)
                return false;
        if (!g1. nVertices)
                return true;
        for (int_t i = 0; i < g1. nVertices; ++ i)
        {
                if (g1. edgeEnds [i] != g2. edgeEnds [i])
                        return false;
        }
        int_t nEdges = g1. edgeEnds [g1. nVertices - 1];
        for (int_t i = 0; i < nEdges; ++ i)
        {
                if (g1. edges [i] != g2. edges [i])
                        return false;
        }
        return true;
} /* operator == */

template <class wType1, class wType2>
inline bool operator != (const diGraph<wType1> &g1,
        const diGraph<wType2> &g2)
{
        return !(g1 == g2);
} /* operator != */

// --------------------------------------------------

template <class wType>
inline void diGraph<wType>::swap (diGraph<wType> &g)
{
        my_swap (nVertices, g. nVertices);
        edgeEnds. swap (g. edgeEnds);
        edges. swap (g. edges);
        weights. swap (g. weights);
        return;
} /* diGraph::swap */

template <class wType>
inline void diGraph<wType>::addVertex (void)
{
        edgeEnds [nVertices] = nVertices ? edgeEnds [nVertices - 1] :
                static_cast<int_t> (0);
        ++ nVertices;
        return;
} /* diGraph::addVertex */

template <class wType>
inline void diGraph<wType>::addEdge (int_t target)
{
        if (!nVertices)
                throw "Trying to add an edge to an empty graph.";
//      if (target >= nVertices)
//              throw "Trying to add an edge to a nonexistent vertex.";
        if (target < 0)
                throw "Trying to add an edge to a negative vertex.";
        int_t &offset = edgeEnds [nVertices - 1];
        if (offset + 1 <= 0)
                throw "Too many edges in a diGraph (limit = 2,147,483,647).";
        edges [offset] = target;
        ++ offset;
        return;
} /* diGraph::addEdge */

template <class wType>
inline void diGraph<wType>::addEdge (int_t target, const wType &weight)
{
        if (!nVertices)
                throw "Trying to add an edge to an empty graph.";
//      if (target >= nVertices)
//              throw "Trying to add an edge to a nonexistent vertex.";
        if (target < 0)
                throw "Trying to add an edge to a negative vertex.";
        int_t &offset = edgeEnds [nVertices - 1];
        if (offset + 1 <= 0)
                throw "Too many edges in a diGraph (limit = 2,147,483,647).";
        edges [offset] = target;
        weights [offset] = weight;
        ++ offset;
        return;
} /* diGraph::addEdge */

template <class wType>
inline void diGraph<wType>::setWeight (int_t vertex, int_t i,
        const wType &weight)
{
        if (!vertex)
                weights [i] = weight;
        else
                weights [edgeEnds [vertex - 1] + i] = weight;
        return;
} /* diGraph::setWeight */

template <class wType>
inline void diGraph<wType>::setWeight (int_t edge, const wType &weight)
{
        weights [edge] = weight;
        return;
} /* diGraph::setWeight */

template <class wType> template <class Table>
inline void diGraph<wType>::setWeights (const Table &tab)
{
        if (!nVertices)
                return;
        int_t nEdges = edgeEnds [nVertices - 1];
        for (int_t i = 0; i < nEdges; ++ i)
                weights [i] = tab [i];
        return;
} /* diGraph::setWeights */

template <class wType>
inline void diGraph<wType>::removeVertex (void)
{
        if (!nVertices)
                throw "Trying to remove a vertex from the empty graph.";
        -- nVertices;
        return;
} /* diGraph::removeVertex */

template <class wType>
inline void diGraph<wType>::removeVertex (int_t vertex, bool updateweights)
{
        // make sure that the vertex number is within the scope
        if ((vertex < 0) || (vertex >= nVertices))
                throw "Trying to remove a vertex that is not in the graph.";

        // remove edges that begin or end at the given vertex
        int_t curEdge = 0;
        int_t newEdge = 0;
        int_t skipped = 0;
        for (int_t v = 0; v < nVertices; ++ v)
        {
                // skip the edges that begin at the given vertex
                if (!skipped && (v == vertex))
                {
                        curEdge = edgeEnds [v];
                        ++ skipped;
                        continue;
                }

                // skip the edges that point to the given vertex
                int_t maxEdge = edgeEnds [v];
                for (; curEdge < maxEdge; ++ curEdge)
                {
                        if (edges [curEdge] == vertex)
                                continue;
                        int_t target = edges [curEdge];
                        edges [newEdge] = (target < vertex) ? target :
                                (target - 1);
                        if (updateweights)
                                weights [newEdge] = weights [curEdge];
                        ++ newEdge;
                }
                edgeEnds [v - skipped] = newEdge;
        }

        // decrease the number of vertices
        nVertices -= skipped;

        return;
} /* diGraph::removeVertex */

template <class wType>
inline int_t diGraph<wType>::countVertices (void) const
{
        return nVertices;
} /* diGraph::countVertices */

template <class wType>
inline int_t diGraph<wType>::countEdges (void) const
{
        if (!nVertices)
                return 0;
        else
                return edgeEnds [nVertices - 1];
} /* diGraph::countEdges */

template <class wType>
inline int_t diGraph<wType>::countEdges (int_t vertex) const
{
        if (!vertex)
                return edgeEnds [0];
        else
                return edgeEnds [vertex] - edgeEnds [vertex - 1];
} /* diGraph::countEdges */

template <class wType>
inline int_t diGraph<wType>::getEdge (int_t vertex, int_t i) const
{
        if (!vertex)
                return edges [i];
        else
                return edges [edgeEnds [vertex - 1] + i];
} /* diGraph::getEdge */

template <class wType>
inline const wType &diGraph<wType>::getWeight (int_t vertex, int_t i) const
{
        if (!vertex)
                return weights [i];
        else
                return weights [edgeEnds [vertex - 1] + i];
} /* diGraph::getWeight */

template <class wType>
inline const wType &diGraph<wType>::getWeight (int_t edge) const
{
        return weights [edge];
} /* diGraph::getWeight */

template <class wType> template <class Table>
inline void diGraph<wType>::getWeights (Table &tab) const
{
        if (!nVertices)
                return;
        int_t nEdges = edgeEnds [nVertices - 1];
        for (int_t i = 0; i < nEdges; ++ i)
                tab [i] = weights [i];
        return;
} /* diGraph::getWeights */

template <class wType> template <class Table>
inline void diGraph<wType>::writeEdges (Table &tab) const
{
        int_t curEdge = 0;
        for (int_t curVertex = 0; curVertex < nVertices; ++ curVertex)
        {
                int_t maxEdge = edgeEnds [curVertex];
                while (curEdge < maxEdge)
                {
                        tab [curEdge] [0] = curVertex;
                        tab [curEdge] [1] = edges [curEdge];
                        ++ curEdge;
                }
        }
        return;
} /* diGraph::writeEdges */

template <class wType> template <class wType1>
inline void diGraph<wType>::transpose (diGraph<wType1> &result,
        bool copyweights) const
{
        // prepare the resulting graph
        result. nVertices = nVertices;
        if (!nVertices)
                return;

        // compute the initial offsets for the edge numbers
        for (int_t i = 0; i < nVertices; ++ i)
                result. edgeEnds [i] = 0;
        int_t nEdges = edgeEnds [nVertices - 1];
        for (int_t i = 0; i < nEdges; ++ i)
        {
                if (edges [i] < nVertices - 1)
                        ++ result. edgeEnds [edges [i] + 1];
        }
        for (int_t i = 2; i < nVertices; ++ i)
                result. edgeEnds [i] += result. edgeEnds [i - 1];

        // compute the reversed edges
        int_t curEdge = 0;
        for (int_t i = 0; i < nVertices; ++ i)
        {
                for (; curEdge < edgeEnds [i]; ++ curEdge)
                {
                        int_t j = edges [curEdge];
                        int_t &offset = result. edgeEnds [j];
                        result. edges [offset] = i;
                        if (copyweights)
                                result. weights [offset] = weights [curEdge];
                        ++ offset;
                }
        }
        return;
} /* diGraph::transpose */

template <class wType> template <class Table, class wType1>
inline void diGraph<wType>::subgraph (diGraph<wType1> &g,
        const Table &tab, bool copyweights) const
{
        // compute the new numbers of vertices that remain in the graph
        int_t *numbers = new int_t [nVertices];
        int_t curNumber = 0;
        for (int_t i = 0; i < nVertices; ++ i)
        {
                if (tab [i])
                        numbers [i] = curNumber ++;
                else
                        numbers [i] = -1;
        }

        // copy the edges from the old graph to the new one
        for (int_t i = 0; i < nVertices; ++ i)
        {
                if (numbers [i] < 0)
                        continue;
                g. addVertex ();
                int_t firstEdge = i ? edgeEnds [i - 1] :
                        static_cast<int_t> (0);
                int_t endEdge = edgeEnds [i];
                for (int_t j = firstEdge; j < endEdge; ++ j)
                {
                        int_t edgeEnd = edges [j];
                        if (numbers [edgeEnd] >= 0)
                        {
                                if (copyweights)
                                        g. addEdge (numbers [edgeEnd],
                                                weights [j]);
                                else
                                        g. addEdge (numbers [edgeEnd]);
                        }
                }
        }

        // clean up memory and exit
        delete [] numbers;
        return;
} /* diGraph::subgraph */

// --------------------------------------------------

template <class wType> template <class Table, class Color>
inline void diGraph<wType>::DFScolorRecurrent (Table &tab,
        const Color &color, int_t vertex) const
{
        tab [vertex] = color;
        int_t maxEdge = edgeEnds [vertex];
        for (int_t i = vertex ? edgeEnds [vertex - 1] :
                static_cast<int_t> (0); i < maxEdge; ++ i)
        {
                int_t next = edges [i];
                if (tab [next] != color)
                        DFScolorRecurrent (tab, color, next);
        }
        return;
} /* diGraph::DFScolorRecurrent */

template <class wType> template <class Table, class Color>
inline void diGraph<wType>::DFScolorStack (Table &tab, const Color &color,
        int_t vertex) const
{
        // prepare stacks for the recursion
        std::stack<int_t> s_vertex;
        std::stack<int_t> s_edge;
        std::stack<int_t> s_maxedge;

        // mark the current vertex as visited
        tab [vertex] = color;

        // determine the edges to be visited
        int_t edge = vertex ? edgeEnds [vertex - 1] :
                static_cast<int_t> (0);
        int_t maxedge = edgeEnds [vertex];

        while (1)
        {
                // return to the previous recursion level
                // if all the edges have been followed
                if (edge >= maxedge)
                {
                        // return if this is the initial recursion level
                        if (s_vertex. empty ())
                                return;

                        // restore the variables from the previous level
                        vertex = s_vertex. top ();
                        s_vertex. pop ();
                        edge = s_edge. top ();
                        s_edge. pop ();
                        maxedge = s_maxedge. top ();
                        s_maxedge. pop ();
                        continue;
                }

                // go to the deeper recursion level if possible
                int_t next = edges [edge ++];
                if (tab [next] != color)
                {
                        // store the previous variables at the stacks
                        s_vertex. push (vertex);
                        s_edge. push (edge);
                        s_maxedge. push (maxedge);

                        // set the new vertex
                        vertex = next;
                        
                        // mark the new vertex as visited
                        tab [vertex] = color;
                        
                        // determine the edges to be visited
                        edge = vertex ? edgeEnds [vertex - 1] :
                                static_cast<int_t> (0);
                        maxedge = edgeEnds [vertex];
                }
        }
} /* diGraph::DFScolorStack */

template <class wType> template <class Table, class Color>
inline void diGraph<wType>::DFScolor (Table &tab, const Color &color,
        int_t vertex) const
{
        DFScolorStack (tab, color, vertex);
        return;
} /* diGraph::DFScolor */

// --------------------------------------------------

template <class wType> template <class Table>
inline void diGraph<wType>::DFSfinishTimeRecurrent (Table &tab,
        int_t vertex, int_t &counter) const
{
        // mark the current vertex as visited
        tab [vertex] = -1;

        // call DFS for the other vertices
        for (int_t edge = vertex ? edgeEnds [vertex - 1] :
                static_cast<int_t> (0);
                edge < edgeEnds [vertex]; ++ edge)
        {
                int_t next = edges [edge];
                if (!tab [next])
                        DFSfinishTimeRecurrent (tab, next, counter);
        }

        // record the finishing time for the current vertex and return
        tab [vertex] = ++ counter;
        return;
} /* diGraph::DFSfinishTimeRecurrent */

template <class wType> template <class Table>
inline void diGraph<wType>::DFSfinishTimeStack (Table &tab, int_t vertex,
        int_t &counter) const
{
        // prepare stacks for the recursion
        std::stack<int_t> s_vertex;
        std::stack<int_t> s_edge;
        std::stack<int_t> s_maxedge;

        // mark the current vertex as visited
        tab [vertex] = -1;

        // determine the edges to be visited
        int_t edge = vertex ? edgeEnds [vertex - 1] :
                static_cast<int_t> (0);
        int_t maxedge = edgeEnds [vertex];

        while (1)
        {
                // return to the previous recursion level
                // if all the edges have been followed
                if (edge >= maxedge)
                {
                        // record the finishing time
                        // for the current vertex
                        tab [vertex] = ++ counter;

                        // return if this is the initial recursion level
                        if (s_vertex. empty ())
                                return;

                        // restore the variables from the previous level
                        vertex = s_vertex. top ();
                        s_vertex. pop ();
                        edge = s_edge. top ();
                        s_edge. pop ();
                        maxedge = s_maxedge. top ();
                        s_maxedge. pop ();
                        continue;
                }

                // go to the deeper recursion level if possible
                int_t next = edges [edge ++];
                if (!tab [next])
                {
                        // store the previous variables at the stacks
                        s_vertex. push (vertex);
                        s_edge. push (edge);
                        s_maxedge. push (maxedge);

                        // set the new vertex
                        vertex = next;
                        
                        // mark the new vertex as visited
                        tab [vertex] = -1;
                        
                        // determine the edges to be visited
                        edge = vertex ? edgeEnds [vertex - 1] :
                                static_cast<int_t> (0);
                        maxedge = edgeEnds [vertex];
                }
        }

        return;
} /* diGraph::DFSfinishTimeStack */

template <class wType> template <class Table>
inline void diGraph<wType>::DFSfinishTime (Table &tab) const
{
        // initialize the table and the counter
        for (int_t i = 0; i < nVertices; ++ i)
                tab [i] = 0;
        int_t counter = 0;

        // compute the finishing time for each tree in the DFS forest
        for (int_t i = 0; i < nVertices; ++ i)
        {
                if (!tab [i])
                        DFSfinishTimeStack (tab, i, counter);
        }

        #ifdef DIGRAPH_DEBUG
        int_t *tabdebug = new int_t [nVertices];
        for (int_t i = 0; i < nVertices; ++ i)
                tabdebug [i] = 0;
        int_t counterdebug = 0;
        for (int_t i = 0; i < nVertices; ++ i)
                if (!tabdebug [i])
                        DFSfinishTimeRecurrent (tabdebug, i, counterdebug);
        if (counter != counterdebug)
                throw "DFSfinishTime: Wrong counter.";
        for (int_t i = 0; i < nVertices; ++ i)
        {
                if (tab [i] != tabdebug [i])
                {
                        sbug << "\nDFSfinishTime error: tabRec [" << i <<
                                "] = " << tab [i] << ", tabStack [" << i <<
                                "] = " << tabdebug [i] << ".\n";
                        throw "DFSfinishTime: Wrong numbers.";
                }
        }
        sbug << "DEBUG: DFSfinishTime OK. ";
        #endif // DIGRAPH_DEBUG
        return;
} /* diGraph::DFSfinishTime */

// --------------------------------------------------

template <class wType> template <class Table1, class Table2>
inline bool diGraph<wType>::DFSforestRecurrent (Table1 &tab, Table1 &ntab,
        int_t vertex, int_t treeNumber, int_t countTrees,
        Table2 &compVertices, int_t &curVertex,
        diGraph<wType> *sccGraph, int_t *sccEdgeAdded) const
{
        // add the vertex to the tree
        compVertices [curVertex ++] = vertex;

        // mark the vertex as belonging to the current tree
        tab [vertex] = treeNumber;
//      if (sccGraph)
//              ntab [treeNumber - 1] = countTrees;

        // build the tree recursively or record connections
        bool loop = false;
        for (int_t edge = vertex ? edgeEnds [vertex - 1] :
                static_cast<int_t> (0);
                edge < edgeEnds [vertex]; ++ edge)
        {
                int_t next = edges [edge];
                if (!tab [next])
                        loop |= DFSforestRecurrent (tab, ntab, next,
                                treeNumber, countTrees, compVertices,
                                curVertex, sccGraph, sccEdgeAdded);
                else if (tab [next] == treeNumber)
                {
                        if (sccGraph)
                        {
                                int_t target = ntab [treeNumber - 1];
                                if (sccEdgeAdded [target] != treeNumber)
                                {
                                        sccGraph -> addEdge (target);
                                        sccEdgeAdded [target] = treeNumber;
                                }
                        }
                        loop = true;
                }
                else if (sccGraph)
                {
                        int_t target = ntab [tab [next] - 1];
                        if ((target >= 0) &&
                                (sccEdgeAdded [target] != treeNumber))
                        {
                                sccGraph -> addEdge (target);
                                sccEdgeAdded [target] = treeNumber;
                        }
                }
        }

        return loop;
} /* diGraph::DFSforestRecurrent */

template <class wType> template <class Table1, class Table2>
inline bool diGraph<wType>::DFSforestStack (Table1 &tab, Table1 &ntab,
        int_t vertex, int_t treeNumber, int_t countTrees,
        Table2 &compVertices, int_t &curVertex,
        diGraph<wType> *sccGraph, int_t *sccEdgeAdded) const
{
        // prepare stacks for the recursion
        std::stack<int_t> s_vertex;
        std::stack<int_t> s_edge;
        std::stack<int_t> s_maxedge;
        std::stack<bool> s_loop;

        // add the vertex to the tree
        compVertices [curVertex ++] = vertex;

        // mark the vertex as belonging to the current tree
        tab [vertex] = treeNumber;
//      if (sccGraph)
//              ntab [vertex] = countTrees;

        // build the tree recursively or record connections
        bool loop = false;
        int_t edge = vertex ? edgeEnds [vertex - 1] :
                static_cast<int_t> (0);
        int_t maxedge = edgeEnds [vertex];
        while (1)
        {
                // return to the previous recursion level
                // if all the edges have been followed
                if (edge >= maxedge)
                {
                        // return if this is the initial recursion level
                        if (s_vertex. empty ())
                                return loop;

                        // restore the variables from the previous level
                        vertex = s_vertex. top ();
                        s_vertex. pop ();
                        edge = s_edge. top ();
                        s_edge. pop ();
                        maxedge = s_maxedge. top ();
                        s_maxedge. pop ();
                        loop |= s_loop. top ();
                        s_loop. pop ();
                        continue;
                }

                // go to the deeper recursion level if possible
                int_t next = edges [edge ++];
                if (!tab [next])
                {
                        // store the previous variables at the stacks
                        s_vertex. push (vertex);
                        s_edge. push (edge);
                        s_maxedge. push (maxedge);
                        s_loop. push (loop);

                        // set the new vertex
                        vertex = next;
                        
                        // add the vertex to the tree
                        compVertices [curVertex ++] = vertex;

                        // mark the vertex as belonging to the current tree
                        tab [vertex] = treeNumber;

                        // determine the edges to be visited
                        loop = false;
                        edge = vertex ? edgeEnds [vertex - 1] :
                                static_cast<int_t> (0);
                        maxedge = edgeEnds [vertex];
                }
                else if (tab [next] == treeNumber)
                {
                        if (sccGraph)
                        {
                                int_t target = ntab [treeNumber - 1];
                                if (sccEdgeAdded [target] != treeNumber)
                                {
                                        sccGraph -> addEdge (target);
                                        sccEdgeAdded [target] = treeNumber;
                                }
                        }
                        loop = true;
                }
                else if (sccGraph)
                {
                        int_t target = ntab [tab [next] - 1];
                        if ((target >= 0) &&
                                (sccEdgeAdded [target] != treeNumber))
                        {
                                sccGraph -> addEdge (target);
                                sccEdgeAdded [target] = treeNumber;
                        }
                }
        }

        return loop;
} /* diGraph::DFSforestStack */

template <class wType> template <class Table1, class Table2, class Table3>
inline int_t diGraph<wType>::DFSforest (const Table1 &ordered,
        Table2 &compVertices, Table3 &compEnds, bool nontrivial,
        diGraph<wType> *sccGraph) const
{
        // prepare a table to record the numbers of DFS trees
        // to which the vertices belong (the tree numbers begin with 1)
        int_t *tab = new int_t [nVertices];
        for (int_t i = 0; i < nVertices; ++ i)
                tab [i] = 0;

        // prepare a table to record the numbers of nontrivial trees
        // that correspond to the trees in 'tab' (these numbers begin with 0)
        int_t *ntab = new int_t [nVertices];

        // prepare a table to record the numbers of edges already in the
        // scc graph; "sccEdgeAdded [n] = m" indicates that the edge
        // m -> n has been added to the scc graph
        int_t *sccEdgeAdded = sccGraph ? new int_t [nVertices] :
                static_cast<int_t *> (0);
        if (sccGraph)
        {
                for (int_t n = 0; n < nVertices; ++ n)
                        sccEdgeAdded [n] = -1;
        }

        // prepare the official DFS tree number
        int_t treeNumber = 0;

        // prepare the data for keeping the nontrivial trees information
        int_t countTrees = 0;
        int_t curVertex = 0;

        // compute the DFS trees and connections between them
        for (int_t i = 0; i < nVertices; ++ i)
        {
                // take the next vertex
                int_t vertex = ordered [i];

                // if the vertex already belongs to some tree, skip it
                if (tab [vertex])
                        continue;

                // add a vertex corresponding to the component
                if (sccGraph)
                        sccGraph -> addVertex ();

                // remember the previous vertex number
                int_t prevVertex = curVertex;

                // mark the entire component and record connections graph
                if (sccGraph)
                        ntab [treeNumber] = countTrees;
                ++ treeNumber;
                bool loop = DFSforestStack (tab, ntab, vertex,
                        treeNumber, countTrees, compVertices,
                        curVertex, sccGraph, sccEdgeAdded);

                // update the index bound for the vertex list
                compEnds [countTrees ++] = curVertex;

                // remove the component if it is trivial
                if (nontrivial && !loop)
                {
                        -- countTrees;
                        curVertex = prevVertex;
                        if (sccGraph)
                        {
                                ntab [treeNumber - 1] = -1;
                                sccGraph -> removeVertex ();
                        }
                }
        }

        #ifdef DIGRAPH_DEBUG
        diGraph<wType> *sccGraphdebug = 0;
        if (sccGraph)
                sccGraphdebug = new diGraph<wType>;
        // prepare a table to record the numbers of DFS trees
        // to which the vertices belong (the tree numbers begin with 1)
        int_t *tabdebug = new int_t [nVertices];
        for (int_t i = 0; i < nVertices; ++ i)
                tabdebug [i] = 0;

        // prepare a table to record the numbers of nontrivial trees
        // to which the vertices belong (the tree numbers begin with 0)
        int_t *ntabdebug = new int_t [nVertices];

        // prepare a table to record the numbers of vertices from which
        // edges were added to the scc graph
        int_t *sccEdgeAddeddebug = sccGraph ? new int_t [nVertices] :
                static_cast<int_t> (0);
        if (sccGraph)
        {
                for (int_t n = 0; n < nVertices; ++ n)
                        sccEdgeAddeddebug [n] = -1;
        }
        // prepare the official DFS tree number
        int_t treeNumberdebug = 0;

        // prepare the data for keeping the nontrivial trees information
        int_t countTreesdebug = 0;
        int_t curVertexdebug = 0;

        int_t *compVerticesdebug = new int_t [nVertices];
        int_t *compEndsdebug = new int_t [nVertices];
        
        // compute the DFS trees and connections between them
        for (int_t i = 0; i < nVertices; ++ i)
        {
                // take the next vertex
                int_t vertex = ordered [i];

                // if the vertex already belongs to some tree, skip it
                if (tabdebug [vertex])
                        continue;

                // add a vertex corresponding to the component
                if (sccGraphdebug)
                        sccGraphdebug -> addVertex ();

                // remember the previous vertex number
                int_t prevVertex = curVertexdebug;

                // mark the entire component and record connections graph
                if (sccGraphdebug)
                        ntabdebug [treeNumberdebug] = countTreesdebug;
                ++ treeNumberdebug;
                bool loop = DFSforestRecurrent (tabdebug, ntabdebug, vertex,
                        treeNumberdebug, countTreesdebug, compVerticesdebug,
                        curVertexdebug, sccGraphdebug, sccEdgeAddeddebug);

                // update the index bound for the vertex list
                compEndsdebug [countTreesdebug ++] = curVertexdebug;

                // remove the component if it is trivial
                if (nontrivial && !loop)
                {
                        -- countTreesdebug;
                        curVertexdebug = prevVertex;
                        if (sccGraphdebug)
                        {
                                ntabdebug [treeNumberdebug - 1] = -1;
                                sccGraphdebug -> removeVertex ();
                        }
                }
        }
        if (countTrees != countTreesdebug)
                throw "DFSforest: Wrong countTrees.";
        for (int_t i = 0; i < countTrees; ++ i)
                if (compEnds [i] != compEndsdebug [i])
                        throw "DFSforest: Wrong compEnds.";
        for (int_t i = 0; i < compEndsdebug [countTrees - 1]; ++ i)
                if (compVertices [i] != compVerticesdebug [i])
                        throw "DFSforest: Wrong vertices.";
        if (curVertex != curVertexdebug)
                throw "DFSforest: Wrong curVertex.";
        for (int_t i = 0; i < nVertices; ++ i)
                if (tab [i] != tabdebug [i])
                        throw "DFSforest: Wrong tab.";
        if (sccGraph)
        {
                for (int_t i = 0; i < nVertices; ++ i)
                        if (ntab [i] != ntabdebug [i])
                                throw "DFSforest: Wrong ntab.";
                if (*sccGraph != *sccGraphdebug)
                        throw "DFSforest: Wrong graph.";
        }
        if (sccEdgeAdded)
        {
                for (int_t i = 0; i < nVertices; ++ i)
                        if (sccEdgeAdded [i] != sccEdgeAddeddebug [i])
                                throw "DFSforest: Wrong sccEdgeAdded.";
        }
        if (treeNumber != treeNumberdebug)
                throw "DFSforest: Wrong treeNumber.";
        sbug << "DEBUG: DFSforest OK. ";
        if (!sccGraph)
                sbug << "(Graphs not compared.) ";
        delete [] compVerticesdebug;
        delete [] compEndsdebug;
        if (sccGraphdebug)
                delete sccGraphdebug;
        delete [] ntabdebug;
        delete [] tabdebug;
        if (sccEdgeAddeddebug)
                delete [] sccEdgeAddeddebug;
        #endif // DIGRAPH_DEBUG

        if (sccEdgeAdded)
                delete [] sccEdgeAdded;
        delete [] ntab;
        delete [] tab;
        return countTrees;
} /* diGraph::DFSforest */

template <class wType>
inline int_t diGraph<wType>::shortestPath (int_t source, int_t destination)
        const
{
        // make sure that the given vertex is present in the graph
        if ((source < 0) || (source >= nVertices) ||
                (destination < 0) || (destination >= nVertices))
        {
                throw "diGraph - shortest path: Wrong vertex number.";
        }

        // prepare an array of bits to store the information
        // on whether the given vertices have been visited or not
        BitField visited;
        visited. allocate (nVertices);
        visited. clearall (nVertices);

        // prepare queues for the BFS algorithm
        std::queue<int_t> q_vertex;
        std::queue<int_t> q_depth;

        // set the initial vertex
        int_t vertex = source;
        int_t depth = 0;

        while (1)
        {
                // mark the current vertex as visited
                visited. set (vertex);

                // determine the depth of the vertices that will be analyzed
                ++ depth;

                // determine the edges to be checked
                int_t firstedge = vertex ? edgeEnds [vertex - 1] :
                        static_cast<int_t> (0);
                int_t maxedge = edgeEnds [vertex];

                // put all the unvisited destination vertices on the stack
                for (int_t edge = firstedge; edge < maxedge; ++ edge)
                {
                        // determine the vertex pointed to by this edge
                        int_t next = edges [edge];

                        // if this is the destination vertex then return
                        // the shortest path length; note: this vertex
                        // might be visited if checking a loop, so it is
                        // important to check the destination first
                        if (next == destination)
                        {
                                visited. free ();
                                return depth;
                        }

                        // if the vertex was already visited then skip it
                        if (visited. test (next))
                                continue;

                        // add the vertex to the queue
                        q_vertex. push (next);
                        q_depth. push (depth);
                }

                // if there are no vertices whose neighbors are to be
                // analyzed and the destination vertex was not found
                // then there is no path to that vertex
                if (q_vertex. empty ())
                {
                        visited. free ();
                        return 0;
                }

                // pick up a vertex stored in the queue
                vertex = q_vertex. front ();
                q_vertex. pop ();
                depth = q_depth. front ();
                q_depth. pop ();
        }
} /* diGraph::shortestPath */

template <class wType>
inline int_t diGraph<wType>::shortestLoop (int_t origin) const
{
        return shortestPath (origin, origin);
} /* diGraph::shortestLoop */

template <class wType>
template <class lenTable, class weightsType, class roundType>
inline void diGraph<wType>::Dijkstra (const roundType &rounding,
        int_t source, lenTable &len, weightsType &edgeWeights) const
{
        // use the Fibonacci heap as a priority queue
        FibonacciHeap<posWeight> Q (nVertices);

        // add the vertices to the heap with the initial shortest path
        // lengths: 0 to the source, plus infinity to all the others
        for (int_t v = 0; v < nVertices; ++ v)
        {
                posWeight w (0);
                if (v != source)
                        w. setInfinity ();
                int_t number = Q. Insert (w);
                if (number != v)
                {
                        throw "Wrong implementation of Fibonacci heap "
                                "for this version of Dijkstra.";
                }
        }

        // pick up vertices from the priority queue, record the length
        // of the shortest path to them, and modify the remaining paths
        for (int_t i = 0; i < nVertices; ++ i)
        {
                // extract the minimal vertex from the queue
                int_t minVertex = Q. Minimum ();
                posWeight minWeight = Q. ExtractMinimum ();

                if (minWeight. isInfinity ())
                {
                        len [minVertex] = -1;
                        continue;
                }
                wType minValue = minWeight. getValue ();
                len [minVertex] = minValue;

                // go through all the edges emanating from this vertex
                // and update the path lengths for the target vertices
                int_t edge = minVertex ? edgeEnds [minVertex - 1] :
                        static_cast<int_t> (0);
                int_t maxEdge = edgeEnds [minVertex];
                for (; edge < maxEdge; ++ edge)
                {
                        // determine the vertex at the other end of the edge
                        int_t nextVertex = edges [edge];

                        // if the path that runs through the extracted
                        // vertex is shorter, then make a correction
                        const posWeight &nextWeight = Q. Value (nextVertex);
                        wType newWeight = rounding. add_down (minValue,
                                edgeWeights [edge]);
                        if (newWeight < 0)
                                newWeight = 0;
                        if (nextWeight. isInfinity () ||
                                (newWeight < nextWeight. getValue ()))
                        {
                                Q. DecreaseKey (nextVertex,
                                        posWeight (newWeight));
                        }
                }
        }
        return;
} /* diGraph::Dijkstra */

template <class wType>
template <class lenTable, class roundType>
inline void diGraph<wType>::Dijkstra (const roundType &rounding,
        int_t source, lenTable &len) const
{
        this -> Dijkstra (rounding, source, len, this -> weights);
        return;
} /* diGraph::Dijkstra */

template <class wType>
template <class lenTable>
inline void diGraph<wType>::Dijkstra (int_t source, lenTable &len) const
{
        const dummyRounding<wType> rounding = dummyRounding<wType> ();
        this -> Dijkstra (rounding, source, len);
        return;
} /* diGraph::Dijkstra */

template <class wType> template <class lenTable, class predTable,
        class roundType>
inline bool diGraph<wType>::BellmanFord (const roundType &rounding,
        int_t source, lenTable &len, wType *infinity, predTable pred) const
{
        // make sure the source vertex number is correct
        if ((source < 0) || (source >= nVertices))
                throw "Bellman-Ford: Wrong source vertex number.";

        // prepare marks to indicate finite values (not "infinity")
        BitField finite;
        finite. allocate (nVertices);
        finite. clearall (nVertices);
        finite. set (source);
        len [source] = 0;

        // count the negative vertices
        int_t countNegative = 0;

        // set the initial predecessors
        for (int_t i = 0; i < nVertices; ++ i)
                pred [i] = -1;

        // update the lenghts of the paths repeatedly (max nVertices times)
        bool noNegativeLoop = false;
        int_t counter = 0;
        for (; counter <= nVertices; ++ counter)
        {
                bool modified = false;
                int_t curEdge = 0;
                for (int_t vertex = 0; vertex < nVertices; ++ vertex)
                {
                        int_t maxEdge = edgeEnds [vertex];
                        if (!finite. test (vertex))
                        {
                                curEdge = maxEdge;
                                continue;
                        }
                        for (; curEdge < maxEdge; ++ curEdge)
                        {
                                int_t next = edges [curEdge];
                                wType newlen = rounding. add_down
                                        (len [vertex], weights [curEdge]);
                                if (!finite. test (next))
                                {
                                        finite. set (next);
                                        modified = true;
                                        len [next] = newlen;
                                        pred [next] = vertex;
                                        if (newlen < 0)
                                                ++ countNegative;
                                }
                                else if (newlen < len [next])
                                {
                                        modified = true;
                                        if (!(len [next] < 0) &&
                                                (newlen < 0))
                                        {
                                                ++ countNegative;
                                        }
                                        len [next] = newlen;
                                        pred [next] = vertex;
                                }
                        }
                }
                if (countNegative == nVertices)
                {
                        noNegativeLoop = false;
                        ++ counter;
                        break;
                }
                if (!modified)
                {
                        noNegativeLoop = true;
                        ++ counter;
                        break;
                }
        }

        // show a message on how many loops have been done
        if (false && chomp::homology::sbug. show)
        {
                chomp::homology::sbug << "Bellman-Ford: " <<
                        counter << ((counter > 1) ? " loops (" : " loop (") <<
                        nVertices << " vertices, " << countNegative <<
                        " negative). " <<
                        (noNegativeLoop ? "No negative loops.\n" :
                        "A negative loop found.\n");
        }

        // compute the value for the infinity and set the undefined distances
        if (infinity && noNegativeLoop)
        {
                wType infty (0);
                bool first = true;
                for (int_t i = 0; i < nVertices; ++ i)
                {
                        if (!finite. test (i))
                                continue;
                        if (first)
                        {
                                infty = len [i];
                                first = false;
                        }
                        else if (infty < len [i])
                        {
                                infty = len [i];
                        }
                }
                infty = infty + 1;
                for (int_t i = 0; i < nVertices; ++ i)
                {
                        if (!finite. test (i))
                                len [i] = infty;
                }
                *infinity = infty;
        }

        finite. free ();
        return noNegativeLoop;
} /* diGraph::BellmanFord */

template <class wType> template <class lenTable, class predTable>
inline bool diGraph<wType>::BellmanFord (int_t source, lenTable &len,
        wType *infinity, predTable pred) const
{
        const dummyRounding<wType> rounding = dummyRounding<wType> ();
        return this -> BellmanFord (rounding, source, len, infinity, pred);
} /* diGraph::BellmanFord */

template <class wType> template <class roundType>
inline bool diGraph<wType>::BellmanFord (const roundType &rounding,
        int_t source) const
{
        std::auto_ptr<wType> len_ptr (new wType [nVertices]);
        wType *len = len_ptr. get ();
        wType *infinity = 0;
        dummyArray tab;
        return BellmanFord (rounding, source, len, infinity, tab);
} /* diGraph::BellmanFord */

template <class wType>
inline bool diGraph<wType>::BellmanFord (int_t source) const
{
        const dummyRounding<wType> rounding = dummyRounding<wType> ();
        return BellmanFord (rounding, source);
} /* diGraph::BellmanFord */

template <class wType>
inline wType diGraph<wType>::Edmonds () const
{
        // create a list of edges with weights and sort this list
        std::vector<edgeTriple> edgeTriples (countEdges ());
        int_t prevEdge = 0;
        int_t curEdge = 0;
        for (int_t vert = 0; vert < nVertices; ++ vert)
        {
                while (curEdge < edgeEnds [vert])
                {
                        edgeTriple &e = edgeTriples [curEdge];
                        e. vert1 = vert;
                        e. vert2 = edges [curEdge];
                        e. weight = weights [curEdge];
                        ++ curEdge;
                }
                prevEdge = curEdge;
        }
        std::sort (edgeTriples. begin (), edgeTriples. end ());

        // create a forest which initially contains single vertices
        std::auto_ptr<int_t> root_ptr (new int_t [nVertices]);
        int_t *root = root_ptr. get ();
        for (int_t vert = 0; vert < nVertices; ++ vert)
        {
                root [vert] = -1;
        }

        // go through the edges and join the trees, but avoid loops
        wType totalWeight = 0;
        int_t nEdges = countEdges ();
        for (int_t curEdge = 0; curEdge < nEdges; ++ curEdge)
        {
                // determine the roots of both vertices of the edge
                // and compress the paths
                edgeTriple &e = edgeTriples [curEdge];
                int_t root1 = e. vert1;
                while (root [root1] >= 0)
                {
                        root1 = root [root1];
                }
                int_t vert1 = e. vert1;
                while (root [vert1] >= 0)
                {
                        int_t next = root [vert1];
                        root [vert1] = root1;
                        vert1 = next;
                }
                int_t root2 = e. vert2;
                while (root [root2] >= 0)
                {
                        root2 = root [root2];
                }
                int_t vert2 = e. vert2;
                while (root [vert2] >= 0)
                {
                        int_t next = root [vert2];
                        root [vert2] = root2;
                        vert2 = next;
                }

                // skip the edge if it closes a loop
                if (root1 == root2)
                        continue;

                // add the weight
                totalWeight += e. weight;

                // join the trees
                root [root1] = root2;
        }
        return totalWeight;
} /* diGraph::Edmonds */

template <class wType>
inline wType diGraph<wType>::EdmondsOld () const
{
        // create a list of edges with weights and sort this list
        std::vector<edgeTriple> edgeTriples (countEdges ());
        int_t prevEdge = 0;
        int_t curEdge = 0;
        for (int_t vert = 0; vert < nVertices; ++ vert)
        {
                while (curEdge < edgeEnds [vert])
                {
                        edgeTriple &e = edgeTriples [curEdge];
                        e. vert1 = vert;
                        e. vert2 = edges [curEdge];
                        e. weight = weights [curEdge];
                        ++ curEdge;
                }
                prevEdge = curEdge;
        }
        std::sort (edgeTriples. begin (), edgeTriples. end ());

        // create a forest which initially contains single vertices
        std::auto_ptr<int_t> forest_ptr (new int_t [nVertices]);
        int_t *forest = forest_ptr. get ();
        std::auto_ptr<int_t> next_ptr (new int_t [nVertices]);
        int_t *next = next_ptr. get ();
        std::auto_ptr<int_t> prev_ptr (new int_t [nVertices]);
        int_t *prev = prev_ptr. get ();
        for (int_t vert = 0; vert < nVertices; ++ vert)
        {
                forest [vert] = vert;
                next [vert] = -1;
                prev [vert] = -1;
        }

        // go through the edges and join the trees, but avoid loops
        wType totalWeight = 0;
        int_t nEdges = countEdges ();
        for (int_t curEdge = 0; curEdge < nEdges; ++ curEdge)
        {
                // check the edge and skip it if it closes a loop
                edgeTriple &e = edgeTriples [curEdge];
                if (forest [e. vert1] == forest [e. vert2])
                        continue;

                // add the weight
                totalWeight += e. weight;

                // go to the end of the first tree
                int_t vert1 = e. vert1;
                while (next [vert1] >= 0)
                {
                        vert1 = next [vert1];
                }
                
                // go to the beginning of the second tree
                int_t vert2 = e. vert2;
                while (prev [vert2] >= 0)
                {
                        vert2 = prev [vert2];
                }

                // join the trees and modify the numbers
                next [vert1] = vert2;
                prev [vert2] = vert1;
                int_t treeNumber = forest [vert1];
                while (vert2 >= 0)
                {
                        forest [vert2] = treeNumber;
                        vert2 = next [vert2];
                }
        }
        return totalWeight;
} /* diGraph::EdmondsOld */

template <class wType>
template <class arrayType, class roundType>
inline wType diGraph<wType>::FloydWarshall (const roundType &rounding,
        arrayType &arr, bool setInfinity, bool ignoreNegLoop) const
{
        // do nothing if the graph is empty
        if (!nVertices)
                return 0;

        // prepare marks to indicate finite values (not "infinity")
        BitField *finite = new BitField [nVertices];
        for (int_t i = 0; i < nVertices; ++ i)
        {
                finite [i]. allocate (nVertices);
                finite [i]. clearall (nVertices);
        }

        // create the initial values of the array based on the edge weights
        int_t curEdge = 0;
        for (int_t source = 0; source < nVertices; ++ source)
        {
                bool diagset = false;
                while (curEdge < edgeEnds [source])
                {
                        int_t target = edges [curEdge];
                        const wType &weight = weights [curEdge];
                        if (source == target)
                        {
                                if (weight < 0)
                                {
                                        arr [source] [target] = weight;
                                        diagset = true;
                                }
                        }
                        else
                        {
                                arr [source] [target] = weight;
                                finite [source]. set (target);
                        }
                        ++ curEdge;
                }
                if (!diagset)
                        arr [source] [source] = 0;
                finite [source]. set (source);
        }

        // find the shortest paths between the vertices (dynamic programming)
        for (int_t k = 0; k < nVertices; ++ k)
        {
                for (int_t i = 0; i < nVertices; ++ i)
                {
                        if (!finite [i]. test (k))
                                continue;
                        for (int_t j = 0; j < nVertices; ++ j)
                        {
                                if (!finite [k]. test (j))
                                        continue;
                                const wType sum = rounding. add_down
                                        (arr [i] [k], arr [k] [j]);
                                if (finite [i]. test (j))
                                {
                                        if (sum < arr [i] [j])
                                                arr [i] [j] = sum;
                                }
                                else
                                {
                                        arr [i] [j] = sum;
                                        finite [i]. set (j);
                                }
                        }
                }
        }

        // verify if a negative loop exists by checking for a negative
        // result in the diagonal
        if (!ignoreNegLoop)
        {
                for (int_t i = 0; i < nVertices; ++ i)
                {
                        if (arr [i] [i] < 0)
                                throw "Negative loop in Floyd-Warshall.";
                }
        }

        // prepare a variable to store the returned result
        wType result = 0;

        // compute the value for the infinity and fill in the array
        // if requested to do so
        if (setInfinity)
        {
                wType &infinity = result;
                for (int_t i = 0; i < nVertices; ++ i)
                {
                        for (int_t j = 0; j < nVertices; ++ j)
                        {
                                if (finite [i]. test (j) &&
                                        (infinity <= arr [i] [j]))
                                {
                                        infinity = rounding. add_up
                                                (arr [i] [j], 1);
                                }
                        }
                }
                for (int_t i = 0; i < nVertices; ++ i)
                {
                        for (int_t j = 0; j < nVertices; ++ j)
                        {
                                if (!finite [i]. test (j))
                                        arr [i] [j] = infinity;
                        }
                }
        }

        // otherwise, only compute the minimum path weight
        else
        {
                wType &minWeight = result;
                bool firstTime = true;
                for (int_t i = 0; i < nVertices; ++ i)
                {
                        for (int_t j = 0; j < nVertices; ++ j)
                        {
                                if (finite [i]. test (j))
                                {
                                        if (firstTime)
                                        {
                                                minWeight = arr [i] [j];
                                                firstTime = false;
                                        }
                                        else if (arr [i] [j] < minWeight)
                                        {
                                                minWeight = arr [i] [j];
                                        }
                                }
                        }
                }
        }

        // release the 'finite' bitfields
        for (int_t i = 0; i < nVertices; ++ i)
                finite [i]. free ();
        delete [] finite;

        return result;
} /* diGraph::FloydWarshall */

template <class wType>
template <class arrayType>
inline wType diGraph<wType>::FloydWarshall (arrayType &arr,
        bool setInfinity, bool ignoreNegLoop) const
{
        const dummyRounding<wType> rounding = dummyRounding<wType> ();
        return FloydWarshall (rounding, arr, setInfinity, ignoreNegLoop);
} /* diGraph::FloydWarshall */

template <class wType>
template <class arrayType, class roundType>
inline wType diGraph<wType>::Johnson (const roundType &rounding,
        arrayType &arr, bool setInfinity, bool ignoreNegLoop) const
{
        // DEBUG VERIFICATION
        if (false && sbug. show)
        {
                timeused stopWatch;
                wType res = FloydWarshall (rounding,
                        arr, setInfinity, ignoreNegLoop);
                chomp::homology::sbug << "\n[Floyd-Warshall: " << res <<
                        ", " << (double) stopWatch << " sec]\n";
        }
        // debug time measurement (see below)
//      timeused stopWatch;

        // do nothing if the graph is empty
        if (!nVertices)
                return 0;

        // STEP 1: Compute the shortest paths to any vertex from an
        // artificial extra vertex connected to every other vertex in the
        // graph by an edge of weight zero. Use Bellman-Ford for this.
        wType *len = new wType [nVertices];
        for (int_t i = 0; i < nVertices; ++ i)
                len [i] = 0;

        // update the lenghts of the paths repeatedly (max nVertices times)
        bool noNegativeLoop = false;
        int_t counter = 0;
        for (; counter <= nVertices; ++ counter)
        {
                bool modified = false;
                int_t curEdge = 0;
                for (int_t vertex = 0; vertex < nVertices; ++ vertex)
                {
                        int_t maxEdge = edgeEnds [vertex];
                        for (; curEdge < maxEdge; ++ curEdge)
                        {
                                int_t next = edges [curEdge];
                                wType newlen = rounding. add_down
                                        (len [vertex], weights [curEdge]);
                                if (newlen < len [next])
                                {
                                        // this "if" statement is necessary
                                        // because of a bug in GCC 3.4.2...
                                        if (counter > nVertices)
                                        {
                                                std::cout << vertex;
                                        }
                                        modified = true;
                                        len [next] = newlen;
                                }
                        }
                }
                if (!modified)
                {
                        noNegativeLoop = true;
                        ++ counter;
                        break;
                }
        }
        if (!ignoreNegLoop && !noNegativeLoop)
                throw "Negative loop found in Johnson's algorithm.";

        // STEP 2: Re-weigh the edges using the computed path lengths.
        wType *newWeights = new wType [edgeEnds [nVertices - 1]];
        int_t edge = 0;
        for (int_t vertex = 0; vertex < nVertices; ++ vertex)
        {
                int_t maxEdge = edgeEnds [vertex];
                for (; edge < maxEdge; ++ edge)
                {
                        wType w = weights [edge];
                        w = rounding. add_down (w, len [vertex]);
                        w = rounding. sub_down (w, len [edges [edge]]);
                        newWeights [edge] = (w < 0) ?
                                static_cast<wType> (0) : w;
                }
        }

        // STEP 3: Run the Dijkstra algorithm for every vertex to compute
        // the shortest paths to all the other vertices.
        // Negative entries indicate no path existence.
        for (int_t u = 0; u < nVertices; ++ u)
        {
                this -> Dijkstra (rounding, u, arr [u], newWeights);
        }
        delete [] newWeights;

        // STEP 4: Make a correction to the computed path lengths.
        // Compute the value for infinity if requested to.
        // Otherwise compute the minimum of path lengths.
        wType result = 0;
        if (setInfinity)
        {
                wType &infinity = result;
                for (int_t u = 0; u < nVertices; ++ u)
                {
                        for (int_t v = 0; v < nVertices; ++ v)
                        {
                                wType w = arr [u] [v];
                                if (w < 0)
                                        continue;
                                w = rounding. add_down (w, len [v]);
                                w = rounding. sub_down (w, len [u]);
                                if (w < infinity)
                                        continue;
                                infinity = rounding. add_up (w, 1);
                        }
                }
                for (int_t u = 0; u < nVertices; ++ u)
                {
                        for (int_t v = 0; v < nVertices; ++ v)
                        {
                                wType w = arr [u] [v];
                                if (w < 0)
                                {
                                        arr [u] [v] = infinity;
                                        continue;
                                }
                                w = rounding. add_down (w, len [v]);
                                arr [u] [v] =
                                        rounding. sub_down (w, len [u]);
                        }
                }
        }
        else
        {
                wType &minWeight = result;
                bool firstTime = true;
                for (int_t u = 0; u < nVertices; ++ u)
                {
                        for (int_t v = 0; v < nVertices; ++ v)
                        {
                                wType w = arr [u] [v];
                                if (w < 0)
                                        continue;
                                w = rounding. add_down (w, len [v]);
                                w = rounding. sub_down (w, len [u]);
                                if (firstTime)
                                {
                                        minWeight = w;
                                        firstTime = false;
                                }
                                else if (w < minWeight)
                                        minWeight = w;
                        }
                }
        }
        delete [] len;

        // DEBUG VERIFICATION
        if (false && sbug. show)
        {
//              chomp::homology::sbug << "[Johnson: " << result <<
//                      ", " << (double) stopWatch << " sec]\n";
        }

        return result;
} /* diGraph::Johnson */

template <class wType>
template <class arrayType>
inline wType diGraph<wType>::Johnson (arrayType &arr,
        bool setInfinity, bool ignoreNegLoop) const
{
        const dummyRounding<wType> rounding = dummyRounding<wType> ();
        return Johnson (rounding, arr, setInfinity, ignoreNegLoop);
} /* diGraph::Johnson */

template <class wType>
template <class roundType>
wType diGraph<wType>::minPathWeight (const roundType &rounding,
        bool ignoreNegLoop, int sparseGraph) const
{
        // if the graph is empty, return 0 as the minimum path weight
        if (nVertices <= 0)
                return 0;

        // allocate memory for an array of arrays to store min path weights
        wType **arr = new wType * [nVertices];
        for (int_t i = 0; i < nVertices; ++ i)
                arr [i] = new wType [nVertices];

        // determine whether to run the Floyd-Warshall algorithm
        // or Johnson's algorithm
        bool sparse = false;
        if (sparseGraph == 1)
                sparse = true;
        else if (sparseGraph != 0)
        {
                double nEdgesD = this -> countEdges ();
                double nVerticesD = this -> countVertices ();
                if ((nVerticesD > 3000) && (nEdgesD < nVerticesD * 1000) &&
                        (nEdgesD < nVerticesD * nVerticesD / 50))
                {
                        sparse = true;
                }
        }

        // run the Johnson's or Floyd-Warshall algorithm
        wType minWeight = sparse ?
                this -> Johnson (rounding, arr, false, ignoreNegLoop) :
                this -> FloydWarshall (rounding, arr, false, ignoreNegLoop);

/*      // compute the minimum of all the paths
        wType minWeight = arr [0] [0];
        for (int_t i = 0; i < nVertices; ++ i)
        {
                for (int_t j = 0; j < nVertices; ++ j)
                {
                        const wType &weight = arr [i] [j];
                        if (weight < minWeight)
                                minWeight = weight;
                }
        }
*/
        // release the memory
        for (int_t i = 0; i < nVertices; ++ i)
                delete [] (arr [i]);
        delete [] arr;

        return minWeight;
} /* diGraph::minPathWeight */

template <class wType>
wType diGraph<wType>::minPathWeight (bool ignoreNegLoop, int sparseGraph)
        const
{
        const dummyRounding<wType> rounding = dummyRounding<wType> ();
        return this -> minPathWeight (rounding, ignoreNegLoop, sparseGraph);
} /* diGraph::minPathWeight */

template <class wType> template <class outType>
inline outType &diGraph<wType>::show (outType &out, bool showWeights) const
{
        out << "; Directed graph: " << nVertices << " vertices.\n";
        int_t curEdge = 0;
        for (int_t i = 0; i < nVertices; ++ i)
        {
                for (; curEdge < edgeEnds [i]; ++ curEdge)
                {
                        out << i << " -> " << edges [curEdge];
                        if (showWeights)
                                out << " [" << weights [curEdge] << "]\n";
                        else
                                out << "\n";
                }
        }
        out << "; This is the end of the graph.\n";
        return out;
} /* diGraph::show */

        
// --------------------------------------------------

template <class wType>
inline std::ostream &operator << (std::ostream &out, const diGraph<wType> &g)
{
        return g. show (out, false);
} /* operator << */

// --------------------------------------------------

template <class wType, class Table1, class Table2>
inline int_t SCC (const diGraph<wType> &g, Table1 &compVertices,
        Table2 &compEnds, diGraph<wType> *scc = 0,
        diGraph<wType> *transposed = 0, bool copyweights = false)
{
        // prepare two tables
        int_t nVert = g. countVertices ();
        int_t *ordered = new int_t [nVert];
        int_t *tab = new int_t [nVert];

        // compute the list of vertices in the descending finishing time
        g. DFSfinishTime (tab);
        for (int_t i = 0; i < nVert; ++ i)
                ordered [nVert - tab [i]] = i;
        delete [] tab;

        // create the transposed graph
        diGraph<wType> gT;
        if (!transposed)
                transposed = &gT;
        g. transpose (*transposed, copyweights);

        // extract the DFS forest of gT in the given order of vertices
        int_t n = transposed -> DFSforest (ordered, compVertices, compEnds,
                true, scc);

        // cleanup memory and return the number of components
        delete [] ordered;
        return n;
} /* SCC */

// --------------------------------------------------

template <class wType, class Table1, class Table2>
inline int_t SCC_Tarjan (const diGraph<wType> &g, Table1 &compVertices,
        Table2 &compEnds)
{
        // return the obvious result if the graph is empty
        int_t nVertices = g. countVertices ();
        if (!nVertices)
                return 0;

        // prepare an array of discovery times for all the vertices
        // (zero == not yet discovered)
        std::vector<int_t> dfsIndex (nVertices, 0);

        // prepare an array of the minimal index of a node reachable
        // from each of the vertices
        std::vector<int_t> lowLink (nVertices, 0);

        // prepare an empty stack of nodes
        std::stack<int_t> s_nodes;

        // prepare an array of bits indicating whether the vertices are
        // in the stack of nodes or not
        BitField inTheStack;
        inTheStack. allocate (nVertices);
        inTheStack. clearall (nVertices);

        // prepare the number of strongly connected components
        int_t nComponents = 0;

        // prepare the current position in the array 'compVertices'
        int_t posVertices = 0;

        // remember the next vertex number in the graph to scan
        // whether this vertex has already been visited or not
        int_t vertexToScan = 0;

        // prepare a variable for storing the discovery time in the DFS
        int_t discoveryTime = 0;

        // prepare stacks for the DFS recursion
        std::stack<int_t> s_vertex;
        std::stack<int_t> s_edge;
        std::stack<int_t> s_maxedge;

        // initialize the number of the currently processed vertex
        int_t vertex = -1;

        // initialize the range of edges to be visited
        int_t edge = 0;
        int_t maxedge = 0;

        while (1)
        {
                // return to the previous recursion level
                // if all the edges have been checked
                if (edge >= maxedge)
                {
                        // extract a strongly connected component
                        if ((vertex >= 0) && (lowLink [vertex] ==
                                dfsIndex [vertex]))
                        {
                                int_t v = 0;
                                do
                                {
                                        v = s_nodes. top ();
                                        s_nodes. pop ();
                                        inTheStack. clear (v);
                                        compVertices [posVertices ++] = v;
                                } while (v != vertex);
                                compEnds [nComponents ++] = posVertices;
                        }

                        // if this is the top level of the recursion
                        // then find another unvisited vertex or return
                        if (s_vertex. empty ())
                        {
                                // find an unvisited vertex in the graph
                                while ((vertexToScan < nVertices) &&
                                        (dfsIndex [vertexToScan] != 0))
                                {
                                        ++ vertexToScan;
                                }

                                // return the result if all visited
                                if (vertexToScan == nVertices)
                                {
                                        inTheStack. free ();
                                        return nComponents;
                                }

                                // set the new vertex
                                vertex = vertexToScan ++;

                                // mark the current vertex as visited
                                // and initialize its low link
                                dfsIndex [vertex] = ++ discoveryTime;
                                lowLink [vertex] = discoveryTime;

                                // push this vertex on the stack
                                s_nodes. push (vertex);
                                inTheStack. set (vertex);

                                // determine the edges to be visited
                                edge = 0;
                                maxedge = g. countEdges (vertex);
                        }

                        // otherwise trace back to the previous level
                        else
                        {
                                // remember the current low link index
                                int_t lowLink2 = lowLink [vertex];

                                // restore the variables
                                vertex = s_vertex. top ();
                                s_vertex. pop ();
                                edge = s_edge. top ();
                                s_edge. pop ();
                                maxedge = s_maxedge. top ();
                                s_maxedge. pop ();

                                // update the current low link index
                                if (lowLink [vertex] > lowLink2)
                                        lowLink [vertex] = lowLink2;
                        }
                }

                // analyse the next edge coming out from the current vertex
                else
                {
                        // determine the next vertex
                        int_t next = g. getEdge (vertex, edge ++);

                        // go to a deeper recursion level if unvisited
                        if (dfsIndex [next] == 0)
                        {
                                // store the variables at the stacks
                                s_vertex. push (vertex);
                                s_edge. push (edge);
                                s_maxedge. push (maxedge);

                                // set the new vertex
                                vertex = next;
                        
                                // mark the new vertex as visited
                                dfsIndex [vertex] = ++ discoveryTime;
                                lowLink [vertex] = discoveryTime;
                        
                                // push this vertex on the stack
                                s_nodes. push (vertex);
                                inTheStack. set (vertex);

                                // determine the edges to be visited
                                edge = 0;
                                maxedge = g. countEdges (vertex);
                        }

                        // update the low link index if the vertex has been
                        // visited and is currently in the stack of nodes
                        else if (inTheStack. test (next))
                        {
                                if (lowLink [vertex] > dfsIndex [next])
                                        lowLink [vertex] = dfsIndex [next];
                        }
                }
        }

        // finalize and return the number of strongly connected components
        inTheStack. free ();
        return nComponents;
} /* SCC_Tarjan */

// --------------------------------------------------

template <class wType>
wType diGraph<wType>::minMeanCycleWeight (diGraph<wType> *transposed) const
{
        // find the strongly connected components of the graph
        multitable<int_t> compVertices, compEnds;
        bool copyweights = !!transposed;
        int_t countSCC = SCC (*this, compVertices, compEnds,
                static_cast<diGraph<wType> *> (0), transposed, copyweights);
        if (countSCC <= 0)
                return 0;

        // compute the maximum size of each strongly connected component
        int_t maxCompSize = compEnds [0];
        for (int_t comp = 1; comp < countSCC; ++ comp)
        {
                int_t compSize = compEnds [comp] - compEnds [comp - 1];
                if (maxCompSize < compSize)
                        maxCompSize = compSize;
        }

        // allocate arrays for weights and bit fields
        wType **F = new wType * [maxCompSize + 1];
        BitField *finite = new BitField [maxCompSize + 1];
        for (int_t i = 0; i <= maxCompSize; ++ i)
        {
                F [i] = new wType [maxCompSize];
                finite [i]. allocate (maxCompSize);
        }

        // compute the numbers of vertices in each component
        int_t *numbers = new int_t [nVertices];
        int_t *components = new int_t [nVertices];
        for (int_t i = 0; i < nVertices; ++ i)
                components [i] = -1;
        int_t offset = 0;
        for (int_t comp = 0; comp < countSCC; ++ comp)
        {
                int_t maxOffset = compEnds [comp];
                int_t pos = offset;
                for (; pos < maxOffset; ++ pos)
                {
                        numbers [compVertices [pos]] = pos - offset;
                        components [compVertices [pos]] = comp;
                }
                offset = pos;
        }

        // compute the minimum mean cycle weight for each component
        wType minWeight (0);
        for (int_t comp = 0; comp < countSCC; ++ comp)
        {
                int_t offset = comp ? compEnds [comp - 1] :
                        static_cast<int_t> (0);
                int_t compSize = compEnds [comp] - offset;
                for (int_t i = 0; i <= compSize; ++ i)
                        finite [i]. clearall (compSize);
                F [0] [0] = 0;
                finite [0]. set (0);
                // compute path progressions of given length
                for (int_t len = 1; len <= compSize; ++ len)
                {
                        // process source vertices
                        for (int_t i = 0; i < compSize; ++ i)
                        {
                                if (!finite [len - 1]. test (i))
                                        continue;
                                wType prevWeight = F [len - 1] [i];
                                int_t source = compVertices [offset + i];

                                // process target vertices (and edges)
                                int_t edgeOffset = source ?
                                        edgeEnds [source - 1] :
                                        static_cast<int_t> (0);
                                int_t edgeEnd = edgeEnds [source];
                                for (; edgeOffset < edgeEnd; ++ edgeOffset)
                                {
                                        int_t target = edges [edgeOffset];
                                        if (components [target] != comp)
                                                continue;
                                        int_t j = numbers [target];
                                        wType newWeight = prevWeight +
                                                weights [edgeOffset];
                                        if (!finite [len]. test (j))
                                        {
                                                finite [len]. set (j);
                                                F [len] [j] = newWeight;
                                        }
                                        else if (newWeight < F [len] [j])
                                                F [len] [j] = newWeight;
                                }
                        }
                }

                // compute the minimum mean cycle weight for this component
                wType minCompWeight (0);
                bool firstMin = true;
                for (int_t vert = 0; vert < compSize; ++ vert)
                {
                        if (!finite [compSize]. test (vert))
                                continue;
                        bool firstAverage = true;
                        wType maxAverage = 0;
                        for (int_t first = 0; first < compSize; ++ first)
                        {
                                if (!finite [first]. test (vert))
                                        continue;
                                wType average = (F [compSize] [vert] -
                                        F [first] [vert]) /
                                        (compSize - first);
                                if (firstAverage)
                                {
                                        maxAverage = average;
                                        firstAverage = false;
                                }
                                else if (maxAverage < average)
                                        maxAverage = average;
                        }
                        if (firstMin || (maxAverage < minCompWeight))
                        {
                                if (firstAverage)
                                        throw "Bug in Karp's algorithm";
                                minCompWeight = maxAverage;
                                firstMin = false;
                        }
                }

                // make a correction to the total minimum if necessary
                if (!comp || (minCompWeight < minWeight))
                        minWeight = minCompWeight;
        }

        // release allocated memory
        delete [] components;
        delete [] numbers;
        for (int_t i = 0; i < maxCompSize; ++ i)
        {
                finite [i]. free ();
                delete [] (F [i]);
        }
        delete [] finite;
        delete [] F;

        // return the computed minimum
        return minWeight;
} /* diGraph::minMeanCycleWeight */

template <class wType>
template <class roundType>
wType diGraph<wType>::minMeanCycleWeight (const roundType &rounding,
        diGraph<wType> *transposed) const
{
        // find the strongly connected components of the graph
        multitable<int_t> compVertices, compEnds;
        bool copyweights = !!transposed;
        int_t countSCC = SCC (*this, compVertices, compEnds,
                static_cast<diGraph<wType> *> (0), transposed, copyweights);
        if (countSCC <= 0)
                return 0;

        // compute the maximum size of each strongly connected component
        int_t maxCompSize = compEnds [0];
        for (int_t comp = 1; comp < countSCC; ++ comp)
        {
                int_t compSize = compEnds [comp] - compEnds [comp - 1];
                if (maxCompSize < compSize)
                        maxCompSize = compSize;
        }

        // allocate arrays for weights and bit fields
        wType **FUpper = new wType * [maxCompSize + 1];
        wType **FLower = new wType * [maxCompSize + 1];
        BitField *finite = new BitField [maxCompSize + 1];
        for (int_t i = 0; i <= maxCompSize; ++ i)
        {
                FUpper [i] = new wType [maxCompSize];
                FLower [i] = new wType [maxCompSize];
                finite [i]. allocate (maxCompSize);
        }

        // compute the numbers of vertices in each component
        int_t *numbers = new int_t [nVertices];
        int_t *components = new int_t [nVertices];
        for (int_t i = 0; i < nVertices; ++ i)
                components [i] = -1;
        int_t offset = 0;
        for (int_t comp = 0; comp < countSCC; ++ comp)
        {
                int_t maxOffset = compEnds [comp];
                int_t pos = offset;
                for (; pos < maxOffset; ++ pos)
                {
                        numbers [compVertices [pos]] = pos - offset;
                        components [compVertices [pos]] = comp;
                }
                offset = pos;
        }

        // compute the minimum mean cycle weight for each component
        wType minWeight (0);
        for (int_t comp = 0; comp < countSCC; ++ comp)
        {
                int_t offset = comp ? compEnds [comp - 1] :
                        static_cast<int_t> (0);
                int_t compSize = compEnds [comp] - offset;
                for (int_t i = 0; i <= compSize; ++ i)
                        finite [i]. clearall (compSize);
                FUpper [0] [0] = 0;
                FLower [0] [0] = 0;
                finite [0]. set (0);
                // compute path progressions of given length
                for (int_t len = 1; len <= compSize; ++ len)
                {
                        // process source vertices
                        for (int_t i = 0; i < compSize; ++ i)
                        {
                                if (!finite [len - 1]. test (i))
                                        continue;
                                wType prevUpper = FUpper [len - 1] [i];
                                wType prevLower = FLower [len - 1] [i];
                                int_t source = compVertices [offset + i];

                                // process target vertices (and edges)
                                int_t edgeOffset = source ?
                                        edgeEnds [source - 1] :
                                        static_cast<int_t> (0);
                                int_t edgeEnd = edgeEnds [source];
                                for (; edgeOffset < edgeEnd; ++ edgeOffset)
                                {
                                        int_t target = edges [edgeOffset];
                                        if (components [target] != comp)
                                                continue;
                                        int_t j = numbers [target];
                                        wType newUpper = rounding. add_up
                                                (prevUpper,
                                                weights [edgeOffset]);
                                        wType newLower = rounding. add_down
                                                (prevLower,
                                                weights [edgeOffset]);
                                        if (!finite [len]. test (j))
                                        {
                                                finite [len]. set (j);
                                                FUpper [len] [j] = newUpper;
                                                FLower [len] [j] = newLower;
                                        }
                                        else
                                        {
                                                wType &curUpper =
                                                        FUpper [len] [j];
                                                if (newUpper < curUpper)
                                                        curUpper = newUpper;
                                                wType &curLower =
                                                        FLower [len] [j];
                                                if (newLower < curLower)
                                                        curLower = newLower;
                                        }
                                }
                        }
                }

                // compute the minimum mean cycle weight for this component
                wType minCompWeight (0);
                bool firstMin = true;
                for (int_t vert = 0; vert < compSize; ++ vert)
                {
                        if (!finite [compSize]. test (vert))
                                continue;
                        bool firstAverage = true;
                        wType maxAverage = 0;
                        for (int_t first = 0; first < compSize; ++ first)
                        {
                                if (!finite [first]. test (vert))
                                        continue;
                                const wType diff = rounding. sub_down
                                        (FLower [compSize] [vert],
                                        FUpper [first] [vert]);
                                wType average = rounding. div_down
                                        (diff, compSize - first);
                                if (firstAverage)
                                {
                                        maxAverage = average;
                                        firstAverage = false;
                                }
                                else if (maxAverage < average)
                                        maxAverage = average;
                        }
                        if (firstMin || (maxAverage < minCompWeight))
                        {
                                if (firstAverage)
                                        throw "Bug in Karp's algorithm";
                                minCompWeight = maxAverage;
                                firstMin = false;
                        }
                }

                // make a correction to the total minimum if necessary
                if (!comp || (minCompWeight < minWeight))
                        minWeight = minCompWeight;
        }

        // release allocated memory
        delete [] components;
        delete [] numbers;
        for (int_t i = 0; i < maxCompSize; ++ i)
        {
                finite [i]. free ();
                delete [] (FUpper [i]);
                delete [] (FLower [i]);
        }
        delete [] finite;
        delete [] FUpper;
        delete [] FLower;

        // return the computed minimum
        return minWeight;
} /* diGraph::minMeanCycleWeight_intv */

template <class wType>
template <class arrayType, class roundType>
wType diGraph<wType>::minMeanPathWeight (const roundType &rounding,
        const arrayType &starting, int_t n) const
{
        // allocate arrays for weights and bit fields
        const int nIndices = 2;
        wType **F = new wType * [nIndices];
        BitField *finite = new BitField [nIndices];
        for (int i = 0; i < nIndices; ++ i)
        {
                F [i] = new wType [nVertices];
                finite [i]. allocate (nVertices);
        }

        // set the zero path lengths from the initial vertices
        for (int_t i = 0; i < n; ++ i)
        {
                int_t v = starting [i];
                if ((v < 0) || (v >= nVertices))
                        throw "Starting vertex out of range.";
                F [0] [v] = 0;
                finite [0]. set (v);
        }

        // compute path progressions of given length and average weights
        wType minWeight (0);
        bool firstWeight = true;
        for (int_t len = 1; len <= nVertices; ++ len)
        {
                // determine the indices for previous and current paths
                int_t prevIndex = (len - 1) & 1;
                int_t curIndex = len & 1;
                finite [curIndex]. clearall (nVertices);

                // process source vertices
                for (int_t source = 0; source < nVertices; ++ source)
                {
                        if (!finite [prevIndex]. test (source))
                                continue;
                        wType prevWeight = F [prevIndex] [source];

                        // process target vertices (and edges)
                        int_t edgeOffset = source ?
                                edgeEnds [source - 1] :
                                static_cast<int_t> (0);
                        int_t edgeEnd = edgeEnds [source];
                        for (; edgeOffset < edgeEnd; ++ edgeOffset)
                        {
                                int_t target = edges [edgeOffset];
                                wType newWeight = rounding. add_down
                                        (prevWeight, weights [edgeOffset]);
                                if (!finite [curIndex]. test (target))
                                {
                                        finite [curIndex]. set (target);
                                        F [curIndex] [target] = newWeight;
                                }
                                else if (newWeight < F [curIndex] [target])
                                        F [curIndex] [target] = newWeight;
                        }
                }

                // update the minimum mean path weight
                for (int_t vert = 0; vert < nVertices; ++ vert)
                {
                        if (!finite [curIndex]. test (vert))
                                continue;
                        wType average = rounding. div_down
                                (F [curIndex] [vert], len);
                        if (firstWeight)
                        {
                                minWeight = average;
                                firstWeight = false;
                        }
                        else if (average < minWeight)
                                minWeight = average;
                }
        }

        // release allocated memory
        for (int i = 0; i < nIndices; ++ i)
        {
                finite [i]. free ();
                delete [] (F [i]);
        }
        delete [] finite;
        delete [] F;

        // return the computed minimum
        return minWeight;
} /* diGraph::minMeanPathWeight */

template <class wType>
template <class arrayType>
wType diGraph<wType>::minMeanPathWeight (const arrayType &starting, int_t n)
        const
{
        const dummyRounding<wType> rounding = dummyRounding<wType> ();
        return minMeanPathWeight (rounding, starting, n);
} /* diGraph::minMeanPathWeight */


// --------------------------------------------------
// ------------- OTHER GRAPH ALGORITHMS -------------
// --------------------------------------------------

template <class wType>
inline int_t addRenumEdges (const diGraph<wType> &g, int_t vertex,
        const int_t *newNum, int_t cur, int_t *srcVert,
        diGraph<wType> &result)
{
        int_t nEdges = g. countEdges (vertex);
        // add all the edges that start at the component
        for (int_t edge = 0; edge < nEdges; ++ edge)
        {
                // determine the dest. vertex of the edge
                int_t dest = newNum [g. getEdge (vertex, edge)];
                // if this is an edge to self, then ignore it
                if (dest == cur)
                        continue;
                // if the edge has already been added,
                // then skip it
                if (srcVert [dest] == cur)
                        continue;
                // remember that the dest. vertex has an edge
                // pointing out from the current vertex
                srcVert [dest] = cur;
                // add the edge to the result graph
                result. addEdge (dest);
        }
        return 0;
} /* addRenumEdges */

template <class wType, class Table1, class Table2>
inline int_t collapseVertices (const diGraph<wType> &g,
        int_t compNum, Table1 &compVertices, Table2 &compEnds,
        diGraph<wType> &result, int_t *newNumbers = 0)
{
        // do nothing if the input graph is empty
        int_t nVert = g. countVertices ();
        if (!nVert)
                return 0;

        // compute the new numbers of vertices: newNum [i] is the
        // number of vertex 'i' from 'g' in the resulting graph
        int_t *newNum = newNumbers ? newNumbers : new int_t [nVert];
        for (int_t i = 0; i < nVert; ++ i)
                newNum [i] = -1;
        int_t cur = 0; // current vertex number in the new graph
        int_t pos = 0; // position in compVertices
        while (cur < compNum)
        {
                int_t posEnd = compEnds [cur];
                while (pos < posEnd)
                        newNum [compVertices [pos ++]] = cur;
                ++ cur;
        }
        for (int_t i = 0; i < nVert; ++ i)
        {
                if (newNum [i] < 0)
                        newNum [i] = cur ++;
        }

        // prepare a table to mark the most recent source vertex for edges:
        // srcVert [j] contains i if the edge i -> j has just been added
        int_t newVert = nVert - (compNum ? compEnds [compNum - 1] :
                static_cast<int_t> (0)) + compNum;
        // debug:
        if (cur != newVert)
                throw "DEBUG: Wrong numbers.";
        int_t *srcVert = new int_t [newVert];
        for (int_t i = 0; i < newVert; ++ i)
                srcVert [i] = -1;

        // scan the input graph and create the resulting graph: begin with
        // the vertices in the collapsed groups
        cur = 0, pos = 0;
        while (cur < compNum)
        {
                // add a new vertex to the result graph
                result. addVertex ();
                // for all the vertices in this component...
                int_t posEnd = compEnds [cur];
                while (pos < posEnd)
                {
                        // take the next vertex from the component
                        int_t vertex = compVertices [pos ++];
                        // add the appropriate edges to the result graph
                        addRenumEdges (g, vertex, newNum, cur,
                                srcVert, result);
                }
                ++ cur;
        }
        // process the remaining vertices of the graph
        for (int_t vertex = 0; vertex < nVert; ++ vertex)
        {
                // debug
                if (newNum [vertex] > cur)
                        throw "DEBUG: Wrong order.";
                // if the vertex has already been processed, skip it
                if (newNum [vertex] != cur)
                        continue;
                // add the appropriate edges to the result graph
                result. addVertex ();
                addRenumEdges (g, vertex, newNum, cur, srcVert, result);
                ++ cur;
        }

        // free memory and exit
        delete [] srcVert;
        if (!newNumbers)
                delete [] newNum;
        return 0;
} /* collapseVertices */

template <class wType, class TabSets>
inline int_t addRenumEdges2 (const diGraph<wType> &g, int_t vertex,
        const int_t *newNum, TabSets &numSets,
        int_t cur, int_t *srcVert, diGraph<wType> &result)
{
        int_t nEdges = g. countEdges (vertex);

        // add all the edges that start at this vertex
        for (int_t edge = 0; edge < nEdges; ++ edge)
        {
                // determine the dest. vertex of the edge
                int_t destNumber = newNum [g. getEdge (vertex, edge)];

                // consider all the destination vertices
                int_t destSize = (destNumber < 0) ?
                        numSets [-destNumber]. size () :
                        static_cast<int_t> (1);
                for (int_t i = 0; i < destSize; ++ i)
                {
                        // determine the consecutive destination vertex
                        int_t dest = (destNumber < 0) ?
                                numSets [-destNumber] [i] : destNumber;

                        // if this is an edge to self, then ignore it
                        if (dest == cur)
                                continue;

                        // if the edge has already been added,
                        // then skip it
                        if (srcVert [dest] == cur)
                                continue;

                        // add the edge to the result graph
                        result. addEdge (dest);
                }
        }
        return 0;
} /* addRenumEdges2 */

template <class wType, class Table1, class Table2>
inline int_t collapseVertices2 (const diGraph<wType> &g,
        int_t compNum, Table1 &compVertices, Table2 &compEnds,
        diGraph<wType> &result)
{
        // do nothing if the input graph is empty
        int_t nVert = g. countVertices ();
        if (!nVert)
                return 0;

        // compute the new numbers of vertices: newNum [i] is the
        // number of vertex 'i' from 'g' in the resulting graph,
        // unless negative; then it points to a set of numbers,
        // with -1 meaning "not defined, yet"
        int_t *newNum = new int_t [nVert];
        for (int_t i = 0; i < nVert; ++ i)
                newNum [i] = -1;

        // sets of numbers of vertices; the numbers refering to the sets
        // begin with 2, thus the first two sets are skipped
        multitable<hashedset<int_t> > numSets;
        // the number of the set waiting to be used next
        int_t numSetCur = 2;

        int_t cur = 0; // current vertex number in the new graph
        int_t pos = 0; // position in compVertices
        while (cur < compNum)
        {
                int_t posEnd = compEnds [cur];
                while (pos < posEnd)
                {
                        int_t number = compVertices [pos ++];
                        if (newNum [number] == -1)
                                newNum [number] = cur;
                        else if (newNum [number] < 0)
                                numSets [-newNum [number]]. add (cur);
                        else
                        {
                                numSets [numSetCur]. add (newNum [number]);
                                numSets [numSetCur]. add (cur);
                                newNum [number] = -(numSetCur ++);
                        }
                }
                ++ cur;
        }
        for (int_t i = 0; i < nVert; ++ i)
        {
                if (newNum [i] == -1)
                        newNum [i] = cur ++;
        }

        // prepare a table to mark the most recent source vertex for edges:
        // srcVert [j] contains i if the edge i -> j has just been added
        int_t newVert = cur;
        int_t *srcVert = new int_t [newVert];
        for (int_t i = 0; i < newVert; ++ i)
                srcVert [i] = -1;

        // scan the input graph and create the resulting graph: begin with
        // the vertices in the collapsed groups
        cur = 0;
        pos = 0;
        while (cur < compNum)
        {
                // add a new vertex to the result graph
                result. addVertex ();

                // for all the vertices in this component...
                int_t posEnd = compEnds [cur];
                while (pos < posEnd)
                {
                        // take the next vertex from the component
                        int_t vertex = compVertices [pos ++];

                        // add the appropriate edges to the result graph
                        addRenumEdges2 (g, vertex, newNum, numSets, cur,
                                srcVert, result);
                }
                ++ cur;
        }

        // process the remaining vertices of the graph
        for (int_t vertex = 0; vertex < nVert; ++ vertex)
        {
                // debug
                if (newNum [vertex] > cur)
                        throw "DEBUG: Wrong order 2.";

                // if the vertex has already been processed, skip it
                if (newNum [vertex] != cur)
                        continue;

                // add the appropriate edges to the result graph
                result. addVertex ();
                addRenumEdges2 (g, vertex, newNum, numSets, cur,
                        srcVert, result);
                ++ cur;
        }

        // free memory and exit
        delete [] srcVert;
        delete [] newNum;
        return 0;
} /* collapseVertices2 */

template <class wType>
inline int_t connectionGraph (const diGraph<wType> &g, int_t nVert,
        diGraph<wType> &result)
{
        // remember the number of vertices in the input graph
        int_t nVertG = g. countVertices ();
        if (!nVertG)
                return 0;

        // prepare a bitfield in which visited vertices will be marked
        BitField visited;
        visited. allocate (nVertG);
        visited. clearall (nVertG);

        // run DFS for each of the starting vertices
        for (int_t startVertex = 0; startVertex < nVert; ++ startVertex)
        {
                // add this vertex to the resulting graph
                result. addVertex ();

                // prepare a counter and a stack of visited vertices
                int_t nVisited = 0;
                std::stack<int_t> s_visited;

                // prepare stacks for the DFS algorithm
                std::stack<int_t> s_vertex;
                std::stack<int_t> s_edge;

                // mark the starting vertex as visited
                visited. set (startVertex);
                s_visited. push (startVertex);
                ++ nVisited;

                // use DFS to visit vertices reachable from that vertex
                int_t vertex = startVertex;
                int_t edge = 0;
                while (1)
                {
                        // go back with the recursion
                        // if all the edges have been processed
                        if (edge >= g. countEdges (vertex))
                        {
                                // if this was the top recursion level
                                // then quit
                                if (s_vertex. empty ())
                                        break;

                                // return from the recursion
                                vertex = s_vertex. top ();
                                s_vertex. pop ();
                                edge = s_edge. top ();
                                s_edge. pop ();
                                continue;
                        }

                        // take the next vertex
                        int_t next = g. getEdge (vertex, edge ++);

                        // go to the deeper recursion level if not visited
                        if (!visited. test (next))
                        {
                                // add an edge to the result if necessary
                                if (next < nVert)
                                        result. addEdge (next);

                                // store the previous variables at the stacks
                                s_vertex. push (vertex);
                                s_edge. push (edge);

                                // take the new vertex
                                vertex = next;
                                edge = 0;

                                // mark the new vertex as visited
                                visited. set (vertex);
                                s_visited. push (vertex);
                                ++ nVisited;
                        }
                }

                // if this was the last vertex then skip the rest
                if (startVertex == nVert - 1)
                        break;

                // mark each vertex as non-visited
                if (nVisited > (nVertG >> 6))
                {
                        visited. clearall (nVertG);
                }
                else
                {
                        while (!s_visited. empty ())
                        {
                                int_t vertex = s_visited. top ();
                                s_visited. pop ();
                                visited. clear (vertex);
                        }
                }
        }

        // free memory and exit
        visited. free ();
        return 0;
} /* connectionGraph */

template <class GraphType>
inline int_t computePeriod (const GraphType &g)
{
        // remember the number of vertices in the input graph
        int_t nVertG = g. countVertices ();
        if (!nVertG)
                return 0;

        // prepare an array to record the depth of each visited vertex
        std::vector<int_t> visited (nVertG, 0);

        // run DFS starting at the first vertex
        // prepare stacks for the DFS algorithm
        std::stack<int_t> s_vertex;
        std::stack<int_t> s_edge;

        // mark the starting vertex as visited
        visited [0] = 1;

        // use DFS to visit vertices reachable from that vertex
        int_t vertex = 0;
        int_t edge = 0;
        int_t level = 1;
        int_t period = 0;
        while (1)
        {
                // go back with the recursion
                // if all the edges have been processed
                if (edge >= g. countEdges (vertex))
                {
                        // if this was the top recursion level then quit
                        if (s_vertex. empty ())
                                break;

                        // return from the recursion
                        vertex = s_vertex. top ();
                        s_vertex. pop ();
                        edge = s_edge. top ();
                        s_edge. pop ();
                        -- level;
                        continue;
                }

                // take the next vertex
                int_t next = g. getEdge (vertex, edge ++);

                // update the GCD of the cycle lengths
                if (visited [next])
                {
                        // determine the period provided by the cycle
                        int_t newPeriod = visited [next] - level - 1;
                        if (newPeriod < 0)
                                newPeriod = -newPeriod;

                        // compute the GCD of the old and new periods
                        int_t a = newPeriod;
                        int_t b = period;
                        while (b)
                        {
                                int_t t = b;
                                b = a % b;
                                a = t;
                        }
                        period = a;

                        // if the graph turns out to be aperiodic
                        // then immediately return the result
                        if (period == 1)
                                return period;
                }

                // go to the deeper recursion level if not visited
                else
                {
                        // store the previous variables at the stacks
                        s_vertex. push (vertex);
                        s_edge. push (edge);

                        // take the new vertex
                        vertex = next;
                        edge = 0;
                        ++ level;

                        // mark the new vertex as visited
                        visited [vertex] = level;
                }
        }

        // return the computed peirod and exit
        return period;
} /* computePeriod */

template <class wType>
inline int_t inclusionGraph (const diGraph<wType> &g, int_t nVert,
        diGraph<wType> &result)
{
        // remember the number of vertices in the input graph
        int_t nVertG = g. countVertices ();
        if (!nVertG)
                return 0;

        // mark each vertex as non-visited
        BitField visited;
        visited. allocate (nVertG);
        visited. clearall (nVertG);

        // create the sets of vertices reachable from each vertex
        BitSets lists (nVertG, nVert);

        // prepare stacks for the DFS algorithm
        std::stack<int_t> s_vertex;
        std::stack<int_t> s_edge;

        // use DFS to propagate the reachable sets
        int_t startVertex = 0;
        int_t vertex = 0;
        int_t edge = 0;
        visited. set (vertex);
        while (1)
        {
                // go back with the recursion
                // if all the edges have been processed
                if (edge >= g. countEdges (vertex))
                {
                        // if this was the top recursion level,
                        // then find another starting point or quit
                        if (s_vertex. empty ())
                        {
                                while ((startVertex < nVertG) &&
                                        (visited. test (startVertex)))
                                        ++ startVertex;
                                if (startVertex >= nVertG)
                                        break;
                                vertex = startVertex;
                                edge = 0;
                                visited. set (vertex);
                                continue;
                        }

                        // determine the previous vertex (to trace back to)
                        int_t prev = s_vertex. top ();

                        // add the list of the current vertex
                        // to the list of the previous vertex
                        // unless that was one of the first 'nVert' vertices
                        if (vertex >= nVert)
                                lists. addset (prev, vertex);

                        // return from the recursion
                        vertex = prev;
                        s_vertex. pop ();
                        edge = s_edge. top ();
                        s_edge. pop ();
                        continue;
                }

                // take the next vertex
                int_t next = g. getEdge (vertex, edge ++);

                // add the number to the list if appropriate
                if (next < nVert)
                        lists. add (vertex, next);

                // go to the deeper recursion level if vertex not visited
                if (!visited. test (next))
                {
                        // store the previous variables at the stacks
                        s_vertex. push (vertex);
                        s_edge. push (edge);

                        // take the new vertex and mark as visited
                        vertex = next;
                        edge = 0;
                        visited. set (vertex);
                }

                // add the list of the next vertex to the current one
                // if the vertex has already been visited before
                // and the vertex is not one of the first 'nVert' ones
                else if (next >= nVert)
                {
                        lists. addset (vertex, next);
                }
        }

        // create the result graph based on the lists of edges
        for (int_t vertex = 0; vertex < nVert; ++ vertex)
        {
                result. addVertex ();
                int_t edge = lists. get (vertex);
                while (edge >= 0)
                {
                        result. addEdge (edge);
                        edge = lists. get (vertex, edge + 1);
                }
        }

        // free memory and exit
        visited. free ();
        return 0;
} /* inclusionGraph */

// --------------------------------------------------

template<class wType, class TConn>
inline int_t inclusionGraph (const diGraph<wType> &g, int_t nVert,
        diGraph<wType> &result, TConn &connections)
{
        // remember the number of vertices in the input graph
        int_t nVertG = g. countVertices ();
        if (!nVertG)
                return 0;

        // mark each vertex as non-visited
        BitField visited;
        visited. allocate (nVertG);
        visited. clearall (nVertG);

        // create the sets of vertices reachable from each vertex
        BitSets forwardlists (nVertG, nVert);

        // prepare stacks for the DFS algorithm
        std::stack<int_t> s_vertex;
        std::stack<int_t> s_edge;

        // use DFS to propagate the reachable sets
        int_t startVertex = 0;
        int_t vertex = 0;
        int_t edge = 0;
        visited. set (vertex);
        while (1)
        {
                // go back with the recursion
                // if all the edges have been processed
                if (edge >= g. countEdges (vertex))
                {
                        // if this was the top recursion level,
                        // then find another starting point or quit
                        if (s_vertex. empty ())
                        {
                                while ((startVertex < nVertG) &&
                                        (visited. test (startVertex)))
                                        ++ startVertex;
                                if (startVertex >= nVertG)
                                        break;
                                vertex = startVertex;
                                edge = 0;
                                visited. set (vertex);
                                continue;
                        }

                        // determine the previous vertex (to trace back to)
                        int_t prev = s_vertex. top ();

                        // add the list of the current vertex
                        // to the list of the previous vertex
                        // unless that was one of the first 'nVert' vertices
                        if (vertex >= nVert)
                                forwardlists. addset (prev, vertex);

                        // return from the recursion
                        vertex = prev;
                        s_vertex. pop ();
                        edge = s_edge. top ();
                        s_edge. pop ();
                        continue;
                }

                // take the next vertex
                int_t next = g. getEdge (vertex, edge ++);

                // add the number to the list if appropriate
                if (next < nVert)
                        forwardlists. add (vertex, next);

                // go to the deeper recursion level if vertex not visited
                if (!visited. test (next))
                {
                        // store the previous variables at the stacks
                        s_vertex. push (vertex);
                        s_edge. push (edge);

                        // take the new vertex and mark as visited
                        vertex = next;
                        edge = 0;
                        visited. set (vertex);
                }

                // add the list of the next vertex to the current one
                // if the vertex has already been visited before
                // and the vertex is not one of the first 'nVert' ones
                else if (next >= nVert)
                {
                        forwardlists. addset (vertex, next);
                }
        }

        // ------------------------------------------

        // create the sets of vertices that can reach each vertex
        BitSets backlists (nVertG, nVert);

        diGraph<wType> gT;
        g. transpose (gT);

        // do a simple test for debugging purposes
        if (!s_vertex. empty () || !s_edge. empty ())
                throw "DEBUG: Nonempty stacks in 'inclusionGraph'.";

        // mark all the vertices as unvisited for the backward run
        visited. clearall (nVertG);

        // use DFS to propagate the reachable sets
        startVertex = 0;
        vertex = 0;
        edge = 0;
        visited. set (vertex);
        while (1)
        {
                // go back with the recursion
                // if all the edges have been processed
                if (edge >= gT. countEdges (vertex))
                {
                        // if this was the top recursion level,
                        // then find another starting point or quit
                        if (s_vertex. empty ())
                        {
                                while ((startVertex < nVertG) &&
                                        (visited. test (startVertex)))
                                        ++ startVertex;
                                if (startVertex >= nVertG)
                                        break;
                                vertex = startVertex;
                                edge = 0;
                                visited. set (vertex);
                                continue;
                        }

                        // determine the previous vertex (to trace back to)
                        int_t prev = s_vertex. top ();

                        // add the list of the current vertex
                        // to the list of the previous vertex
                        // unless that was one of the first 'nVert' vertices
                        if (vertex >= nVert)
                                backlists. addset (prev, vertex);

                        // return from the recursion
                        vertex = prev;
                        s_vertex. pop ();
                        edge = s_edge. top ();
                        s_edge. pop ();
                        continue;
                }

                // take the next vertex
                int_t next = gT. getEdge (vertex, edge ++);

                // add the number to the list if appropriate
                if (next < nVert)
                        backlists. add (vertex, next);

                // go to the deeper recursion level if vertex not visited
                if (!visited. test (next))
                {
                        // store the previous variables at the stacks
                        s_vertex. push (vertex);
                        s_edge. push (edge);

                        // take the new vertex and mark as visited
                        vertex = next;
                        edge = 0;
                        visited. set (vertex);
                }

                // add the list of the next vertex to the current one
                // if the vertex has already been visited before
                // and the vertex is not one of the first 'nVert' ones
                else if (next >= nVert)
                {
                        backlists. addset (vertex, next);
                }
        }

        // ------------------------------------------

        // record the connections to the obtained data structure
        for (int_t v = nVert; v < nVertG; ++ v)
        {
                int_t bvertex = backlists. get (v);
                while (bvertex >= 0)
                {
                        int_t fvertex = forwardlists. get (v);
                        while (fvertex >= 0)
                        {
                                connections (bvertex, fvertex, v);
                                fvertex = forwardlists. get (v, fvertex + 1);
                        }
                        bvertex = backlists. get (v, bvertex + 1);
                }
        }

        // ------------------------------------------
        
        // create the result graph based on the lists of edges
        for (int_t vertex = 0; vertex < nVert; ++ vertex)
        {
                result. addVertex ();
                int_t edge = forwardlists. get (vertex);
                while (edge >= 0)
                {
                        result. addEdge (edge);
                        edge = forwardlists. get (vertex, edge + 1);
                }
        }

        // free memory and exit
        visited. free ();
        return 0;
} /* inclusionGraph */

// --------------------------------------------------

template <class wType, class matrix>
inline void graph2matrix (const diGraph<wType> &g, matrix &m)
{
        int_t nVert = g. countVertices ();
        for (int_t v = 0; v < nVert; ++ v)
        {
                int_t nEdges = g. countEdges (v);
                for (int_t e = 0; e < nEdges; ++ e)
                {
                        int_t w = g. getEdge (v, e);
                        m [v] [w] = 1;
                }
        }
        return;
} /* graph2matrix */

template <class wType, class matrix>
inline void matrix2graph (const matrix &m, int_t n, diGraph<wType> &g)
{
        for (int_t v = 0; v < n; ++ v)
        {
                g. addVertex ();
                for (int_t w = 0; w < n; ++ w)
                        if (m [v] [w])
                                g. addEdge (w);
        }
        return;
} /* matrix2graph */

// --------------------------------------------------

template <class matrix>
inline void transitiveClosure (matrix &m, int_t n)
{
        for (int_t k = 0; k < n; ++ k)
        {
                for (int_t i = 0; i < n; ++ i)
                {
                        for (int_t j = 0; j < n; ++ j)
                        {
                                if (m [i] [k] && m [k] [j])
                                        m [i] [j] = 1;
                        }
                }
        }
        return;
} /* transitiveClosure */

template <class matrix>
inline void transitiveReduction (matrix &m, int_t n)
{
        for (int_t k = n - 1; k >= 0; -- k)
        {
                for (int_t i = 0; i < n; ++ i)
                {
                        for (int_t j = 0; j < n; ++ j)
                        {
                                if (m [i] [k] && m [k] [j])
                                        m [i] [j] = 0;
                        }
                }
        }
        return;
} /* transitiveReduction */

template <class wType>
inline void transitiveReduction (const diGraph<wType> &g,
        diGraph<wType> &gRed)
{
        int_t nVert = g. countVertices ();
        if (nVert <= 0)
                return;
        flatMatrix<char> m (nVert);
        m. clear (0);
        graph2matrix (g, m);
        transitiveClosure (m, nVert);
        transitiveReduction (m, nVert);
        matrix2graph (m, nVert, gRed);
        return;
} /* transitiveReduction */


} // namespace homology
} // namespace chomp

#endif // _CHOMP_STRUCT_DIGRAPH_H_