%% ``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 via the world wide web 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.
%%
-module(rand).
-author("Paulo Sergio Almeida ").
-export([start/0, start/1, stop/0]).
-export([uniform/0, uniform/1, uniform_vector/1]).
-export([permutation/1, weighted_permutation/1]).
-compile({inline, [{uniform_s,1}, {uniform_s,2}]}).
%%% Utilities for randomness.
%%% Random number generator server.
%%% Used to obtain a single random number sequence in an application.
%%% Provides uniform/0 and uniform/1 as the random module.
%%% Also provides functions to generate a random vector in a single
%%% invocation (minimizing msg passing), and to generate random
%%% permutations and weighted random permutations of lists.
%%%
%%% Uses same method as random module from stdlib:
%%% The method is attributed to B.A. Wichmann and I.D.Hill
%%% See "An efficient and portable pseudo-random number generator",
%%% Journal of Applied Statistics. AS183. 1982. Also Byte March 1987.
%%% interface functions
start() ->
Seed = {3172, 9814, 20125}, % random:seed0
start(Seed).
start(Seed) ->
Pid = spawn(fun() -> loop(Seed) end),
case catch register(rand, Pid) of
true -> ok;
_ -> exit(Pid, kill), ok
end.
stop() ->
rand ! {stop, self()}.
uniform() ->
rand ! {uniform, self()},
receive
{rand, Val} -> Val
end.
uniform(N) ->
rand ! {uniform, N, self()},
receive
{rand, Val} -> Val
end.
uniform_vector(N) ->
rand ! {uniform_vector, N, self()},
receive
{rand, Val} -> Val
end.
%% permutation(List) ->
%% Creates a random permutation of the list;
permutation(List) ->
Rands = uniform_vector(length(List)),
Perm = lists:sort(rzip(Rands, List, [])),
[X || {_W, X} <- Perm].
rzip([X | T1], [Y | T2], Acc) ->
rzip(T1, T2, [{X, Y} | Acc]);
rzip([], [], Acc) -> Acc.
%% weighted_permutation(List) ->
%% Receives a list of {Elem, Weight} and
%% creates a random permutation of the list according to weights;
%%
%% Based on "Weighted random sampling with a reservoir"
%% By P. S. Efraimidis and P. G. Spirakis
%% Information Processing Letters Volume 97, Issue 5 (March 2006)
%%
%% Elements with more weight appear with higher probability first in the
%% list.
weighted_permutation(List) ->
Rands = uniform_vector(length(List)),
Perm = lists:sort(weigh(List, Rands, [])),
[X || {_W, X} <- Perm].
weigh([{X, W} | T1], [R | T2], Acc) when W /= 0 ->
weigh(T1, T2, [{-math:pow(R, 1/W), X} | Acc]);
weigh([{X, _W} | T1], [_R | T2], Acc) ->
weigh(T1, T2, [{0.0, X} | Acc]);
weigh([], [], Acc) -> Acc.
%%% server loop
loop(Seed) ->
{Val, NewSeed} =
receive
{uniform, From} -> uniform_s(Seed);
{uniform, N, From} -> uniform_s(N, Seed);
{uniform_vector, N, From} -> univect(N, Seed, []);
{stop, From} -> exit(normal)
end,
From ! {rand, Val},
loop(NewSeed).
%%% implementation
univect(0, Seed, Acc) -> {Acc, Seed};
univect(N, Seed, Acc) ->
{Val, NewSeed} = uniform_s(Seed),
univect(N-1, NewSeed, [Val | Acc]).
%% same as random module in stdlib; copied here to avoid inter-module calls
%% and be able to inline
uniform_s({A1, A2, A3}) ->
B1 = (A1*171) rem 30269,
B2 = (A2*172) rem 30307,
B3 = (A3*170) rem 30323,
R = A1/30269 + A2/30307 + A3/30323,
{R - trunc(R), {B1,B2,B3}}.
uniform_s(N, State0) when N >= 1 ->
{F, State1} = uniform_s(State0),
{trunc(F * N) + 1, State1}.