Try dedicated shuffle_exsss

Author: RaimoNiskanen

lib/stdlib/src/rand.erl

@@ -1832,6 +1832,9 @@ and the [`NewState`](`t:state/0`).
               when
       List         :: list(),
       State        :: state().
+shuffle_s(List, {#{type:=exsss} = AlgHandler, S}) ->
+    {ShuffledList, _P1, NewS} = shuffle_s_exsss(List, S),
+    {ShuffledList, {AlgHandler, NewS}};
 shuffle_s(List, {#{type:=mwc59} = AlgHandler, CX0}) ->
     {ShuffledList, _P1, CX1} = shuffle_s_mwc59(List, CX0),
     {ShuffledList, {AlgHandler, CX1}};
@@ -2384,6 +2387,86 @@ exsplus_jump(S, [AS0|AS1], J, N) ->
             exsplus_jump(NS, [AS0|AS1], J bsr 1, N-1)
     end.
 
+%% -------
+
+shuffle_s_exsss(List, S) ->
+    shuffle_exsss_r(List, [], 1, S).
+
+%% Leaf cases - random permutations for 0..3 elements
+shuffle_exsss_r([], Acc, P, S) ->
+    {Acc, P, S};
+shuffle_exsss_r([X], Acc, P, S) ->
+    {[X | Acc], P, S};
+shuffle_exsss_r([X, Y], Acc, P, S) ->
+    shuffle_exsss_r_2(X, Acc, P, S, Y);
+shuffle_exsss_r([X, Y, Z], Acc, P, S) ->
+    shuffle_exsss_r_3(X, Acc, P, S, Y, Z);
+%% General case - split and recursive shuffle
+shuffle_exsss_r([_, _, _ | _] = List, Acc, P, S) ->
+    %% P and S is bitstream cache and state
+    shuffle_exsss_r(List, Acc, P, S, [], [], [], []).
+%%
+%% Split L into 4 random subsets
+%%
+shuffle_exsss_r([], Acc0, P0, S0, Zero, One, Two, Three) ->
+    %% Split done, recursively shuffle the splitted lists onto Acc
+    {Acc1, P1, S1} = shuffle_exsss_r(Zero, Acc0, P0, S0),
+    {Acc2, P2, S2} = shuffle_exsss_r(One, Acc1, P1, S1),
+    {Acc3, P3, S3} = shuffle_exsss_r(Two, Acc2, P2, S2),
+    shuffle_exsss_r(Three, Acc3, P3, S3);
+shuffle_exsss_r([X | L], Acc, P0, S, Zero, One, Two, Three)
+  when is_integer(P0), 3 < P0, P0 =< ?MASK(59) ->
+    P1 = P0 bsr 2,
+    case P0 band 3 of
+        0 -> shuffle_exsss_r(L, Acc, P1, S, [X | Zero], One, Two, Three);
+        1 -> shuffle_exsss_r(L, Acc, P1, S, Zero, [X | One], Two, Three);
+        2 -> shuffle_exsss_r(L, Acc, P1, S, Zero, One, [X | Two], Three);
+        3 -> shuffle_exsss_r(L, Acc, P1, S, Zero, One, Two, [X | Three])
+    end;
+shuffle_exsss_r([_ | _] = L, Acc, _P, [S1|S0], Zero, One, Two, Three) ->
+    S0_1 = ?MASK(58, S0),
+    S1_1 = ?exs_next(S0_1, S1, Tmp1),
+    P    = ?scramble_starstar(S0_1, Tmp2, Tmp3) bor ?BIT(58),
+    shuffle_exsss_r(L, Acc, P, [S0_1|S1_1], Zero, One, Two, Three).
+
+%% Permute 2 elements
+shuffle_exsss_r_2(X, Acc, P, S, Y)
+  when is_integer(P), 1 < P, P =< ?MASK(59) ->
+    {case P band 1 of
+         0 -> [Y, X | Acc];
+         1 -> [X, Y | Acc]
+     end, P bsr 1, S};
+shuffle_exsss_r_2(X, Acc, _P, [S1|S0], Y) ->
+    S0_1 = ?MASK(58, S0),
+    S1_1 = ?exs_next(S0_1, S1, Tmp1),
+    P    = ?scramble_starstar(S0_1, Tmp2, Tmp3) bor ?BIT(58),
+    shuffle_exsss_r_2(X, Acc, P, [S0_1|S1_1], Y).
+
+%% Permute 3 elements
+%%
+%% Uses 3 random bits per iteration with a probability of 1/4
+%% to reject and retry, which on average is 3 * 4/3
+%% (infinite sum of (1/4)^k) = 4 bits per permutation
+shuffle_exsss_r_3(X, Acc, P0, S, Y, Z)
+  when is_integer(P0), 7 < P0, P0 =< ?MASK(59) ->
+    P1 = P0 bsr 3,
+    case P0 band 7 of
+        0 -> {[Z, Y, X | Acc], P1, S};
+        1 -> {[Y, Z, X | Acc], P1, S};
+        2 -> {[Z, X, Y | Acc], P1, S};
+        3 -> {[X, Z, Y | Acc], P1, S};
+        4 -> {[Y, X, Z | Acc], P1, S};
+        5 -> {[X, Y, Z | Acc], P1, S};
+        _ -> % Reject and retry
+            shuffle_exsss_r_3(X, Acc, P1, S, Y, Z)
+    end;
+shuffle_exsss_r_3(X, Acc, _P, [S1|S0], Y, Z) ->
+    S0_1 = ?MASK(58, S0),
+    S1_1 = ?exs_next(S0_1, S1, Tmp1),
+    P    = ?scramble_starstar(S0_1, Tmp2, Tmp3) bor ?BIT(58),
+    shuffle_exsss_r_3(X, Acc, P, [S0_1|S1_1], Y, Z).
+
+
 %% =====================================================================
 %% exs1024 PRNG: Xorshift1024*
 %% Algorithm by Sebastiano Vigna