lemon/unionfind.h
author Peter Kovacs <kpeter@inf.elte.hu>
Sat, 08 Jan 2011 22:49:09 +0100
changeset 1032 62ba43576f85
parent 875 07ec2b52e53d
child 1092 dceba191c00d
permissions -rw-r--r--
Doc group for TSP algorithms (#386)
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2010
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 #ifndef LEMON_UNION_FIND_H
    20 #define LEMON_UNION_FIND_H
    21 
    22 //!\ingroup auxdat
    23 //!\file
    24 //!\brief Union-Find data structures.
    25 //!
    26 
    27 #include <vector>
    28 #include <list>
    29 #include <utility>
    30 #include <algorithm>
    31 #include <functional>
    32 
    33 #include <lemon/core.h>
    34 
    35 namespace lemon {
    36 
    37   /// \ingroup auxdat
    38   ///
    39   /// \brief A \e Union-Find data structure implementation
    40   ///
    41   /// The class implements the \e Union-Find data structure.
    42   /// The union operation uses rank heuristic, while
    43   /// the find operation uses path compression.
    44   /// This is a very simple but efficient implementation, providing
    45   /// only four methods: join (union), find, insert and size.
    46   /// For more features, see the \ref UnionFindEnum class.
    47   ///
    48   /// It is primarily used in Kruskal algorithm for finding minimal
    49   /// cost spanning tree in a graph.
    50   /// \sa kruskal()
    51   ///
    52   /// \pre You need to add all the elements by the \ref insert()
    53   /// method.
    54   template <typename IM>
    55   class UnionFind {
    56   public:
    57 
    58     ///\e
    59     typedef IM ItemIntMap;
    60     ///\e
    61     typedef typename ItemIntMap::Key Item;
    62 
    63   private:
    64     // If the items vector stores negative value for an item then
    65     // that item is root item and it has -items[it] component size.
    66     // Else the items[it] contains the index of the parent.
    67     std::vector<int> items;
    68     ItemIntMap& index;
    69 
    70     bool rep(int idx) const {
    71       return items[idx] < 0;
    72     }
    73 
    74     int repIndex(int idx) const {
    75       int k = idx;
    76       while (!rep(k)) {
    77         k = items[k] ;
    78       }
    79       while (idx != k) {
    80         int next = items[idx];
    81         const_cast<int&>(items[idx]) = k;
    82         idx = next;
    83       }
    84       return k;
    85     }
    86 
    87   public:
    88 
    89     /// \brief Constructor
    90     ///
    91     /// Constructor of the UnionFind class. You should give an item to
    92     /// integer map which will be used from the data structure. If you
    93     /// modify directly this map that may cause segmentation fault,
    94     /// invalid data structure, or infinite loop when you use again
    95     /// the union-find.
    96     UnionFind(ItemIntMap& m) : index(m) {}
    97 
    98     /// \brief Returns the index of the element's component.
    99     ///
   100     /// The method returns the index of the element's component.
   101     /// This is an integer between zero and the number of inserted elements.
   102     ///
   103     int find(const Item& a) {
   104       return repIndex(index[a]);
   105     }
   106 
   107     /// \brief Clears the union-find data structure
   108     ///
   109     /// Erase each item from the data structure.
   110     void clear() {
   111       items.clear();
   112     }
   113 
   114     /// \brief Inserts a new element into the structure.
   115     ///
   116     /// This method inserts a new element into the data structure.
   117     ///
   118     /// The method returns the index of the new component.
   119     int insert(const Item& a) {
   120       int n = items.size();
   121       items.push_back(-1);
   122       index.set(a,n);
   123       return n;
   124     }
   125 
   126     /// \brief Joining the components of element \e a and element \e b.
   127     ///
   128     /// This is the \e union operation of the Union-Find structure.
   129     /// Joins the component of element \e a and component of
   130     /// element \e b. If \e a and \e b are in the same component then
   131     /// it returns false otherwise it returns true.
   132     bool join(const Item& a, const Item& b) {
   133       int ka = repIndex(index[a]);
   134       int kb = repIndex(index[b]);
   135 
   136       if ( ka == kb )
   137         return false;
   138 
   139       if (items[ka] < items[kb]) {
   140         items[ka] += items[kb];
   141         items[kb] = ka;
   142       } else {
   143         items[kb] += items[ka];
   144         items[ka] = kb;
   145       }
   146       return true;
   147     }
   148 
   149     /// \brief Returns the size of the component of element \e a.
   150     ///
   151     /// Returns the size of the component of element \e a.
   152     int size(const Item& a) {
   153       int k = repIndex(index[a]);
   154       return - items[k];
   155     }
   156 
   157   };
   158 
   159   /// \ingroup auxdat
   160   ///
   161   /// \brief A \e Union-Find data structure implementation which
   162   /// is able to enumerate the components.
   163   ///
   164   /// The class implements a \e Union-Find data structure
   165   /// which is able to enumerate the components and the items in
   166   /// a component. If you don't need this feature then perhaps it's
   167   /// better to use the \ref UnionFind class which is more efficient.
   168   ///
   169   /// The union operation uses rank heuristic, while
   170   /// the find operation uses path compression.
   171   ///
   172   /// \pre You need to add all the elements by the \ref insert()
   173   /// method.
   174   ///
   175   template <typename IM>
   176   class UnionFindEnum {
   177   public:
   178 
   179     ///\e
   180     typedef IM ItemIntMap;
   181     ///\e
   182     typedef typename ItemIntMap::Key Item;
   183 
   184   private:
   185 
   186     ItemIntMap& index;
   187 
   188     // If the parent stores negative value for an item then that item
   189     // is root item and it has ~(items[it].parent) component id.  Else
   190     // the items[it].parent contains the index of the parent.
   191     //
   192     // The \c next and \c prev provides the double-linked
   193     // cyclic list of one component's items.
   194     struct ItemT {
   195       int parent;
   196       Item item;
   197 
   198       int next, prev;
   199     };
   200 
   201     std::vector<ItemT> items;
   202     int firstFreeItem;
   203 
   204     struct ClassT {
   205       int size;
   206       int firstItem;
   207       int next, prev;
   208     };
   209 
   210     std::vector<ClassT> classes;
   211     int firstClass, firstFreeClass;
   212 
   213     int newClass() {
   214       if (firstFreeClass == -1) {
   215         int cdx = classes.size();
   216         classes.push_back(ClassT());
   217         return cdx;
   218       } else {
   219         int cdx = firstFreeClass;
   220         firstFreeClass = classes[firstFreeClass].next;
   221         return cdx;
   222       }
   223     }
   224 
   225     int newItem() {
   226       if (firstFreeItem == -1) {
   227         int idx = items.size();
   228         items.push_back(ItemT());
   229         return idx;
   230       } else {
   231         int idx = firstFreeItem;
   232         firstFreeItem = items[firstFreeItem].next;
   233         return idx;
   234       }
   235     }
   236 
   237 
   238     bool rep(int idx) const {
   239       return items[idx].parent < 0;
   240     }
   241 
   242     int repIndex(int idx) const {
   243       int k = idx;
   244       while (!rep(k)) {
   245         k = items[k].parent;
   246       }
   247       while (idx != k) {
   248         int next = items[idx].parent;
   249         const_cast<int&>(items[idx].parent) = k;
   250         idx = next;
   251       }
   252       return k;
   253     }
   254 
   255     int classIndex(int idx) const {
   256       return ~(items[repIndex(idx)].parent);
   257     }
   258 
   259     void singletonItem(int idx) {
   260       items[idx].next = idx;
   261       items[idx].prev = idx;
   262     }
   263 
   264     void laceItem(int idx, int rdx) {
   265       items[idx].prev = rdx;
   266       items[idx].next = items[rdx].next;
   267       items[items[rdx].next].prev = idx;
   268       items[rdx].next = idx;
   269     }
   270 
   271     void unlaceItem(int idx) {
   272       items[items[idx].prev].next = items[idx].next;
   273       items[items[idx].next].prev = items[idx].prev;
   274 
   275       items[idx].next = firstFreeItem;
   276       firstFreeItem = idx;
   277     }
   278 
   279     void spliceItems(int ak, int bk) {
   280       items[items[ak].prev].next = bk;
   281       items[items[bk].prev].next = ak;
   282       int tmp = items[ak].prev;
   283       items[ak].prev = items[bk].prev;
   284       items[bk].prev = tmp;
   285 
   286     }
   287 
   288     void laceClass(int cls) {
   289       if (firstClass != -1) {
   290         classes[firstClass].prev = cls;
   291       }
   292       classes[cls].next = firstClass;
   293       classes[cls].prev = -1;
   294       firstClass = cls;
   295     }
   296 
   297     void unlaceClass(int cls) {
   298       if (classes[cls].prev != -1) {
   299         classes[classes[cls].prev].next = classes[cls].next;
   300       } else {
   301         firstClass = classes[cls].next;
   302       }
   303       if (classes[cls].next != -1) {
   304         classes[classes[cls].next].prev = classes[cls].prev;
   305       }
   306 
   307       classes[cls].next = firstFreeClass;
   308       firstFreeClass = cls;
   309     }
   310 
   311   public:
   312 
   313     UnionFindEnum(ItemIntMap& _index)
   314       : index(_index), items(), firstFreeItem(-1),
   315         firstClass(-1), firstFreeClass(-1) {}
   316 
   317     /// \brief Inserts the given element into a new component.
   318     ///
   319     /// This method creates a new component consisting only of the
   320     /// given element.
   321     ///
   322     int insert(const Item& item) {
   323       int idx = newItem();
   324 
   325       index.set(item, idx);
   326 
   327       singletonItem(idx);
   328       items[idx].item = item;
   329 
   330       int cdx = newClass();
   331 
   332       items[idx].parent = ~cdx;
   333 
   334       laceClass(cdx);
   335       classes[cdx].size = 1;
   336       classes[cdx].firstItem = idx;
   337 
   338       firstClass = cdx;
   339 
   340       return cdx;
   341     }
   342 
   343     /// \brief Inserts the given element into the component of the others.
   344     ///
   345     /// This methods inserts the element \e a into the component of the
   346     /// element \e comp.
   347     void insert(const Item& item, int cls) {
   348       int rdx = classes[cls].firstItem;
   349       int idx = newItem();
   350 
   351       index.set(item, idx);
   352 
   353       laceItem(idx, rdx);
   354 
   355       items[idx].item = item;
   356       items[idx].parent = rdx;
   357 
   358       ++classes[~(items[rdx].parent)].size;
   359     }
   360 
   361     /// \brief Clears the union-find data structure
   362     ///
   363     /// Erase each item from the data structure.
   364     void clear() {
   365       items.clear();
   366       firstClass = -1;
   367       firstFreeItem = -1;
   368     }
   369 
   370     /// \brief Finds the component of the given element.
   371     ///
   372     /// The method returns the component id of the given element.
   373     int find(const Item &item) const {
   374       return ~(items[repIndex(index[item])].parent);
   375     }
   376 
   377     /// \brief Joining the component of element \e a and element \e b.
   378     ///
   379     /// This is the \e union operation of the Union-Find structure.
   380     /// Joins the component of element \e a and component of
   381     /// element \e b. If \e a and \e b are in the same component then
   382     /// returns -1 else returns the remaining class.
   383     int join(const Item& a, const Item& b) {
   384 
   385       int ak = repIndex(index[a]);
   386       int bk = repIndex(index[b]);
   387 
   388       if (ak == bk) {
   389         return -1;
   390       }
   391 
   392       int acx = ~(items[ak].parent);
   393       int bcx = ~(items[bk].parent);
   394 
   395       int rcx;
   396 
   397       if (classes[acx].size > classes[bcx].size) {
   398         classes[acx].size += classes[bcx].size;
   399         items[bk].parent = ak;
   400         unlaceClass(bcx);
   401         rcx = acx;
   402       } else {
   403         classes[bcx].size += classes[acx].size;
   404         items[ak].parent = bk;
   405         unlaceClass(acx);
   406         rcx = bcx;
   407       }
   408       spliceItems(ak, bk);
   409 
   410       return rcx;
   411     }
   412 
   413     /// \brief Returns the size of the class.
   414     ///
   415     /// Returns the size of the class.
   416     int size(int cls) const {
   417       return classes[cls].size;
   418     }
   419 
   420     /// \brief Splits up the component.
   421     ///
   422     /// Splitting the component into singleton components (component
   423     /// of size one).
   424     void split(int cls) {
   425       int fdx = classes[cls].firstItem;
   426       int idx = items[fdx].next;
   427       while (idx != fdx) {
   428         int next = items[idx].next;
   429 
   430         singletonItem(idx);
   431 
   432         int cdx = newClass();
   433         items[idx].parent = ~cdx;
   434 
   435         laceClass(cdx);
   436         classes[cdx].size = 1;
   437         classes[cdx].firstItem = idx;
   438 
   439         idx = next;
   440       }
   441 
   442       items[idx].prev = idx;
   443       items[idx].next = idx;
   444 
   445       classes[~(items[idx].parent)].size = 1;
   446 
   447     }
   448 
   449     /// \brief Removes the given element from the structure.
   450     ///
   451     /// Removes the element from its component and if the component becomes
   452     /// empty then removes that component from the component list.
   453     ///
   454     /// \warning It is an error to remove an element which is not in
   455     /// the structure.
   456     /// \warning This running time of this operation is proportional to the
   457     /// number of the items in this class.
   458     void erase(const Item& item) {
   459       int idx = index[item];
   460       int fdx = items[idx].next;
   461 
   462       int cdx = classIndex(idx);
   463       if (idx == fdx) {
   464         unlaceClass(cdx);
   465         items[idx].next = firstFreeItem;
   466         firstFreeItem = idx;
   467         return;
   468       } else {
   469         classes[cdx].firstItem = fdx;
   470         --classes[cdx].size;
   471         items[fdx].parent = ~cdx;
   472 
   473         unlaceItem(idx);
   474         idx = items[fdx].next;
   475         while (idx != fdx) {
   476           items[idx].parent = fdx;
   477           idx = items[idx].next;
   478         }
   479 
   480       }
   481 
   482     }
   483 
   484     /// \brief Gives back a representant item of the component.
   485     ///
   486     /// Gives back a representant item of the component.
   487     Item item(int cls) const {
   488       return items[classes[cls].firstItem].item;
   489     }
   490 
   491     /// \brief Removes the component of the given element from the structure.
   492     ///
   493     /// Removes the component of the given element from the structure.
   494     ///
   495     /// \warning It is an error to give an element which is not in the
   496     /// structure.
   497     void eraseClass(int cls) {
   498       int fdx = classes[cls].firstItem;
   499       unlaceClass(cls);
   500       items[items[fdx].prev].next = firstFreeItem;
   501       firstFreeItem = fdx;
   502     }
   503 
   504     /// \brief LEMON style iterator for the representant items.
   505     ///
   506     /// ClassIt is a lemon style iterator for the components. It iterates
   507     /// on the ids of the classes.
   508     class ClassIt {
   509     public:
   510       /// \brief Constructor of the iterator
   511       ///
   512       /// Constructor of the iterator
   513       ClassIt(const UnionFindEnum& ufe) : unionFind(&ufe) {
   514         cdx = unionFind->firstClass;
   515       }
   516 
   517       /// \brief Constructor to get invalid iterator
   518       ///
   519       /// Constructor to get invalid iterator
   520       ClassIt(Invalid) : unionFind(0), cdx(-1) {}
   521 
   522       /// \brief Increment operator
   523       ///
   524       /// It steps to the next representant item.
   525       ClassIt& operator++() {
   526         cdx = unionFind->classes[cdx].next;
   527         return *this;
   528       }
   529 
   530       /// \brief Conversion operator
   531       ///
   532       /// It converts the iterator to the current representant item.
   533       operator int() const {
   534         return cdx;
   535       }
   536 
   537       /// \brief Equality operator
   538       ///
   539       /// Equality operator
   540       bool operator==(const ClassIt& i) {
   541         return i.cdx == cdx;
   542       }
   543 
   544       /// \brief Inequality operator
   545       ///
   546       /// Inequality operator
   547       bool operator!=(const ClassIt& i) {
   548         return i.cdx != cdx;
   549       }
   550 
   551     private:
   552       const UnionFindEnum* unionFind;
   553       int cdx;
   554     };
   555 
   556     /// \brief LEMON style iterator for the items of a component.
   557     ///
   558     /// ClassIt is a lemon style iterator for the components. It iterates
   559     /// on the items of a class. By example if you want to iterate on
   560     /// each items of each classes then you may write the next code.
   561     ///\code
   562     /// for (ClassIt cit(ufe); cit != INVALID; ++cit) {
   563     ///   std::cout << "Class: ";
   564     ///   for (ItemIt iit(ufe, cit); iit != INVALID; ++iit) {
   565     ///     std::cout << toString(iit) << ' ' << std::endl;
   566     ///   }
   567     ///   std::cout << std::endl;
   568     /// }
   569     ///\endcode
   570     class ItemIt {
   571     public:
   572       /// \brief Constructor of the iterator
   573       ///
   574       /// Constructor of the iterator. The iterator iterates
   575       /// on the class of the \c item.
   576       ItemIt(const UnionFindEnum& ufe, int cls) : unionFind(&ufe) {
   577         fdx = idx = unionFind->classes[cls].firstItem;
   578       }
   579 
   580       /// \brief Constructor to get invalid iterator
   581       ///
   582       /// Constructor to get invalid iterator
   583       ItemIt(Invalid) : unionFind(0), idx(-1) {}
   584 
   585       /// \brief Increment operator
   586       ///
   587       /// It steps to the next item in the class.
   588       ItemIt& operator++() {
   589         idx = unionFind->items[idx].next;
   590         if (idx == fdx) idx = -1;
   591         return *this;
   592       }
   593 
   594       /// \brief Conversion operator
   595       ///
   596       /// It converts the iterator to the current item.
   597       operator const Item&() const {
   598         return unionFind->items[idx].item;
   599       }
   600 
   601       /// \brief Equality operator
   602       ///
   603       /// Equality operator
   604       bool operator==(const ItemIt& i) {
   605         return i.idx == idx;
   606       }
   607 
   608       /// \brief Inequality operator
   609       ///
   610       /// Inequality operator
   611       bool operator!=(const ItemIt& i) {
   612         return i.idx != idx;
   613       }
   614 
   615     private:
   616       const UnionFindEnum* unionFind;
   617       int idx, fdx;
   618     };
   619 
   620   };
   621 
   622   /// \ingroup auxdat
   623   ///
   624   /// \brief A \e Extend-Find data structure implementation which
   625   /// is able to enumerate the components.
   626   ///
   627   /// The class implements an \e Extend-Find data structure which is
   628   /// able to enumerate the components and the items in a
   629   /// component. The data structure is a simplification of the
   630   /// Union-Find structure, and it does not allow to merge two components.
   631   ///
   632   /// \pre You need to add all the elements by the \ref insert()
   633   /// method.
   634   template <typename IM>
   635   class ExtendFindEnum {
   636   public:
   637 
   638     ///\e
   639     typedef IM ItemIntMap;
   640     ///\e
   641     typedef typename ItemIntMap::Key Item;
   642 
   643   private:
   644 
   645     ItemIntMap& index;
   646 
   647     struct ItemT {
   648       int cls;
   649       Item item;
   650       int next, prev;
   651     };
   652 
   653     std::vector<ItemT> items;
   654     int firstFreeItem;
   655 
   656     struct ClassT {
   657       int firstItem;
   658       int next, prev;
   659     };
   660 
   661     std::vector<ClassT> classes;
   662 
   663     int firstClass, firstFreeClass;
   664 
   665     int newClass() {
   666       if (firstFreeClass != -1) {
   667         int cdx = firstFreeClass;
   668         firstFreeClass = classes[cdx].next;
   669         return cdx;
   670       } else {
   671         classes.push_back(ClassT());
   672         return classes.size() - 1;
   673       }
   674     }
   675 
   676     int newItem() {
   677       if (firstFreeItem != -1) {
   678         int idx = firstFreeItem;
   679         firstFreeItem = items[idx].next;
   680         return idx;
   681       } else {
   682         items.push_back(ItemT());
   683         return items.size() - 1;
   684       }
   685     }
   686 
   687   public:
   688 
   689     /// \brief Constructor
   690     ExtendFindEnum(ItemIntMap& _index)
   691       : index(_index), items(), firstFreeItem(-1),
   692         classes(), firstClass(-1), firstFreeClass(-1) {}
   693 
   694     /// \brief Inserts the given element into a new component.
   695     ///
   696     /// This method creates a new component consisting only of the
   697     /// given element.
   698     int insert(const Item& item) {
   699       int cdx = newClass();
   700       classes[cdx].prev = -1;
   701       classes[cdx].next = firstClass;
   702       if (firstClass != -1) {
   703         classes[firstClass].prev = cdx;
   704       }
   705       firstClass = cdx;
   706 
   707       int idx = newItem();
   708       items[idx].item = item;
   709       items[idx].cls = cdx;
   710       items[idx].prev = idx;
   711       items[idx].next = idx;
   712 
   713       classes[cdx].firstItem = idx;
   714 
   715       index.set(item, idx);
   716 
   717       return cdx;
   718     }
   719 
   720     /// \brief Inserts the given element into the given component.
   721     ///
   722     /// This methods inserts the element \e item a into the \e cls class.
   723     void insert(const Item& item, int cls) {
   724       int idx = newItem();
   725       int rdx = classes[cls].firstItem;
   726       items[idx].item = item;
   727       items[idx].cls = cls;
   728 
   729       items[idx].prev = rdx;
   730       items[idx].next = items[rdx].next;
   731       items[items[rdx].next].prev = idx;
   732       items[rdx].next = idx;
   733 
   734       index.set(item, idx);
   735     }
   736 
   737     /// \brief Clears the union-find data structure
   738     ///
   739     /// Erase each item from the data structure.
   740     void clear() {
   741       items.clear();
   742       classes.clear();
   743       firstClass = firstFreeClass = firstFreeItem = -1;
   744     }
   745 
   746     /// \brief Gives back the class of the \e item.
   747     ///
   748     /// Gives back the class of the \e item.
   749     int find(const Item &item) const {
   750       return items[index[item]].cls;
   751     }
   752 
   753     /// \brief Gives back a representant item of the component.
   754     ///
   755     /// Gives back a representant item of the component.
   756     Item item(int cls) const {
   757       return items[classes[cls].firstItem].item;
   758     }
   759 
   760     /// \brief Removes the given element from the structure.
   761     ///
   762     /// Removes the element from its component and if the component becomes
   763     /// empty then removes that component from the component list.
   764     ///
   765     /// \warning It is an error to remove an element which is not in
   766     /// the structure.
   767     void erase(const Item &item) {
   768       int idx = index[item];
   769       int cdx = items[idx].cls;
   770 
   771       if (idx == items[idx].next) {
   772         if (classes[cdx].prev != -1) {
   773           classes[classes[cdx].prev].next = classes[cdx].next;
   774         } else {
   775           firstClass = classes[cdx].next;
   776         }
   777         if (classes[cdx].next != -1) {
   778           classes[classes[cdx].next].prev = classes[cdx].prev;
   779         }
   780         classes[cdx].next = firstFreeClass;
   781         firstFreeClass = cdx;
   782       } else {
   783         classes[cdx].firstItem = items[idx].next;
   784         items[items[idx].next].prev = items[idx].prev;
   785         items[items[idx].prev].next = items[idx].next;
   786       }
   787       items[idx].next = firstFreeItem;
   788       firstFreeItem = idx;
   789 
   790     }
   791 
   792 
   793     /// \brief Removes the component of the given element from the structure.
   794     ///
   795     /// Removes the component of the given element from the structure.
   796     ///
   797     /// \warning It is an error to give an element which is not in the
   798     /// structure.
   799     void eraseClass(int cdx) {
   800       int idx = classes[cdx].firstItem;
   801       items[items[idx].prev].next = firstFreeItem;
   802       firstFreeItem = idx;
   803 
   804       if (classes[cdx].prev != -1) {
   805         classes[classes[cdx].prev].next = classes[cdx].next;
   806       } else {
   807         firstClass = classes[cdx].next;
   808       }
   809       if (classes[cdx].next != -1) {
   810         classes[classes[cdx].next].prev = classes[cdx].prev;
   811       }
   812       classes[cdx].next = firstFreeClass;
   813       firstFreeClass = cdx;
   814     }
   815 
   816     /// \brief LEMON style iterator for the classes.
   817     ///
   818     /// ClassIt is a lemon style iterator for the components. It iterates
   819     /// on the ids of classes.
   820     class ClassIt {
   821     public:
   822       /// \brief Constructor of the iterator
   823       ///
   824       /// Constructor of the iterator
   825       ClassIt(const ExtendFindEnum& ufe) : extendFind(&ufe) {
   826         cdx = extendFind->firstClass;
   827       }
   828 
   829       /// \brief Constructor to get invalid iterator
   830       ///
   831       /// Constructor to get invalid iterator
   832       ClassIt(Invalid) : extendFind(0), cdx(-1) {}
   833 
   834       /// \brief Increment operator
   835       ///
   836       /// It steps to the next representant item.
   837       ClassIt& operator++() {
   838         cdx = extendFind->classes[cdx].next;
   839         return *this;
   840       }
   841 
   842       /// \brief Conversion operator
   843       ///
   844       /// It converts the iterator to the current class id.
   845       operator int() const {
   846         return cdx;
   847       }
   848 
   849       /// \brief Equality operator
   850       ///
   851       /// Equality operator
   852       bool operator==(const ClassIt& i) {
   853         return i.cdx == cdx;
   854       }
   855 
   856       /// \brief Inequality operator
   857       ///
   858       /// Inequality operator
   859       bool operator!=(const ClassIt& i) {
   860         return i.cdx != cdx;
   861       }
   862 
   863     private:
   864       const ExtendFindEnum* extendFind;
   865       int cdx;
   866     };
   867 
   868     /// \brief LEMON style iterator for the items of a component.
   869     ///
   870     /// ClassIt is a lemon style iterator for the components. It iterates
   871     /// on the items of a class. By example if you want to iterate on
   872     /// each items of each classes then you may write the next code.
   873     ///\code
   874     /// for (ClassIt cit(ufe); cit != INVALID; ++cit) {
   875     ///   std::cout << "Class: ";
   876     ///   for (ItemIt iit(ufe, cit); iit != INVALID; ++iit) {
   877     ///     std::cout << toString(iit) << ' ' << std::endl;
   878     ///   }
   879     ///   std::cout << std::endl;
   880     /// }
   881     ///\endcode
   882     class ItemIt {
   883     public:
   884       /// \brief Constructor of the iterator
   885       ///
   886       /// Constructor of the iterator. The iterator iterates
   887       /// on the class of the \c item.
   888       ItemIt(const ExtendFindEnum& ufe, int cls) : extendFind(&ufe) {
   889         fdx = idx = extendFind->classes[cls].firstItem;
   890       }
   891 
   892       /// \brief Constructor to get invalid iterator
   893       ///
   894       /// Constructor to get invalid iterator
   895       ItemIt(Invalid) : extendFind(0), idx(-1) {}
   896 
   897       /// \brief Increment operator
   898       ///
   899       /// It steps to the next item in the class.
   900       ItemIt& operator++() {
   901         idx = extendFind->items[idx].next;
   902         if (fdx == idx) idx = -1;
   903         return *this;
   904       }
   905 
   906       /// \brief Conversion operator
   907       ///
   908       /// It converts the iterator to the current item.
   909       operator const Item&() const {
   910         return extendFind->items[idx].item;
   911       }
   912 
   913       /// \brief Equality operator
   914       ///
   915       /// Equality operator
   916       bool operator==(const ItemIt& i) {
   917         return i.idx == idx;
   918       }
   919 
   920       /// \brief Inequality operator
   921       ///
   922       /// Inequality operator
   923       bool operator!=(const ItemIt& i) {
   924         return i.idx != idx;
   925       }
   926 
   927     private:
   928       const ExtendFindEnum* extendFind;
   929       int idx, fdx;
   930     };
   931 
   932   };
   933 
   934   /// \ingroup auxdat
   935   ///
   936   /// \brief A \e Union-Find data structure implementation which
   937   /// is able to store a priority for each item and retrieve the minimum of
   938   /// each class.
   939   ///
   940   /// A \e Union-Find data structure implementation which is able to
   941   /// store a priority for each item and retrieve the minimum of each
   942   /// class. In addition, it supports the joining and splitting the
   943   /// components. If you don't need this feature then you makes
   944   /// better to use the \ref UnionFind class which is more efficient.
   945   ///
   946   /// The union-find data strcuture based on a (2, 16)-tree with a
   947   /// tournament minimum selection on the internal nodes. The insert
   948   /// operation takes O(1), the find, set, decrease and increase takes
   949   /// O(log(n)), where n is the number of nodes in the current
   950   /// component.  The complexity of join and split is O(log(n)*k),
   951   /// where n is the sum of the number of the nodes and k is the
   952   /// number of joined components or the number of the components
   953   /// after the split.
   954   ///
   955   /// \pre You need to add all the elements by the \ref insert()
   956   /// method.
   957   template <typename V, typename IM, typename Comp = std::less<V> >
   958   class HeapUnionFind {
   959   public:
   960 
   961     ///\e
   962     typedef V Value;
   963     ///\e
   964     typedef typename IM::Key Item;
   965     ///\e
   966     typedef IM ItemIntMap;
   967     ///\e
   968     typedef Comp Compare;
   969 
   970   private:
   971 
   972     static const int cmax = 16;
   973 
   974     ItemIntMap& index;
   975 
   976     struct ClassNode {
   977       int parent;
   978       int depth;
   979 
   980       int left, right;
   981       int next, prev;
   982     };
   983 
   984     int first_class;
   985     int first_free_class;
   986     std::vector<ClassNode> classes;
   987 
   988     int newClass() {
   989       if (first_free_class < 0) {
   990         int id = classes.size();
   991         classes.push_back(ClassNode());
   992         return id;
   993       } else {
   994         int id = first_free_class;
   995         first_free_class = classes[id].next;
   996         return id;
   997       }
   998     }
   999 
  1000     void deleteClass(int id) {
  1001       classes[id].next = first_free_class;
  1002       first_free_class = id;
  1003     }
  1004 
  1005     struct ItemNode {
  1006       int parent;
  1007       Item item;
  1008       Value prio;
  1009       int next, prev;
  1010       int left, right;
  1011       int size;
  1012     };
  1013 
  1014     int first_free_node;
  1015     std::vector<ItemNode> nodes;
  1016 
  1017     int newNode() {
  1018       if (first_free_node < 0) {
  1019         int id = nodes.size();
  1020         nodes.push_back(ItemNode());
  1021         return id;
  1022       } else {
  1023         int id = first_free_node;
  1024         first_free_node = nodes[id].next;
  1025         return id;
  1026       }
  1027     }
  1028 
  1029     void deleteNode(int id) {
  1030       nodes[id].next = first_free_node;
  1031       first_free_node = id;
  1032     }
  1033 
  1034     Comp comp;
  1035 
  1036     int findClass(int id) const {
  1037       int kd = id;
  1038       while (kd >= 0) {
  1039         kd = nodes[kd].parent;
  1040       }
  1041       return ~kd;
  1042     }
  1043 
  1044     int leftNode(int id) const {
  1045       int kd = ~(classes[id].parent);
  1046       for (int i = 0; i < classes[id].depth; ++i) {
  1047         kd = nodes[kd].left;
  1048       }
  1049       return kd;
  1050     }
  1051 
  1052     int nextNode(int id) const {
  1053       int depth = 0;
  1054       while (id >= 0 && nodes[id].next == -1) {
  1055         id = nodes[id].parent;
  1056         ++depth;
  1057       }
  1058       if (id < 0) {
  1059         return -1;
  1060       }
  1061       id = nodes[id].next;
  1062       while (depth--) {
  1063         id = nodes[id].left;
  1064       }
  1065       return id;
  1066     }
  1067 
  1068 
  1069     void setPrio(int id) {
  1070       int jd = nodes[id].left;
  1071       nodes[id].prio = nodes[jd].prio;
  1072       nodes[id].item = nodes[jd].item;
  1073       jd = nodes[jd].next;
  1074       while (jd != -1) {
  1075         if (comp(nodes[jd].prio, nodes[id].prio)) {
  1076           nodes[id].prio = nodes[jd].prio;
  1077           nodes[id].item = nodes[jd].item;
  1078         }
  1079         jd = nodes[jd].next;
  1080       }
  1081     }
  1082 
  1083     void push(int id, int jd) {
  1084       nodes[id].size = 1;
  1085       nodes[id].left = nodes[id].right = jd;
  1086       nodes[jd].next = nodes[jd].prev = -1;
  1087       nodes[jd].parent = id;
  1088     }
  1089 
  1090     void pushAfter(int id, int jd) {
  1091       int kd = nodes[id].parent;
  1092       if (nodes[id].next != -1) {
  1093         nodes[nodes[id].next].prev = jd;
  1094         if (kd >= 0) {
  1095           nodes[kd].size += 1;
  1096         }
  1097       } else {
  1098         if (kd >= 0) {
  1099           nodes[kd].right = jd;
  1100           nodes[kd].size += 1;
  1101         }
  1102       }
  1103       nodes[jd].next = nodes[id].next;
  1104       nodes[jd].prev = id;
  1105       nodes[id].next = jd;
  1106       nodes[jd].parent = kd;
  1107     }
  1108 
  1109     void pushRight(int id, int jd) {
  1110       nodes[id].size += 1;
  1111       nodes[jd].prev = nodes[id].right;
  1112       nodes[jd].next = -1;
  1113       nodes[nodes[id].right].next = jd;
  1114       nodes[id].right = jd;
  1115       nodes[jd].parent = id;
  1116     }
  1117 
  1118     void popRight(int id) {
  1119       nodes[id].size -= 1;
  1120       int jd = nodes[id].right;
  1121       nodes[nodes[jd].prev].next = -1;
  1122       nodes[id].right = nodes[jd].prev;
  1123     }
  1124 
  1125     void splice(int id, int jd) {
  1126       nodes[id].size += nodes[jd].size;
  1127       nodes[nodes[id].right].next = nodes[jd].left;
  1128       nodes[nodes[jd].left].prev = nodes[id].right;
  1129       int kd = nodes[jd].left;
  1130       while (kd != -1) {
  1131         nodes[kd].parent = id;
  1132         kd = nodes[kd].next;
  1133       }
  1134       nodes[id].right = nodes[jd].right;
  1135     }
  1136 
  1137     void split(int id, int jd) {
  1138       int kd = nodes[id].parent;
  1139       nodes[kd].right = nodes[id].prev;
  1140       nodes[nodes[id].prev].next = -1;
  1141 
  1142       nodes[jd].left = id;
  1143       nodes[id].prev = -1;
  1144       int num = 0;
  1145       while (id != -1) {
  1146         nodes[id].parent = jd;
  1147         nodes[jd].right = id;
  1148         id = nodes[id].next;
  1149         ++num;
  1150       }
  1151       nodes[kd].size -= num;
  1152       nodes[jd].size = num;
  1153     }
  1154 
  1155     void pushLeft(int id, int jd) {
  1156       nodes[id].size += 1;
  1157       nodes[jd].next = nodes[id].left;
  1158       nodes[jd].prev = -1;
  1159       nodes[nodes[id].left].prev = jd;
  1160       nodes[id].left = jd;
  1161       nodes[jd].parent = id;
  1162     }
  1163 
  1164     void popLeft(int id) {
  1165       nodes[id].size -= 1;
  1166       int jd = nodes[id].left;
  1167       nodes[nodes[jd].next].prev = -1;
  1168       nodes[id].left = nodes[jd].next;
  1169     }
  1170 
  1171     void repairLeft(int id) {
  1172       int jd = ~(classes[id].parent);
  1173       while (nodes[jd].left != -1) {
  1174         int kd = nodes[jd].left;
  1175         if (nodes[jd].size == 1) {
  1176           if (nodes[jd].parent < 0) {
  1177             classes[id].parent = ~kd;
  1178             classes[id].depth -= 1;
  1179             nodes[kd].parent = ~id;
  1180             deleteNode(jd);
  1181             jd = kd;
  1182           } else {
  1183             int pd = nodes[jd].parent;
  1184             if (nodes[nodes[jd].next].size < cmax) {
  1185               pushLeft(nodes[jd].next, nodes[jd].left);
  1186               if (less(jd, nodes[jd].next) ||
  1187                   nodes[jd].item == nodes[pd].item) {
  1188                 nodes[nodes[jd].next].prio = nodes[jd].prio;
  1189                 nodes[nodes[jd].next].item = nodes[jd].item;
  1190               }
  1191               popLeft(pd);
  1192               deleteNode(jd);
  1193               jd = pd;
  1194             } else {
  1195               int ld = nodes[nodes[jd].next].left;
  1196               popLeft(nodes[jd].next);
  1197               pushRight(jd, ld);
  1198               if (less(ld, nodes[jd].left) ||
  1199                   nodes[ld].item == nodes[pd].item) {
  1200                 nodes[jd].item = nodes[ld].item;
  1201                 nodes[jd].prio = nodes[ld].prio;
  1202               }
  1203               if (nodes[nodes[jd].next].item == nodes[ld].item) {
  1204                 setPrio(nodes[jd].next);
  1205               }
  1206               jd = nodes[jd].left;
  1207             }
  1208           }
  1209         } else {
  1210           jd = nodes[jd].left;
  1211         }
  1212       }
  1213     }
  1214 
  1215     void repairRight(int id) {
  1216       int jd = ~(classes[id].parent);
  1217       while (nodes[jd].right != -1) {
  1218         int kd = nodes[jd].right;
  1219         if (nodes[jd].size == 1) {
  1220           if (nodes[jd].parent < 0) {
  1221             classes[id].parent = ~kd;
  1222             classes[id].depth -= 1;
  1223             nodes[kd].parent = ~id;
  1224             deleteNode(jd);
  1225             jd = kd;
  1226           } else {
  1227             int pd = nodes[jd].parent;
  1228             if (nodes[nodes[jd].prev].size < cmax) {
  1229               pushRight(nodes[jd].prev, nodes[jd].right);
  1230               if (less(jd, nodes[jd].prev) ||
  1231                   nodes[jd].item == nodes[pd].item) {
  1232                 nodes[nodes[jd].prev].prio = nodes[jd].prio;
  1233                 nodes[nodes[jd].prev].item = nodes[jd].item;
  1234               }
  1235               popRight(pd);
  1236               deleteNode(jd);
  1237               jd = pd;
  1238             } else {
  1239               int ld = nodes[nodes[jd].prev].right;
  1240               popRight(nodes[jd].prev);
  1241               pushLeft(jd, ld);
  1242               if (less(ld, nodes[jd].right) ||
  1243                   nodes[ld].item == nodes[pd].item) {
  1244                 nodes[jd].item = nodes[ld].item;
  1245                 nodes[jd].prio = nodes[ld].prio;
  1246               }
  1247               if (nodes[nodes[jd].prev].item == nodes[ld].item) {
  1248                 setPrio(nodes[jd].prev);
  1249               }
  1250               jd = nodes[jd].right;
  1251             }
  1252           }
  1253         } else {
  1254           jd = nodes[jd].right;
  1255         }
  1256       }
  1257     }
  1258 
  1259 
  1260     bool less(int id, int jd) const {
  1261       return comp(nodes[id].prio, nodes[jd].prio);
  1262     }
  1263 
  1264   public:
  1265 
  1266     /// \brief Returns true when the given class is alive.
  1267     ///
  1268     /// Returns true when the given class is alive, ie. the class is
  1269     /// not nested into other class.
  1270     bool alive(int cls) const {
  1271       return classes[cls].parent < 0;
  1272     }
  1273 
  1274     /// \brief Returns true when the given class is trivial.
  1275     ///
  1276     /// Returns true when the given class is trivial, ie. the class
  1277     /// contains just one item directly.
  1278     bool trivial(int cls) const {
  1279       return classes[cls].left == -1;
  1280     }
  1281 
  1282     /// \brief Constructs the union-find.
  1283     ///
  1284     /// Constructs the union-find.
  1285     /// \brief _index The index map of the union-find. The data
  1286     /// structure uses internally for store references.
  1287     HeapUnionFind(ItemIntMap& _index)
  1288       : index(_index), first_class(-1),
  1289         first_free_class(-1), first_free_node(-1) {}
  1290 
  1291     /// \brief Clears the union-find data structure
  1292     ///
  1293     /// Erase each item from the data structure.
  1294     void clear() {
  1295       nodes.clear();
  1296       classes.clear();
  1297       first_free_node = first_free_class = first_class = -1;
  1298     }
  1299 
  1300     /// \brief Insert a new node into a new component.
  1301     ///
  1302     /// Insert a new node into a new component.
  1303     /// \param item The item of the new node.
  1304     /// \param prio The priority of the new node.
  1305     /// \return The class id of the one-item-heap.
  1306     int insert(const Item& item, const Value& prio) {
  1307       int id = newNode();
  1308       nodes[id].item = item;
  1309       nodes[id].prio = prio;
  1310       nodes[id].size = 0;
  1311 
  1312       nodes[id].prev = -1;
  1313       nodes[id].next = -1;
  1314 
  1315       nodes[id].left = -1;
  1316       nodes[id].right = -1;
  1317 
  1318       nodes[id].item = item;
  1319       index[item] = id;
  1320 
  1321       int class_id = newClass();
  1322       classes[class_id].parent = ~id;
  1323       classes[class_id].depth = 0;
  1324 
  1325       classes[class_id].left = -1;
  1326       classes[class_id].right = -1;
  1327 
  1328       if (first_class != -1) {
  1329         classes[first_class].prev = class_id;
  1330       }
  1331       classes[class_id].next = first_class;
  1332       classes[class_id].prev = -1;
  1333       first_class = class_id;
  1334 
  1335       nodes[id].parent = ~class_id;
  1336 
  1337       return class_id;
  1338     }
  1339 
  1340     /// \brief The class of the item.
  1341     ///
  1342     /// \return The alive class id of the item, which is not nested into
  1343     /// other classes.
  1344     ///
  1345     /// The time complexity is O(log(n)).
  1346     int find(const Item& item) const {
  1347       return findClass(index[item]);
  1348     }
  1349 
  1350     /// \brief Joins the classes.
  1351     ///
  1352     /// The current function joins the given classes. The parameter is
  1353     /// an STL range which should be contains valid class ids. The
  1354     /// time complexity is O(log(n)*k) where n is the overall number
  1355     /// of the joined nodes and k is the number of classes.
  1356     /// \return The class of the joined classes.
  1357     /// \pre The range should contain at least two class ids.
  1358     template <typename Iterator>
  1359     int join(Iterator begin, Iterator end) {
  1360       std::vector<int> cs;
  1361       for (Iterator it = begin; it != end; ++it) {
  1362         cs.push_back(*it);
  1363       }
  1364 
  1365       int class_id = newClass();
  1366       { // creation union-find
  1367 
  1368         if (first_class != -1) {
  1369           classes[first_class].prev = class_id;
  1370         }
  1371         classes[class_id].next = first_class;
  1372         classes[class_id].prev = -1;
  1373         first_class = class_id;
  1374 
  1375         classes[class_id].depth = classes[cs[0]].depth;
  1376         classes[class_id].parent = classes[cs[0]].parent;
  1377         nodes[~(classes[class_id].parent)].parent = ~class_id;
  1378 
  1379         int l = cs[0];
  1380 
  1381         classes[class_id].left = l;
  1382         classes[class_id].right = l;
  1383 
  1384         if (classes[l].next != -1) {
  1385           classes[classes[l].next].prev = classes[l].prev;
  1386         }
  1387         classes[classes[l].prev].next = classes[l].next;
  1388 
  1389         classes[l].prev = -1;
  1390         classes[l].next = -1;
  1391 
  1392         classes[l].depth = leftNode(l);
  1393         classes[l].parent = class_id;
  1394 
  1395       }
  1396 
  1397       { // merging of heap
  1398         int l = class_id;
  1399         for (int ci = 1; ci < int(cs.size()); ++ci) {
  1400           int r = cs[ci];
  1401           int rln = leftNode(r);
  1402           if (classes[l].depth > classes[r].depth) {
  1403             int id = ~(classes[l].parent);
  1404             for (int i = classes[r].depth + 1; i < classes[l].depth; ++i) {
  1405               id = nodes[id].right;
  1406             }
  1407             while (id >= 0 && nodes[id].size == cmax) {
  1408               int new_id = newNode();
  1409               int right_id = nodes[id].right;
  1410 
  1411               popRight(id);
  1412               if (nodes[id].item == nodes[right_id].item) {
  1413                 setPrio(id);
  1414               }
  1415               push(new_id, right_id);
  1416               pushRight(new_id, ~(classes[r].parent));
  1417 
  1418               if (less(~classes[r].parent, right_id)) {
  1419                 nodes[new_id].item = nodes[~classes[r].parent].item;
  1420                 nodes[new_id].prio = nodes[~classes[r].parent].prio;
  1421               } else {
  1422                 nodes[new_id].item = nodes[right_id].item;
  1423                 nodes[new_id].prio = nodes[right_id].prio;
  1424               }
  1425 
  1426               id = nodes[id].parent;
  1427               classes[r].parent = ~new_id;
  1428             }
  1429             if (id < 0) {
  1430               int new_parent = newNode();
  1431               nodes[new_parent].next = -1;
  1432               nodes[new_parent].prev = -1;
  1433               nodes[new_parent].parent = ~l;
  1434 
  1435               push(new_parent, ~(classes[l].parent));
  1436               pushRight(new_parent, ~(classes[r].parent));
  1437               setPrio(new_parent);
  1438 
  1439               classes[l].parent = ~new_parent;
  1440               classes[l].depth += 1;
  1441             } else {
  1442               pushRight(id, ~(classes[r].parent));
  1443               while (id >= 0 && less(~(classes[r].parent), id)) {
  1444                 nodes[id].prio = nodes[~(classes[r].parent)].prio;
  1445                 nodes[id].item = nodes[~(classes[r].parent)].item;
  1446                 id = nodes[id].parent;
  1447               }
  1448             }
  1449           } else if (classes[r].depth > classes[l].depth) {
  1450             int id = ~(classes[r].parent);
  1451             for (int i = classes[l].depth + 1; i < classes[r].depth; ++i) {
  1452               id = nodes[id].left;
  1453             }
  1454             while (id >= 0 && nodes[id].size == cmax) {
  1455               int new_id = newNode();
  1456               int left_id = nodes[id].left;
  1457 
  1458               popLeft(id);
  1459               if (nodes[id].prio == nodes[left_id].prio) {
  1460                 setPrio(id);
  1461               }
  1462               push(new_id, left_id);
  1463               pushLeft(new_id, ~(classes[l].parent));
  1464 
  1465               if (less(~classes[l].parent, left_id)) {
  1466                 nodes[new_id].item = nodes[~classes[l].parent].item;
  1467                 nodes[new_id].prio = nodes[~classes[l].parent].prio;
  1468               } else {
  1469                 nodes[new_id].item = nodes[left_id].item;
  1470                 nodes[new_id].prio = nodes[left_id].prio;
  1471               }
  1472 
  1473               id = nodes[id].parent;
  1474               classes[l].parent = ~new_id;
  1475 
  1476             }
  1477             if (id < 0) {
  1478               int new_parent = newNode();
  1479               nodes[new_parent].next = -1;
  1480               nodes[new_parent].prev = -1;
  1481               nodes[new_parent].parent = ~l;
  1482 
  1483               push(new_parent, ~(classes[r].parent));
  1484               pushLeft(new_parent, ~(classes[l].parent));
  1485               setPrio(new_parent);
  1486 
  1487               classes[r].parent = ~new_parent;
  1488               classes[r].depth += 1;
  1489             } else {
  1490               pushLeft(id, ~(classes[l].parent));
  1491               while (id >= 0 && less(~(classes[l].parent), id)) {
  1492                 nodes[id].prio = nodes[~(classes[l].parent)].prio;
  1493                 nodes[id].item = nodes[~(classes[l].parent)].item;
  1494                 id = nodes[id].parent;
  1495               }
  1496             }
  1497             nodes[~(classes[r].parent)].parent = ~l;
  1498             classes[l].parent = classes[r].parent;
  1499             classes[l].depth = classes[r].depth;
  1500           } else {
  1501             if (classes[l].depth != 0 &&
  1502                 nodes[~(classes[l].parent)].size +
  1503                 nodes[~(classes[r].parent)].size <= cmax) {
  1504               splice(~(classes[l].parent), ~(classes[r].parent));
  1505               deleteNode(~(classes[r].parent));
  1506               if (less(~(classes[r].parent), ~(classes[l].parent))) {
  1507                 nodes[~(classes[l].parent)].prio =
  1508                   nodes[~(classes[r].parent)].prio;
  1509                 nodes[~(classes[l].parent)].item =
  1510                   nodes[~(classes[r].parent)].item;
  1511               }
  1512             } else {
  1513               int new_parent = newNode();
  1514               nodes[new_parent].next = nodes[new_parent].prev = -1;
  1515               push(new_parent, ~(classes[l].parent));
  1516               pushRight(new_parent, ~(classes[r].parent));
  1517               setPrio(new_parent);
  1518 
  1519               classes[l].parent = ~new_parent;
  1520               classes[l].depth += 1;
  1521               nodes[new_parent].parent = ~l;
  1522             }
  1523           }
  1524           if (classes[r].next != -1) {
  1525             classes[classes[r].next].prev = classes[r].prev;
  1526           }
  1527           classes[classes[r].prev].next = classes[r].next;
  1528 
  1529           classes[r].prev = classes[l].right;
  1530           classes[classes[l].right].next = r;
  1531           classes[l].right = r;
  1532           classes[r].parent = l;
  1533 
  1534           classes[r].next = -1;
  1535           classes[r].depth = rln;
  1536         }
  1537       }
  1538       return class_id;
  1539     }
  1540 
  1541     /// \brief Split the class to subclasses.
  1542     ///
  1543     /// The current function splits the given class. The join, which
  1544     /// made the current class, stored a reference to the
  1545     /// subclasses. The \c splitClass() member restores the classes
  1546     /// and creates the heaps. The parameter is an STL output iterator
  1547     /// which will be filled with the subclass ids. The time
  1548     /// complexity is O(log(n)*k) where n is the overall number of
  1549     /// nodes in the splitted classes and k is the number of the
  1550     /// classes.
  1551     template <typename Iterator>
  1552     void split(int cls, Iterator out) {
  1553       std::vector<int> cs;
  1554       { // splitting union-find
  1555         int id = cls;
  1556         int l = classes[id].left;
  1557 
  1558         classes[l].parent = classes[id].parent;
  1559         classes[l].depth = classes[id].depth;
  1560 
  1561         nodes[~(classes[l].parent)].parent = ~l;
  1562 
  1563         *out++ = l;
  1564 
  1565         while (l != -1) {
  1566           cs.push_back(l);
  1567           l = classes[l].next;
  1568         }
  1569 
  1570         classes[classes[id].right].next = first_class;
  1571         classes[first_class].prev = classes[id].right;
  1572         first_class = classes[id].left;
  1573 
  1574         if (classes[id].next != -1) {
  1575           classes[classes[id].next].prev = classes[id].prev;
  1576         }
  1577         classes[classes[id].prev].next = classes[id].next;
  1578 
  1579         deleteClass(id);
  1580       }
  1581 
  1582       {
  1583         for (int i = 1; i < int(cs.size()); ++i) {
  1584           int l = classes[cs[i]].depth;
  1585           while (nodes[nodes[l].parent].left == l) {
  1586             l = nodes[l].parent;
  1587           }
  1588           int r = l;
  1589           while (nodes[l].parent >= 0) {
  1590             l = nodes[l].parent;
  1591             int new_node = newNode();
  1592 
  1593             nodes[new_node].prev = -1;
  1594             nodes[new_node].next = -1;
  1595 
  1596             split(r, new_node);
  1597             pushAfter(l, new_node);
  1598             setPrio(l);
  1599             setPrio(new_node);
  1600             r = new_node;
  1601           }
  1602           classes[cs[i]].parent = ~r;
  1603           classes[cs[i]].depth = classes[~(nodes[l].parent)].depth;
  1604           nodes[r].parent = ~cs[i];
  1605 
  1606           nodes[l].next = -1;
  1607           nodes[r].prev = -1;
  1608 
  1609           repairRight(~(nodes[l].parent));
  1610           repairLeft(cs[i]);
  1611 
  1612           *out++ = cs[i];
  1613         }
  1614       }
  1615     }
  1616 
  1617     /// \brief Gives back the priority of the current item.
  1618     ///
  1619     /// Gives back the priority of the current item.
  1620     const Value& operator[](const Item& item) const {
  1621       return nodes[index[item]].prio;
  1622     }
  1623 
  1624     /// \brief Sets the priority of the current item.
  1625     ///
  1626     /// Sets the priority of the current item.
  1627     void set(const Item& item, const Value& prio) {
  1628       if (comp(prio, nodes[index[item]].prio)) {
  1629         decrease(item, prio);
  1630       } else if (!comp(prio, nodes[index[item]].prio)) {
  1631         increase(item, prio);
  1632       }
  1633     }
  1634 
  1635     /// \brief Increase the priority of the current item.
  1636     ///
  1637     /// Increase the priority of the current item.
  1638     void increase(const Item& item, const Value& prio) {
  1639       int id = index[item];
  1640       int kd = nodes[id].parent;
  1641       nodes[id].prio = prio;
  1642       while (kd >= 0 && nodes[kd].item == item) {
  1643         setPrio(kd);
  1644         kd = nodes[kd].parent;
  1645       }
  1646     }
  1647 
  1648     /// \brief Increase the priority of the current item.
  1649     ///
  1650     /// Increase the priority of the current item.
  1651     void decrease(const Item& item, const Value& prio) {
  1652       int id = index[item];
  1653       int kd = nodes[id].parent;
  1654       nodes[id].prio = prio;
  1655       while (kd >= 0 && less(id, kd)) {
  1656         nodes[kd].prio = prio;
  1657         nodes[kd].item = item;
  1658         kd = nodes[kd].parent;
  1659       }
  1660     }
  1661 
  1662     /// \brief Gives back the minimum priority of the class.
  1663     ///
  1664     /// Gives back the minimum priority of the class.
  1665     const Value& classPrio(int cls) const {
  1666       return nodes[~(classes[cls].parent)].prio;
  1667     }
  1668 
  1669     /// \brief Gives back the minimum priority item of the class.
  1670     ///
  1671     /// \return Gives back the minimum priority item of the class.
  1672     const Item& classTop(int cls) const {
  1673       return nodes[~(classes[cls].parent)].item;
  1674     }
  1675 
  1676     /// \brief Gives back a representant item of the class.
  1677     ///
  1678     /// Gives back a representant item of the class.
  1679     /// The representant is indpendent from the priorities of the
  1680     /// items.
  1681     const Item& classRep(int id) const {
  1682       int parent = classes[id].parent;
  1683       return nodes[parent >= 0 ? classes[id].depth : leftNode(id)].item;
  1684     }
  1685 
  1686     /// \brief LEMON style iterator for the items of a class.
  1687     ///
  1688     /// ClassIt is a lemon style iterator for the components. It iterates
  1689     /// on the items of a class. By example if you want to iterate on
  1690     /// each items of each classes then you may write the next code.
  1691     ///\code
  1692     /// for (ClassIt cit(huf); cit != INVALID; ++cit) {
  1693     ///   std::cout << "Class: ";
  1694     ///   for (ItemIt iit(huf, cit); iit != INVALID; ++iit) {
  1695     ///     std::cout << toString(iit) << ' ' << std::endl;
  1696     ///   }
  1697     ///   std::cout << std::endl;
  1698     /// }
  1699     ///\endcode
  1700     class ItemIt {
  1701     private:
  1702 
  1703       const HeapUnionFind* _huf;
  1704       int _id, _lid;
  1705 
  1706     public:
  1707 
  1708       /// \brief Default constructor
  1709       ///
  1710       /// Default constructor
  1711       ItemIt() {}
  1712 
  1713       ItemIt(const HeapUnionFind& huf, int cls) : _huf(&huf) {
  1714         int id = cls;
  1715         int parent = _huf->classes[id].parent;
  1716         if (parent >= 0) {
  1717           _id = _huf->classes[id].depth;
  1718           if (_huf->classes[id].next != -1) {
  1719             _lid = _huf->classes[_huf->classes[id].next].depth;
  1720           } else {
  1721             _lid = -1;
  1722           }
  1723         } else {
  1724           _id = _huf->leftNode(id);
  1725           _lid = -1;
  1726         }
  1727       }
  1728 
  1729       /// \brief Increment operator
  1730       ///
  1731       /// It steps to the next item in the class.
  1732       ItemIt& operator++() {
  1733         _id = _huf->nextNode(_id);
  1734         return *this;
  1735       }
  1736 
  1737       /// \brief Conversion operator
  1738       ///
  1739       /// It converts the iterator to the current item.
  1740       operator const Item&() const {
  1741         return _huf->nodes[_id].item;
  1742       }
  1743 
  1744       /// \brief Equality operator
  1745       ///
  1746       /// Equality operator
  1747       bool operator==(const ItemIt& i) {
  1748         return i._id == _id;
  1749       }
  1750 
  1751       /// \brief Inequality operator
  1752       ///
  1753       /// Inequality operator
  1754       bool operator!=(const ItemIt& i) {
  1755         return i._id != _id;
  1756       }
  1757 
  1758       /// \brief Equality operator
  1759       ///
  1760       /// Equality operator
  1761       bool operator==(Invalid) {
  1762         return _id == _lid;
  1763       }
  1764 
  1765       /// \brief Inequality operator
  1766       ///
  1767       /// Inequality operator
  1768       bool operator!=(Invalid) {
  1769         return _id != _lid;
  1770       }
  1771 
  1772     };
  1773 
  1774     /// \brief Class iterator
  1775     ///
  1776     /// The iterator stores
  1777     class ClassIt {
  1778     private:
  1779 
  1780       const HeapUnionFind* _huf;
  1781       int _id;
  1782 
  1783     public:
  1784 
  1785       ClassIt(const HeapUnionFind& huf)
  1786         : _huf(&huf), _id(huf.first_class) {}
  1787 
  1788       ClassIt(const HeapUnionFind& huf, int cls)
  1789         : _huf(&huf), _id(huf.classes[cls].left) {}
  1790 
  1791       ClassIt(Invalid) : _huf(0), _id(-1) {}
  1792 
  1793       const ClassIt& operator++() {
  1794         _id = _huf->classes[_id].next;
  1795         return *this;
  1796       }
  1797 
  1798       /// \brief Equality operator
  1799       ///
  1800       /// Equality operator
  1801       bool operator==(const ClassIt& i) {
  1802         return i._id == _id;
  1803       }
  1804 
  1805       /// \brief Inequality operator
  1806       ///
  1807       /// Inequality operator
  1808       bool operator!=(const ClassIt& i) {
  1809         return i._id != _id;
  1810       }
  1811 
  1812       operator int() const {
  1813         return _id;
  1814       }
  1815 
  1816     };
  1817 
  1818   };
  1819 
  1820   //! @}
  1821 
  1822 } //namespace lemon
  1823 
  1824 #endif //LEMON_UNION_FIND_H