diff options
Diffstat (limited to 'lib/stdlib/src/rand.erl')
| -rw-r--r-- | lib/stdlib/src/rand.erl | 220 |
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 %% ===================================================================== |