# Fast, stable, almost in-place radix and merge sorts

I've developed LSD radix sort algorithm that is

• stable,
• require only $O(\sqrt{RN})$ extra space when we sort into R buckets.

The same technique also allows to perform stable, efficient, almost in-place merge & quick sorts; as well as stable almost inplace merge/partition operations; and also applicable to some other list processing algorithms.

I describe it in my answer to this question. I believe that the proposed modifications to classic algorithms has important practical value. Other known techniques for (almost) in-place stable sort:

• WikiSort and GrailSort that have even smaller memory requirements but seems to perform more data moves than the simple mergesort
• highly theoretical in-place stable sorts that doesn't even pretend to have practical usefulness

EDIT: I've found this technique described 4 years ago in section 3.2.3 of paper A comprehensive study of main-memory partitioning and its application to large-scale comparison- and radix-sort

• Looks like you need to do some literature search. Use a search engine. – Yuval Filmus Jun 27 '18 at 14:16
• How can I perform this search? Googling for "stable in-place LSD radix sort" and its variants doesn't help. The best I can find is papers about in-place stable radix sort that is too slow for any practical usage. – Bulat Jun 28 '18 at 0:19
• In order to answer your question, I would need to do a literature search. It seems more appropriate for you to do it instead. You might be lucky and have an expert on this particular topic answer your question, but there are many topics and only few experts around. – Yuval Filmus Jun 28 '18 at 7:07
• I answered in the question body. The best I had found is TAOCP assignment w/o any description how it should be implemented. Never ever mentioned in other books. – Bulat Jun 28 '18 at 20:41
• TAoCP assignment 5.2.5 13, rated "term project"? – greybeard Jun 29 '18 at 4:07

Here is description of the proposed algorithm.

The algorithm is further improvement of the one proposed in Parallel radix sort with virtual memory and write-combining. They implemented fastest parallel LSD radix sort that avoids the counting pass by using "indefinite" space for output buckets. Each bucket allocs space on-demand using hardware paged memory support. A comprehensive study of main-memory partitioning and its application to large-scale comparison- and radix-sort took this idea further, proposing software implementation of the paged memory, that allows to perform stable radix/merge/quick sort using only $O(\sqrt{NR})$ extra memory.

Basic operations of radix, merge and quick sort have something in common:

• radix sort pass: read from 1 stream and write into R streams
• merge: read from R streams and write into 1 stream
• partition: read from 1 stream and write into R streams

Each of these operations process its input and output streams in strictly sequential manner. On every step, input plus output streams contain exactly N elements total - we neither produce nor consume any data items, but only move them between streams. This allows to implement in-place list-processing versions of these algorithms that don't move any data at all and only modify next list pointers. And this brings us to idea that we can perform these operations on array in-place, allocating only a little more space to perform home-made memory management of the array space.

To accomplish that, we split input and output streams into fixed-size pages and maintain list of pages in the each stream. When an input stream consumed all data on a page - this page is moved into the freepage list. When output stream needs more space - it allocs a page from the freepage list. Each stream may have incomplete first and last page, so in the worst case we may have 2*R*PAGESIZE unused items at any given moment - it's the extra space we need to alloc in addition to the page management structures. It's easy to find that the optimal PAGESIZE is $O(\sqrt{N/R})$ and the entire extra memory usage is $O(\sqrt{NR})$.

A few concrete examples:

• 2-way sorting of 1 million 8-byte items, i.e. 8 MB: we need 4 extra pages and can fit page number into 16 bits. 2 KB pages are optimal, and we need 16 KB of extra space total.
• 256-way sorting of the same array: we need 512 extra pages, while page number still fits into 16 bits. 256 byte pages are optimal, and we need 192 KB of extra space total. It's 2% of 8 MB, which is the largest proportion among these examples.
• 2-way sorting of 128 millions 8-byte items, i.e. 1 GB: we need 4 extra pages, and page number still fits into 16 bits. 32 KB pages are optimal, and we need 192 KB of extra space total. You may note that the two later cases has the same R*N and hence the same memory requirements.
• 256-way sorting of the same array: we need 512 extra pages, and page number finally requires full 32 bits. 4 KB pages are optimal, and we need 3 MB of extra space total.

TL;DR: sorting more than 1 MB of data usually requires less than 10% extra space, and the proportion reduces as datasize grows.

# Algorithm: utility functions

Let's define data structures for our memory manager:

// Number of items per page
const int PAGESIZE = 1024; // just for example

const int PAGE_COUNT = N/PAGESIZE;

// next_page[page] is either -1 (meaning EOF)
// or index of the next page in the same stream
int next_page[PAGE_COUNT];

// Input and output streams
StreamDescriptor in_streams[R];
StreamDescriptor out_streams[R];

struct StreamDescriptor
{
int first_page;
int first_element;  // index of first element on the first page
int last_page;
int last_element;  // index of last element on the last page + 1
// all pages between the first and last ones are fully populated
};


Here is the "memory manager" which allows to alloc new pages for output streams and release pages consumed from input streams:

// head of the available pages LIFO queue
int first_available_page;

int AllocPage()
{
int page = first_available_page;
assert(page != -1);
first_available_page = next_page[first_available_page];
return page;
}

void ReleasePage(int page)
{
next_page[page] = first_available_page;
first_available_page = page;
}

// Main data storage (holding N elements) provided by the caller
T* MainArray;

// Extra memory allocated for temporary data storage
T* ExtraPages;

// Return address of given memory page
T* Page (int page)
{
if (page < PAGE_COUNT)
return MainArray + PAGESIZE*page;
else
return ExtraPages + PAGESIZE*(page - PAGE_COUNT);
}


And these operations read item from given input stream and write item into given output stream:

T get_next_item (int input_stream)
{
StreamDescriptor s = in_streams[input_stream];
if (s.first_page == s.last_page  &&  s.first_element == s.last_element)
{
ReleasePage(s.first_page);
return EOF;
}
if (s.first_element == PAGESIZE)
{
int next = next_page[s.first_page];
ReleasePage(s.first_page);
if (next == -1)
return EOF;
s.first_page = next;
s.first_element = 0;
}
return Page(s.first_page)[s.first_element++];
}

void put_item (int output_stream, T item)
{
StreamDescriptor s = out_streams[output_stream];
if (s.last_element == PAGESIZE)
{
int next = AllocPage();
next_page[s.last_page] = next;
s.last_page = next;
s.last_element = 0;
}
Page(s.last_page)[s.last_element++] = item;
}


Initialization depends on the amount of input and output streams. For example, prior to first pass of LSD sort we do the following:

void LsdRadixSort (T* data_to_sort, int N)
{
MainArray = data_to_sort;
ExtraPages = new T[2*R*PAGESIZE];

// Setup first input stream as containing the entire MainArray contents
in_streams[0].first_page = 0;
in_streams[0].first_element = 0;
in_streams[0].last_page = PAGE_COUNT-1;
in_streams[0].last_element = N % PAGESIZE;

for (int i=0; i<PAGE_COUNT-1; i++)
next_page[i] = i+1;
next_page[PAGE_COUNT-1] = -1;

// Alloc one initial page for each outstream
for (int i=0; i<R; i++)
{
int page = PAGE_COUNT+i;
out_streams[i].first_page = page;
out_streams[i].first_element = 0;
out_streams[i].last_page = page;
out_streams[i].last_element = 0;
}

first_available_page = -1;

...
}


Prior to second and following passes of LSD sort, we perform this:

// Make old out_streams new in_streams
for (int i=0; i<R; i++)
{
in_streams[i] = out_streams[i];
}

cur_input_stream = 0;

// Alloc one initial page for each outstream
for (int i=0; i<R; i++)
{
int page = AllocPage();
out_streams[i].first_page = page;
out_streams[i].first_element = 0;
out_streams[i].last_page = page;
out_streams[i].last_element = 0;
}


In the first pass, data should be read with get_next_item(0) calls until it returns EOF. On the second and following passes, data should be read using the following procedure, also up to EOF:

// Enumerate contents of all input streams in order
T get_next_item_2nd()
{
if (cur_input_stream >= R)
return EOF;

T item = get_next_item(cur_input_stream);
if (item==EOF)
{
++cur_input_stream;
return get_next_item_2nd();
}

return item;
}


One pass of LSD radix sort:

void one_pass(int pass)
{
for(;;)
{
T item  =  pass==1? get_next_item(0) : get_next_item_2nd();
if (item==EOF)
break;
int bucket = calc_bucket(item,pass);
put_item(bucket,item);
}
}


# Algorithm: reordering data

Once sorting is finished, the caller may either read data in sorted order using the same get_next_item_2nd() procedure, or data may be placed in proper order in original array.

Reordering requires two steps:

1. compact data into single stream, removing "holes"
2. move each page of the compacted stream into right place

For the first step, we can just append to the first stream contents of remaining streams:

// Make old out_streams new in_streams (except for the first one)
for (int i=1; i<R; i++)
{
in_streams[i] = out_streams[i];
}

cur_input_stream = 1;

// Move data from all streams into the first one
for(;;)
{
T item = get_next_item_2nd();
if (item==EOF)
break;
put_item(0,item);
}


For the second step, we start with enumerating pages in the single stream remaining, and marking all other (i.e. free) pages with -1:

// final_place and next_page are the same array used in different roles
int* const final_place = next_page;

int page_i = out_streams[0].first_page;

for (int i=0; i<PAGE_COUNT; i++)
{
int next = next_page[page_i];
final_place[page_i] = i;
page_i = next;
}

while (first_available_page != -1)
{
int next = next_page[first_available_page];
final_place[first_available_page] = -1;
first_available_page = next;
}


From this point on, final_place[pg] is a page where current contents of page pg should be finally placed, or -1 for a page whose contents is garbage and so can be overwritten. Now all that remains is the classic "move them all to their final places" algorithm:

for (int i=0; i<PAGE_COUNT+2*R; i++)
{
int page = i;

while (final_place[page] != -1  &&  final_place[page] != page)
{
mem_exchange( Page(page), Page(final_place[page]), PAGESIZE*sizeof(T));
std::swap( final_place[page], final_place[final_place[page]]);
}
}

• What are the benefits of using a software implementation of paged memory instead of using hardware support for paging? – D.W. Jun 29 '18 at 2:52
• Each stream may have incomplete first and last page isn't that first for input streams or last for output streams? extra space we need to alloc in addition to the page management structures I think neglecting space needed for management is where you follow Wassenberg&Sanders in cheating about linear additional space. – greybeard Jun 29 '18 at 4:14
• @greybeard input stream may easily have incomplete last page just because it initially had uneven amount of items. With 2-way sorting and PAGESIZE=sqrt(N) we need just 4 extra pages for "overflowed" data and sqrt(N) entries in next_page array, so O(sqrt(N)) extra memory total. With R-way sorting and PAGESIZE=sqrt(N/R) we need 2*R extra pages of total size 2*sqrt(NR), plus sqrt(NR) entries in next_page array. – Bulat Jun 29 '18 at 4:33
• @D.W. I interested in inclusion of this technique into CUDA Thrust and Boost.Sort libraries. So 0) sorting even 20 MB array by 256-way radix sort will require 5GB memory space, i.e. impossible on 32-bit CPUs; 1) handling of OS-specific functionality (duplicated for each specific OS) in otherwise 100% Portable C library will be boring; 2) some modern environments (GPU, WebAssembly) probably don't have such functionality at all; 3) sometimes optimal page size will be less than 4 KB (see examples in the answer); 4) performing separate system call for every 4-32 KB map/unmap may be too expensive – Bulat Jun 29 '18 at 8:49