File 2961-Change-mwc59-to-returns-59-bits-and-use-the-high.patch of Package erlang
From 41b986f5bb43c8e9a7463bbffbea1f4cc3044fcb Mon Sep 17 00:00:00 2001
From: Raimo Niskanen <raimo@erlang.org>
Date: Tue, 10 May 2022 00:29:49 +0200
Subject: [PATCH 1/8] Change mwc59 to returns 59 bits and use the high
The high bits are more to be trusted, and to get it efficient
all 59 bits are returned so the application can shift down
for a power of 2 range.
---
lib/stdlib/doc/src/rand.xml | 39 ++++++++++++++++++++++--------
lib/stdlib/src/rand.erl | 44 ++++++++++++++++------------------
lib/stdlib/test/rand_SUITE.erl | 41 +++++++++++++++----------------
3 files changed, 69 insertions(+), 55 deletions(-)
diff --git a/lib/stdlib/doc/src/rand.xml b/lib/stdlib/doc/src/rand.xml
index 25dec14a51..1cf2cca6fd 100644
--- a/lib/stdlib/doc/src/rand.xml
+++ b/lib/stdlib/doc/src/rand.xml
@@ -926,14 +926,14 @@ end.</pre>
<seemfa marker="#mwc59_fast_value/1">
<c>mwc59_fast_value</c>
</seemfa>
- is a fast scrambler that gives 32 decent bits, but still
- some problems in 2- and 3-dimensional collision tests
- show through.
+ is a fast scrambler that returns a 59-bit number with
+ 32 decent bits, but still some problems in 2- and
+ 3-dimensional collision tests show through.
The slightly slower
<seemfa marker="#mwc59_value/1">
<c>mwc59_value</c>
</seemfa>
- scrambler returns 58 bits of very good quality, and
+ scrambler returns 59 bits of very good quality, and
<seemfa marker="#mwc59_float/1"><c>mwc59_float</c></seemfa>
returns a <c>float()</c> of very good quality.
</p>
@@ -978,13 +978,22 @@ end.</pre>
<fsummary>Return the generator value.</fsummary>
<desc>
<p>
- Returns a 58-bit value <c><anno>V</anno></c>
+ Returns a 59-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 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 decent quality.
+ enough that the 32 highest bits are of decent quality.
+ </p>
+ <p>
+ To extract a power of two number it is recommended
+ to shift down the high bits, but for an arbitrary
+ range a `rem` operation on the whole value
+ can be used. Be careful to not accidentaly create
+ a bignum when handling the value <c><anno>V</anno></c>.
+ </p>
+ <p>
It is not recommended to generate numbers
in a range > 2^32 with this function.
</p>
@@ -995,15 +1004,25 @@ end.</pre>
<fsummary>Return the generator value.</fsummary>
<desc>
<p>
- Returns a 58-bit value <c><anno>V</anno></c>
+ Returns a 59-bit value <c><anno>V</anno></c>
from a generator state <c><anno>CX</anno></c>.
The generator state is scrambled using
an 4-bit followed by a 27-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 very good quality.
- Beware of bias in the generated numbers if
- generating on a non power of 2 range too close to 2^58.
+ enough that all 59 bits are of very good quality.
+ </p>
+ <p>
+ To extract a power of two number it is recommended
+ to shift down the high bits, but for an arbitrary
+ range a `rem` operation on the whole value
+ can be used. Be careful to not accidentaly create
+ a bignum when handling the value <c><anno>V</anno></c>.
+ </p>
+ <p>
+ As usual: beware of bias in the generated numbers if
+ generating on range too close to 2^59 that is not
+ a power of 2.
</p>
</desc>
</func>
diff --git a/lib/stdlib/src/rand.erl b/lib/stdlib/src/rand.erl
index 6c85a79323..adefa8ab89 100644
--- a/lib/stdlib/src/rand.erl
+++ b/lib/stdlib/src/rand.erl
@@ -53,6 +53,7 @@
exs1024_next/1, exs1024_calc/2,
exro928_next_state/4,
exrop_next/1, exrop_next_s/2,
+ mwc59_value/1,
get_52/1, normal_kiwi/1]}).
-define(DEFAULT_ALG_HANDLER, exsss).
@@ -1507,32 +1508,37 @@ mwc59(CX0) -> % when is_integer(CX0), 1 =< CX0, CX0 < ?MWC59_P ->
X = ?MASK(?MWC59_B, CX),
?MWC59_A * X + C.
--spec mwc59_fast_value(CX :: mwc59_state()) -> V :: 0..?MASK(58).
+-spec mwc59_fast_value(CX :: mwc59_state()) -> V :: 0..?MASK(59).
mwc59_fast_value(CX1) -> % when is_integer(CX1), 1 =< CX1, CX1 < ?MWC59_P ->
- CX = ?MASK(58, CX1),
- CX bxor ?BSL(58, CX, ?MWC59_XS).
+ CX = ?MASK(59, CX1),
+ CX bxor ?BSL(59, CX, ?MWC59_XS).
--spec mwc59_value(CX :: mwc59_state()) -> V :: 0..?MASK(58).
+-spec mwc59_value(CX :: mwc59_state()) -> V :: 0..?MASK(59).
mwc59_value(CX1) -> % when is_integer(CX1), 1 =< CX1, CX1 < ?MWC59_P ->
- CX = ?MASK(58, CX1),
- CX2 = CX bxor ?BSL(58, CX, ?MWC59_XS1),
- CX2 bxor ?BSL(58, CX2, ?MWC59_XS2).
+ CX = ?MASK(59, CX1),
+ CX2 = CX bxor ?BSL(59, CX, ?MWC59_XS1),
+ CX2 bxor ?BSL(59, 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),
- CX3 = CX2 bxor ?BSL(53, CX2, ?MWC59_XS2),
- CX3 * ?TWO_POW_MINUS53.
+ mwc59_value(CX1) * (1/(1 bsl 59)).
-spec mwc59_seed() -> CX :: mwc59_state().
mwc59_seed() ->
{A1, A2, A3} = default_seed(),
- mwc59_seed(seed64_int([A1, A2, A3])).
+ {X1,_} = splitmix64_next(A1),
+ {X2,_} = splitmix64_next(A2),
+ {X3,_} = splitmix64_next(A3),
+ mwc59_seed_state(X1 bxor X2 bxor X3).
-spec mwc59_seed(S :: integer()) -> CX :: mwc59_state().
mwc59_seed(S) when is_integer(S) ->
- mod(?MWC59_P - 1, S) + 1.
+ {X,_} = splitmix64_next(mod(?BIT(64), S)),
+ mwc59_seed_state(X).
+
+mwc59_seed_state(X) ->
+ %% Ensure non-zero carry and within state range
+ ?BIT(?MWC59_B) bor ?MASK(58, X).
@@ -1612,16 +1618,6 @@ 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:
%%
@@ -1988,6 +1984,8 @@ bc(V, B, N) -> bc(V, B bsr 1, N - 1).
%% Non-negative rem
+mod(Q, X) when 0 =< X, X < Q ->
+ X;
mod(Q, X) ->
Y = X rem Q,
if
diff --git a/lib/stdlib/test/rand_SUITE.erl b/lib/stdlib/test/rand_SUITE.erl
index c03a687b3a..04449a3d29 100644
--- a/lib/stdlib/test/rand_SUITE.erl
+++ b/lib/stdlib/test/rand_SUITE.erl
@@ -220,10 +220,10 @@ mwc59_api(CX0, 0) ->
V = 262716604851324112,
{V, V} = {V0, V},
W0 = rand:mwc59_value(CX0),
- W = 240206843777063376,
+ W = 528437219928775120,
{W, W} = {W0, W},
F0 = rand:mwc59_float(CX0),
- F = (W band ((1 bsl 53) - 1)) * (1 / (1 bsl 53)),
+ F = (W bsr (59-53)) * (1 / (1 bsl 53)),
{F, F} = {F0, F},
ok;
mwc59_api(CX, N)
@@ -232,9 +232,9 @@ mwc59_api(CX, N)
W = rand:mwc59_value(CX),
F = rand:mwc59_float(CX),
true = 0 =< V,
- true = V < 1 bsl 58,
+ true = V < 1 bsl 59,
true = 0 =< W,
- true = W < 1 bsl 58,
+ true = W < 1 bsl 59,
true = 0.0 =< F,
true = F < 1.0,
mwc59_api(rand:mwc59(CX), N - 1).
@@ -1231,13 +1231,13 @@ do_measure(Iterations) ->
Range = 1 bsl 32,
fun (St0) ->
St1 = rand:mwc59(St0),
- case rand:mwc59_fast_value(St1) band ((1 bsl 32)-1) of
+ case rand:mwc59_fast_value(St1) bsr (59-32) of
R when is_integer(R), 0 =< R, R < Range ->
St1
end
end
end,
- {mwc59,fast_mask}, Iterations,
+ {mwc59,fast_shift}, Iterations,
TMarkUniform32Bit, OverheadUniform32Bit),
_ =
measure_1(
@@ -1245,16 +1245,13 @@ do_measure(Iterations) ->
Range = 1 bsl 32,
fun (St0) ->
St1 = rand:mwc59(St0),
- case
- rand:mwc59_value(St1) band
- ((1 bsl 32)-1)
- of
+ case rand:mwc59_value(St1) bsr (59-32) of
R when is_integer(R), 0 =< R, R < Range ->
St1
end
end
end,
- {mwc59,value_mask}, Iterations,
+ {mwc59,value_shift}, Iterations,
TMarkUniform32Bit, OverheadUniform32Bit),
_ =
measure_1(
@@ -1369,12 +1366,12 @@ do_measure(Iterations) ->
_ =
measure_1(
fun (_Mod, _State) ->
- Range = 1 bsl 58,
+ Range = (1 bsl 59) - 1,
fun (St0) ->
St1 = rand:mwc59(St0),
V = rand:mwc59_fast_value(St1),
if
- is_integer(V), 0 =< V, V < Range ->
+ is_integer(V), 0 =< V, V =< Range ->
St1
end
end
@@ -1384,12 +1381,12 @@ do_measure(Iterations) ->
_ =
measure_1(
fun (_Mod, _State) ->
- Range = 1 bsl 58,
+ Range = (1 bsl 59) - 1,
fun (St0) ->
St1 = rand:mwc59(St0),
V = rand:mwc59_value(St1),
if
- is_integer(V), 0 =< V, V < Range ->
+ is_integer(V), 0 =< V, V =< Range ->
St1
end
end
@@ -1600,7 +1597,6 @@ do_measure(Iterations) ->
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
@@ -1791,19 +1787,20 @@ mwc59_bytes(N, R0, Bin) when is_integer(N), 7*4 =< N ->
R2 = rand:mwc59(R1),
R3 = rand:mwc59(R2),
R4 = rand:mwc59(R3),
- V1 = rand:mwc59_fast_value(R1),
- V2 = rand:mwc59_fast_value(R2),
- V3 = rand:mwc59_fast_value(R3),
- V4 = rand:mwc59_fast_value(R4),
+ Shift = 59 - 56,
+ V1 = rand:mwc59_fast_value(R1) bsr Shift,
+ V2 = rand:mwc59_fast_value(R2) bsr Shift,
+ V3 = rand:mwc59_fast_value(R3) bsr Shift,
+ V4 = rand:mwc59_fast_value(R4) bsr Shift,
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_fast_value(R1),
+ V = rand:mwc59_fast_value(R1) bsr (59-56),
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_fast_value(R1),
+ V = rand:mwc59_fast_value(R1) bsr (59-Bits),
{<<Bin/binary, V:Bits>>, R1};
mwc59_bytes(0, R0, Bin) ->
{Bin, R0}.
--
2.35.3
