Heap

The heap is a fundamental data structure. The commonly used priority queue is heap-based, and heap sort is also built on the heap. We previously implemented a heap from first principles; now let’s analyze the STL heap implementation. While the STL version is naturally more elaborate, the underlying principle is the same.

Heap Implementation

The STL provides a suite of heap algorithms. Since a heap is stored as an array, how do we check whether a given array is a valid heap? The implementation is straightforward: verify that every parent is greater than or equal to its children (max-heap property). is_heap has other overloads not shown here.

  /**
   *  @brief  Determines whether a range is a heap using comparison functor.
   *  @param  __first  Start of range.
   *  @param  __last   End of range.
   *  @param  __comp   Comparison functor to use.
   *  @return  True if range is a heap, false otherwise.
   *  @ingroup heap_algorithms
  */
template<typename _RandomAccessIterator, typename _Compare>
inline bool is_heap(_RandomAccessIterator __first, _RandomAccessIterator __last, _Compare __comp) { // check if the array satisfies the heap property
    // concept requirements
    __glibcxx_function_requires(_RandomAccessIteratorConcept<_RandomAccessIterator>)
    __glibcxx_requires_valid_range(__first, __last);
    __glibcxx_requires_irreflexive_pred(__first, __last, __comp);

    const auto __dist = std::distance(__first, __last);
    typedef __decltype(__comp) _Cmp;
    __gnu_cxx::__ops::_Iter_comp_iter<_Cmp> __cmp(_GLIBCXX_MOVE(__comp));
    return std::__is_heap_until(__first, __dist, __cmp) == __dist;
}

template<typename _RandomAccessIterator, typename _Distance, typename _Compare>
_Distance __is_heap_until(_RandomAccessIterator __first, _Distance __n, _Compare& __comp) {     
    _Distance __parent = 0;
    for (_Distance __child = 1; __child < __n; ++__child) {     // straightforward: for a max-heap, parent >= child
	    if (__comp(__first + __parent, __first + __child))
	        return __child;
	    if ((__child & 1) == 0)
	        ++__parent;
	}
    return __n;
}

If an array does not satisfy the heap property, how do we transform it into a heap? make_heap does this. It starts from the last non-leaf node and sifts down, checking and swapping parent-child relationships all the way to the root.

  /**
   *  @brief  Construct a heap over a range.
   *  @param  __first  Start of heap.
   *  @param  __last   End of heap.
   *  @ingroup heap_algorithms
   *
   *  This operation makes the elements in [__first,__last) into a heap.
  */
template<typename _RandomAccessIterator>
inline void make_heap(_RandomAccessIterator __first, _RandomAccessIterator __last) {
    // concept requirements
    __glibcxx_function_requires(_Mutable_RandomAccessIteratorConcept<_RandomAccessIterator>)
    __glibcxx_function_requires(_LessThanComparableConcept<typename iterator_traits<_RandomAccessIterator>::value_type>)
    __glibcxx_requires_valid_range(__first, __last);
    __glibcxx_requires_irreflexive(__first, __last);

    __gnu_cxx::__ops::_Iter_less_iter __comp;
    std::__make_heap(__first, __last, __comp);
}

template<typename _RandomAccessIterator, typename _Compare>
void __make_heap(_RandomAccessIterator __first, _RandomAccessIterator __last, _Compare& __comp) {
    typedef typename iterator_traits<_RandomAccessIterator>::value_type _ValueType;
    typedef typename iterator_traits<_RandomAccessIterator>::difference_type _DistanceType;

    if (__last - __first < 2)   
	    return;

    const _DistanceType __len = __last - __first;
    _DistanceType __parent = (__len - 2) / 2;
    while (true) {
	    _ValueType __value = _GLIBCXX_MOVE(*(__first + __parent));
	    std::__adjust_heap(__first, __parent, __len, _GLIBCXX_MOVE(__value), __comp);
	    if (__parent == 0)
	        return;
	    __parent--;
	}
}

template<typename _RandomAccessIterator, typename _Distance, typename _Tp, typename _Compare>
void __adjust_heap(_RandomAccessIterator __first, _Distance __holeIndex, _Distance __len, _Tp __value, _Compare __comp) {
    const _Distance __topIndex = __holeIndex;
    _Distance __secondChild = __holeIndex;
    while (__secondChild < (__len - 1) / 2) {
	    __secondChild = 2 * (__secondChild + 1);
	    if (__comp(__first + __secondChild, __first + (__secondChild - 1)))
	        __secondChild--;
	    *(__first + __holeIndex) = _GLIBCXX_MOVE(*(__first + __secondChild));
	    __holeIndex = __secondChild;
	}
    
    if ((__len & 1) == 0 && __secondChild == (__len - 2) / 2) {
	    __secondChild = 2 * (__secondChild + 1);
	    *(__first + __holeIndex) = _GLIBCXX_MOVE(*(__first + (__secondChild - 1)));
	    __holeIndex = __secondChild - 1;
	}
    __decltype(__gnu_cxx::__ops::__iter_comp_val(_GLIBCXX_MOVE(__comp)))
	__cmp(_GLIBCXX_MOVE(__comp));
    std::__push_heap(__first, __holeIndex, __topIndex, _GLIBCXX_MOVE(__value), __cmp);
}

template<typename _RandomAccessIterator, typename _Distance, typename _Tp, typename _Compare>
void __push_heap(_RandomAccessIterator __first, _Distance __holeIndex, _Distance __topIndex, _Tp __value, _Compare& __comp) {
    _Distance __parent = (__holeIndex - 1) / 2;
    while (__holeIndex > __topIndex && __comp(__first + __parent, __value)) {
	    *(__first + __holeIndex) = _GLIBCXX_MOVE(*(__first + __parent));
	    __holeIndex = __parent;
	    __parent = (__holeIndex - 1) / 2;
	}
    *(__first + __holeIndex) = _GLIBCXX_MOVE(__value);
}

After adding an element to the heap, push_heap restores the heap property by sifting the new element upward until the heap invariant is satisfied.

/**
   *  @brief  Push an element onto a heap.
   *  @param  __first  Start of heap.
   *  @param  __last   End of heap + element.
   *  @ingroup heap_algorithms
   *
   *  This operation pushes the element at last-1 onto the valid heap
   *  over the range [__first,__last-1).  After completion,
   *  [__first,__last) is a valid heap.
*/
template<typename _RandomAccessIterator>
inline void push_heap(_RandomAccessIterator __first, _RandomAccessIterator __last) {
    typedef typename iterator_traits<_RandomAccessIterator>::value_type _ValueType;
    typedef typename iterator_traits<_RandomAccessIterator>::difference_type _DistanceType;

    // concept requirements
    __glibcxx_function_requires(_Mutable_RandomAccessIteratorConcept<_RandomAccessIterator>)
    __glibcxx_function_requires(_LessThanComparableConcept<_ValueType>)
    __glibcxx_requires_valid_range(__first, __last);
    __glibcxx_requires_irreflexive(__first, __last);
    __glibcxx_requires_heap(__first, __last - 1);

    __gnu_cxx::__ops::_Iter_less_val __comp;
    _ValueType __value = _GLIBCXX_MOVE(*(__last - 1));
    std::__push_heap(__first, _DistanceType((__last - __first) - 1), _DistanceType(0), _GLIBCXX_MOVE(__value), __comp);
}

Popping the top element moves it to the end of the container, then sifts the replacement root downward to restore the heap property.

  /**
   *  @brief  Pop an element off a heap using comparison functor.
   *  @param  __first  Start of heap.
   *  @param  __last   End of heap.
   *  @param  __comp   Comparison functor to use.
   *  @ingroup heap_algorithms
   *
   *  This operation pops the top of the heap.  The elements __first
   *  and __last-1 are swapped and [__first,__last-1) is made into a
   *  heap.  Comparisons are made using comp.
  */
template<typename _RandomAccessIterator, typename _Compare>
inline void pop_heap(_RandomAccessIterator __first, _RandomAccessIterator __last, _Compare __comp) {
    // concept requirements
    __glibcxx_function_requires(_Mutable_RandomAccessIteratorConcept<_RandomAccessIterator>)
    __glibcxx_requires_valid_range(__first, __last);
    __glibcxx_requires_irreflexive_pred(__first, __last, __comp);
    __glibcxx_requires_non_empty_range(__first, __last);
    __glibcxx_requires_heap_pred(__first, __last, __comp);

    if (__last - __first > 1) {
	    typedef __decltype(__comp) _Cmp;
	    __gnu_cxx::__ops::_Iter_comp_iter<_Cmp> __cmp(_GLIBCXX_MOVE(__comp));
	    --__last;
	    std::__pop_heap(__first, __last, __last, __cmp);
	}
}

template<typename _RandomAccessIterator, typename _Compare>
inline void __pop_heap(_RandomAccessIterator __first, _RandomAccessIterator __last, _RandomAccessIterator __result, _Compare& __comp) {
    typedef typename iterator_traits<_RandomAccessIterator>::value_type _ValueType;
    typedef typename iterator_traits<_RandomAccessIterator>::difference_type _DistanceType;

    _ValueType __value = _GLIBCXX_MOVE(*__result);
    *__result = _GLIBCXX_MOVE(*__first);
    std::__adjust_heap(__first, _DistanceType(0), _DistanceType(__last - __first), _GLIBCXX_MOVE(__value), __comp);
}

The STL also provides a heap sort implementation:

 /**
   *  @brief  Sort a heap using comparison functor.
   *  @param  __first  Start of heap.
   *  @param  __last   End of heap.
   *  @param  __comp   Comparison functor to use.
   *  @ingroup heap_algorithms
   *
   *  This operation sorts the valid heap in the range [__first,__last).
   *  Comparisons are made using __comp.
  */
template<typename _RandomAccessIterator, typename _Compare>
inline void sort_heap(_RandomAccessIterator __first, _RandomAccessIterator __last, _Compare __comp) {
    // concept requirements
    __glibcxx_function_requires(_Mutable_RandomAccessIteratorConcept<_RandomAccessIterator>)
    __glibcxx_requires_valid_range(__first, __last);
    __glibcxx_requires_irreflexive_pred(__first, __last, __comp);
    __glibcxx_requires_heap_pred(__first, __last, __comp);

    typedef __decltype(__comp) _Cmp;
    __gnu_cxx::__ops::_Iter_comp_iter<_Cmp> __cmp(_GLIBCXX_MOVE(__comp));
    std::__sort_heap(__first, __last, __cmp);
}

template<typename _RandomAccessIterator, typename _Compare>
void __sort_heap(_RandomAccessIterator __first, _RandomAccessIterator __last, _Compare& __comp) {
    while (__last - __first > 1) {
	    --__last;
	    std::__pop_heap(__first, __last, __last, __comp);
	}
}