1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
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}}.
|
