The important difference between sort(rand(...)) and this algorithm is that the sort(.) procedure is at best O(n log(n)) time complexity, while this algorithm is linear time, i.e. O(n). For short length vectors the sort(.) approach will be quicker (due to a more efficient implementation of the sort(.) function), but for large n the differences can be significant. As an example:

>> tic;sort(rand(1,10^7),2);toc
Elapsed time is 1.757046 seconds.

vs.

>> tic;genrandsorted(1,10^7);toc
Elapsed time is 0.531082 seconds.

I am guessing that if one implemented the Bentley & Saxe algorithm in C it would probably be comparable to sort(.) for smaller 'n' (and of course, faster again for large 'n').

the program is well written, but i believe the state space representation is wrong.
Though the implemented part is as mentioned in the refernce paper still , I argue about the 'R' in arma_ConvertToSS

Taking the MA parameters directly is not the proper representation.
instead R=[ c1-a1;c2- a2-( c1-a1);.....];

' SPECTRAL ESTIMATION FOR NOISY SIGNALS OBSERVED THROUGH A LINEAR SYSTEM' Check this paper

The important difference between sort(rand(...)) and this algorithm is that the sort(.) procedure is at best O(n log(n)) time complexity, while this algorithm is linear time, i.e. O(n). For short length vectors the sort(.) approach will be quicker (due to a more efficient implementation of the sort(.) function), but for large n the differences can be significant. As an example:

>> tic;sort(rand(1,10^7),2);toc
Elapsed time is 1.757046 seconds.

vs.

>> tic;genrandsorted(1,10^7);toc
Elapsed time is 0.531082 seconds.

I am guessing that if one implemented the Bentley & Saxe algorithm in C it would probably be comparable to sort(.) for smaller 'n' (and of course, faster again for large 'n').

I acknowledge that this is indeed a clever algorithm but in matlab
sort(rand(..)), as already suggested by Wolgang Schwanghart, is much to be prefered. Therefore this submission has only some educational value, for which the internal comments should be expanded.

What makes the function different to the command
x = sort(rand(n,p),2);
? I don't have the cited paper at hand, so it might be good to explain the difference, if there are any.