File 2951-Implement-the-custom-made-mwc59-PRNG.patch of Package erlang

Overview Repositories Revisions Requests Users Attributes Meta

File 2951-Implement-the-custom-made-mwc59-PRNG.patch of Package erlang (Revision 0b3a6d6e89afeaba4dcdf8fb3d620103)

Currently displaying revision 0b3a6d6e89afeaba4dcdf8fb3d620103 , Show latest

From 46dd97f5f62f3d01461a9a63770ef80149e0b745 Mon Sep 17 00:00:00 2001
From: Raimo Niskanen <raimo@erlang.org>
Date: Thu, 7 Apr 2022 16:22:32 +0200
Subject: [PATCH 1/6] Implement the custom made `mwc59` PRNG

The `mwc59`, Multiply With Carry 59-bit, PRNG uses 32-bit "digits"
and a multiplier to keep intermediate results within 59 bit unsigned.
On that two different scramblers can be applied to cover the statistical
flaws of a power of 2 multiplier generator.

Polish the test suite.
---
 lib/stdlib/doc/src/rand.xml    | 174 ++++++++++++++++++++
 lib/stdlib/src/rand.erl        | 172 ++++++++++++++++++-
 lib/stdlib/test/rand_SUITE.erl | 292 +++++++++++++++++++++++++++++----
 3 files changed, 595 insertions(+), 43 deletions(-)

diff --git a/lib/stdlib/doc/src/rand.xml b/lib/stdlib/doc/src/rand.xml
index 8bd6e56364..884251ab02 100644
--- a/lib/stdlib/doc/src/rand.xml
+++ b/lib/stdlib/doc/src/rand.xml
@@ -416,6 +416,10 @@ tests. We suggest to use a sign test to extract a random Boolean value.</pre>
 <name name="lcg35_state"/>
 <desc>0 .. (2^35 - 1)</desc>
 </datatype>
+ <datatype>
+ <name name="mwc59_state"/>
+ <desc>1 .. ((16#1ffb072 * 2^29 - 1) - 1)</desc>
+ </datatype>
 </datatypes>
 
 
@@ -940,6 +944,23 @@ end.</pre>
 </note>
 </desc>
 </func>
+ <func>
+ <name name="mcg35_seed" arity="0" since="OTP 25.0"/>
+ <name name="mcg35_seed" arity="1" since="OTP 25.0"/>
+ <fsummary>Create a generator state.</fsummary>
+ <desc>
+ 
+ Returns a generator state <c><anno>X</anno></c>.
+ It is set it to <c><anno>S</anno></c>, after folding it
+ into the state's range, if out of range.
+ 
+ 
+ Without <c><anno>S</anno></c>,
+ the generator state is created as for
+ <seemfa marker="#seed_s/1"><c>seed_s(atom())</c></seemfa>.
+ 
+ </desc>
+ </func>
 <func>
 <name name="lcg35" arity="1" since="OTP 25.0"/>
 <fsummary>Return a new generator state / number.</fsummary>
@@ -998,5 +1019,158 @@ end.</pre>
 </note>
 </desc>
 </func>
+ <func>
+ <name name="lcg35_seed" arity="0" since="OTP 25.0"/>
+ <name name="lcg35_seed" arity="1" since="OTP 25.0"/>
+ <fsummary>Create a generator state.</fsummary>
+ <desc>
+ 
+ Returns a generator state <c><anno>X</anno></c>.
+ It is set it to <c><anno>S</anno></c>, after folding it
+ into the state's range, if out of range.
+ 
+ 
+ Without <c><anno>S</anno></c>,
+ the generator state is created as for
+ <seemfa marker="#seed_s/1"><c>seed_s(atom())</c></seemfa>.
+ 
+ </desc>
+ </func>
+ <func>
+ <name name="mwc59" arity="1" since="OTP 25.0"/>
+ <fsummary>Return a new generator state.</fsummary>
+ <desc>
+ 
+ Returns a new generator state <c><anno>CX1</anno></c>,
+ according to a Multiply With Carry generator,
+ which is an efficient implementation of a
+ Multiplicative Congruential Generator with
+ a power of 2 multiplier and a prime modulus.
+ 
+ 
+ This generator uses the multiplier 2^32 and the modulus
+ 16#7f17555&nbsp;*&nbsp;2^32&nbsp;-&nbsp;1,
+ which have been selected,
+ in collaboration with Sebastiano Vigna,
+ to avoid bignum operations
+ and still get good statistical quality.
+ It can be written as: 
+ <c>C&nbsp;=&nbsp;<anno>CX0</anno>&nbsp;bsr&nbsp;32</c> 
+ <c>X&nbsp;=&nbsp;<anno>CX0</anno>&nbsp;band&nbsp;((1&nbsp;bsl&nbsp;32)-1))</c> 
+ <c><anno>CX1</anno>&nbsp;=&nbsp;16#7f17555&nbsp;*&nbsp;X&nbsp;+&nbsp;C</c>
+ 
+ 
+ The quality of the generated bits is best in the low end;
+ the low 16 bits of the generator state
+ <c><anno>CX1</anno></c> can be used without any scrambling.
+ <seemfa marker="#mwc59_value/1"><c>mwc59_value</c></seemfa>
+ is a fast scrambler that gives 32 high quality bits.
+ With the slightly slower
+ <seemfa marker="#mwc59_full_value/1">
+ <c>mwc59_full_value</c>
+ </seemfa>
+ all 58 returned bits are of good quality, and
+ <seemfa marker="#mwc59_float/1"><c>mwc59_float</c></seemfa>
+ returns a good quality <c>float()</c>.
+ 
+ 
+ On a typical 64 bit Erlang VM this generator executes
+ in below 6% (1/15) of the time
+ for the default algorithm in this module.
+ Adding the
+ <seemfa marker="#mwc59_value/1"><c>mwc59_value</c></seemfa>
+ scrambler the total time becomes 10% (1/10),
+ and with
+ <seemfa marker="#mwc59_full_value/1">
+ <c>mwc59_full_value</c>
+ </seemfa>
+ it becomes 12% (1/8) of the time for the default algorithm.
+ With
+ <seemfa marker="#mwc59_float/1"><c>mwc59_float</c></seemfa>
+ the total time is 40% of the time for the default
+ algorithm generating a <c>float()</c>.
+ 
+ <note>
+ 
+ This generator is a niche generator for high speed
+ applications. It has a much shorter period
+ than the default generator, which in itself
+ is a quality concern, although when used with the
+ value scramblers it passes strict PRNG tests
+ (TestU01 BigCrush and PractRand 2 TB).
+ The generator is significantly faster than
+ <seemfa marker="#exsp_next/1"><c>exsp_next/1</c></seemfa>
+ but with a bit lower quality.
+ 
+ </note>
+ </desc>
+ </func>
+ <func>
+ <name name="mwc59_value" arity="1" since="OTP 25.0"/>
+ <fsummary>Return the generator value.</fsummary>
+ <desc>
+ 
+ Returns a 58-bit value <c><anno>V</anno></c>
+ from a generator state <c><anno>CX</anno></c>.
+ The generator state is scrambled using
+ a 10-bit xorshift which masks
+ the statistical imperfecions of the base generator
+ <seemfa marker="#mwc59/1"><c>mwc59</c></seemfa>
+ enough that the 32 lowest bits are of good quality.
+ It is not recommended to generate numbers
+ in a range > 2^32 with this function.
+ 
+ </desc>
+ </func>
+ <func>
+ <name name="mwc59_full_value" arity="1" since="OTP 25.0"/>
+ <fsummary>Return the generator value.</fsummary>
+ <desc>
+ 
+ Returns a 58-bit value <c><anno>V</anno></c>
+ from a generator state <c><anno>CX</anno></c>.
+ The generator state is scrambled using
+ an 8-bit followed by a 16-bit xorshift which masks
+ the statistical imperfecions of the base generator
+ <seemfa marker="#mwc59/1"><c>mwc59</c></seemfa>
+ enough that all 58 bits are of good quality.
+ Beware of bias in the generated numbers if
+ generating on a non power of 2 range too close to 2^58.
+ 
+ </desc>
+ </func>
+ <func>
+ <name name="mwc59_float" arity="1" since="OTP 25.0"/>
+ <fsummary>Return a generated float.</fsummary>
+ <desc>
+ 
+ Returns the generator value <c><anno>V</anno></c>
+ from a generator state <c><anno>CX</anno></c>,
+ as a <c>float()</c>.
+ The generator state is scrambled as with
+ <seemfa marker="#mwc59_value/1">
+ <c>mwc59_full_value/1</c>
+ </seemfa>
+ before converted to a <c>float()</c>.
+ 
+ </desc>
+ </func>
+ <func>
+ <name name="mwc59_seed" arity="0" since="OTP 25.0"/>
+ <name name="mwc59_seed" arity="1" since="OTP 25.0"/>
+ <fsummary>Create a generator state.</fsummary>
+ <desc>
+ 
+ Returns a generator state <c><anno>CX</anno></c>.
+ It is set it to <c><anno>S</anno></c>, after folding it
+ into the state's range, if out of range.
+ 
+ 
+ Without <c><anno>S</anno></c>,
+ the generator state is created as for
+ <seemfa marker="#seed_s/1"><c>seed_s(atom())</c></seemfa>.
+ 
+ </desc>
+ </func>
 </funcs>
 </erlref>
diff --git a/lib/stdlib/src/rand.erl b/lib/stdlib/src/rand.erl
index f92d4212b4..aa228d6eaa 100644
--- a/lib/stdlib/src/rand.erl
+++ b/lib/stdlib/src/rand.erl
@@ -37,7 +37,11 @@
 	]).
 
 %% Utilities
--export([exsp_next/1, exsp_jump/1, splitmix64_next/1, mcg35/1, lcg35/1]).
+-export([exsp_next/1, exsp_jump/1, splitmix64_next/1,
+ mcg35/1, mcg35_seed/0, mcg35_seed/1,
+ lcg35/1, lcg35_seed/0, lcg35_seed/1,
+ mwc59/1, mwc59_value/1, mwc59_full_value/1, mwc59_float/1,
+ mwc59_seed/0, mwc59_seed/1]).
 
 %% Test, dev and internal
 -export([exro928_jump_2pow512/1, exro928_jump_2pow20/1,
@@ -145,7 +149,7 @@
 state/0, export_state/0, seed/0]).
 -export_type(
 [exsplus_state/0, exro928_state/0, exrop_state/0, exs1024_state/0,
- exs64_state/0, dummy_state/0]).
+ exs64_state/0, mwc59_state/0, dummy_state/0]).
 -export_type(
 [uint58/0, uint64/0, splitmix64_state/0, mcg35_state/0, lcg35_state/0]).
 
@@ -272,9 +276,12 @@ seed_s({Alg, AlgState}) when is_atom(Alg) ->
 {AlgHandler,_SeedFun} = mk_alg(Alg),
 {AlgHandler,AlgState};
 seed_s(Alg) ->
- seed_s(Alg, {erlang:phash2([{node(),self()}]),
-		 erlang:system_time(),
-		 erlang:unique_integer()}).
+ seed_s(Alg, default_seed()).
+
+default_seed() ->
+ {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.
@@ -866,7 +873,6 @@ exsss_seed({A1, A2, A3}) ->
 ?MASK(58, V_b + ?BSL(58, V_b, 3)) % * 9
 end).
 
-
 %% Advance state and generate 58bit unsigned integer
 %%
 -dialyzer({no_improper_lists, exsp_next/1}).
@@ -1425,6 +1431,8 @@ dummy_uniform({AlgHandler,R}) ->
 %%
 dummy_seed(L) when is_list(L) ->
 case L of
+ [] ->
+ erlang:error(zero_seed);
 [X] when is_integer(X) ->
 ?MASK(58, X);
 [X|_] when is_integer(X) ->
@@ -1487,6 +1495,14 @@ mcg35(X0) -> % when is_integer(X0), 1 =< X0, X0 < ?MCG35_M ->
 true -> X1 - ?MCG35_M
 end.
 
+-spec mcg35_seed() -> X :: mcg35_state().
+mcg35_seed() ->
+ {A1, A2, A3} = default_seed(),
+ mcg35_seed(seed64_int([A1, A2, A3])).
+
+-spec mcg35_seed(S :: integer()) -> X :: mcg35_state().
+mcg35_seed(S) when is_integer(S) ->
+ mod(?MCG35_M - 1, S) + 1.
 
 %% =====================================================================
 %% lcg35 PRNG: Multiplicative Linear Congruential Generator
@@ -1519,8 +1535,7 @@ mcg35(X0) -> % when is_integer(X0), 1 =< X0, X0 < ?MCG35_M ->
 -type lcg35_state() :: 0..?MASK(?LCG35_B).
 
 -spec lcg35(X0 :: lcg35_state()) -> X1 :: lcg35_state().
-%%lcg35(X0) when is_integer(X0), 0 =< X0, X0 =< ?MASK(?LCG35_B) ->
-lcg35(X0) ->
+lcg35(X0) -> % when is_integer(X0), 0 =< X0, X0 =< ?MASK(?LCG35_B) ->
 %% The mask operation on the input tricks the JIT into
 %% realizing that all following operations does not
 %% need bignum handling. The suggested guard test above
@@ -1528,6 +1543,123 @@ lcg35(X0) ->
 %% needs much more native code to execute than a 'band'.
 ?MASK(?LCG35_B, ?LCG35_A * ?MASK(?LCG35_B, X0) + ?LCG35_C).
 
+-spec lcg35_seed() -> X :: lcg35_state().
+lcg35_seed() ->
+ {A1, A2, A3} = default_seed(),
+ lcg35_seed(seed64_int([A1, A2, A3])).
+
+-spec lcg35_seed(S :: integer()) -> X :: lcg35_state().
+lcg35_seed(S) when is_integer(S) ->
+ ?MASK(?LCG35_B, S).
+
+%% =====================================================================
+%% mcg58 PRNG: Multiply With Carry generator
+%%
+%% Parameters deduced in collaboration with
+%% Prof. Sebastiano Vigna of the University of Milano.
+%%
+%% X = CX0 & (2^B - 1) % Low B bits - digit
+%% C = CX0 >> B % High bits - carry
+%% CX1 = A * X0 + C0
+%%
+%% An MWC generator is an efficient alternative implementation of
+%% a Multiplicative Congruential Generator, that is, the generator
+%% CX1 = (CX0 * 2^B) rem P
+%% where P is the safe prime (A * 2^B - 1), that generates
+%% the same sequence in the reverse order. The generator
+%% CX1 = (A * CX0) rem P
+%% that uses the multiplicative inverse mod P is, indeed,
+%% an exact equevalent to the corresponding MWC generator.
+%%
+%% An MWC generator has, due to the power of two multiplier
+%% in the corresponding MCG, got known statistical weaknesses
+%% in the spectral score for 3 dimensions, so it should be used
+%% with a scrambler that hides the flaws. The scramblers
+%% have been tried out in the PractRand and TestU01 frameworks
+%% and settled for a single Xorshift to get B good bits,
+%% and a double Xorshift to get all bits good enough.
+%%
+%% The chosen parameters are:
+%% A = 16#7f17555
+%% B = 32
+%% Single Xorshift: 10
+%% Double Xorshift: 8, 16
+%%
+%% These parameters gives the MWC "digit" size 32 bits
+%% which gives them theoretical statistical guarantees,
+%% and keeps the state in 59 bits.
+%%
+%% The state should only be used to mask or rem out low bits.
+%% The scramblers return 58 bits from which a number should
+%% be masked or rem:ed out.
+%%
+%% =====================================================================
+%%% -define(MWC_A, (6)).
+%%% -define(MWC_B, (3)).
+
+%%% -define(MWC59_A, (16#20075dc0)).
+%%% -define(MWC59_A, (16#1ffb0729)).
+%%% -define(MWC59_B, (29)).
+
+%%% -define(MWC59_A, (16#7fa6502)).
+-define(MWC59_A, (16#7f17555)).
+%%% -define(MWC59_A, (16#3f35301)).
+-define(MWC59_B, (32)).
+-define(MWC59_P, ((?MWC59_A bsl ?MWC59_B) - 1)).
+
+-define(MWC59_XS, 10).
+-define(MWC59_XS1, 8).
+-define(MWC59_XS2, 16).
+
+-type mwc59_state() :: 1..?MWC59_P-1.
+
+-spec mwc59(CX0 :: mwc59_state()) -> CX1 :: mwc59_state().
+mwc59(CX0) -> % when is_integer(CX0), 1 =< CX0, CX0 < ?MWC59_P ->
+ CX = ?MASK(58, CX0),
+ C = CX bsr ?MWC59_B,
+ X = ?MASK(?MWC59_B, CX),
+ ?MWC59_A * X + C.
+
+-spec mwc59_value(CX :: mwc59_state()) -> V :: 0..?MASK(58).
+mwc59_value(CX1) -> % when is_integer(CX0), 1 =< CX0, CX0 < ?MWC59_P ->
+ CX = ?MASK(58, CX1),
+ CX bxor ?BSL(58, CX, ?MWC59_XS).
+
+-spec mwc59_full_value(CX :: mwc59_state()) -> V :: 0..?MASK(58).
+mwc59_full_value(CX1) -> % when is_integer(CX0), 1 =< CX0, CX0 < ?MWC59_P ->
+ CX = ?MASK(58, CX1),
+ CX2 = CX bxor ?BSL(58, CX, ?MWC59_XS1),
+ CX2 bxor ?BSL(58, CX2, ?MWC59_XS2).
+
+-spec mwc59_float(CX :: mwc59_state()) -> V :: float().
+mwc59_float(CX1) ->
+ CX = ?MASK(53, CX1),
+ CX2 = CX bxor ?BSL(53, CX, ?MWC59_XS1),
+ (CX2 bxor ?BSL(53, CX2, ?MWC59_XS2)) * ?TWO_POW_MINUS53.
+
+-spec mwc59_seed() -> CX :: mwc59_state().
+mwc59_seed() ->
+ {A1, A2, A3} = default_seed(),
+ mwc59_seed(seed64_int([A1, A2, A3])).
+
+-spec mwc59_seed(S :: integer()) -> CX :: mwc59_state().
+mwc59_seed(S) when is_integer(S) ->
+ mod(?MWC59_P - 1, S) + 1.
+
+
+
+%%% %% Verification by equivalent MCG generator
+%%% mwc59_r(CX1) ->
+%%% (CX1 bsl ?MWC59_B) rem ?MWC59_P. % Reverse
+%%% %%% (CX1 * ?MWC59_A) rem ?MWC59_P. % Forward
+%%%
+%%% mwc59(CX0, 0) ->
+%%% CX0;
+%%% mwc59(CX0, N) ->
+%%% CX1 = mwc59(CX0),
+%%% CX0 = mwc59_r(CX1),
+%%% mwc59(CX1, N - 1).
+
 
 %% =====================================================================
 %% Mask and fill state list, ensure not all zeros
@@ -1592,6 +1724,16 @@ seed64(X_0) ->
 	 ZX
 end.
 
+%% Create a splitmixed (big) integer from a list of integers
+seed64_int(As) ->
+ seed64_int(As, 0, 0).
+%%
+seed64_int([], _X, Y) ->
+ Y;
+seed64_int([A | As], X0, Y) ->
+ {Z, X1} = splitmix64_next(A bxor X0),
+ seed64_int(As, X1, (Y bsl 64) bor Z).
+
 %% =====================================================================
 %% The SplitMix64 generator:
 %%
@@ -1955,7 +2097,19 @@ bc64(V) -> ?BC(V, 64).
 %% Linear from high bit - higher probability first gives faster execution
 bc(V, B, N) when B =< V -> N;
 bc(V, B, N) -> bc(V, B bsr 1, N - 1).
- 
+
+
+%% Non-negative rem
+mod(Q, X) ->
+ Y = X rem Q,
+ if
+ Y < 0 ->
+ Y + Q;
+ true ->
+ Y
+ end.
+
+
 make_float(S, E, M) ->
 <<F/float>> = <<S:1, E:11, M:52>>,
 F.
diff --git a/lib/stdlib/test/rand_SUITE.erl b/lib/stdlib/test/rand_SUITE.erl
index 7ce68b349d..93d7b2914a 100644
--- a/lib/stdlib/test/rand_SUITE.erl
+++ b/lib/stdlib/test/rand_SUITE.erl
@@ -35,7 +35,7 @@ all() ->
 [seed, interval_int, interval_float,
 bytes_count,
 api_eq,
- mcg35_api, mcg35_rem, lcg35_api,
+ mcg35_api, mcg35_rem, lcg35_api, mwc59_api,
 exsp_next_api, exsp_jump_api, splitmix64_next_api,
 reference,
 {group, basic_stats},
@@ -254,6 +254,39 @@ lcg35_api(X0, N)
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+%% Verify mwc59 behaviour
+%%
+mwc59_api(Config) when is_list(Config) ->
+ mwc59_api(1, 1000000).
+
+mwc59_api(CX0, 0) ->
+ CX = 187860517065527182,
+ {CX, CX} = {CX0, CX},
+ V0 = rand:mwc59_value(CX0),
+ V = 230807595801982862,
+ {V, V} = {V0, V},
+ W0 = rand:mwc59_full_value(CX0),
+ W = 202476383090409870,
+ {W, W} = {W0, W},
+ F0 = rand:mwc59_float(CX0),
+ F = (W band ((1 bsl 53) - 1)) * (1 / (1 bsl 53)),
+ {F, F} = {F0, F},
+ ok;
+mwc59_api(CX, N)
+ when is_integer(CX), 1 =< CX, CX < (16#7f17555 bsl 32) - 1 ->
+ V = rand:mwc59_value(CX),
+ W = rand:mwc59_full_value(CX),
+ F = rand:mwc59_float(CX),
+ true = 0 =< V,
+ true = V < 1 bsl 58,
+ true = 0 =< W,
+ true = W < 1 bsl 58,
+ true = 0.0 =< F,
+ true = F < 1.0,
+ mwc59_api(rand:mwc59(CX), N - 1).
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 %% Verify exsp_next behaviour
 %%
 exsp_next_api(Config) when is_list(Config) ->
@@ -1132,7 +1165,7 @@ do_measure(Iterations) ->
 St1
 end
 end
- end, mcg35_inline, Iterations,
+ end, {mcg35,mod}, Iterations,
 TMarkUniformRange10000, OverheadUniformRange1000),
 _ =
 measure_1(
@@ -1147,7 +1180,52 @@ do_measure(Iterations) ->
 end
 end
 end,
- lcg35_inline, Iterations,
+ {lcg35,mod}, Iterations,
+ TMarkUniformRange10000, OverheadUniformRange1000),
+ _ =
+ measure_1(
+ fun (_Mod, _State) ->
+ Range = 10000,
+ fun (St0) ->
+ St1 = rand:mwc59(St0),
+ %% Just a 'rem' with slightly skewed distribution
+ case (St1 rem Range) + 1 of
+ R when is_integer(R), 0 < R, R =< Range ->
+ St1
+ end
+ end
+ end,
+ {mwc59,raw_mod}, Iterations,
+ TMarkUniformRange10000, OverheadUniformRange1000),
+ _ =
+ measure_1(
+ fun (_Mod, _State) ->
+ Range = 10000,
+ fun (St0) ->
+ St1 = rand:mwc59(St0),
+ %% Just a 'rem' with slightly skewed distribution
+ case (rand:mwc59_value(St1) rem Range) + 1 of
+ R when is_integer(R), 0 < R, R =< Range ->
+ St1
+ end
+ end
+ end,
+ {mwc59,value_mod}, Iterations,
+ TMarkUniformRange10000, OverheadUniformRange1000),
+ _ =
+ measure_1(
+ fun (_Mod, _State) ->
+ Range = 10000,
+ fun (St0) ->
+ St1 = rand:mwc59(St0),
+ %% Just a 'rem' with slightly skewed distribution
+ case (rand:mwc59_full_value(St1) rem Range) + 1 of
+ R when is_integer(R), 0 < R, R =< Range ->
+ St1
+ end
+ end
+ end,
+ {mwc59,full_mod}, Iterations,
 TMarkUniformRange10000, OverheadUniformRange1000),
 _ =
 measure_1(
@@ -1162,7 +1240,7 @@ do_measure(Iterations) ->
 end
 end
 end,
- exsp_inline, Iterations,
+ {exsp,next}, Iterations,
 TMarkUniformRange10000, OverheadUniformRange1000),
 _ =
 measure_1(
@@ -1220,7 +1298,52 @@ do_measure(Iterations) ->
 end
 end
 end,
- lcg35_inline, Iterations,
+ {lcg35,shift}, Iterations,
+ TMarkUniform32Bit, OverheadUniform32Bit),
+ _ =
+ measure_1(
+ fun (_Mod, _State) ->
+ Range = 1 bsl 32,
+ fun (St0) ->
+ St1 = rand:mwc59(St0),
+ case St1 band ((1 bsl 32)-1) of
+ R when is_integer(R), 0 =< R, R < Range ->
+ St1
+ end
+ end
+ end,
+ {mwc59,raw_mask}, Iterations,
+ TMarkUniform32Bit, OverheadUniform32Bit),
+ _ =
+ measure_1(
+ fun (_Mod, _State) ->
+ Range = 1 bsl 32,
+ fun (St0) ->
+ St1 = rand:mwc59(St0),
+ case rand:mwc59_value(St1) band ((1 bsl 32)-1) of
+ R when is_integer(R), 0 =< R, R < Range ->
+ St1
+ end
+ end
+ end,
+ {mwc59,value_mask}, Iterations,
+ TMarkUniform32Bit, OverheadUniform32Bit),
+ _ =
+ measure_1(
+ fun (_Mod, _State) ->
+ Range = 1 bsl 32,
+ fun (St0) ->
+ St1 = rand:mwc59(St0),
+ case
+ rand:mwc59_full_value(St1) band
+ ((1 bsl 32)-1)
+ of
+ R when is_integer(R), 0 =< R, R < Range ->
+ St1
+ end
+ end
+ end,
+ {mwc59,full_mask}, Iterations,
 TMarkUniform32Bit, OverheadUniform32Bit),
 _ =
 measure_1(
@@ -1234,7 +1357,7 @@ do_measure(Iterations) ->
 end
 end
 end,
- exsp_inline, Iterations,
+ {exsp,shift}, Iterations,
 TMarkUniform32Bit, OverheadUniform32Bit),
 _ =
 measure_1(
@@ -1328,7 +1451,7 @@ do_measure(Iterations) ->
 end
 end
 end,
- mcg35_inline, Iterations,
+ {mcg35,raw}, Iterations,
 TMarkUniformFullRange, OverheadUniformFullRange),
 _ =
 measure_1(
@@ -1341,21 +1464,52 @@ do_measure(Iterations) ->
 end
 end
 end,
- lcg35_inline, Iterations,
+ {lcg35,raw}, Iterations,
 TMarkUniformFullRange, OverheadUniformFullRange),
 _ =
 measure_1(
 fun (_Mod, _State) ->
- Range = 1 bsl 64,
+ Range = (16#7f17555 bsl 32) - 1,
 fun (St0) ->
- {V, St1} = rand:splitmix64_next(St0),
+ St1 = rand:mwc59(St0),
+ V = St1,
+ if
+ is_integer(V), 1 =< V, V < Range ->
+ St1
+ end
+ end
+ end,
+ {mwc59,raw}, Iterations,
+ TMarkUniformFullRange, OverheadUniformFullRange),
+ _ =
+ measure_1(
+ fun (_Mod, _State) ->
+ Range = 1 bsl 58,
+ fun (St0) ->
+ St1 = rand:mwc59(St0),
+ V = rand:mwc59_value(St1),
 if
 is_integer(V), 0 =< V, V < Range ->
 St1
 end
 end
 end,
- splitmix64_inline, Iterations,
+ {mwc59,value}, Iterations,
+ TMarkUniformFullRange, OverheadUniformFullRange),
+ _ =
+ measure_1(
+ fun (_Mod, _State) ->
+ Range = 1 bsl 58,
+ fun (St0) ->
+ St1 = rand:mwc59(St0),
+ V = rand:mwc59_full_value(St1),
+ if
+ is_integer(V), 0 =< V, V < Range ->
+ St1
+ end
+ end
+ end,
+ {mwc59,full_value}, Iterations,
 TMarkUniformFullRange, OverheadUniformFullRange),
 _ =
 measure_1(
@@ -1369,7 +1523,21 @@ do_measure(Iterations) ->
 end
 end
 end,
- exsp_inline, Iterations,
+ {exsp,next}, Iterations,
+ TMarkUniformFullRange, OverheadUniformFullRange),
+ _ =
+ measure_1(
+ fun (_Mod, _State) ->
+ Range = 1 bsl 64,
+ fun (St0) ->
+ {V, St1} = rand:splitmix64_next(St0),
+ if
+ is_integer(V), 0 =< V, V < Range ->
+ St1
+ end
+ end
+ end,
+ splitmix64_inline, Iterations,
 TMarkUniformFullRange, OverheadUniformFullRange),
 _ =
 measure_1(
@@ -1413,7 +1581,7 @@ do_measure(Iterations) ->
 end
 end
 end,
- lcg35_procdict, Iterations,
+ {lcg35,procdict}, Iterations,
 TMarkUniformFullRange, OverheadUniformFullRange),
 %%
 ct:pal("~nRNG uniform integer full range + 1 performance~n",[]),
@@ -1502,7 +1670,16 @@ do_measure(Iterations) ->
 ?CHECK_BYTE_SIZE(
 lcg35_bytes(ByteSize, St0), ByteSize, Bin, St1)
 end
- end, lcg35_bytes, Iterations,
+ end, {lcg35,bytes}, Iterations,
+ TMarkBytes1, OverheadBytes1),
+ _ =
+ measure_1(
+ fun (_Mod, _State) ->
+ fun (St0) ->
+ ?CHECK_BYTE_SIZE(
+ mwc59_bytes(ByteSize, St0), ByteSize, Bin, St1)
+ end
+ end, {mwc59,bytes}, Iterations,
 TMarkBytes1, OverheadBytes1),
 %%
 ByteSize2 = 1000, % At about 100 bytes crypto_bytes breaks even to exsss
@@ -1529,7 +1706,16 @@ do_measure(Iterations) ->
 ?CHECK_BYTE_SIZE(
 lcg35_bytes(ByteSize2, St0), ByteSize2, Bin, St1)
 end
- end, lcg35_bytes, Iterations div 50,
+ end, {lcg35,bytes}, Iterations div 50,
+ TMarkBytes2, OverheadBytes2),
+ _ =
+ measure_1(
+ fun (_Mod, _State) ->
+ fun (St0) ->
+ ?CHECK_BYTE_SIZE(
+ mwc59_bytes(ByteSize2, St0), ByteSize2, Bin, St1)
+ end
+ end, {mwc59,bytes}, Iterations div 50,
 TMarkBytes2, OverheadBytes2),
 %%
 ct:pal("~nRNG uniform float performance~n",[]),
@@ -1554,7 +1740,7 @@ do_measure(Iterations) ->
 end
 end
 end,
- mcg35_inline, Iterations,
+ {mcg35,float35}, Iterations,
 TMarkUniformFloat, OverheadUniformFloat),
 _ =
 measure_1(
@@ -1568,7 +1754,21 @@ do_measure(Iterations) ->
 end
 end
 end,
- lcg35_inline, Iterations,
+ {lcg35,float35}, Iterations,
+ TMarkUniformFloat, OverheadUniformFloat),
+ _ =
+ measure_1(
+ fun (_Mod, _State) ->
+ fun (St0) ->
+ St1 = rand:mwc59(St0),
+ %% Too few bits for full mantissa
+ case rand:mwc59_float(St1) of
+ R when is_float(R), 0.0 =< R, R < 1.0 ->
+ St1
+ end
+ end
+ end,
+ {mwc59,float}, Iterations,
 TMarkUniformFloat, OverheadUniformFloat),
 _ =
 measure_1(
@@ -1581,7 +1781,7 @@ do_measure(Iterations) ->
 end
 end
 end,
- exsp_inline, Iterations,
+ {exsp,float}, Iterations,
 TMarkUniformFloat, OverheadUniformFloat),
 %%
 ct:pal("~nRNG uniform_real float performance~n",[]),
@@ -1710,26 +1910,25 @@ measure_init(Alg) ->
 {?MODULE, undefined};
 system_time ->
 {?MODULE, undefined};
- mcg35_inline ->
- {_, S} = rand:seed_s(dummy),
- {rand, (S rem ((1 bsl 35)-31 - 1)) + 1};
- lcg35_inline ->
- {_, S} = rand:seed_s(dummy),
- {rand, S bsr (58-35)};
- lcg35_bytes ->
- {_, S} = rand:seed_s(dummy),
- {rand, S bsr (58-35)};
 splitmix64_inline ->
 {rand, erlang:unique_integer()};
- lcg35_procdict ->
- {_, S} = rand:seed_s(dummy),
- _ = put(Alg, S bsr (58-35)),
- {rand, undefined};
- exsp_inline ->
- {_, S} = rand:seed_s(exsp),
- {rand, S};
 procdict ->
 {rand, rand:seed(exsss)};
+ {Name, Tag} ->
+ case Name of
+ lcg35 when Tag =:= procdict ->
+ _ = put(lcg35_procdict, rand:lcg35_seed()),
+ {rand, undefined};
+ lcg35 ->
+ {rand, rand:lcg35_seed()};
+ mcg35 ->
+ {rand, rand:mcg35_seed()};
+ mwc59 ->
+ {rand, rand:mwc59_seed()};
+ exsp ->
+ {_, S} = rand:seed_s(exsp),
+ {rand, S}
+ end;
 _ ->
 {rand, rand:seed_s(Alg)}
 end.
@@ -1762,6 +1961,31 @@ lcg35_bytes(N, R0, Bin) ->
 V = R1 bsr (35 - 32),
 lcg35_bytes(N - 4, R1, <<Bin/binary, V:32>>).
 
+mwc59_bytes(N, R0) ->
+ mwc59_bytes(N, R0, <<>>).
+%%
+mwc59_bytes(N, R0, Bin) when is_integer(N), 7*4 =< N ->
+ R1 = rand:mwc59(R0),
+ R2 = rand:mwc59(R1),
+ R3 = rand:mwc59(R2),
+ R4 = rand:mwc59(R3),
+ V1 = rand:mwc59_value(R1),
+ V2 = rand:mwc59_value(R2),
+ V3 = rand:mwc59_value(R3),
+ V4 = rand:mwc59_value(R4),
+ mwc59_bytes(N-7*4, R4, <<Bin/binary, V1:56, V2:56, V3:56, V4:56>>);
+mwc59_bytes(N, R0, Bin) when is_integer(N), 7 =< N ->
+ R1 = rand:mwc59(R0),
+ V = rand:mwc59_value(R1),
+ mwc59_bytes(N-7, R1, <<Bin/binary, V:56>>);
+mwc59_bytes(N, R0, Bin) when is_integer(N), 0 < N ->
+ R1 = rand:mwc59(R0),
+ Bits = N bsl 3,
+ V = rand:mwc59_value(R1),
+ {<<Bin/binary, V:Bits>>, R1};
+mwc59_bytes(0, R0, Bin) ->
+ {Bin, R0}.
+
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% The jump sequence tests has two parts
-- 
2.35.3

Places

File 2951-Implement-the-custom-made-mwc59-PRNG.patch of Package erlang (Revision 0b3a6d6e89afeaba4dcdf8fb3d620103)

Places