Radix sort function in R
Putting these comments into an answer...

As pointed out, it will be a tall order to beat or even approach the computation times of the native sort function, as it is compiled from C. Often with R, the trick to make your code faster consists in composing with some of these fast compiled building blocks that are native functions. In particular, we would look to substitute for loops like the one you have with vectorized functions. Unfortunately, your use of a for loop here is particular in the sense that each iteration has a side-effect to the x vector. This means that the order of the operations is important and we cannot run the loop iterations independently in parallel or via vectorized functions.

If we cannot get rid of the for loop, we are left trying to optimize the code within the body of the loop. Let's start with a profile of your code:

x &lt;- sample(10000)
Rprof(tmp &lt;- tempfile())
for (i in 1:10) z &lt;- radix_sort(x)
Rprof()
summaryRprof(tmp)$by.total
#                         total.time total.pct self.time self.pct
# "radix_sort"                  8.26     99.76      0.72     8.70
# "split"                       7.34     88.65      0.06     0.72
# "split.default"               7.28     87.92      0.54     6.52
# "as.factor"                   6.74     81.40      0.08     0.97
# "factor"                      6.64     80.19      1.72    20.77
# "as.character"                4.34     52.42      4.34    52.42
# "unique"                      0.42      5.07      0.04     0.48
# "unique.default"              0.38      4.59      0.38     4.59
# "%%"                          0.14      1.69      0.14     1.69
# "get_digit"                   0.14      1.69      0.00     0.00
# "sort.list"                   0.12      1.45      0.02     0.24
# "order"                       0.08      0.97      0.06     0.72
# "unlist"                      0.06      0.72      0.06     0.72
# [...]


We see that the main culprit here is the use of split. Another surprise here is the use of order (it seems to be applied when figuring out the levels of x_digit_i) which you could consider like cheating given your objective.

So what alternative do we have to your use of split/unsplit? Essentially, you have a vector x that you want to reorder based on a vector of digits x_digit_i. One way is to use outer to create a matrix of TRUE/FALSE where each column locates a different digit (for a total of ten columns):

z &lt;- outer(x_digit_i, 0:9, "==")


Then, you want to turn this matrix into a vector of indices, such that the first few indices will locate the zeroes, then the ones, etc (the equivalent of idx &lt;- order(x_digit_i). You can do so by doing:

idx &lt;- row(z)[z]


or (harder to understand but a bit faster)

idx &lt;- 1L + (which(z) - 1L) %% length(x)  


Finally, you just have to do:

x &lt;- x[idx]


Also note that since you are dealing with integers, your code might be a little faster (and more robust by avoiding floating point errors) if you make sure to use integers everywhere possible. In particular, your get_digit function could be rewritten as follows:

get_digit &lt;- function(x, d) (x %% as.integer(10^d)) %/% as.integer(10^(d-1))


The benchmarks below show decent progress (where my suggestions are implemented under the name sort_radix2) bridging the gap with the native sort. I hope it helps!

x &lt;- sample(100)
microbenchmark(radix_sort(x), radix_sort2(x), sort(x))
# Unit: microseconds
#            expr     min       lq       mean   median       uq      max neval
#   radix_sort(x) 964.692 972.3675 1025.35180 984.3775 1012.178 2233.397   100
#  radix_sort2(x) 250.642 256.5720  282.58952 261.2910  282.449 1266.061   100
#         sort(x)  82.270  86.1605   92.22669  88.0230   90.943  223.249   100

x &lt;- sample(10000)
microbenchmark(radix_sort(x), radix_sort2(x), sort(x))
# Unit: microseconds
#            expr       min         lq       mean     median         uq        max neval
#   radix_sort(x) 71939.706 76147.1715 80028.7541 78389.8140 81512.4140 144632.484   100
#  radix_sort2(x) 24218.810 27613.3190 34841.8724 29477.7115 31772.9415 143283.337   100
#         sort(x)   411.691   454.4015   563.4825   492.6165   558.0925   3412.719   100