%%
%% %CopyrightBegin%
%%
%% Copyright Ericsson AB 2015. All Rights Reserved.
%%
%% The contents of this file are subject to the Erlang Public License,
%% Version 1.1, (the "License"); you may not use this file except in
%% compliance with the License. You should have received a copy of the
%% Erlang Public License along with this software. If not, it can be
%% retrieved online at http://www.erlang.org/.
%%
%% Software distributed under the License is distributed on an "AS IS"
%% basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See
%% the License for the specific language governing rights and limitations
%% under the License.
%%
%% %CopyrightEnd%
%%
%% =====================================================================
%% Multiple PRNG module for Erlang/OTP
%% Copyright (c) 2015 Kenji Rikitake
%% =====================================================================
-module(rand).
-export([seed_s/1, seed_s/2, seed/1, seed/2,
export_seed/0, export_seed_s/1,
uniform/0, uniform/1, uniform_s/1, uniform_s/2]).
-compile({inline, [exs64_next/1, exsplus_next/1,
exs1024_next/1, exs1024_calc/2]}).
-define(DEFAULT_ALG_HANDLER, exsplus).
-define(SEED_DICT, rand_seed).
%% =====================================================================
%% Types
%% =====================================================================
%% This depends on the algorithm handler function
-opaque alg_seed() :: exs64_state() | exsplus_state() | exs1024_state().
%% This is the algorithm handler function within this module
-type alg_handler() :: #{type => alg(),
max => integer(),
uniform => fun(),
uniform_n => fun()}.
%% Internal state
-type state() :: {alg_handler(), alg_seed()}.
-type alg() :: exs64 | exsplus | exs1024.
-export_type([alg/0, alg_handler/0, state/0, alg_seed/0]).
%% =====================================================================
%% API
%% =====================================================================
%% Return algorithm and seed so that RNG state can be recreated with seed/1
-spec export_seed() -> undefined | {alg(), alg_seed()}.
export_seed() ->
case seed_get() of
{#{type:=Alg}, Seed} -> {Alg, Seed};
_ -> undefined
end.
-spec export_seed_s(state()) -> {alg(), alg_seed()}.
export_seed_s({#{type:=Alg}, Seed}) -> {Alg, Seed}.
%% seed(Alg) seeds RNG with runtime dependent values
%% and return the NEW state
%% seed({Alg,Seed}) setup RNG with a previously exported seed
%% and return the NEW state
-spec seed(alg() | {alg(), alg_seed()}) -> state().
seed(Alg) ->
R = seed_s(Alg),
_ = seed_put(R),
R.
-spec seed_s(alg() | {alg(), alg_seed()}) -> state().
seed_s(Alg) when is_atom(Alg) ->
seed_s(Alg, {erlang:phash2([{node(),self()}]),
erlang:system_time(),
erlang:unique_integer()});
seed_s({Alg0, Seed}) ->
{Alg,_SeedFun} = mk_alg(Alg0),
{Alg, Seed}.
%% seed/2: seeds RNG with the algorithm and given values
%% and returns the NEW state.
-spec seed(Alg :: alg(), {integer(), integer(), integer()}) -> state().
seed(Alg0, S0) ->
State = seed_s(Alg0, S0),
_ = seed_put(State),
State.
-spec seed_s(Alg :: alg(), {integer(), integer(), integer()}) -> state().
seed_s(Alg0, S0 = {_, _, _}) ->
{Alg, Seed} = mk_alg(Alg0),
AS = Seed(S0),
{Alg, AS}.
%%% uniform/0, uniform/1, uniform_s/1, uniform_s/2 are all
%%% uniformly distributed random numbers.
%% uniform/0: returns a random float X where 0.0 < X < 1.0,
%% updating the state in the process dictionary.
-spec uniform() -> float().
uniform() ->
{X, Seed} = uniform_s(seed_get()),
_ = seed_put(Seed),
X.
%% uniform/1: given an integer N >= 1,
%% uniform/1 returns a random integer X where 1 =< X =< N,
%% updating the state in the process dictionary.
-spec uniform(N :: pos_integer()) -> pos_integer().
uniform(N) ->
{X, Seed} = uniform_s(N, seed_get()),
_ = seed_put(Seed),
X.
%% uniform_s/1: given a state, uniform_s/1
%% returns a random float X where 0.0 < X < 1.0,
%% and a new state.
-spec uniform_s(state()) -> {float(), NewS :: state()}.
uniform_s(State = {#{uniform:=Uniform}, _}) ->
Uniform(State).
%% uniform_s/2: given an integer N >= 1 and a state, uniform_s/2
%% uniform_s/2 returns a random integer X where 1 =< X =< N,
%% and a new state.
-spec uniform_s(N::pos_integer(), state()) -> {pos_integer(), NewS::state()}.
uniform_s(N, State = {#{uniform_n:=Uniform, max:=Max}, _})
when 0 < N, N =< Max ->
Uniform(N, State);
uniform_s(N, State0 = {#{uniform:=Uniform}, _})
when is_integer(N), 0 < N ->
{F, State} = Uniform(State0),
{trunc(F * N) + 1, State}.
%% =====================================================================
%% Internal functions
-define(UINT21MASK, 16#00000000001fffff).
-define(UINT32MASK, 16#00000000ffffffff).
-define(UINT33MASK, 16#00000001ffffffff).
-define(UINT39MASK, 16#0000007fffffffff).
-define(UINT58MASK, 16#03ffffffffffffff).
-define(UINT64MASK, 16#ffffffffffffffff).
-type uint64() :: 0..16#ffffffffffffffff.
-type uint58() :: 0..16#03ffffffffffffff.
-spec seed_put(state()) -> undefined | state().
seed_put(Seed) ->
put(?SEED_DICT, Seed).
seed_get() ->
case get(?SEED_DICT) of
undefined -> seed(?DEFAULT_ALG_HANDLER);
Old -> Old % no type checking here
end.
%% Setup alg record
mk_alg(exs64) ->
{#{type=>exs64, max=>?UINT64MASK,
uniform=>fun exs64_uniform/1, uniform_n=>fun exs64_uniform/2},
fun exs64_seed/1};
mk_alg(exsplus) ->
{#{type=>exsplus, max=>?UINT58MASK,
uniform=>fun exsplus_uniform/1, uniform_n=>fun exsplus_uniform/2},
fun exsplus_seed/1};
mk_alg(exs1024) ->
{#{type=>exs1024, max=>?UINT64MASK,
uniform=>fun exs1024_uniform/1, uniform_n=>fun exs1024_uniform/2},
fun exs1024_seed/1}.
%% =====================================================================
%% exs64 PRNG: Xorshift64*
%% Algorithm by Sebastiano Vigna
%% Reference URL: http://xorshift.di.unimi.it/
%% =====================================================================
-type exs64_state() :: uint64().
exs64_seed({A1, A2, A3}) ->
{V1, _} = exs64_next(((A1 band ?UINT32MASK) * 4294967197 + 1)),
{V2, _} = exs64_next(((A2 band ?UINT32MASK) * 4294967231 + 1)),
{V3, _} = exs64_next(((A3 band ?UINT32MASK) * 4294967279 + 1)),
((V1 * V2 * V3) rem (?UINT64MASK - 1)) + 1.
%% Advance xorshift64* state for one step and generate 64bit unsigned integer
-spec exs64_next(exs64_state()) -> {uint64(), exs64_state()}.
exs64_next(R) ->
R1 = R bxor (R bsr 12),
R2 = R1 bxor ((R1 band ?UINT39MASK) bsl 25),
R3 = R2 bxor (R2 bsr 27),
{(R3 * 2685821657736338717) band ?UINT64MASK, R3}.
exs64_uniform({Alg, R0}) ->
{V, R1} = exs64_next(R0),
{V / 18446744073709551616, {Alg, R1}}.
exs64_uniform(Max, {Alg, R}) ->
{V, R1} = exs64_next(R),
{(V rem Max) + 1, {Alg, R1}}.
%% =====================================================================
%% exsplus PRNG: Xorshift116+
%% Algorithm by Sebastiano Vigna
%% Reference URL: http://xorshift.di.unimi.it/
%% 58 bits fits into an immediate on 64bits erlang and is thus much faster.
%% Modification of the original Xorshift128+ algorithm to 116
%% by Sebastiano Vigna, a lot of thanks for his help and work.
%% =====================================================================
-type exsplus_state() :: [uint58()|uint58()].
exsplus_seed({A1, A2, A3}) ->
{_, R1} = exsplus_next([(((A1 * 4294967197) + 1) band ?UINT58MASK)|
(((A2 * 4294967231) + 1) band ?UINT58MASK)]),
{_, R2} = exsplus_next([(((A3 * 4294967279) + 1) band ?UINT58MASK)|
tl(R1)]),
R2.
%% Advance xorshift116+ state for one step and generate 58bit unsigned integer
-spec exsplus_next(exsplus_state()) -> {uint58(), exsplus_state()}.
exsplus_next([S1|S0]) ->
%% Note: members s0 and s1 are swapped here
S11 = (S1 bxor (S1 bsl 24)) band ?UINT58MASK,
S12 = S11 bxor S0 bxor (S11 bsr 11) bxor (S0 bsr 41),
{(S0 + S12) band ?UINT58MASK, [S0|S12]}.
exsplus_uniform({Alg, R0}) ->
{I, R1} = exsplus_next(R0),
{I / (?UINT58MASK+1), {Alg, R1}}.
exsplus_uniform(Max, {Alg, R}) ->
{V, R1} = exsplus_next(R),
{(V rem Max) + 1, {Alg, R1}}.
%% =====================================================================
%% exs1024 PRNG: Xorshift1024*
%% Algorithm by Sebastiano Vigna
%% Reference URL: http://xorshift.di.unimi.it/
%% =====================================================================
-type exs1024_state() :: {list(uint64()), list(uint64())}.
exs1024_seed({A1, A2, A3}) ->
B1 = (((A1 band ?UINT21MASK) + 1) * 2097131) band ?UINT21MASK,
B2 = (((A2 band ?UINT21MASK) + 1) * 2097133) band ?UINT21MASK,
B3 = (((A3 band ?UINT21MASK) + 1) * 2097143) band ?UINT21MASK,
{exs1024_gen1024((B1 bsl 43) bor (B2 bsl 22) bor (B3 bsl 1) bor 1),
[]}.
%% Generate a list of 16 64-bit element list
%% of the xorshift64* random sequence
%% from a given 64-bit seed.
%% Note: dependent on exs64_next/1
-spec exs1024_gen1024(uint64()) -> list(uint64()).
exs1024_gen1024(R) ->
exs1024_gen1024(16, R, []).
exs1024_gen1024(0, _, L) ->
L;
exs1024_gen1024(N, R, L) ->
{X, R2} = exs64_next(R),
exs1024_gen1024(N - 1, R2, [X|L]).
%% Calculation of xorshift1024*.
%% exs1024_calc(S0, S1) -> {X, NS1}.
%% X: random number output
-spec exs1024_calc(uint64(), uint64()) -> {uint64(), uint64()}.
exs1024_calc(S0, S1) ->
S11 = S1 bxor ((S1 band ?UINT33MASK) bsl 31),
S12 = S11 bxor (S11 bsr 11),
S01 = S0 bxor (S0 bsr 30),
NS1 = S01 bxor S12,
{(NS1 * 1181783497276652981) band ?UINT64MASK, NS1}.
%% Advance xorshift1024* state for one step and generate 64bit unsigned integer
-spec exs1024_next(exs1024_state()) -> {uint64(), exs1024_state()}.
exs1024_next({[S0,S1|L3], RL}) ->
{X, NS1} = exs1024_calc(S0, S1),
{X, {[NS1|L3], [S0|RL]}};
exs1024_next({[H], RL}) ->
NL = [H|lists:reverse(RL)],
exs1024_next({NL, []}).
exs1024_uniform({Alg, R0}) ->
{V, R1} = exs1024_next(R0),
{V / 18446744073709551616, {Alg, R1}}.
exs1024_uniform(Max, {Alg, R}) ->
{V, R1} = exs1024_next(R),
{(V rem Max) + 1, {Alg, R1}}.
