Fastest way to combine pre-sorted lists?

Question

I need to combine two long, already-sorted lists into one sorted list. My effort below takes advantage of the fact that the lists are already sorted; but it's far slower than just using Join and Sort, which ignore the fact that the lists were sorted.

a = Sort[RandomInteger[10^8, 10^6]];
b = Sort[RandomInteger[10^8, 10^7]];
Timing[c1 = Sort[Flatten[{a, b}]];]
Timing[c2 = Sort[Join[a, b]];]

Timing[ 
 aIndex = 1;
 bIndex = 1;
 c3 = ConstantArray[0, Length[a] + Length[b]];
 cIndex = 1; 
 While[aIndex <= Length[a] && bIndex <= Length[b],c3[[cIndex++]] = If[a[[aIndex]] <= b[[bIndex]], a[[aIndex++]], b[[bIndex++]]]];

 c3[[cIndex ;; Length[c3]]] = If[aIndex > Length[a], b[[bIndex ;; Length[b]]], a[[aIndex ;; Length[a]]]];]

(* {1.14063, Null} {0.24063, Null} {70.4063, Null} *)

Can we do better?

Answer 1

I rewrote your procedural loop for a compiled function:

mergeSort = Compile[{{a, _Integer, 1}, {b, _Integer, 1}},
  Block[{aIndex = 1, la = Length[a], lb = Length[b], 
    lc = Length[a] + Length[b],
    bIndex = 1, cIndex = 1,
    c3 = Table[0, {Length[a] + Length[b]}]},
   For[cIndex = 1, aIndex <= la && bIndex <= lb, cIndex++, 
    c3[[cIndex]] = 
     If[a[[aIndex]] <= b[[bIndex]], a[[aIndex++]], b[[bIndex++]]]];
   c3[[cIndex ;; lc]] = 
    If[aIndex > la, b[[bIndex ;; lb]], a[[aIndex ;; la]]];
   c3]
  , CompilationTarget -> "C", Parallelization -> True, 
  "RuntimeOptions" -> "Speed"]

I also went with corey's suggestion:

mergeSort2 = Compile[{{a, _Integer, 1}, {b, _Integer, 1}},
  Sort[Join[a, b]]
  , CompilationTarget -> "C", Parallelization -> True, 
  "RuntimeOptions" -> "Speed"]

The timings are very close.

RepeatedTiming[mergeSort[a, b];]
RepeatedTiming[Sort[Join[a, b]];]
RepeatedTiming[mergeSort2[a, b];]

{0.299, Null}
{0.33, Null}
{0.330, Null}

The list order doesn't make a difference:

RepeatedTiming[mergeSort[b, a];]
RepeatedTiming[Sort[Join[b, a]];]
RepeatedTiming[mergeSort2[b, a];]

{0.302, Null}
{0.271, Null}
{0.308, Null}

As you have seen in your tests, Sort@Flatten[{a,b}] is significantly slower. Turns out, this is unsurprising:

mergeSort3 = Compile[{{a, _Integer, 1}, {b, _Integer, 1}},
  Sort[Flatten[{a, b}]]
  , CompilationTarget -> "C", Parallelization -> True, 
  "RuntimeOptions" -> "Speed"]

Calling mergeSort3[a,b] returns errors about non-tensor objects being generated: indeed, {a, b} is a list of two lists of different lengths.

A concatenation of two lists is likely very favorable for a real merge-sort (nothing to do with the names of my functions). Sort is (almost certainly) implemented in low-level code which would be very hard to beat, though I agree, that dropping the assumption that the two lists are pre-sorted introduces some overhead. However, this overhead seems to be on the order of the error-bar of timings.

Update

On request of kglr:

RepeatedTiming[#[[Ordering@#]] &@Join[a, b];]
RepeatedTiming[#[[Ordering@#]] &@Join[b, a];]

{0.33, Null}
{0.276, Null}

I've run this a few times and Join[b,a] seems consistently slightly faster than Join[a,b], however otherwise it's more or less the same as the other functions.

Answer 2

This question is related to: Complement on pre-sorted lists

Giving basically the same answer I did there I don't think Sort is being "wasteful" in this application and I doubt you will be able to substantially improve upon it. By definition sorting must still take place outside of a trivial case like all elements of one list being larger than any of those in the other.

As Simon Woods already observed the Sort algorithm is more efficient on lists that are partially ordered. Although the algorithm is not the same take as illustration How can I collect data for visualization of quick sort? Using SeedRandom[0]; s = RandomSample @ Range @ 100; the sort takes 55 steps. However if we pre-sort each half of the list with s = Join @@ Sort /@ Partition[s, 50]; the remaining sort only takes 22 steps.

Trying this test with Sort itself shows that it may be even more efficient in this regard.

SeedRandom[0];

s = RandomSample @ Range @ 2*^6;
s2 = Join @@ Sort /@ Partition[s, 1*^6];

Sort[s];  // RepeatedTiming // First
Sort[s2]; // RepeatedTiming // First

0.27
0.037