Ford-Johnson merge-insertion sort

Question

The Ford-Johnson algorithm, also known as merge-insertion sort (the name was probably given by Knuth) is an in-place sorting algorithm designed to perform as few comparisons as possible to sort a collection. Unfortunately, the algorithm needs some specific data structures to keep track of the order of the elements and is generally too slow to be practical. Anyway, it is an interesting algorithm from a computer science point of view; while not performing an optimal number of comparisons, it's still a reference and one of the best known comparison sorts when it comes to reducing the number of comparisons.

The algorithm is not the simplest in the world and I couldn't find a proper explanation online, so I will try to explain it as I can based on the descriptions in The Art of Computer Programming, Volume 3 by Donald Knuth.

Making the best of binary search

To perform a minimal number of comparisons, we need to take into account the following observation about binary search: the maximal number of comparisons needed to perform a binary search on a sorted sequence is the same when the number of elements is \$2^n\$ and when it is \$2^{n+1}-1\$. For example, looking for an element in a sorted sequence of \$8\$ or \$15\$ elements requires the same number of comparisons.

Many insertion-based sorting algorithms perform binary searches to find where to insert elements, but most of them don't take that property of binary search into account.

The merge-insertion sort

The Ford-Johnson merge-insertion sort is a three-step algorithm, let \$n\$ be the number of elements to sort:

1. Split the collection in \$n/2\$ pairs of \$2\$ elements and order these elements pairwise. If the number of elements is odd, the last element of the collection isn't paired with any element.

2. Recursively sort the pairs of elements by their highest value. If the number of elements is odd, the last element is not considered a highest value and is left at the end of the collection. Consider that the highest values form a sorted list that we will call the main chain while the rest of the elements will be known as pend elements. Tag the elements of the main chain with the names \$a_1, a_2, a_3, ..., a_{n/2}\$ then tag the pend elements with the names \$b_1, b_2, b_3, ..., b_{n/2}\$. For every \$k\$, we have the relation \$b_k \le a_k\$.

3. Insert the pend elements into the main chain. We know that the first pend element \$b_1\$ is lesser than \$a_1\$; we consider it to be part of the main chain which then becomes \$\{b_1, a_1, a_2, a_3, ..., a_{n/2}\}\$. Now, we need to insert the other pend elements into the main chain in a way that ensures that the size of the insertion area is a power of \$2\$ minus \$1\$ as often as possible. Basically, we will insert \$b_3\$ in \$\{b_1, a_1, a_2\}\$ (we know that it is less than \$a_3\$), then we will insert \$b_2\$ into \$\{b_1, a_1, b_3\}\$ no matter where \$b_3\$ was inserted. Note that during these insertions, the size of the insertion area is always at most 3.

   The value of the next pend element \$b_k\$ to insert into the main chain while minimizing the number of comparisons actually corresponds to the next Jacobsthal number. We inserted the element \$3\$ first, so the next will be \$5\$ then \$11\$ then \$21\$, etc...

   To sum up, the insertion order of the first pend elements into the main chain is as follows: \$b_1, b_3, b_2, b_5, b_4, b_{11}, b_{10}, b_9, b_8, b_7, b_6, ...\$.

To be honest, this explanation is probably not clear enough, and I really recommend that you look for other descriptions of the algorithm from whatever resource is available to you; they generally come with pictures to make it more obvious what is happening. I can at least give you links to a step-by-step description of the unrolled algorithm for 5 elements; The Art of Computer Programming, Volume 3 contains descriptions of the algorithm for 5 and 21 elements in section 5.3.1 as well as a general description of the algorithm (if you have access to a copy of the book).

Gimme teh codez

I did what I could to describe the algorithm. Now, let's give you the code. Most optimizations are described in the comments near the relevant parts. Also, some parts are ugly, such as the full double move of the collection to and from a cache in the end, but it is hard to get rid of this specific part. To perform the recursive sort smoothly, the algorithm uses a group_iterator class that "views" a contiguous group of elements from the collection of a given size but only uses the last element of the group when performing a comparison, even though it is able to swap the whole viewed area with another part of the collection.

#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <functional>
#include <list>
#include <iterator>
#include <type_traits>
#include <vector>

////////////////////////////////////////////////////////////
// Iterator used to sort views of the collection

template<typename Iterator>
class group_iterator
{
    private:

        Iterator _it;
        std::size_t _size;

    public:

        ////////////////////////////////////////////////////////////
        // Public types

        using iterator_category = std::random_access_iterator_tag;
        using iterator_type     = Iterator;
        using value_type        = typename std::iterator_traits<Iterator>::value_type;
        using difference_type   = typename std::iterator_traits<Iterator>::difference_type;
        using pointer           = typename std::iterator_traits<Iterator>::pointer;
        using reference         = typename std::iterator_traits<Iterator>::reference;

        ////////////////////////////////////////////////////////////
        // Constructors

        group_iterator() = default;

        group_iterator(Iterator it, std::size_t size):
            _it(it),
            _size(size)
        {}

        ////////////////////////////////////////////////////////////
        // Members access

        auto base() const
            -> iterator_type
        {
            return _it;
        }

        auto size() const
            -> std::size_t
        {
            return _size;
        }

        ////////////////////////////////////////////////////////////
        // Element access

        auto operator*() const
            -> reference
        {
            return _it[_size - 1];
        }

        auto operator->() const
            -> pointer
        {
            return &(operator*());
        }

        ////////////////////////////////////////////////////////////
        // Increment/decrement operators

        auto operator++()
            -> group_iterator&
        {
            _it += _size;
            return *this;
        }

        auto operator++(int)
            -> group_iterator
        {
            auto tmp = *this;
            operator++();
            return tmp;
        }

        auto operator--()
            -> group_iterator&
        {
            _it -= _size;
            return *this;
        }

        auto operator--(int)
            -> group_iterator
        {
            auto tmp = *this;
            operator--();
            return tmp;
        }

        auto operator+=(std::size_t increment)
            -> group_iterator&
        {
            _it += _size * increment;
            return *this;
        }

        auto operator-=(std::size_t increment)
            -> group_iterator&
        {
            _it -= _size * increment;
            return *this;
        }

        ////////////////////////////////////////////////////////////
        // Elements access operators

        auto operator[](std::size_t pos)
            -> decltype(_it[pos * _size + _size - 1])
        {
            return _it[pos * _size + _size - 1];
        }

        auto operator[](std::size_t pos) const
            -> decltype(_it[pos * _size + _size - 1])
        {
            return _it[pos * _size + _size - 1];
        }
};

template<typename Iterator1, typename Iterator2>
auto iter_swap(group_iterator<Iterator1> lhs, group_iterator<Iterator2> rhs)
    -> void
{
    std::swap_ranges(lhs.base(), lhs.base() + lhs.size(), rhs.base());
}

////////////////////////////////////////////////////////////
// Comparison operators

template<typename Iterator1, typename Iterator2>
auto operator==(const group_iterator<Iterator1>& lhs,
                const group_iterator<Iterator2>& rhs)
    -> bool
{
    return lhs.base() == rhs.base();
}

template<typename Iterator1, typename Iterator2>
auto operator!=(const group_iterator<Iterator1>& lhs,
                const group_iterator<Iterator2>& rhs)
    -> bool
{
    return lhs.base() != rhs.base();
}

////////////////////////////////////////////////////////////
// Relational operators

template<typename Iterator1, typename Iterator2>
auto operator<(const group_iterator<Iterator1>& lhs,
               const group_iterator<Iterator2>& rhs)
    -> bool
{
    return lhs.base() < rhs.base();
}

template<typename Iterator1, typename Iterator2>
auto operator<=(const group_iterator<Iterator1>& lhs,
                const group_iterator<Iterator2>& rhs)
    -> bool
{
    return lhs.base() <= rhs.base();
}

template<typename Iterator1, typename Iterator2>
auto operator>(const group_iterator<Iterator1>& lhs,
               const group_iterator<Iterator2>& rhs)
    -> bool
{
    return lhs.base() > rhs.base();
}

template<typename Iterator1, typename Iterator2>
auto operator>=(const group_iterator<Iterator1>& lhs,
                const group_iterator<Iterator2>& rhs)
    -> bool
{
    return lhs.base >= rhs.base();
}

////////////////////////////////////////////////////////////
// Arithmetic operators

template<typename Iterator>
auto operator+(group_iterator<Iterator> it, std::size_t size)
    -> group_iterator<Iterator>
{
    return it += size;
}

template<typename Iterator>
auto operator+(std::size_t size, group_iterator<Iterator> it)
    -> group_iterator<Iterator>
{
    return it += size;
}

template<typename Iterator>
auto operator-(group_iterator<Iterator> it, std::size_t size)
    -> group_iterator<Iterator>
{
    return it -= size;
}

template<typename Iterator>
auto operator-(const group_iterator<Iterator>& lhs, const group_iterator<Iterator>& rhs)
    -> typename group_iterator<Iterator>::difference_type
{
    return (lhs.base() - rhs.base()) / lhs.size();
}

////////////////////////////////////////////////////////////
// Construction function

template<typename Iterator>
auto make_group_iterator(Iterator it, std::size_t size)
    -> group_iterator<Iterator>
{
    return { it, size };
}

template<typename Iterator>
auto make_group_iterator(group_iterator<Iterator> it, std::size_t size)
    -> group_iterator<Iterator>
{
    return { it.base(), size * it.size() };
}

////////////////////////////////////////////////////////////
// Merge-insertion sort

template<
    typename RandomAccessIterator,
    typename Compare
>
auto merge_insertion_sort_impl(RandomAccessIterator first, RandomAccessIterator last,
                               Compare compare)
{
    // Cache all the differences between a Jacobsthal number and its
    // predecessor that fit in 64 bits, starting with the difference
    // between the Jacobsthal numbers 4 and 3 (the previous ones are
    // unneeded)
    static constexpr std::uint_least64_t jacobsthal_diff[] = {
        2u, 2u, 6u, 10u, 22u, 42u, 86u, 170u, 342u, 682u, 1366u,
        2730u, 5462u, 10922u, 21846u, 43690u, 87382u, 174762u, 349526u, 699050u,
        1398102u, 2796202u, 5592406u, 11184810u, 22369622u, 44739242u, 89478486u,
        178956970u, 357913942u, 715827882u, 1431655766u, 2863311530u, 5726623062u,
        11453246122u, 22906492246u, 45812984490u, 91625968982u, 183251937962u,
        366503875926u, 733007751850u, 1466015503702u, 2932031007402u, 5864062014806u,
        11728124029610u, 23456248059222u, 46912496118442u, 93824992236886u, 187649984473770u,
        375299968947542u, 750599937895082u, 1501199875790165u, 3002399751580331u,
        6004799503160661u, 12009599006321322u, 24019198012642644u, 48038396025285288u,
        96076792050570576u, 192153584101141152u, 384307168202282304u, 768614336404564608u,
        1537228672809129216u, 3074457345618258432u, 6148914691236516864u
    };

    using std::iter_swap;

    auto size = std::distance(first, last);
    if (size < 2) return;

    // Whether there is a stray element not in a pair
    // at the end of the chain
    bool has_stray = (size % 2 != 0);

    ////////////////////////////////////////////////////////////
    // Group elements by pairs

    auto end = has_stray ? std::prev(last) : last;
    for (auto it = first ; it != end ; it += 2)
    {
        if (compare(it[1], it[0]))
        {
            iter_swap(it, it + 1);
        }
    }

    ////////////////////////////////////////////////////////////
    // Recursively sort the pairs by max

    merge_insertion_sort(
        make_group_iterator(first, 2),
        make_group_iterator(end, 2),
        compare
    );

    ////////////////////////////////////////////////////////////
    // Separate main chain and pend elements

    // Small node struct for pend elements
    struct node
    {
        RandomAccessIterator it;
        typename std::list<RandomAccessIterator>::iterator next;
    };

    // The first pend element is always part of the main chain,
    // so we can safely initialize the list with the first two
    // elements of the sequence
    std::list<RandomAccessIterator> chain = { first, std::next(first) };
    std::list<node> pend;

    for (auto it = first + 2 ; it != end ; it += 2)
    {
        auto tmp = chain.insert(chain.end(), std::next(it));
        pend.push_back({it, tmp});
    }

    // Add the last element to pend if it exists, when it
    // exists, it always has to be inserted in the full chain,
    // so giving it chain.end() as end insertion point is ok
    if (has_stray)
    {
        pend.push_back({end, chain.end()});
    }

    ////////////////////////////////////////////////////////////
    // Binary insertion into the main chain

    for (int k = 0 ; ; ++k)
    {
        // Find next index
        auto dist = jacobsthal_diff[k];
        if (dist >= pend.size()) break;
        auto it = pend.begin();
        std::advance(it, dist);

        while (true)
        {
            auto insertion_point = std::upper_bound(
                chain.begin(), it->next, it->it,
                [=](auto lhs, auto rhs) {
                    return compare(*lhs, *rhs);
                }
            );
            chain.insert(insertion_point, it->it);

            it = pend.erase(it);
            if (it == pend.begin()) break;
            --it;
        }
    }

    // If there are elements left, insert them too
    while (not pend.empty())
    {
        auto it = std::prev(pend.end());
        auto insertion_point = std::upper_bound(
            chain.begin(), it->next, it->it,
            [=](auto lhs, auto rhs) {
                return compare(*lhs, *rhs);
            }
        );
        chain.insert(insertion_point, it->it);
        pend.pop_back();
    }

    ////////////////////////////////////////////////////////////
    // Move values in order to a cache then back to origin

    std::vector<typename std::iterator_traits<RandomAccessIterator>::value_type> cache;
    cache.reserve(size);

    for (auto&& it: chain)
    {
        auto begin = it.base();
        auto end = begin + it.size();
        std::move(begin, end, std::back_inserter(cache));
    }
    std::move(cache.begin(), cache.end(), first.base());
}

template<
    typename RandomAccessIterator,
    typename Compare = std::less<>
>
auto merge_insertion_sort(RandomAccessIterator first, RandomAccessIterator last,
                          Compare compare={})
    -> void
{
    merge_insertion_sort_impl(
        make_group_iterator(first, 1),
        make_group_iterator(last, 1),
        compare
    );
}

The algorithm is slow but was never expected to be fast, only to perform fewer comparisons than most sorting algorithms. Now, is it possible to improve it, be it from an elegance or performance (and of course correctness) point of view?

---

As a side note, I still maintain an implementation of the algorithm in a library of mine. Even if most errors have been highlighted in my self-answer, I might have forgotten some of them. Don't hesitate to answer if you think it can still be improved, or to ping me if some link is dead.

Answer 1

First of all, there is an error

Not a fatal error though since the algorithm still sorts the collection properly, but it actually means that the algorithm doesn't perform as few comparisons as it should. The line std::advance(it, dist); advances the iterator one step too far, so the binary insertion is sometimes done in a main chain too big compared to what it should be (more than \$2^n-1\$ elements). The obvious solution is to advance the iterator by dist - 1 instead of dist; however, removing 1 from every element in jacobsthal_diff is also a solution.

We don't need to remove elements from pend

Instead of erasing elements from pend once they have been inserted into chain, we can instead track the first used iterator in pend corresponding to a Jacobsthal diff, add the next Jacobsthal diff to find the next such iterator, and decrease that iterator until it encounters the previous remembered Jacobsthal diff iterator. Not having to remove the elements from pend means that we don't need to store nodes in a container supporting fast deletion from the beginning. Basically, since we don't remove anything, we can switch to an std::vector<node> for pend.

Use the original iterators

Since we only have to add a Jacobsthal diff number to find the next element of pend, it means that we can perform the same operation on the original collection to find the iterator to insert into the main chain (every element with an even index in [first, last) is a pend iterator). It means that we can drop this information from pend and only store an std::vector<typename std::list<RandomAccessIterator>::value_type> instead of an std::vector<node>. The maximal size of the vector should be (size + 1) / 2 - 1 so we can directly reserve that amount of elements.

We can insert the remaining elements in any order

At first I thought that the remaining pend elements had to be inserted in reverse order (the elements left when the farthest pend element whose index corresponds to a Jacobsthal number has been inserted). However, it appears that we can insert them in any order thanks to the properties of binary search. Therefore, inserting them in ascending order shoud probably ease the CPU's task.

Smaller things

• iter_swap isn't right: the overload of iter_swap for group_iterator is designed to swap several group_iterator of different types, which isn't quite right. Not only does it look like it can cause problems, but apparently it can also cause ADL problems: in a more complex case, the compiler found the unqualified call to iter_swap ambiguous. The solution was to make iter_swap work only with group_iterators of the same type:

  template<typename Iterator>
  auto iter_swap(group_iterator<Iterator> lhs, group_iterator<Iterator> rhs)
      -> void
  {
      std::swap_ranges(lhs.base(), lhs.base() + lhs.size(), rhs.base());
  }
  

• The recursion is needlessly complicated: merge_insertion_sort_impl calls merge_insertion_sort which calls... merge_insertion_sort_impl without adding anything significant while it introduces yet another indirection. While it is likely to get elided by the compiler, making a direct recursive call of merge_insertion_sort_impl makes things easier for everyone.

• There is still a small error in operator>= for group_iterator: the parenthesis after lhs.base are missing, which is likely to cause a compilation error if the function is ever called.

Putting it all together

Once we stick all these remarks together, merge_insertion_sort_impl looks like this:

template<
    typename RandomAccessIterator,
    typename Compare
>
auto merge_insertion_sort_impl(RandomAccessIterator first, RandomAccessIterator last,
                               Compare compare)
{
    // Cache all the differences between a Jacobsthal number and its
    // predecessor that fit in 64 bits, starting with the difference
    // between the Jacobsthal numbers 4 and 3 (the previous ones are
    // unneeded)
    static constexpr std::uint_fast64_t jacobsthal_diff[] = {
        2u, 2u, 6u, 10u, 22u, 42u, 86u, 170u, 342u, 682u, 1366u,
        2730u, 5462u, 10922u, 21846u, 43690u, 87382u, 174762u, 349526u, 699050u,
        1398102u, 2796202u, 5592406u, 11184810u, 22369622u, 44739242u, 89478486u,
        178956970u, 357913942u, 715827882u, 1431655766u, 2863311530u, 5726623062u,
        11453246122u, 22906492246u, 45812984490u, 91625968982u, 183251937962u,
        366503875926u, 733007751850u, 1466015503702u, 2932031007402u, 5864062014806u,
        11728124029610u, 23456248059222u, 46912496118442u, 93824992236886u, 187649984473770u,
        375299968947542u, 750599937895082u, 1501199875790165u, 3002399751580331u,
        6004799503160661u, 12009599006321322u, 24019198012642644u, 48038396025285288u,
        96076792050570576u, 192153584101141152u, 384307168202282304u, 768614336404564608u,
        1537228672809129216u, 3074457345618258432u, 6148914691236516864u
    };

    using std::iter_swap;

    auto size = std::distance(first, last);
    if (size < 2) return;

    // Whether there is a stray element not in a pair
    // at the end of the chain
    bool has_stray = (size % 2 != 0);

    ////////////////////////////////////////////////////////////
    // Group elements by pairs

    auto end = has_stray ? std::prev(last) : last;
    for (auto it = first ; it != end ; it += 2)
    {
        if (compare(it[1], it[0]))
        {
            iter_swap(it, it + 1);
        }
    }

    ////////////////////////////////////////////////////////////
    // Recursively sort the pairs by max

    merge_insertion_sort_impl(
        make_group_iterator(first, 2),
        make_group_iterator(end, 2),
        compare
    );

    ////////////////////////////////////////////////////////////
    // Separate main chain and pend elements

    // The first pend element is always part of the main chain,
    // so we can safely initialize the list with the first two
    // elements of the sequence
    std::list<RandomAccessIterator> chain = { first, std::next(first) };

    // Upper bounds for the insertion of pend elements
    std::vector<typename std::list<RandomAccessIterator>::iterator> pend;
    pend.reserve((size + 1) / 2 - 1);

    for (auto it = first + 2 ; it != end ; it += 2)
    {
        auto tmp = chain.insert(std::end(chain), std::next(it));
        pend.push_back(tmp);
    }

    // Add the last element to pend if it exists; when it
    // exists, it always has to be inserted in the full chain,
    // so giving it chain.end() as end insertion point is ok
    if (has_stray)
    {
        pend.push_back(std::end(chain));
    }

    ////////////////////////////////////////////////////////////
    // Binary insertion into the main chain

    auto current_it = first + 2;
    auto current_pend = std::begin(pend);

    for (int k = 0 ; ; ++k)
    {
        // Should be safe: in this code, std::distance should always return
        // a positive number, so there is no risk of comparing funny values
        using size_type = std::common_type_t<
            std::uint_fast64_t,
            typename std::list<RandomAccessIterator>::difference_type
        >;

        // Find next index
        auto dist = jacobsthal_diff[k];
        if (dist > static_cast<size_type>(std::distance(current_pend, std::end(pend)))) break;

        auto it = std::next(current_it, dist * 2);
        auto pe = std::next(current_pend, dist);

        do
        {
            --pe;
            it -= 2;

            auto insertion_point = std::upper_bound(
                std::begin(chain), *pe, *it,
                [=](const auto& lhs, const auto& rhs) {
                    return compare(lhs, *rhs);
                }
            );
            chain.insert(insertion_point, it);
        } while (pe != current_pend);

        std::advance(current_it, dist * 2);
        std::advance(current_pend, dist);
    }

    // If there are pend elements left, insert them into
    // the main chain, the order of insertion does not
    // matter so forward traversal is ok
    while (current_pend != std::end(pend))
    {
        auto insertion_point = std::upper_bound(
            std::begin(chain), *current_pend, *current_it,
            [=](const auto& lhs, const auto& rhs) {
                return compare(lhs, *rhs);
            }
        );
        chain.insert(insertion_point, current_it);
        current_it += 2;
        ++current_pend;
    }

    ////////////////////////////////////////////////////////////
    // Move values in order to a cache then back to origin

    std::vector<typename std::iterator_traits<RandomAccessIterator>::value_type> cache;
    cache.reserve(size);

    for (auto&& it: chain)
    {
        auto begin = it.base();
        auto end = begin + it.size();
        std::move(begin, end, std::back_inserter(cache));
    }
    std::move(std::begin(cache), std::end(cache), first.base());
}

The algorithm remains several orders of magnitude slower than most common sorting algorithms, but it's more correct and a bit faster than the original version in the question.