# Algorithm in RandomChoice[{w1,w2,…}->{e1,e2,…},n]?

I wonder if someone happens to know the actual algorithm implemented in the particular overload of RandomChoice that handles explicitly presented discrete distributions, namely

RandomChoice[{w1, w2, ...}->{e1, e1, ...}]
gives a pseudorandom choice weighted by the w1.


It's blazing fast, as illustrated in my test notebook here: https://www.dropbox.com/s/rwm80fut60v926b/FastNonUniformPseudoRandoms002.cdf

(btw, such things are really useful in, say, particle filters).

I know of at least four algorithms for solving this problem. A greedy-space array of outcomes has O(1) time but O(S) space where S is the sum of the frequency counts (and only works when the weights are integers or can be scaled up to integers). Two algorithms explicitly invert the CDF, either by linear search -- O(N) time, O(N) space, or binary search -- O(log N) time, O(N) space. The best I know is Walker's Alias Method http://en.wikipedia.org/wiki/Alias_method

which is O(N) space, O(1) time. I implemented this in my notebook and it beats the pants off the others (as expected) but is still 30 times slower than the built-in (though still just O(1)).

I wonder if the built-in is just an optimized Walker's or is it something better?