aboutsummaryrefslogtreecommitdiffstats
path: root/lib/stdlib/src/rand.erl
diff options
context:
space:
mode:
Diffstat (limited to 'lib/stdlib/src/rand.erl')
-rw-r--r--lib/stdlib/src/rand.erl302
1 files changed, 302 insertions, 0 deletions
diff --git a/lib/stdlib/src/rand.erl b/lib/stdlib/src/rand.erl
new file mode 100644
index 0000000000..0cafb35dd8
--- /dev/null
+++ b/lib/stdlib/src/rand.erl
@@ -0,0 +1,302 @@
+%%
+%% %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}}.