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.erl220
1 files changed, 178 insertions, 42 deletions
diff --git a/lib/stdlib/src/rand.erl b/lib/stdlib/src/rand.erl
index 93409d95df..dfd102f9ef 100644
--- a/lib/stdlib/src/rand.erl
+++ b/lib/stdlib/src/rand.erl
@@ -19,7 +19,7 @@
%%
%% =====================================================================
%% Multiple PRNG module for Erlang/OTP
-%% Copyright (c) 2015 Kenji Rikitake
+%% Copyright (c) 2015-2016 Kenji Rikitake
%% =====================================================================
-module(rand).
@@ -27,11 +27,14 @@
-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,
+ jump/0, jump/1,
normal/0, normal_s/1
]).
-compile({inline, [exs64_next/1, exsplus_next/1,
+ exsplus_jump/1,
exs1024_next/1, exs1024_calc/2,
+ exs1024_jump/1,
get_52/1, normal_kiwi/1]}).
-define(DEFAULT_ALG_HANDLER, exsplus).
@@ -42,19 +45,31 @@
%% =====================================================================
%% This depends on the algorithm handler function
--type alg_seed() :: exs64_state() | exsplus_state() | exs1024_state().
+-type alg_state() ::
+ exs64_state() | exsplus_state() | exs1024_state() | term().
+
%% This is the algorithm handler function within this module
--type alg_handler() :: #{type := alg(),
- max := integer(),
- next := fun(),
- uniform := fun(),
- uniform_n := fun()}.
-
-%% Internal state
--opaque state() :: {alg_handler(), alg_seed()}.
--type alg() :: exs64 | exsplus | exs1024.
--opaque export_state() :: {alg(), alg_seed()}.
--export_type([alg/0, state/0, export_state/0]).
+-type alg_handler() ::
+ #{type := alg(),
+ max := integer() | infinity,
+ next :=
+ fun((alg_state()) -> {non_neg_integer(), alg_state()}),
+ uniform :=
+ fun((state()) -> {float(), state()}),
+ uniform_n :=
+ fun((pos_integer(), state()) -> {pos_integer(), state()}),
+ jump :=
+ fun((state()) -> state())}.
+
+%% Algorithm state
+-type state() :: {alg_handler(), alg_state()}.
+-type builtin_alg() :: exs64 | exsplus | exs1024.
+-type alg() :: builtin_alg() | atom().
+-type export_state() :: {alg(), alg_state()}.
+-export_type(
+ [builtin_alg/0, alg/0, alg_handler/0, alg_state/0,
+ state/0, export_state/0]).
+-export_type([exs64_state/0, exsplus_state/0, exs1024_state/0]).
%% =====================================================================
%% API
@@ -68,7 +83,7 @@ export_seed() ->
_ -> undefined
end.
--spec export_seed_s(state()) -> export_state().
+-spec export_seed_s(State :: state()) -> export_state().
export_seed_s({#{type:=Alg}, Seed}) -> {Alg, Seed}.
%% seed(Alg) seeds RNG with runtime dependent values
@@ -77,31 +92,37 @@ export_seed_s({#{type:=Alg}, Seed}) -> {Alg, Seed}.
%% seed({Alg,Seed}) setup RNG with a previously exported seed
%% and return the NEW state
--spec seed(AlgOrExpState::alg() | export_state()) -> state().
+-spec seed(
+ AlgOrStateOrExpState :: builtin_alg() | state() | export_state()) ->
+ state().
seed(Alg) ->
- R = seed_s(Alg),
- _ = seed_put(R),
- R.
+ seed_put(seed_s(Alg)).
--spec seed_s(AlgOrExpState::alg() | export_state()) -> state().
-seed_s(Alg) when is_atom(Alg) ->
- seed_s(Alg, {erlang:phash2([{node(),self()}]),
- erlang:system_time(),
- erlang:unique_integer()});
+-spec seed_s(
+ AlgOrStateOrExpState :: builtin_alg() | state() | export_state()) ->
+ state().
+seed_s({AlgHandler, _Seed} = State) when is_map(AlgHandler) ->
+ State;
seed_s({Alg0, Seed}) ->
{Alg,_SeedFun} = mk_alg(Alg0),
- {Alg, Seed}.
+ {Alg, Seed};
+seed_s(Alg) ->
+ seed_s(Alg, {erlang:phash2([{node(),self()}]),
+ erlang:system_time(),
+ erlang:unique_integer()}).
%% seed/2: seeds RNG with the algorithm and given values
%% and returns the NEW state.
--spec seed(Alg :: alg(), {integer(), integer(), integer()}) -> state().
+-spec seed(
+ Alg :: builtin_alg(), Seed :: {integer(), integer(), integer()}) ->
+ state().
seed(Alg0, S0) ->
- State = seed_s(Alg0, S0),
- _ = seed_put(State),
- State.
+ seed_put(seed_s(Alg0, S0)).
--spec seed_s(Alg :: alg(), {integer(), integer(), integer()}) -> state().
+-spec seed_s(
+ Alg :: builtin_alg(), Seed :: {integer(), integer(), integer()}) ->
+ state().
seed_s(Alg0, S0 = {_, _, _}) ->
{Alg, Seed} = mk_alg(Alg0),
AS = Seed(S0),
@@ -113,7 +134,7 @@ seed_s(Alg0, S0 = {_, _, _}) ->
%% uniform/0: returns a random float X where 0.0 < X < 1.0,
%% updating the state in the process dictionary.
--spec uniform() -> X::float().
+-spec uniform() -> X :: float().
uniform() ->
{X, Seed} = uniform_s(seed_get()),
_ = seed_put(Seed),
@@ -123,7 +144,7 @@ uniform() ->
%% uniform/1 returns a random integer X where 1 =< X =< N,
%% updating the state in the process dictionary.
--spec uniform(N :: pos_integer()) -> X::pos_integer().
+-spec uniform(N :: pos_integer()) -> X :: pos_integer().
uniform(N) ->
{X, Seed} = uniform_s(N, seed_get()),
_ = seed_put(Seed),
@@ -133,7 +154,7 @@ uniform(N) ->
%% returns a random float X where 0.0 < X < 1.0,
%% and a new state.
--spec uniform_s(state()) -> {X::float(), NewS :: state()}.
+-spec uniform_s(State :: state()) -> {X :: float(), NewState :: state()}.
uniform_s(State = {#{uniform:=Uniform}, _}) ->
Uniform(State).
@@ -141,7 +162,8 @@ uniform_s(State = {#{uniform:=Uniform}, _}) ->
%% uniform_s/2 returns a random integer X where 1 =< X =< N,
%% and a new state.
--spec uniform_s(N::pos_integer(), state()) -> {X::pos_integer(), NewS::state()}.
+-spec uniform_s(N :: pos_integer(), State :: state()) ->
+ {X :: pos_integer(), NewState :: state()}.
uniform_s(N, State = {#{uniform_n:=Uniform, max:=Max}, _})
when 0 < N, N =< Max ->
Uniform(N, State);
@@ -150,6 +172,25 @@ uniform_s(N, State0 = {#{uniform:=Uniform}, _})
{F, State} = Uniform(State0),
{trunc(F * N) + 1, State}.
+%% jump/1: given a state, jump/1
+%% returns a new state which is equivalent to that
+%% after a large number of call defined for each algorithm.
+%% The large number is algorithm dependent.
+
+-spec jump(state()) -> NewState :: state().
+jump(State = {#{jump:=Jump}, _}) ->
+ Jump(State).
+
+%% jump/0: read the internal state and
+%% apply the jump function for the state as in jump/1
+%% and write back the new value to the internal state,
+%% then returns the new value.
+
+-spec jump() -> NewState :: state().
+
+jump() ->
+ seed_put(jump(seed_get())).
+
%% normal/0: returns a random float with standard normal distribution
%% updating the state in the process dictionary.
@@ -163,7 +204,7 @@ normal() ->
%% The Ziggurat Method for generating random variables - Marsaglia and Tsang
%% Paper and reference code: http://www.jstatsoft.org/v05/i08/
--spec normal_s(state()) -> {float(), NewS :: state()}.
+-spec normal_s(State :: state()) -> {float(), NewState :: state()}.
normal_s(State0) ->
{Sign, R, State} = get_52(State0),
Idx = R band 16#FF,
@@ -192,9 +233,10 @@ normal_s(State0) ->
-type uint64() :: 0..16#ffffffffffffffff.
-type uint58() :: 0..16#03ffffffffffffff.
--spec seed_put(state()) -> undefined | state().
+-spec seed_put(state()) -> state().
seed_put(Seed) ->
- put(?SEED_DICT, Seed).
+ put(?SEED_DICT, Seed),
+ Seed.
seed_get() ->
case get(?SEED_DICT) of
@@ -205,15 +247,18 @@ seed_get() ->
%% Setup alg record
mk_alg(exs64) ->
{#{type=>exs64, max=>?UINT64MASK, next=>fun exs64_next/1,
- uniform=>fun exs64_uniform/1, uniform_n=>fun exs64_uniform/2},
+ uniform=>fun exs64_uniform/1, uniform_n=>fun exs64_uniform/2,
+ jump=>fun exs64_jump/1},
fun exs64_seed/1};
mk_alg(exsplus) ->
{#{type=>exsplus, max=>?UINT58MASK, next=>fun exsplus_next/1,
- uniform=>fun exsplus_uniform/1, uniform_n=>fun exsplus_uniform/2},
+ uniform=>fun exsplus_uniform/1, uniform_n=>fun exsplus_uniform/2,
+ jump=>fun exsplus_jump/1},
fun exsplus_seed/1};
mk_alg(exs1024) ->
{#{type=>exs1024, max=>?UINT64MASK, next=>fun exs1024_next/1,
- uniform=>fun exs1024_uniform/1, uniform_n=>fun exs1024_uniform/2},
+ uniform=>fun exs1024_uniform/1, uniform_n=>fun exs1024_uniform/2,
+ jump=>fun exs1024_jump/1},
fun exs1024_seed/1}.
%% =====================================================================
@@ -222,7 +267,7 @@ mk_alg(exs1024) ->
%% Reference URL: http://xorshift.di.unimi.it/
%% =====================================================================
--type exs64_state() :: uint64().
+-opaque exs64_state() :: uint64().
exs64_seed({A1, A2, A3}) ->
{V1, _} = exs64_next(((A1 band ?UINT32MASK) * 4294967197 + 1)),
@@ -246,6 +291,9 @@ exs64_uniform(Max, {Alg, R}) ->
{V, R1} = exs64_next(R),
{(V rem Max) + 1, {Alg, R1}}.
+exs64_jump(_) ->
+ erlang:error(not_implemented).
+
%% =====================================================================
%% exsplus PRNG: Xorshift116+
%% Algorithm by Sebastiano Vigna
@@ -254,7 +302,7 @@ exs64_uniform(Max, {Alg, R}) ->
%% Modification of the original Xorshift128+ algorithm to 116
%% by Sebastiano Vigna, a lot of thanks for his help and work.
%% =====================================================================
--type exsplus_state() :: nonempty_improper_list(uint58(), uint58()).
+-opaque exsplus_state() :: nonempty_improper_list(uint58(), uint58()).
-dialyzer({no_improper_lists, exsplus_seed/1}).
@@ -283,13 +331,47 @@ exsplus_uniform(Max, {Alg, R}) ->
{V, R1} = exsplus_next(R),
{(V rem Max) + 1, {Alg, R1}}.
+%% This is the jump function for the exsplus generator, equivalent
+%% to 2^64 calls to next/1; it can be used to generate 2^52
+%% non-overlapping subsequences for parallel computations.
+%% Note: the jump function takes 116 times of the execution time of
+%% next/1.
+
+%% -define(JUMPCONST, 16#000d174a83e17de2302f8ea6bc32c797).
+%% split into 58-bit chunks
+%% and two iterative executions
+
+-define(JUMPCONST1, 16#02f8ea6bc32c797).
+-define(JUMPCONST2, 16#345d2a0f85f788c).
+-define(JUMPELEMLEN, 58).
+
+-dialyzer({no_improper_lists, exsplus_jump/1}).
+-spec exsplus_jump(state()) -> state().
+exsplus_jump({Alg, S}) ->
+ {S1, AS1} = exsplus_jump(S, [0|0], ?JUMPCONST1, ?JUMPELEMLEN),
+ {_, AS2} = exsplus_jump(S1, AS1, ?JUMPCONST2, ?JUMPELEMLEN),
+ {Alg, AS2}.
+
+-dialyzer({no_improper_lists, exsplus_jump/4}).
+exsplus_jump(S, AS, _, 0) ->
+ {S, AS};
+exsplus_jump(S, [AS0|AS1], J, N) ->
+ {_, NS} = exsplus_next(S),
+ case (J band 1) of
+ 1 ->
+ [S0|S1] = S,
+ exsplus_jump(NS, [(AS0 bxor S0)|(AS1 bxor S1)], J bsr 1, N-1);
+ 0 ->
+ exsplus_jump(NS, [AS0|AS1], J bsr 1, N-1)
+ end.
+
%% =====================================================================
%% exs1024 PRNG: Xorshift1024*
%% Algorithm by Sebastiano Vigna
%% Reference URL: http://xorshift.di.unimi.it/
%% =====================================================================
--type exs1024_state() :: {list(uint64()), list(uint64())}.
+-opaque exs1024_state() :: {list(uint64()), list(uint64())}.
exs1024_seed({A1, A2, A3}) ->
B1 = (((A1 band ?UINT21MASK) + 1) * 2097131) band ?UINT21MASK,
@@ -340,6 +422,60 @@ exs1024_uniform(Max, {Alg, R}) ->
{V, R1} = exs1024_next(R),
{(V rem Max) + 1, {Alg, R1}}.
+%% This is the jump function for the exs1024 generator, equivalent
+%% to 2^512 calls to next(); it can be used to generate 2^512
+%% non-overlapping subsequences for parallel computations.
+%% Note: the jump function takes ~2000 times of the execution time of
+%% next/1.
+
+%% Jump constant here split into 58 bits for speed
+-define(JUMPCONSTHEAD, 16#00242f96eca9c41d).
+-define(JUMPCONSTTAIL,
+ [16#0196e1ddbe5a1561,
+ 16#0239f070b5837a3c,
+ 16#03f393cc68796cd2,
+ 16#0248316f404489af,
+ 16#039a30088bffbac2,
+ 16#02fea70dc2d9891f,
+ 16#032ae0d9644caec4,
+ 16#0313aac17d8efa43,
+ 16#02f132e055642626,
+ 16#01ee975283d71c93,
+ 16#00552321b06f5501,
+ 16#00c41d10a1e6a569,
+ 16#019158ecf8aa1e44,
+ 16#004e9fc949d0b5fc,
+ 16#0363da172811fdda,
+ 16#030e38c3b99181f2,
+ 16#0000000a118038fc]).
+-define(JUMPTOTALLEN, 1024).
+-define(RINGLEN, 16).
+
+-spec exs1024_jump(state()) -> state().
+
+exs1024_jump({Alg, {L, RL}}) ->
+ P = length(RL),
+ AS = exs1024_jump({L, RL},
+ [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+ ?JUMPCONSTTAIL, ?JUMPCONSTHEAD, ?JUMPELEMLEN, ?JUMPTOTALLEN),
+ {ASL, ASR} = lists:split(?RINGLEN - P, AS),
+ {Alg, {ASL, lists:reverse(ASR)}}.
+
+exs1024_jump(_, AS, _, _, _, 0) ->
+ AS;
+exs1024_jump(S, AS, [H|T], _, 0, TN) ->
+ exs1024_jump(S, AS, T, H, ?JUMPELEMLEN, TN);
+exs1024_jump({L, RL}, AS, JL, J, N, TN) ->
+ {_, NS} = exs1024_next({L, RL}),
+ case (J band 1) of
+ 1 ->
+ AS2 = lists:zipwith(fun(X, Y) -> X bxor Y end,
+ AS, L ++ lists:reverse(RL)),
+ exs1024_jump(NS, AS2, JL, J bsr 1, N-1, TN-1);
+ 0 ->
+ exs1024_jump(NS, AS, JL, J bsr 1, N-1, TN-1)
+ end.
+
%% =====================================================================
%% Ziggurat cont
%% =====================================================================
generated by cgit v1.2.3 (git 2.25.1) at 2025-01-01 16:03:55 +0000