otpbp

-module(otpbp_rand).

-ifndef(HAVE_rand__uniform_real_0).
% OTP 21.0
-export([uniform_real/0]).
-endif.
-ifndef(HAVE_rand__uniform_real_s_1).
% OTP 21.0
-export([uniform_real_s/1]).
-endif.
-ifndef(HAVE_rand__splitmix64_next_1).
% OTP 25.0
-export([splitmix64_next/1]).
-endif.
-ifndef(HAVE_rand__exsp_next_1).
% OTP 25.0
-export([exsp_next/1]).
-endif.
-ifndef(HAVE_rand__exsp_jump_1).
% OTP 25.0
-export([exsp_jump/1]).
-endif.
-ifndef(HAVE_rand__mwc59_1).
% OTP 25.0
-export([mwc59/1]).
-endif.
-ifndef(HAVE_rand__mwc59_value32_1).
% OTP 25.0
-export([mwc59_value32/1]).
-endif.
-ifndef(HAVE_rand__mwc59_value_1).
% OTP 25.0
-export([mwc59_value/1]).
-endif.
-ifndef(HAVE_rand__mwc59_float_1).
% OTP 25.0
-export([mwc59_float/1]).
-endif.
-ifndef(HAVE_rand__mwc59_seed_0).
% OTP 25.0
-export([mwc59_seed/0]).
-endif.
-ifndef(HAVE_rand__mwc59_seed_1).
% OTP 25.0
-export([mwc59_seed/1]).
-endif.
-ifndef(HAVE_rand__bytes_1).
% OTP 24.0
-export([bytes/1]).
-endif.
-ifndef(HAVE_rand__bytes_s_2).
% OTP 24.0
-export([bytes_s/2]).
-endif.
-ifndef(HAVE_rand__bytes_s_2).
% OTP 24.0
-export([bytes_s/2]).
-endif.
-ifndef(HAVE_rand__shuffle_1).
% OTP 29.0
-export([shuffle/1]).
-endif.
-ifndef(HAVE_rand__shuffle_s_2).
% OTP 29.0
-export([shuffle_s/2]).
-endif.

-ifndef(HAVE_rand__uniform_real_0).
-ifdef(HAVE_rand__uniform_real_s_1).
-import(rand, [uniform_real_s/1]).
-endif.
-endif.
-ifndef(HAVE_rand__exsp_jump_1).
-ifdef(HAVE_rand__exsp_next_1).
-import(rand, [exsp_next/1]).
-endif.
-endif.
-ifndef(HAVE_rand__bytes_1).
-ifdef(HAVE_rand__bytes_s_2).
-import(rand, [bytes_s/2]).
-endif.
-endif.
-ifndef(HAVE_rand__shuffle_1).
-ifdef(HAVE_rand__shuffle_s_2).
-import(rand, [shuffle_s/2]).
-endif.
-endif.

-define(SEED_DICT, rand_seed).
-define(BIT(Bits), (1 bsl (Bits))).
-define(MASK(Bits), (?BIT(Bits) - 1)).
-define(MASK(Bits, X), ((X) band ?MASK(Bits))).
-define(BSL(Bits, X, N), (?MASK((Bits) - (N), (X)) bsl (N))).

-define(MWC59_XS1, 4).
-define(MWC59_XS2, 27).

-ifndef(HAVE_rand__bytes_1).
bytes(N) ->
    {Bytes, State} = bytes_s(N, seed_get()),
    seed_put(State),
    Bytes.

-ifndef(NEED_seed_get_0).
-define(NEED_seed_get_0, true).
-endif.
-ifndef(NEED_seed_put_1).
-define(NEED_seed_put_1, true).
-endif.
-endif.

-ifndef(HAVE_rand__bytes_s_2).
bytes_s(N, {#{bits := Bits, next := Next} = AlgHandler, R}) when is_integer(N), 0 =< N ->
    bytes_r(N, AlgHandler, Next, R, Bits, maps:get(weak_low_bits, AlgHandler, 0));
bytes_s(N, {#{max := Mask, next := Next} = AlgHandler, R}) when is_integer(N), 0 =< N, ?MASK(58) =< Mask ->
    %% Old spec - assume 58 bits and 2 weak low bits
    %% giving 56 bits i.e precisely 7 bytes per generated number
    bytes_r(N, AlgHandler, Next, R, 58, 2).

-compile({inline, bytes_r/6}).
bytes_r(N, AlgHandler, Next, R, Bits, WeakLowBits) ->
    %% We use whole bytes from each generator word,
    %% GoodBytes: that number of bytes
    GoodBytes = (Bits - WeakLowBits) bsr 3,
    GoodBits = GoodBytes bsl 3,
    %% Shift: how many bits of each generator word to waste by shifting right - we use the bits from the big end
    bytes_r(N, AlgHandler, Next, R, <<>>, GoodBytes, GoodBits, Bits - GoodBits).

bytes_r(N0, AlgHandler, Next, R0, Bytes0, GoodBytes, GoodBits, Shift) when GoodBytes < N0 ->
    {V, R1} = Next(R0),
    bytes_r(N0 - GoodBytes, AlgHandler, Next, R1, <<Bytes0/binary, (V bsr Shift):GoodBits>>, GoodBytes, GoodBits, Shift);
bytes_r(N, AlgHandler, Next, R0, Bytes, _GoodBytes, GoodBits, _Shift) ->
    {V, R1} = Next(R0),
    Bits = N bsl 3,
    %% Use the big end bits
    {<<Bytes/binary, (V bsr (GoodBits - Bits)):Bits>>, {AlgHandler, R1}}.
-endif.

-ifndef(HAVE_rand__uniform_real_0).
uniform_real() ->
    {X, Seed} = uniform_real_s(seed_get()),
    seed_put(Seed),
    X.

-ifndef(NEED_seed_get_0).
-define(NEED_seed_get_0, true).
-endif.
-ifndef(NEED_seed_put_1).
-define(NEED_seed_put_1, true).
-endif.
-endif.

-ifndef(HAVE_rand__uniform_real_s_1).
-define(BC(V, N), bc((V), ?BIT((N) - 1), N)).

uniform_real_s({#{next := Next} = Alg, R0}) ->
    {V1, R1} = Next(R0),
    uniform_real_s(Alg, V1, R1, Next).

-compile({inline, uniform_real_s/5}).
uniform_real_s(#{bits := Bits} = Alg, V1, R1, Next) ->
    case V1 bsr (Bits - 56) of
        M1 when ?BIT(55) =< M1 -> %% We have 56 bits - waste 3
            {(M1 bsr 3) * math:pow(2.0, -53), {Alg, R1}};
        M1 when ?BIT(54) =< M1 -> %% We have 55 bits - waste 2
            {(M1 bsr 2) * math:pow(2.0, -54), {Alg, R1}};
        M1 when ?BIT(53) =< M1 -> %% We have 54 bits - waste 1
            {(M1 bsr 1) * math:pow(2.0, -55), {Alg, R1}};
        M1 when ?BIT(52) =< M1 -> %% We have 53 bits - use all
            {M1 * math:pow(2.0, -56), {Alg, R1}};
        M1 -> %% Need more bits
            {V2, R2} = Next(R1),
            uniform_real_s(Alg, Next, M1, -56, R2, V2, Bits)
    end;
uniform_real_s(#{max := _} = Alg, V1, R1, Next) ->
    case ?MASK(56, V1) of
        M1 when ?BIT(55) =< M1 -> %% We have 56 bits - waste 3
            {(M1 bsr 3) * math:pow(2.0, -53), {Alg, R1}};
        M1 when ?BIT(54) =< M1 -> %% We have 55 bits - waste 2
            {(M1 bsr 2) * math:pow(2.0, -54), {Alg, R1}};
        M1 when ?BIT(53) =< M1 -> %% We have 54 bits - waste 1
            {(M1 bsr 1) * math:pow(2.0, -55), {Alg, R1}};
        M1 when ?BIT(52) =< M1 -> %% We have 53 bits - use all
            {M1 * math:pow(2.0, -56), {Alg, R1}};
        M1 -> %% Need more bits
            {V2, R2} = Next(R1),
            uniform_real_s(Alg, Next, M1, -56, R2, V2, 56)
    end.

uniform_real_s(Alg, _Next, M0, -1064, R1, V1, Bits) -> % 19*56
    B0 = (53 - ?BC(M0, 52)), % Missing bits
    {(((M0 bsl B0) bor (V1 bsr (Bits - B0))) * math:pow(2.0, -1064 - B0)), {Alg, R1}};
uniform_real_s(Alg, _Next, M0, BitNo, R1, V1, Bits) when ?BIT(51) =< M0 ->
    {(((M0 bsl 1) bor (V1 bsr (Bits - 1))) * math:pow(2.0, BitNo - 1)), {Alg, R1}};
uniform_real_s(Alg, _Next, M0, BitNo, R1, V1, Bits) when ?BIT(50) =< M0 ->
    {(((M0 bsl 2) bor (V1 bsr (Bits - 2))) * math:pow(2.0, BitNo - 2)), {Alg, R1}};
uniform_real_s(Alg, _Next, M0, BitNo, R1, V1, Bits) when ?BIT(49) =< M0 ->
    {(((M0 bsl 3) bor (V1 bsr (Bits - 3))) * math:pow(2.0, BitNo - 3)), {Alg, R1}};
uniform_real_s(Alg, Next, 0, BitNo, R1, V1, Bits) ->
    case V1 bsr (Bits - 56) of
        M1 when ?BIT(55) =< M1 -> %% We have 56 bits - waste 3
            {(M1 bsr 3) * math:pow(2.0, BitNo - 53), {Alg, R1}};
        M1 when ?BIT(54) =< M1 -> %% We have 55 bits - waste 2
            {(M1 bsr 2) * math:pow(2.0, BitNo - 54), {Alg, R1}};
        M1 when ?BIT(53) =< M1 -> %% We have 54 bits - waste 1
            {(M1 bsr 1) * math:pow(2.0, BitNo - 55), {Alg, R1}};
        M1 when ?BIT(52) =< M1 -> %% We have 53 bits - use all
            {M1 * math:pow(2.0, BitNo - 56), {Alg, R1}};
        M1 when BitNo =/= -1008; ?BIT(42) =< M1 -> %% Need more bits; We have 43 bits - get the last bits
            uniform_real_s(Alg, Next, M1, BitNo - 56, R1);
        _ -> uniform_real_s({Alg, R1})
    end;
uniform_real_s(Alg, _Next, M0, BitNo, R1, V1, Bits) ->
    %% Fill up to 53 bits
    B0 = 53 - ?BC(M0, 49), % Number of bits we need to append
    {(((M0 bsl B0) bor (V1 bsr (Bits - B0))) * math:pow(2.0, BitNo - B0)), {Alg, R1}}.

uniform_real_s(Alg, Next, M0, BitNo, R0) ->
    {V1, R1} = Next(R0),
    case Alg of
        #{bits := Bits} -> uniform_real_s(Alg, Next, M0, BitNo, R1, V1, Bits);
        #{max := _} -> uniform_real_s(Alg, Next, M0, BitNo, R1, ?MASK(56, V1), 56)
    end.

bc(V, B, N) when B =< V -> N;
bc(V, B, N) -> bc(V, B bsr 1, N - 1).
-endif.

-ifndef(HAVE_rand__splitmix64_next_1).
splitmix64_next(X_0) ->
    X = ?MASK(64, X_0 + 16#9e3779b97f4a7c15),
    Z_0 = ?MASK(64, (X bxor (X bsr 30)) * 16#bf58476d1ce4e5b9),
    Z_1 = ?MASK(64, (Z_0 bxor (Z_0 bsr 27)) * 16#94d049bb133111eb),
    {?MASK(64, Z_1 bxor (Z_1 bsr 31)), X}.
-endif.

-ifndef(HAVE_rand__exsp_next_1).
-define(exs_next(S0, S1, S1_b),
        begin
        S1_b = ?MASK(58, S1) bxor ?BSL(58, S1, 24),
        S1_b bxor S0 bxor (S1_b bsr 11) bxor (S0 bsr 41)
        end).

exsp_next([S1|S0]) ->
    S0_1 = ?MASK(58, S0),
    NewS1 = ?exs_next(S0_1, S1, S1_b),
    {?MASK(58, S0_1 + NewS1), [S0_1|NewS1]}.
-endif.

-ifndef(HAVE_rand__exsp_jump_1).
-define(JUMPCONST1, 16#02f8ea6bc32c797).
-define(JUMPCONST2, 16#345d2a0f85f788c).
-define(JUMPELEMLEN, 58).

exsp_jump(S) ->
    {S1, AS1} = exsplus_jump(S, [0|0], ?JUMPCONST1, ?JUMPELEMLEN),
    {_,  AS2} = exsplus_jump(S1, AS1,  ?JUMPCONST2, ?JUMPELEMLEN),
    AS2.

-dialyzer({no_improper_lists, exsplus_jump/4}).
exsplus_jump(S, AS, _, 0) -> {S, AS};
exsplus_jump(S, [AS0|AS1], J, N) ->
    {_, NS} = exsp_next(S),
    exsplus_jump(NS,
                 case ?MASK(1, J) of
                     1 ->
                         [S0|S1] = S,
                         [AS0 bxor S0|AS1 bxor S1];
                     0 -> [AS0|AS1]
                 end,
                 J bsr 1,
                 N - 1).
-endif.

-ifndef(HAVE_rand__mwc59_1).
-define(MWC59_A, 16#7fa6502).
-define(MWC59_B, 32).
mwc59(CX0) ->
    CX = ?MASK(59, CX0),
    ?MWC59_A * ?MASK(?MWC59_B, CX) + CX bsr ?MWC59_B.
-endif.

-ifndef(HAVE_rand__mwc59_value32_1).
-define(MWC59_XS, 8).
mwc59_value32(CX1) ->
    CX = ?MASK(32, CX1),
    CX bxor ?BSL(32, CX, ?MWC59_XS).
-endif.

-ifndef(HAVE_rand__mwc59_value_1).
mwc59_value(CX1) ->
    CX = ?MASK(59, CX1),
    CX2 = CX bxor ?BSL(59, CX, ?MWC59_XS1),
    CX2 bxor ?BSL(59, CX2, ?MWC59_XS2).
-endif.

-ifndef(HAVE_rand__mwc59_float_1).
-define(TWO_POW_MINUS53, 1.11022302462515657e-16).
mwc59_float(CX1) ->
    CX = ?MASK(53, CX1),
    CX2 = CX bxor ?BSL(53, CX, ?MWC59_XS1),
    (CX2 bxor ?BSL(53, CX2, ?MWC59_XS2)) * ?TWO_POW_MINUS53.
-endif.

-ifndef(HAVE_rand__mwc59_seed_0).
mwc59_seed() ->
    {A1, A2, A3} = default_seed(),
    hash58(A1) bxor hash58(A2) bxor hash58(A3) + 1.

-compile({inline, default_seed/0}).
default_seed() -> {erlang:phash2([{node(), self()}]), erlang:system_time(), erlang:unique_integer()}.

-ifndef(NEED_hash58_1).
-define(NEED_hash58_1, true).
-endif.
-endif.

-ifndef(HAVE_rand__mwc59_seed_1).
mwc59_seed(S) when is_integer(S), 0 =< S, S =< ?MASK(58) -> hash58(S) + 1.

-ifndef(NEED_hash58_1).
-define(NEED_hash58_1, true).
-endif.
-endif.

-ifndef(HAVE_rand__shuffle_1).
shuffle(List) ->
    {ShuffledList, State} = shuffle_s(List, seed_get()),
    seed_put(State),
    ShuffledList.

-ifndef(NEED_seed_get_0).
-define(NEED_seed_get_0, true).
-endif.
-ifndef(NEED_seed_put_1).
-define(NEED_seed_put_1, true).
-endif.
-endif.

-ifndef(HAVE_rand__shuffle_s_2).
shuffle_s(List, {#{next := Next} = AlgHandler, R0}) when is_list(List) ->
    [P0|S0] = shuffle_init_bitstream(R0, Next,
                                     case AlgHandler of
                                         #{bits := _} -> maps:get(weak_low_bits, AlgHandler, 0);
                                         #{max := _} -> 2 % Old spec - assume 2 weak low bits
                                     end),
    {ShuffledList, _P1, [R1|_]} = shuffle_r(List, [], P0, S0),
    {ShuffledList, {AlgHandler, R1}}.

-dialyzer({no_improper_lists, shuffle_init_bitstream/3}).
-compile({inline, shuffle_init_bitstream/3}).
shuffle_init_bitstream(R, Next, WeakLowBits) -> [1|[R|{Next, WeakLowBits}]].

-dialyzer({no_improper_lists, shuffle_new_bits/1}).
shuffle_new_bits([R0|{Next, WeakLowBits} = W]) ->
    {V, R1} = Next(R0),
    %% Setting the top bit M here provides the marker for when we are out of random bits: P =:= 1
    M = 1 bsl 56,
    [((V bsr WeakLowBits) band (M - 1)) bor M|[R1|W]].

%% Leaf cases - random permutations for 0..3 elements
shuffle_r([], Acc, P, S) -> {Acc, P, S};
shuffle_r([X], Acc, P, S) -> {[X|Acc], P, S};
shuffle_r([X, Y], Acc, P, S) -> shuffle_r(X, Acc, P, S, Y);
shuffle_r([X, Y, Z], Acc, P, S) -> shuffle_r(X, Acc, P, S, Y, Z);
%% General case - split and recursive shuffle
shuffle_r([_, _, _|_] = List, Acc, P, S) ->
    %% P and S is bitstream cache and state
    shuffle_r(List, Acc, P, S, [], [], [], []).
%% Split L into 4 random subsets
shuffle_r([], Acc0, P0, S0, Zero, One, Two, Three) ->
    %% Split done, recursively shuffle the splitted lists onto Acc
    {Acc1, P1, S1} = shuffle_r(Zero, Acc0, P0, S0),
    {Acc2, P2, S2} = shuffle_r(One, Acc1, P1, S1),
    {Acc3, P3, S3} = shuffle_r(Two, Acc2, P2, S2),
    shuffle_r(Three, Acc3, P3, S3);
shuffle_r([X|L], Acc, P0, S, Zero, One, Two, Three) when is_integer(P0), 3 < P0, P0 < 1 bsl 57 ->
    P1 = P0 bsr 2,
    case P0 band 3 of
        0 -> shuffle_r(L, Acc, P1, S, [X|Zero], One, Two, Three);
        1 -> shuffle_r(L, Acc, P1, S, Zero, [X|One], Two, Three);
        2 -> shuffle_r(L, Acc, P1, S, Zero, One, [X|Two], Three);
        3 -> shuffle_r(L, Acc, P1, S, Zero, One, Two, [X|Three])
    end;
shuffle_r([_|_] = L, Acc, _P, S0, Zero, One, Two, Three) ->
    [P|S1] = shuffle_new_bits(S0),
    shuffle_r(L, Acc, P, S1, Zero, One, Two, Three).

%% Permute 2 elements
shuffle_r(X, Acc, P, S, Y) when is_integer(P), 1 < P, P < 1 bsl 57 ->
    {case P band 1 of
         0 -> [Y, X|Acc];
         1 -> [X, Y|Acc]
     end,
     P bsr 1,
     S};
shuffle_r(X, Acc, _P, S0, Y) ->
    [P|S1] = shuffle_new_bits(S0),
    shuffle_r(X, Acc, P, S1, 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_r(X, Acc, P0, S, Y, Z) when is_integer(P0), 7 < P0, P0 < 1 bsl 57 ->
    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};
        _ -> shuffle_r(X, Acc, P1, S, Y, Z) % Reject and retry
    end;
shuffle_r(X, Acc, _P, S0, Y, Z) ->
    [P|S1] = shuffle_new_bits(S0),
    shuffle_r(X, Acc, P, S1, Y, Z).
-endif.

-ifdef(NEED_seed_get_0).
-define(DEFAULT_ALG_HANDLER, exrop).
seed_get() ->
    case get(?SEED_DICT) of
        undefined -> rand:seed(?DEFAULT_ALG_HANDLER);
        Old -> Old  % no type checking here
    end.
-endif.

-ifdef(NEED_seed_put_1).
seed_put(Seed) ->
    put(?SEED_DICT, Seed),
    Seed.
-endif.

-ifdef(NEED_hash58_1).
hash58(X) ->
    X0 = ?MASK(58, X),
    X1 = ?MASK(58, (X0 bxor (X0 bsr 29)) * 16#351afd7ed558ccd),
    X2 = ?MASK(58, (X1 bxor (X1 bsr 29)) * 16#0ceb9fe1a85ec53),
    X2 bxor (X2 bsr 29).
-endif.