/* -*- C++ -*-

  *

  * This file is a part of LEMON, a generic C++ optimization library

  *

  * Copyright (C) 2003-2008

  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

  * (Egervary Research Group on Combinatorial Optimization, EGRES).

  *

  * Permission to use, modify and distribute this software is granted

  * provided that this copyright notice appears in all copies. For

  * precise terms see the accompanying LICENSE file.

  *

  * This software is provided "AS IS" with no warranty of any kind,

  * express or implied, and with no claim as to its suitability for any

  * purpose.

  *

  */

 #ifndef LEMON_UNION_FIND_H

 #define LEMON_UNION_FIND_H

 //!\ingroup auxdat

 //!\file

 //!\brief Union-Find data structures.

 //!

 #include <vector>

 #include <list>

 #include <utility>

 #include <algorithm>

 #include <functional>

 #include <lemon/bits/invalid.h>

 namespace lemon {

   /// \ingroup auxdat

   ///

   /// \brief A \e Union-Find data structure implementation

   ///

   /// The class implements the \e Union-Find data structure. 

   /// The union operation uses rank heuristic, while

   /// the find operation uses path compression.

   /// This is a very simple but efficient implementation, providing 

   /// only four methods: join (union), find, insert and size.

   /// For more features see the \ref UnionFindEnum class.

   ///

   /// It is primarily used in Kruskal algorithm for finding minimal

   /// cost spanning tree in a graph.

   /// \sa kruskal()

   ///

   /// \pre You need to add all the elements by the \ref insert()

   /// method.  

   template <typename _ItemIntMap>

   class UnionFind {

   public:

     typedef _ItemIntMap ItemIntMap;

     typedef typename ItemIntMap::Key Item;

   private:

     // If the items vector stores negative value for an item then

     // that item is root item and it has -items[it] component size.

     // Else the items[it] contains the index of the parent.

     std::vector<int> items;

     ItemIntMap& index;

     bool rep(int idx) const {

       return items[idx] < 0;

     }

     int repIndex(int idx) const {

       int k = idx;

       while (!rep(k)) {

         k = items[k] ;

       }

       while (idx != k) {

         int next = items[idx];

         const_cast<int&>(items[idx]) = k;

         idx = next;

       }

       return k;

     }

   public:

     /// \brief Constructor

     ///

     /// Constructor of the UnionFind class. You should give an item to

     /// integer map which will be used from the data structure. If you

     /// modify directly this map that may cause segmentation fault,

     /// invalid data structure, or infinite loop when you use again

     /// the union-find.

     UnionFind(ItemIntMap& m) : index(m) {}

     /// \brief Returns the index of the element's component.

     ///

     /// The method returns the index of the element's component.

     /// This is an integer between zero and the number of inserted elements.

     ///

     int find(const Item& a) {

       return repIndex(index[a]);

     }

     /// \brief Clears the union-find data structure

     ///

     /// Erase each item from the data structure.

     void clear() {

       items.clear();

     }

     /// \brief Inserts a new element into the structure.

     ///

     /// This method inserts a new element into the data structure. 

     ///

     /// The method returns the index of the new component.

     int insert(const Item& a) {

       int n = items.size();

       items.push_back(-1);

       index.set(a,n);

       return n;

     }

     /// \brief Joining the components of element \e a and element \e b.

     ///

     /// This is the \e union operation of the Union-Find structure. 

     /// Joins the component of element \e a and component of

     /// element \e b. If \e a and \e b are in the same component then

     /// it returns false otherwise it returns true.

     bool join(const Item& a, const Item& b) {

       int ka = repIndex(index[a]);

       int kb = repIndex(index[b]);

       if ( ka == kb ) 

 	return false;

       if (items[ka] < items[kb]) {

 	items[ka] += items[kb];

 	items[kb] = ka;

       } else {

 	items[kb] += items[ka];

 	items[ka] = kb;

       }

       return true;

     }

     /// \brief Returns the size of the component of element \e a.

     ///

     /// Returns the size of the component of element \e a.

     int size(const Item& a) {

       int k = repIndex(index[a]);

       return - items[k];

     }

   };

   /// \ingroup auxdat

   ///

   /// \brief A \e Union-Find data structure implementation which

   /// is able to enumerate the components.

   ///

   /// The class implements a \e Union-Find data structure

   /// which is able to enumerate the components and the items in

   /// a component. If you don't need this feature then perhaps it's

   /// better to use the \ref UnionFind class which is more efficient.

   ///

   /// The union operation uses rank heuristic, while

   /// the find operation uses path compression.

   ///

   /// \pre You need to add all the elements by the \ref insert()

   /// method.

   ///

   template <typename _ItemIntMap>

   class UnionFindEnum {

   public:

     typedef _ItemIntMap ItemIntMap;

     typedef typename ItemIntMap::Key Item;

   private:

     ItemIntMap& index;

     // If the parent stores negative value for an item then that item

     // is root item and it has ~(items[it].parent) component id.  Else

     // the items[it].parent contains the index of the parent.

     //

     // The \c next and \c prev provides the double-linked

     // cyclic list of one component's items.

     struct ItemT {

       int parent;

       Item item;

       int next, prev;

     };

     std::vector<ItemT> items;

     int firstFreeItem;

     struct ClassT {

       int size;

       int firstItem;

       int next, prev;

     };

     std::vector<ClassT> classes;

     int firstClass, firstFreeClass;

     int newClass() {

       if (firstFreeClass == -1) {

 	int cdx = classes.size();

 	classes.push_back(ClassT());

 	return cdx;

       } else {

 	int cdx = firstFreeClass;

 	firstFreeClass = classes[firstFreeClass].next;

 	return cdx;

       }

     }

     int newItem() {

       if (firstFreeItem == -1) {

 	int idx = items.size();

 	items.push_back(ItemT());

 	return idx;

       } else {

 	int idx = firstFreeItem;

 	firstFreeItem = items[firstFreeItem].next;

 	return idx;

       }

     }

     bool rep(int idx) const {

       return items[idx].parent < 0;

     }

     int repIndex(int idx) const {

       int k = idx;

       while (!rep(k)) {

         k = items[k].parent;

       }

       while (idx != k) {

         int next = items[idx].parent;

         const_cast<int&>(items[idx].parent) = k;

         idx = next;

       }

       return k;

     }

     int classIndex(int idx) const {

       return ~(items[repIndex(idx)].parent);

     }

     void singletonItem(int idx) {

       items[idx].next = idx;

       items[idx].prev = idx;

     }

     void laceItem(int idx, int rdx) {

       items[idx].prev = rdx;

       items[idx].next = items[rdx].next;

       items[items[rdx].next].prev = idx;

       items[rdx].next = idx;

     }

     void unlaceItem(int idx) {

       items[items[idx].prev].next = items[idx].next;

       items[items[idx].next].prev = items[idx].prev;

       items[idx].next = firstFreeItem;

       firstFreeItem = idx;

     }

     void spliceItems(int ak, int bk) {

       items[items[ak].prev].next = bk;

       items[items[bk].prev].next = ak;

       int tmp = items[ak].prev;

       items[ak].prev = items[bk].prev;

       items[bk].prev = tmp;

     }

     void laceClass(int cls) {

       if (firstClass != -1) {

         classes[firstClass].prev = cls;

       }

       classes[cls].next = firstClass;

       classes[cls].prev = -1;

       firstClass = cls;

     } 

     void unlaceClass(int cls) {

       if (classes[cls].prev != -1) {

         classes[classes[cls].prev].next = classes[cls].next;

       } else {

         firstClass = classes[cls].next;

       }

       if (classes[cls].next != -1) {

         classes[classes[cls].next].prev = classes[cls].prev;

       }

       classes[cls].next = firstFreeClass;

       firstFreeClass = cls;

     } 

   public:

     UnionFindEnum(ItemIntMap& _index) 

       : index(_index), items(), firstFreeItem(-1), 

 	firstClass(-1), firstFreeClass(-1) {}

     /// \brief Inserts the given element into a new component.

     ///

     /// This method creates a new component consisting only of the

     /// given element.

     ///

     int insert(const Item& item) {

       int idx = newItem();

       index.set(item, idx);

       singletonItem(idx);

       items[idx].item = item;

       int cdx = newClass();

       items[idx].parent = ~cdx;

       laceClass(cdx);

       classes[cdx].size = 1;

       classes[cdx].firstItem = idx;

       firstClass = cdx;

       return cdx;

     }

     /// \brief Inserts the given element into the component of the others.

     ///

     /// This methods inserts the element \e a into the component of the

     /// element \e comp. 

     void insert(const Item& item, int cls) {

       int rdx = classes[cls].firstItem;

       int idx = newItem();

       index.set(item, idx);

       laceItem(idx, rdx);

       items[idx].item = item;

       items[idx].parent = rdx;

       ++classes[~(items[rdx].parent)].size;

     }

     /// \brief Clears the union-find data structure

     ///

     /// Erase each item from the data structure.

     void clear() {

       items.clear();

       firstClass = -1;

       firstFreeItem = -1;

     }

     /// \brief Finds the component of the given element.

     ///

     /// The method returns the component id of the given element.

     int find(const Item &item) const {

       return ~(items[repIndex(index[item])].parent);

     }

     /// \brief Joining the component of element \e a and element \e b.

     ///

     /// This is the \e union operation of the Union-Find structure. 

     /// Joins the component of element \e a and component of

     /// element \e b. If \e a and \e b are in the same component then

     /// returns -1 else returns the remaining class.

     int join(const Item& a, const Item& b) {

       int ak = repIndex(index[a]);

       int bk = repIndex(index[b]);

       if (ak == bk) {

 	return -1;

       }

       int acx = ~(items[ak].parent);

       int bcx = ~(items[bk].parent);

       int rcx;

       if (classes[acx].size > classes[bcx].size) {

 	classes[acx].size += classes[bcx].size;

 	items[bk].parent = ak;

         unlaceClass(bcx);

 	rcx = acx;

       } else {

 	classes[bcx].size += classes[acx].size;

 	items[ak].parent = bk;

         unlaceClass(acx);

 	rcx = bcx;

       }

       spliceItems(ak, bk);

       return rcx;

     }

     /// \brief Returns the size of the class.

     ///

     /// Returns the size of the class.

     int size(int cls) const {

       return classes[cls].size;

     }

     /// \brief Splits up the component. 

     ///

     /// Splitting the component into singleton components (component

     /// of size one).

     void split(int cls) {

       int fdx = classes[cls].firstItem;

       int idx = items[fdx].next;

       while (idx != fdx) {

         int next = items[idx].next;

 	singletonItem(idx);

 	int cdx = newClass();        

         items[idx].parent = ~cdx;

 	laceClass(cdx);

 	classes[cdx].size = 1;

 	classes[cdx].firstItem = idx;

         idx = next;

       }

       items[idx].prev = idx;

       items[idx].next = idx;

       classes[~(items[idx].parent)].size = 1;

     }

     /// \brief Removes the given element from the structure.

     ///

     /// Removes the element from its component and if the component becomes

     /// empty then removes that component from the component list.

     ///

     /// \warning It is an error to remove an element which is not in

     /// the structure.

     /// \warning This running time of this operation is proportional to the

     /// number of the items in this class.

     void erase(const Item& item) {

       int idx = index[item];

       int fdx = items[idx].next;

       int cdx = classIndex(idx);

       if (idx == fdx) {

 	unlaceClass(cdx);

 	items[idx].next = firstFreeItem;

 	firstFreeItem = idx;

 	return;

       } else {

 	classes[cdx].firstItem = fdx;

 	--classes[cdx].size;

 	items[fdx].parent = ~cdx;

 	unlaceItem(idx);

 	idx = items[fdx].next;

 	while (idx != fdx) {

 	  items[idx].parent = fdx;

 	  idx = items[idx].next;

 	}

       }

     }

     /// \brief Gives back a representant item of the component.

     ///

     /// Gives back a representant item of the component.

     Item item(int cls) const {

       return items[classes[cls].firstItem].item;

     }

     /// \brief Removes the component of the given element from the structure.

     ///

     /// Removes the component of the given element from the structure.

     ///

     /// \warning It is an error to give an element which is not in the

     /// structure.

     void eraseClass(int cls) {

       int fdx = classes[cls].firstItem;

       unlaceClass(cls);

       items[items[fdx].prev].next = firstFreeItem;

       firstFreeItem = fdx;

     }

     /// \brief Lemon style iterator for the representant items.

     ///

     /// ClassIt is a lemon style iterator for the components. It iterates

     /// on the ids of the classes.

     class ClassIt {

     public:

       /// \brief Constructor of the iterator

       ///

       /// Constructor of the iterator

       ClassIt(const UnionFindEnum& ufe) : unionFind(&ufe) {

         cdx = unionFind->firstClass;

       }

       /// \brief Constructor to get invalid iterator

       ///

       /// Constructor to get invalid iterator

       ClassIt(Invalid) : unionFind(0), cdx(-1) {}

       /// \brief Increment operator

       ///

       /// It steps to the next representant item.

       ClassIt& operator++() {

         cdx = unionFind->classes[cdx].next;

         return *this;

       }

       /// \brief Conversion operator

       ///

       /// It converts the iterator to the current representant item.

       operator int() const {

         return cdx;

       }

       /// \brief Equality operator

       ///

       /// Equality operator

       bool operator==(const ClassIt& i) { 

         return i.cdx == cdx;

       }

       /// \brief Inequality operator

       ///

       /// Inequality operator

       bool operator!=(const ClassIt& i) { 

         return i.cdx != cdx;

       }

     private:

       const UnionFindEnum* unionFind;

       int cdx;

     };

     /// \brief Lemon style iterator for the items of a component.

     ///

     /// ClassIt is a lemon style iterator for the components. It iterates

     /// on the items of a class. By example if you want to iterate on

     /// each items of each classes then you may write the next code.

     ///\code

     /// for (ClassIt cit(ufe); cit != INVALID; ++cit) {

     ///   std::cout << "Class: ";

     ///   for (ItemIt iit(ufe, cit); iit != INVALID; ++iit) {

     ///     std::cout << toString(iit) << ' ' << std::endl;

     ///   }

     ///   std::cout << std::endl;

     /// }

     ///\endcode

     class ItemIt {

     public:

       /// \brief Constructor of the iterator

       ///

       /// Constructor of the iterator. The iterator iterates

       /// on the class of the \c item.

       ItemIt(const UnionFindEnum& ufe, int cls) : unionFind(&ufe) {

         fdx = idx = unionFind->classes[cls].firstItem;

       }

       /// \brief Constructor to get invalid iterator

       ///

       /// Constructor to get invalid iterator

       ItemIt(Invalid) : unionFind(0), idx(-1) {}

       /// \brief Increment operator

       ///

       /// It steps to the next item in the class.

       ItemIt& operator++() {

         idx = unionFind->items[idx].next;

         if (idx == fdx) idx = -1;

         return *this;

       }

       /// \brief Conversion operator

       ///

       /// It converts the iterator to the current item.

       operator const Item&() const {

         return unionFind->items[idx].item;

       }

       /// \brief Equality operator

       ///

       /// Equality operator

       bool operator==(const ItemIt& i) { 

         return i.idx == idx;

       }

       /// \brief Inequality operator

       ///

       /// Inequality operator

       bool operator!=(const ItemIt& i) { 

         return i.idx != idx;

       }

     private:

       const UnionFindEnum* unionFind;

       int idx, fdx;

     };

   };

   /// \ingroup auxdat

   ///

   /// \brief A \e Extend-Find data structure implementation which

   /// is able to enumerate the components.

   ///

   /// The class implements an \e Extend-Find data structure which is

   /// able to enumerate the components and the items in a

   /// component. The data structure is a simplification of the

   /// Union-Find structure, and it does not allow to merge two components.

   ///

   /// \pre You need to add all the elements by the \ref insert()

   /// method.

   template <typename _ItemIntMap>

   class ExtendFindEnum {

   public:

     typedef _ItemIntMap ItemIntMap;

     typedef typename ItemIntMap::Key Item;

   private:

     ItemIntMap& index;

     struct ItemT {

       int cls;

       Item item;

       int next, prev;

     };

     std::vector<ItemT> items;

     int firstFreeItem;

     struct ClassT {

       int firstItem;

       int next, prev;

     };

     std::vector<ClassT> classes;

     int firstClass, firstFreeClass;

     int newClass() {

       if (firstFreeClass != -1) {

 	int cdx = firstFreeClass;

 	firstFreeClass = classes[cdx].next;

 	return cdx;

       } else {

 	classes.push_back(ClassT());

 	return classes.size() - 1;

       }

     }

     int newItem() {

       if (firstFreeItem != -1) {

 	int idx = firstFreeItem;

 	firstFreeItem = items[idx].next;

 	return idx;

       } else {

 	items.push_back(ItemT());

 	return items.size() - 1;

       }

     }

   public:

     /// \brief Constructor

     ExtendFindEnum(ItemIntMap& _index) 

       : index(_index), items(), firstFreeItem(-1), 

 	classes(), firstClass(-1), firstFreeClass(-1) {}

     /// \brief Inserts the given element into a new component.

     ///

     /// This method creates a new component consisting only of the

     /// given element.

     int insert(const Item& item) {

       int cdx = newClass();

       classes[cdx].prev = -1;

       classes[cdx].next = firstClass;

       if (firstClass != -1) {

 	classes[firstClass].prev = cdx;

       }

       firstClass = cdx;

       int idx = newItem();

       items[idx].item = item;

       items[idx].cls = cdx;

       items[idx].prev = idx;

       items[idx].next = idx;

       classes[cdx].firstItem = idx;

       index.set(item, idx);

       return cdx;

     }

     /// \brief Inserts the given element into the given component.

     ///

     /// This methods inserts the element \e item a into the \e cls class.

     void insert(const Item& item, int cls) {

       int idx = newItem();

       int rdx = classes[cls].firstItem;

       items[idx].item = item;

       items[idx].cls = cls;

       items[idx].prev = rdx;

       items[idx].next = items[rdx].next;

       items[items[rdx].next].prev = idx;

       items[rdx].next = idx;

       index.set(item, idx);

     }

     /// \brief Clears the union-find data structure

     ///

     /// Erase each item from the data structure.

     void clear() {

       items.clear();

       classes.clear;

       firstClass = firstFreeClass = firstFreeItem = -1;

     }

     /// \brief Gives back the class of the \e item.

     ///

     /// Gives back the class of the \e item.

     int find(const Item &item) const {

       return items[index[item]].cls;

     }

     /// \brief Gives back a representant item of the component.

     ///

     /// Gives back a representant item of the component.

     Item item(int cls) const {

       return items[classes[cls].firstItem].item;

     }

     /// \brief Removes the given element from the structure.

     ///

     /// Removes the element from its component and if the component becomes

     /// empty then removes that component from the component list.

     ///

     /// \warning It is an error to remove an element which is not in

     /// the structure.

     void erase(const Item &item) {

       int idx = index[item];

       int cdx = items[idx].cls;

       if (idx == items[idx].next) {

 	if (classes[cdx].prev != -1) {

 	  classes[classes[cdx].prev].next = classes[cdx].next; 

 	} else {

 	  firstClass = classes[cdx].next;

 	}

 	if (classes[cdx].next != -1) {

 	  classes[classes[cdx].next].prev = classes[cdx].prev; 

 	}

 	classes[cdx].next = firstFreeClass;

 	firstFreeClass = cdx;

       } else {

 	classes[cdx].firstItem = items[idx].next;

 	items[items[idx].next].prev = items[idx].prev;

 	items[items[idx].prev].next = items[idx].next;

       }

       items[idx].next = firstFreeItem;

       firstFreeItem = idx;

     }    

     /// \brief Removes the component of the given element from the structure.

     ///

     /// Removes the component of the given element from the structure.

     ///

     /// \warning It is an error to give an element which is not in the

     /// structure.

     void eraseClass(int cdx) {

       int idx = classes[cdx].firstItem;

       items[items[idx].prev].next = firstFreeItem;

       firstFreeItem = idx;

       if (classes[cdx].prev != -1) {

 	classes[classes[cdx].prev].next = classes[cdx].next; 

       } else {

 	firstClass = classes[cdx].next;

       }

       if (classes[cdx].next != -1) {

 	classes[classes[cdx].next].prev = classes[cdx].prev; 

       }

       classes[cdx].next = firstFreeClass;

       firstFreeClass = cdx;

     }

     /// \brief Lemon style iterator for the classes.

     ///

     /// ClassIt is a lemon style iterator for the components. It iterates

     /// on the ids of classes.

     class ClassIt {

     public:

       /// \brief Constructor of the iterator

       ///

       /// Constructor of the iterator

       ClassIt(const ExtendFindEnum& ufe) : extendFind(&ufe) {

         cdx = extendFind->firstClass;

       }

       /// \brief Constructor to get invalid iterator

       ///

       /// Constructor to get invalid iterator

       ClassIt(Invalid) : extendFind(0), cdx(-1) {}

       /// \brief Increment operator

       ///

       /// It steps to the next representant item.

       ClassIt& operator++() {

         cdx = extendFind->classes[cdx].next;

         return *this;

       }

       /// \brief Conversion operator

       ///

       /// It converts the iterator to the current class id.

       operator int() const {

         return cdx;

       }

       /// \brief Equality operator

       ///

       /// Equality operator

       bool operator==(const ClassIt& i) { 

         return i.cdx == cdx;

       }

       /// \brief Inequality operator

       ///

       /// Inequality operator

       bool operator!=(const ClassIt& i) { 

         return i.cdx != cdx;

       }

     private:

       const ExtendFindEnum* extendFind;

       int cdx;

     };

     /// \brief Lemon style iterator for the items of a component.

     ///

     /// ClassIt is a lemon style iterator for the components. It iterates

     /// on the items of a class. By example if you want to iterate on

     /// each items of each classes then you may write the next code.

     ///\code

     /// for (ClassIt cit(ufe); cit != INVALID; ++cit) {

     ///   std::cout << "Class: ";

     ///   for (ItemIt iit(ufe, cit); iit != INVALID; ++iit) {

     ///     std::cout << toString(iit) << ' ' << std::endl;

     ///   }

     ///   std::cout << std::endl;

     /// }

     ///\endcode

     class ItemIt {

     public:

       /// \brief Constructor of the iterator

       ///

       /// Constructor of the iterator. The iterator iterates

       /// on the class of the \c item.

       ItemIt(const ExtendFindEnum& ufe, int cls) : extendFind(&ufe) {

         fdx = idx = extendFind->classes[cls].firstItem;

       }

       /// \brief Constructor to get invalid iterator

       ///

       /// Constructor to get invalid iterator

       ItemIt(Invalid) : extendFind(0), idx(-1) {}

       /// \brief Increment operator

       ///

       /// It steps to the next item in the class.

       ItemIt& operator++() {

         idx = extendFind->items[idx].next;

 	if (fdx == idx) idx = -1;

         return *this;

       }

       /// \brief Conversion operator

       ///

       /// It converts the iterator to the current item.

       operator const Item&() const {

         return extendFind->items[idx].item;

       }

       /// \brief Equality operator

       ///

       /// Equality operator

       bool operator==(const ItemIt& i) { 

         return i.idx == idx;

       }

       /// \brief Inequality operator

       ///

       /// Inequality operator

       bool operator!=(const ItemIt& i) { 

         return i.idx != idx;

       }

     private:

       const ExtendFindEnum* extendFind;

       int idx, fdx;

     };

   };

   /// \ingroup auxdat

   ///

   /// \brief A \e Union-Find data structure implementation which

   /// is able to store a priority for each item and retrieve the minimum of

   /// each class.

   ///

   /// A \e Union-Find data structure implementation which is able to

   /// store a priority for each item and retrieve the minimum of each

   /// class. In addition, it supports the joining and splitting the

   /// components. If you don't need this feature then you makes

   /// better to use the \ref UnionFind class which is more efficient.

   ///

   /// The union-find data strcuture based on a (2, 16)-tree with a

   /// tournament minimum selection on the internal nodes. The insert

   /// operation takes O(1), the find, set, decrease and increase takes

   /// O(log(n)), where n is the number of nodes in the current

   /// component.  The complexity of join and split is O(log(n)*k),

   /// where n is the sum of the number of the nodes and k is the

   /// number of joined components or the number of the components

   /// after the split.

   ///

   /// \pre You need to add all the elements by the \ref insert()

   /// method.

   ///

   template <typename _Value, typename _ItemIntMap, 

             typename _Comp = std::less<_Value> >

   class HeapUnionFind {

   public:

     typedef _Value Value;

     typedef typename _ItemIntMap::Key Item;

     typedef _ItemIntMap ItemIntMap;

     typedef _Comp Comp;

   private:

     static const int cmax = 16;

     ItemIntMap& index;

     struct ClassNode {

       int parent;

       int depth;

       int left, right;

       int next, prev;

     };

     int first_class;

     int first_free_class;

     std::vector<ClassNode> classes;

     int newClass() {

       if (first_free_class < 0) {

         int id = classes.size();

         classes.push_back(ClassNode());

         return id;

       } else {

         int id = first_free_class;

         first_free_class = classes[id].next;

         return id;

       }

     }

     void deleteClass(int id) {

       classes[id].next = first_free_class;

       first_free_class = id;

     }

     struct ItemNode {

       int parent;

       Item item;

       Value prio;

       int next, prev;

       int left, right;

       int size;

     };

     int first_free_node;

     std::vector<ItemNode> nodes;

     int newNode() {

       if (first_free_node < 0) {

         int id = nodes.size();

         nodes.push_back(ItemNode());

         return id;

       } else {

         int id = first_free_node;

         first_free_node = nodes[id].next;

         return id;

       }

     }

     void deleteNode(int id) {

       nodes[id].next = first_free_node;

       first_free_node = id;

     }

     Comp comp;

     int findClass(int id) const {

       int kd = id;

       while (kd >= 0) {

         kd = nodes[kd].parent;

       }

       return ~kd;

     }

     int leftNode(int id) const {

       int kd = ~(classes[id].parent);

       for (int i = 0; i < classes[id].depth; ++i) {

         kd = nodes[kd].left;

       }

       return kd;

     }

     int nextNode(int id) const {

       int depth = 0;

       while (id >= 0 && nodes[id].next == -1) {

         id = nodes[id].parent;

         ++depth;

       }

       if (id < 0) {

         return -1;

       }

       id = nodes[id].next;

       while (depth--) {

         id = nodes[id].left;      

       }

       return id;

     }

     void setPrio(int id) {

       int jd = nodes[id].left;

       nodes[id].prio = nodes[jd].prio;

       nodes[id].item = nodes[jd].item;

       jd = nodes[jd].next;

       while (jd != -1) {

         if (comp(nodes[jd].prio, nodes[id].prio)) {

           nodes[id].prio = nodes[jd].prio;

           nodes[id].item = nodes[jd].item;

         }

         jd = nodes[jd].next;

       }

     }

     void push(int id, int jd) {

       nodes[id].size = 1;

       nodes[id].left = nodes[id].right = jd;

       nodes[jd].next = nodes[jd].prev = -1;

       nodes[jd].parent = id;

     }

     void pushAfter(int id, int jd) {

       int kd = nodes[id].parent;

       if (nodes[id].next != -1) {

         nodes[nodes[id].next].prev = jd;

         if (kd >= 0) {

           nodes[kd].size += 1;

         }

       } else {

         if (kd >= 0) {

           nodes[kd].right = jd;

           nodes[kd].size += 1;

         }

       }

       nodes[jd].next = nodes[id].next;

       nodes[jd].prev = id;

       nodes[id].next = jd;

       nodes[jd].parent = kd;

     }

     void pushRight(int id, int jd) {

       nodes[id].size += 1;

       nodes[jd].prev = nodes[id].right;

       nodes[jd].next = -1;

       nodes[nodes[id].right].next = jd;

       nodes[id].right = jd;

       nodes[jd].parent = id;

     }

     void popRight(int id) {

       nodes[id].size -= 1;

       int jd = nodes[id].right;

       nodes[nodes[jd].prev].next = -1;

       nodes[id].right = nodes[jd].prev;

     }

     void splice(int id, int jd) {

       nodes[id].size += nodes[jd].size;

       nodes[nodes[id].right].next = nodes[jd].left;

       nodes[nodes[jd].left].prev = nodes[id].right;

       int kd = nodes[jd].left;

       while (kd != -1) {

         nodes[kd].parent = id;

         kd = nodes[kd].next;

       }

       nodes[id].right = nodes[jd].right;

     }

     void split(int id, int jd) {

       int kd = nodes[id].parent;

       nodes[kd].right = nodes[id].prev;

       nodes[nodes[id].prev].next = -1;

       nodes[jd].left = id;

       nodes[id].prev = -1;

       int num = 0;

       while (id != -1) {

         nodes[id].parent = jd;

         nodes[jd].right = id;

         id = nodes[id].next;

         ++num;

       }      

       nodes[kd].size -= num;

       nodes[jd].size = num;

     }

     void pushLeft(int id, int jd) {

       nodes[id].size += 1;

       nodes[jd].next = nodes[id].left;

       nodes[jd].prev = -1;

       nodes[nodes[id].left].prev = jd;

       nodes[id].left = jd;

       nodes[jd].parent = id;

     }

     void popLeft(int id) {

       nodes[id].size -= 1;

       int jd = nodes[id].left;

       nodes[nodes[jd].next].prev = -1;

       nodes[id].left = nodes[jd].next;

     }

     void repairLeft(int id) {

       int jd = ~(classes[id].parent);

       while (nodes[jd].left != -1) {

 	int kd = nodes[jd].left;

 	if (nodes[jd].size == 1) {

 	  if (nodes[jd].parent < 0) {

 	    classes[id].parent = ~kd;

 	    classes[id].depth -= 1;

 	    nodes[kd].parent = ~id;

 	    deleteNode(jd);

 	    jd = kd;

 	  } else {

 	    int pd = nodes[jd].parent;

 	    if (nodes[nodes[jd].next].size < cmax) {

 	      pushLeft(nodes[jd].next, nodes[jd].left);

 	      if (less(nodes[jd].left, nodes[jd].next)) {

 		nodes[nodes[jd].next].prio = nodes[nodes[jd].left].prio;

 		nodes[nodes[jd].next].item = nodes[nodes[jd].left].item;

 	      } 

 	      popLeft(pd);

 	      deleteNode(jd);

 	      jd = pd;

 	    } else {

 	      int ld = nodes[nodes[jd].next].left;

 	      popLeft(nodes[jd].next);

 	      pushRight(jd, ld);

 	      if (less(ld, nodes[jd].left)) {

 		nodes[jd].item = nodes[ld].item;

 		nodes[jd].prio = nodes[jd].prio;

 	      }

 	      if (nodes[nodes[jd].next].item == nodes[ld].item) {

 		setPrio(nodes[jd].next);

 	      }

 	      jd = nodes[jd].left;

 	    }

 	  }

 	} else {

 	  jd = nodes[jd].left;

 	}

       }

     }    

     void repairRight(int id) {

       int jd = ~(classes[id].parent);

       while (nodes[jd].right != -1) {

 	int kd = nodes[jd].right;

 	if (nodes[jd].size == 1) {

 	  if (nodes[jd].parent < 0) {

 	    classes[id].parent = ~kd;

 	    classes[id].depth -= 1;

 	    nodes[kd].parent = ~id;

 	    deleteNode(jd);

 	    jd = kd;

 	  } else {

 	    int pd = nodes[jd].parent;

 	    if (nodes[nodes[jd].prev].size < cmax) {

 	      pushRight(nodes[jd].prev, nodes[jd].right);

 	      if (less(nodes[jd].right, nodes[jd].prev)) {

 		nodes[nodes[jd].prev].prio = nodes[nodes[jd].right].prio;

 		nodes[nodes[jd].prev].item = nodes[nodes[jd].right].item;

 	      } 

 	      popRight(pd);

 	      deleteNode(jd);

 	      jd = pd;

 	    } else {

 	      int ld = nodes[nodes[jd].prev].right;

 	      popRight(nodes[jd].prev);

 	      pushLeft(jd, ld);

 	      if (less(ld, nodes[jd].right)) {

 		nodes[jd].item = nodes[ld].item;

 		nodes[jd].prio = nodes[jd].prio;

 	      }

 	      if (nodes[nodes[jd].prev].item == nodes[ld].item) {

 		setPrio(nodes[jd].prev);

 	      }

 	      jd = nodes[jd].right;

 	    }

 	  }

 	} else {

 	  jd = nodes[jd].right;

 	}

       }

     }

     bool less(int id, int jd) const {

       return comp(nodes[id].prio, nodes[jd].prio);

     }

     bool equal(int id, int jd) const {

       return !less(id, jd) && !less(jd, id);

     }

   public:

     /// \brief Returns true when the given class is alive.

     ///

     /// Returns true when the given class is alive, ie. the class is

     /// not nested into other class.

     bool alive(int cls) const {

       return classes[cls].parent < 0;

     }

     /// \brief Returns true when the given class is trivial.

     ///

     /// Returns true when the given class is trivial, ie. the class

     /// contains just one item directly.

     bool trivial(int cls) const {

       return classes[cls].left == -1;

     }

     /// \brief Constructs the union-find.

     ///

     /// Constructs the union-find.  

     /// \brief _index The index map of the union-find. The data

     /// structure uses internally for store references.

     HeapUnionFind(ItemIntMap& _index) 

       : index(_index), first_class(-1), 

 	first_free_class(-1), first_free_node(-1) {}

     /// \brief Insert a new node into a new component.

     ///

     /// Insert a new node into a new component.

     /// \param item The item of the new node.

     /// \param prio The priority of the new node.

     /// \return The class id of the one-item-heap.

     int insert(const Item& item, const Value& prio) {

       int id = newNode();

       nodes[id].item = item;

       nodes[id].prio = prio;

       nodes[id].size = 0;

       nodes[id].prev = -1;

       nodes[id].next = -1;

       nodes[id].left = -1;

       nodes[id].right = -1;

       nodes[id].item = item;

       index[item] = id;

       int class_id = newClass();

       classes[class_id].parent = ~id;

       classes[class_id].depth = 0;

       classes[class_id].left = -1;

       classes[class_id].right = -1;

       if (first_class != -1) {

         classes[first_class].prev = class_id;

       }

       classes[class_id].next = first_class;

       classes[class_id].prev = -1;

       first_class = class_id;

       nodes[id].parent = ~class_id;

       return class_id;

     }

     /// \brief The class of the item.

     ///

     /// \return The alive class id of the item, which is not nested into

     /// other classes.

     ///

     /// The time complexity is O(log(n)).

     int find(const Item& item) const {

       return findClass(index[item]);

     }

     /// \brief Joins the classes.

     ///

     /// The current function joins the given classes. The parameter is

     /// an STL range which should be contains valid class ids. The

     /// time complexity is O(log(n)*k) where n is the overall number

     /// of the joined nodes and k is the number of classes.

     /// \return The class of the joined classes.

     /// \pre The range should contain at least two class ids.

     template <typename Iterator>

     int join(Iterator begin, Iterator end) {

       std::vector<int> cs;

       for (Iterator it = begin; it != end; ++it) {

         cs.push_back(*it);

       }

       int class_id = newClass();

       { // creation union-find

         if (first_class != -1) {

           classes[first_class].prev = class_id;

         }

         classes[class_id].next = first_class;

         classes[class_id].prev = -1;

         first_class = class_id;

         classes[class_id].depth = classes[cs[0]].depth;

         classes[class_id].parent = classes[cs[0]].parent;

         nodes[~(classes[class_id].parent)].parent = ~class_id;

         int l = cs[0];

         classes[class_id].left = l;

         classes[class_id].right = l;

         if (classes[l].next != -1) {

           classes[classes[l].next].prev = classes[l].prev;

         }

         classes[classes[l].prev].next = classes[l].next;

         classes[l].prev = -1;

         classes[l].next = -1;

         classes[l].depth = leftNode(l);

         classes[l].parent = class_id;

       }

       { // merging of heap

         int l = class_id;

         for (int ci = 1; ci < int(cs.size()); ++ci) {

           int r = cs[ci];

           int rln = leftNode(r);

           if (classes[l].depth > classes[r].depth) {

             int id = ~(classes[l].parent);

             for (int i = classes[r].depth + 1; i < classes[l].depth; ++i) {

               id = nodes[id].right;

             }

             while (id >= 0 && nodes[id].size == cmax) {

               int new_id = newNode();

               int right_id = nodes[id].right;

               popRight(id);

               if (nodes[id].item == nodes[right_id].item) {

                 setPrio(id);

               }

               push(new_id, right_id);

               pushRight(new_id, ~(classes[r].parent));

               setPrio(new_id);

               id = nodes[id].parent;

               classes[r].parent = ~new_id;

             }

             if (id < 0) {

               int new_parent = newNode();

               nodes[new_parent].next = -1;

               nodes[new_parent].prev = -1;

               nodes[new_parent].parent = ~l;

               push(new_parent, ~(classes[l].parent));

               pushRight(new_parent, ~(classes[r].parent));

               setPrio(new_parent);

               classes[l].parent = ~new_parent;

               classes[l].depth += 1;

             } else {

               pushRight(id, ~(classes[r].parent));

               while (id >= 0 && less(~(classes[r].parent), id)) {

                 nodes[id].prio = nodes[~(classes[r].parent)].prio;

                 nodes[id].item = nodes[~(classes[r].parent)].item;

                 id = nodes[id].parent;

               }

             }

           } else if (classes[r].depth > classes[l].depth) {

             int id = ~(classes[r].parent);

             for (int i = classes[l].depth + 1; i < classes[r].depth; ++i) {

               id = nodes[id].left;

             }

             while (id >= 0 && nodes[id].size == cmax) {

               int new_id = newNode();

               int left_id = nodes[id].left;

               popLeft(id);

               if (nodes[id].prio == nodes[left_id].prio) {

                 setPrio(id);

               }

               push(new_id, left_id);

               pushLeft(new_id, ~(classes[l].parent));

               setPrio(new_id);

               id = nodes[id].parent;

               classes[l].parent = ~new_id;

             }

             if (id < 0) {

               int new_parent = newNode();

               nodes[new_parent].next = -1;

               nodes[new_parent].prev = -1;

               nodes[new_parent].parent = ~l;

               push(new_parent, ~(classes[r].parent));

               pushLeft(new_parent, ~(classes[l].parent));

               setPrio(new_parent);

               classes[r].parent = ~new_parent;

               classes[r].depth += 1;

             } else {

               pushLeft(id, ~(classes[l].parent));

               while (id >= 0 && less(~(classes[l].parent), id)) {

                 nodes[id].prio = nodes[~(classes[l].parent)].prio;

                 nodes[id].item = nodes[~(classes[l].parent)].item;

                 id = nodes[id].parent;

               }

             }

             nodes[~(classes[r].parent)].parent = ~l;

             classes[l].parent = classes[r].parent;

             classes[l].depth = classes[r].depth;

           } else {

             if (classes[l].depth != 0 && 

                 nodes[~(classes[l].parent)].size + 

                 nodes[~(classes[r].parent)].size <= cmax) {

               splice(~(classes[l].parent), ~(classes[r].parent));

               deleteNode(~(classes[r].parent));

               if (less(~(classes[r].parent), ~(classes[l].parent))) {

                 nodes[~(classes[l].parent)].prio = 

                   nodes[~(classes[r].parent)].prio;

                 nodes[~(classes[l].parent)].item = 

                   nodes[~(classes[r].parent)].item;

               }

             } else {

               int new_parent = newNode();

               nodes[new_parent].next = nodes[new_parent].prev = -1;

               push(new_parent, ~(classes[l].parent));

               pushRight(new_parent, ~(classes[r].parent));

               setPrio(new_parent);

               classes[l].parent = ~new_parent;

               classes[l].depth += 1;

               nodes[new_parent].parent = ~l;

             }

           }

           if (classes[r].next != -1) {

             classes[classes[r].next].prev = classes[r].prev;

           }

           classes[classes[r].prev].next = classes[r].next;

           classes[r].prev = classes[l].right;

           classes[classes[l].right].next = r;

           classes[l].right = r;

           classes[r].parent = l;

           classes[r].next = -1;

           classes[r].depth = rln;

         }

       }

       return class_id;

     }

     /// \brief Split the class to subclasses.

     ///

     /// The current function splits the given class. The join, which

     /// made the current class, stored a reference to the

     /// subclasses. The \c splitClass() member restores the classes

     /// and creates the heaps. The parameter is an STL output iterator

     /// which will be filled with the subclass ids. The time

     /// complexity is O(log(n)*k) where n is the overall number of

     /// nodes in the splitted classes and k is the number of the

     /// classes.

     template <typename Iterator>

     void split(int cls, Iterator out) {

       std::vector<int> cs;

       { // splitting union-find

         int id = cls;

         int l = classes[id].left;

         classes[l].parent = classes[id].parent;

         classes[l].depth = classes[id].depth;

         nodes[~(classes[l].parent)].parent = ~l;

         *out++ = l;

         while (l != -1) {

           cs.push_back(l);

           l = classes[l].next;

         }

         classes[classes[id].right].next = first_class;

         classes[first_class].prev = classes[id].right;

         first_class = classes[id].left;

         if (classes[id].next != -1) {

           classes[classes[id].next].prev = classes[id].prev;

         }

         classes[classes[id].prev].next = classes[id].next;

         deleteClass(id);

       }

       {

         for (int i = 1; i < int(cs.size()); ++i) {

           int l = classes[cs[i]].depth;

           while (nodes[nodes[l].parent].left == l) {

             l = nodes[l].parent;

           }

           int r = l;    

           while (nodes[l].parent >= 0) {

             l = nodes[l].parent;

             int new_node = newNode();

             nodes[new_node].prev = -1;

             nodes[new_node].next = -1;

             split(r, new_node);

             pushAfter(l, new_node);

             setPrio(l);

             setPrio(new_node);

             r = new_node;

           }

           classes[cs[i]].parent = ~r;

           classes[cs[i]].depth = classes[~(nodes[l].parent)].depth;

           nodes[r].parent = ~cs[i];

           nodes[l].next = -1;

           nodes[r].prev = -1;

           repairRight(~(nodes[l].parent));

           repairLeft(cs[i]);

           *out++ = cs[i];

         }

       }

     }

     /// \brief Gives back the priority of the current item.

     ///

     /// \return Gives back the priority of the current item.

     const Value& operator[](const Item& item) const {

       return nodes[index[item]].prio;

     }

     /// \brief Sets the priority of the current item.

     ///

     /// Sets the priority of the current item.

     void set(const Item& item, const Value& prio) {

       if (comp(prio, nodes[index[item]].prio)) {

         decrease(item, prio);

       } else if (!comp(prio, nodes[index[item]].prio)) {

         increase(item, prio);

       }

     }

     /// \brief Increase the priority of the current item.

     ///

     /// Increase the priority of the current item.

     void increase(const Item& item, const Value& prio) {

       int id = index[item];

       int kd = nodes[id].parent;

       nodes[id].prio = prio;

       while (kd >= 0 && nodes[kd].item == item) {

         setPrio(kd);

         kd = nodes[kd].parent;

       }

     }

     /// \brief Increase the priority of the current item.

     ///

     /// Increase the priority of the current item.

     void decrease(const Item& item, const Value& prio) {

       int id = index[item];

       int kd = nodes[id].parent;

       nodes[id].prio = prio;

       while (kd >= 0 && less(id, kd)) {

         nodes[kd].prio = prio;

         nodes[kd].item = item;

         kd = nodes[kd].parent;

       }

     }

     /// \brief Gives back the minimum priority of the class.

     ///

     /// \return Gives back the minimum priority of the class.

     const Value& classPrio(int cls) const {

       return nodes[~(classes[cls].parent)].prio;

     }

     /// \brief Gives back the minimum priority item of the class.

     ///

     /// \return Gives back the minimum priority item of the class.

     const Item& classTop(int cls) const {

       return nodes[~(classes[cls].parent)].item;

     }

     /// \brief Gives back a representant item of the class.

     /// 

     /// The representant is indpendent from the priorities of the

     /// items. 

     /// \return Gives back a representant item of the class.

     const Item& classRep(int id) const {

       int parent = classes[id].parent;

       return nodes[parent >= 0 ? classes[id].depth : leftNode(id)].item;

     }

     /// \brief Lemon style iterator for the items of a class.

     ///

     /// ClassIt is a lemon style iterator for the components. It iterates

     /// on the items of a class. By example if you want to iterate on

     /// each items of each classes then you may write the next code.

     ///\code

     /// for (ClassIt cit(huf); cit != INVALID; ++cit) {

     ///   std::cout << "Class: ";

     ///   for (ItemIt iit(huf, cit); iit != INVALID; ++iit) {

     ///     std::cout << toString(iit) << ' ' << std::endl;

     ///   }

     ///   std::cout << std::endl;

     /// }

     ///\endcode

     class ItemIt {

     private:

       const HeapUnionFind* _huf;

       int _id, _lid;

     public:

       /// \brief Default constructor 

       ///

       /// Default constructor 

       ItemIt() {}

       ItemIt(const HeapUnionFind& huf, int cls) : _huf(&huf) {

         int id = cls;

         int parent = _huf->classes[id].parent;

         if (parent >= 0) {

           _id = _huf->classes[id].depth;

           if (_huf->classes[id].next != -1) {

             _lid = _huf->classes[_huf->classes[id].next].depth;

           } else {

             _lid = -1;

           }

         } else {

           _id = _huf->leftNode(id);

           _lid = -1;

         } 

       }

       /// \brief Increment operator

       ///

       /// It steps to the next item in the class.

       ItemIt& operator++() {

         _id = _huf->nextNode(_id);

         return *this;

       }

       /// \brief Conversion operator

       ///

       /// It converts the iterator to the current item.

       operator const Item&() const {

         return _huf->nodes[_id].item;

       }

       /// \brief Equality operator

       ///

       /// Equality operator

       bool operator==(const ItemIt& i) { 

         return i._id == _id;

       }

       /// \brief Inequality operator

       ///

       /// Inequality operator

       bool operator!=(const ItemIt& i) { 

         return i._id != _id;

       }

       /// \brief Equality operator

       ///

       /// Equality operator

       bool operator==(Invalid) { 

         return _id == _lid;

       }

       /// \brief Inequality operator

       ///

       /// Inequality operator

       bool operator!=(Invalid) { 

         return _id != _lid;

       }

     };

     /// \brief Class iterator

     ///

     /// The iterator stores 

     class ClassIt {

     private:

       const HeapUnionFind* _huf;

       int _id;

     public:

       ClassIt(const HeapUnionFind& huf) 

         : _huf(&huf), _id(huf.first_class) {}

       ClassIt(const HeapUnionFind& huf, int cls) 

         : _huf(&huf), _id(huf.classes[cls].left) {}

       ClassIt(Invalid) : _huf(0), _id(-1) {}

       const ClassIt& operator++() {

         _id = _huf->classes[_id].next;

 	return *this;

       }

       /// \brief Equality operator

       ///

       /// Equality operator

       bool operator==(const ClassIt& i) { 

         return i._id == _id;

       }

       /// \brief Inequality operator

       ///

       /// Inequality operator

       bool operator!=(const ClassIt& i) { 

         return i._id != _id;

       }      

       operator int() const {

 	return _id;

       }

     };

   };

   //! @}

 } //namespace lemon

 #endif //LEMON_UNION_FIND_H