Explain the shuffle algorithm

Author: RaimoNiskanen

lib/stdlib/src/rand.erl

@@ -1328,15 +1328,50 @@ shuffle_s(List, State) when is_list(List) ->
     shuffle_r(List, State, []).
 
 %% Random-split-and-shuffle algorithm suggested by Richard A. O'Keefe
-%% on ErlangForums, as I interpreted it...
+%% on ErlangForums, as I interpreted it...  "basically a randomized
+%% quicksort", shall we name it Quickshuffle?
 %%
-%% Randomly split the list in two lists,
-%% recursively shuffle the two smaller lists,
-%% randomize the order between the lists according to their size.
+%% Randomly split the list in two lists, and recursively shuffle
+%% the two smaller lists.
 %%
-%% This is equivalent to assigning a random number to each
-%% element and sorting, but extending the numbers on demand
-%% while there still are duplicates.
+%% How the algorithm works and why it is correct can be explained like this:
+%%
+%% The objective is, given a list of elements, to return a random
+%% permutation of those elements so that every possible permutation
+%% has the same probability to be returned.
+%%
+%% One of the two correct and bias free algorithms described on the Wikipedia
+%% page for Fisher-Yates shuffling is to assign a random number to each
+%% element in the list and order the elements according to the numbers.
+%% For this to be correct the generated numbers must not have duplicates.
+%%
+%% This algorithm does that, but assigning a number and ordering
+%% the elements is more or less the same step, which is taken
+%% one binary bit at the time.
+%%
+%% It can be seen as, to each element, assign a fixpoint number
+%% of infinite length starting with bit weight 1/2, continuing with 1/4,
+%% and so on..., but reveal it incrementally.
+%%
+%% The list to shuffle is traversed, and a random bit is generated
+%% for each element.  If it is a 0, the element is assigned the zero bit
+%% by moving it to the head of the list Zero, and if it is a 1,
+%% to the head of the list One.
+%%
+%% Then the list Zero is recursively shuffled onto the accumulator list Acc,
+%% after that the list One.  By that all elements in Zero are ordered
+%% before the ones in One, according to the generated numbers.
+%% The order is actually not important as long as it is consistent.
+%%
+%% The algorithm recurses until the Zero or One list has length
+%% 0 or 1, which is when the generated fixpoint number has no duplicate.
+%% The fixpoint number in itself only exists in the guise of the
+%% recursion call stack, that is whether an element belongs to list
+%% Zero or One on each recursion level.
+%%
+%% As an optimization, since the algorithm is equivalent to its objective
+%% to randomly permute a list, we can when reaching a small enough list
+%% as in 4 or less instead do an explicit random permutation of the list.
 
 %% Leaf cases - random permutations for 0..4 elements
 shuffle_r([], State, Acc) ->
@@ -1387,35 +1422,34 @@ shuffle_r([X, Y, Z, Q], State0, Acc) ->
          23 -> [Y, X, Z, Q | Acc];
          24 -> [X, Y, Z, Q | Acc]
      end, State1};
-%%
 %% General case - split and recursive shuffle
 shuffle_r([_, _, _, _ | _] = List, State0, Acc0) ->
-    {Left, Right, State1} = shuffle_split(List, State0),
-    {Acc1, State2} = shuffle_r(Left, State1, Acc0),
-    shuffle_r(Right, State2, Acc1).
+    {Zero, One, State1} = shuffle_split(List, State0),
+    {Acc1, State2} = shuffle_r(Zero, State1, Acc0),
+    shuffle_r(One, State2, Acc1).
 
-%% Split L into two random subsets: Left and Right
+%% Split L into two random subsets: Zero and One
 %%
 shuffle_split(L, State) ->
     shuffle_split(L, State, 1, [], []).
 %%
-shuffle_split([], State, _P, Left, Right) ->
-    {Left, Right, State};
-shuffle_split([_ | _] = L, State0, 1, Left, Right) ->
+shuffle_split([], State, _P, Zero, One) ->
+    {Zero, One, State};
+shuffle_split([_ | _] = L, State0, 1, Zero, One) ->
     M = 1 bsl 56,
-    case rand:uniform_s(M, State0) of
+    case uniform_s(M, State0) of
         {V, State1} when is_integer(V), 1 =< V, V =< M ->
             %% Setting the top bit M here provides the marker
             %% for when we are out of random bits: P =:= 1
-            shuffle_split(L, State1, (V - 1) + M, Left, Right)
+            shuffle_split(L, State1, (V - 1) + M, Zero, One)
     end;
-shuffle_split([X | L], State, P, Left, Right)
-  when is_integer(P), 1 =< P, P < 1 bsl 57 ->
+shuffle_split([X | L], State, P, Zero, One)
+  when is_integer(P), 1 < P, P < 1 bsl 57 ->
     case P band 1 of
         0 ->
-            shuffle_split(L, State, P bsr 1, [X | Left], Right);
+            shuffle_split(L, State, P bsr 1, [X | Zero], One);
         1 ->
-            shuffle_split(L, State, P bsr 1, Left, [X | Right])
+            shuffle_split(L, State, P bsr 1, Zero, [X | One])
     end.
 
 %% =====================================================================