stdlib: Add new random functionality/module

The old random module contains an old algorithm which have flaws
and the api requires the user to invoke seed and or checking if seed have
been invoked, if a non constant seed is to be used.

The api contains the following features:
- The user can invoke rand:unform/[0|1] directly and get a non constant seeding.
- The api is split in functional and non functional functions,
  i.e. usage of _s functions will not affect the process dictionary.
- The api contains several algorithms with different characteristics and
  can be extended with new algorithms in the future.
- Contains state of the art random number generators.
- Default algorithm is taylor made for erlang to be fast on 64bits machines.

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}}.