Quick QuickSort implementation

Question

I am playing with julialang and loving it, but at the same time noticed that they have a benchmark comparison with Mathematica. I have already submited a better version of recursion_fibonacci and now am looking into recursion_quicksort. I would like to find an implementation in Mathematica that is comparable to C-lang.

Currently, the benchmark uses the following code

(* numeric vector sort *)

ClearAll[qsort];
(* qsort[ain_, loin_, hiin_] := Module[
    {a = ain, i = loin, j = hiin, lo = loin, hi = hiin, pivot},
    While[ i < hi,
        pivot = a[[BitShiftRight[lo + hi] ]];
        While[ i <= j,
            While[a[[i]] < pivot, i++];
            While[a[[j]] > pivot, j--];
            If[ i <= j,
                a[[{i,j}]] = a[[{j, i}]];
                i++; j--;
            ];
        ];
        If[ lo < j, a = qsort[a, lo, j] ];
        {lo, j} = {i, hi};
    ];
    a
]; *)
qsort = Compile[
    {{ain, _Real, 1}, {loin, _Integer}, {hiin, _Integer}},
    Module[
        {a = ain, i = loin, j = hiin, lo = loin, hi = hiin, pivot},
        While[ i < hi,
            pivot = a[[ Floor[(lo + hi)/2] ]];
            While[ i <= j,
                While[a[[i]] < pivot, i++];
                While[a[[j]] > pivot, j--];
                If[ i <= j,
                    a[[{i,j}]] = a[[{j, i}]];
                    i++; j--;
                ];
            ];
            If[ lo < j, a[[lo;;j]] = qsort[ a[[lo;;j]], 1, j - lo + 1] ];
            {lo, j} = {i, hi};
        ];
        a
    ]
];


ClearAll[sortperf];
sortperf[n_] := Module[{vec = RandomReal[1, n]}, qsort[vec, 1, n]];

test[OrderedQ[sortperf[5000]] ];
timeit[sortperf[5000], "recursion_quicksort"];

where there is compiled and uncompiled versions of quicksort algorithm. On my laptop the compiled version takes 10.3ms, while uncompiled version takes 137.8ms. I think there is space for improvement since the inbuilt method Sort[] takes only 0.379ms.

How do we speed-up the quicksort algorithm? Bonus points if we don't use Compile[]

Helper functions to run the code above

Needs["CCompilerDriver`"];
If[ Length[CCompilers[]] > 0,
    $CompilationTarget = "C"
];

ClearAll[timeit];
SetAttributes[timeit, HoldFirst];
timeit[ex_, name_String] := Module[
    {t},
    t = Infinity;
    Do[
        t = Min[t, N[First[AbsoluteTiming[ex]]]];
        ,
        {i, 1, 5}
    ];
    If[$printOutput, Print["mathematica,", name, ",", t*1000];  ];
];

ClearAll[test];
SetAttributes[test, HoldFirst];
test[ex_] := Assert[ex];
On[Assert];

Answer 1

I also got angry about those randomly picked and ill-implemented benchmarks by the Julia team. I appreciate their efforts (jit compilers are useful), but the Fibonacci example was straight away ridiculous.

Here is a compiled quick sort implementation that employs a stack in order to avoid recursive calls. The problem with recursion in CompiledFunctions is that it needs to call the main evaluator - and such calls are notoriously slow. I synthesized this implementation from the benchmark implementation and from this guide on geeksforgeeks.org.

Admittedly, one may argue whether this is still a "pure" Mathematica implementation or not (imho, it is rather C), but the Julia team also employed Compile. Hence, once started, one has to do that correctly.

quicksortiterative = Compile[{{ain, _Real, 1}},
   Module[{iter, l, h, a, stack, top, x, i, j, ai, aj},
    a = ain;
    top = 0;
    stack = Table[0, {Length[a] + 1}];
    stack[[++top]] = 1;
    stack[[++top]] = Length[a];
    While[top >= 1,
     h = Compile`GetElement[stack, top--];
     l = Compile`GetElement[stack, top--];
     x = Compile`GetElement[a, l + Quotient[h - l, 2]];
     i = l;
     j = h;
     While[i <= j,
      ai = Compile`GetElement[a, i];
      While[ai < x, ai = Compile`GetElement[a, ++i]];
      aj = Compile`GetElement[a, j];
      While[aj > x, aj = Compile`GetElement[a, --j]];
      If[i <= j,
       a[[i++]] = aj;
       a[[j--]] = ai;
       ];
      ];
     If[j > l, stack[[++top]] = l; stack[[++top]] = j];
     If[i < h, stack[[++top]] = i; stack[[++top]] = h];
     ];
    a
    ],
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"
   ];

To my own surpise, this actually beats Mathematica's built-in Sort function:

n = 1000000;
a = RandomReal[{-1, 1}, n];
r1 = quicksortiterative[a]; // RepeatedTiming // First
r2 = Sort[a]; // RepeatedTiming // First
r1 == r2

0.105

0.127

True