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