Merge pull request #1969 from RaimoNiskanen/raimo/stdlib/rand-xorshift116ss

OTP-14731 Implement 'exsss' (Xorshift116**) as new default 'rand' algorithm

The new algorithm is a combination of the Xorshift116 ('exsp') state update and a new scrambler  "StarStar" from the 2018 paper "Scrambled Linear Pseudorandom Number Generators" by David Blackman and Sebastiano Vigna. This combination should not have the caveat of weak low bits that the previous default algorithm(s) have had, with the cost of about 10% lower speed.

1 files changed, 54 insertions, 27 deletions

diff --git a/lib/stdlib/doc/src/rand.xml b/lib/stdlib/doc/src/rand.xml
index 25eec216ef..8e657698c6 100644
--- a/lib/stdlib/doc/src/rand.xml
+++ b/lib/stdlib/doc/src/rand.xml
@@ -38,34 +38,50 @@
     <p>
       This module provides a pseudo random number generator.
       The module contains a number of algorithms.
-      The uniform distribution algorithms use the
+      The uniform distribution algorithms are based on the
       <url href="http://xorshift.di.unimi.it">
-	xoroshiro116+ and xorshift1024* algorithms by Sebastiano Vigna.
+	Xoroshiro and Xorshift algorithms
       </url>
+      by Sebastiano Vigna.
       The normal distribution algorithm uses the
       <url href="http://www.jstatsoft.org/v05/i08">
 	Ziggurat Method by Marsaglia and Tsang
       </url>
       on top of the uniform distribution algorithm.
     </p>
-    <p>For some algorithms, jump functions are provided for generating
-        non-overlapping sequences for parallel computations.
-        The jump functions perform calculations
-        equivalent to perform a large number of repeated calls
-        for calculating new states. </p>
+    <p>
+      For most algorithms, jump functions are provided for generating
+      non-overlapping sequences for parallel computations.
+      The jump functions perform calculations
+      equivalent to perform a large number of repeated calls
+      for calculating new states.
+    </p>
 
     <p>The following algorithms are provided:</p>
 
     <taglist>
-      <tag><c>exrop</c></tag>
+      <tag><c>exsss</c></tag>
       <item>
-        <p>Xoroshiro116+, 58 bits precision and period of 2^116-1</p>
+        <p>Xorshift116**, 58 bits precision and period of 2^116-1</p>
         <p>Jump function: equivalent to 2^64 calls</p>
-      </item>
-      <tag><c>exs1024s</c></tag>
-      <item>
-        <p>Xorshift1024*, 64 bits precision and a period of 2^1024-1</p>
-        <p>Jump function: equivalent to 2^512 calls</p>
+	<p>
+	  This is the Xorshift116 generator combined with the StarStar scrambler
+	  from the 2018 paper by David Blackman and Sebastiano Vigna:
+	  <url href="http://vigna.di.unimi.it/ftp/papers/ScrambledLinear.pdf">
+	    Scrambled Linear Pseudorandom Number Generators
+	  </url>
+	</p>
+	<p>
+	  The generator does not need 58-bit rotates so it is faster
+	  than the Xoroshiro116 generator, and when combined with
+	  the StarStar scrambler it does not have any weak low bits
+	  like <c>exrop</c> (Xoroshiro116+).
+	</p>
+	<p>
+	  Alas, this combination is about 10% slower than <c>exrop</c>,
+	  but is despite that the default algorithm thanks to its
+	  statistical qualities.
+	</p>
       </item>
       <tag><c>exro928ss</c></tag>
       <item>
@@ -77,8 +93,8 @@
 	  <url href="http://vigna.di.unimi.it/ftp/papers/ScrambledLinear.pdf">
 	    Scrambled Linear Pseudorandom Number Generators
 	  </url>
-	  that on a 64 bit Erlang system executes only about 30% slower than
-	  the default <c>exrop</c> algorithm but with much longer period
+	  that on a 64 bit Erlang system executes only about 40% slower than
+	  the default <c>exsss</c> algorithm but with much longer period
 	  and better statistical properties, and on the flip side
 	  a larger state.
 	</p>
@@ -87,6 +103,16 @@
 	  the 58 bit adaption.
 	</p>
       </item>
+      <tag><c>exrop</c></tag>
+      <item>
+        <p>Xoroshiro116+, 58 bits precision and period of 2^116-1</p>
+        <p>Jump function: equivalent to 2^64 calls</p>
+      </item>
+      <tag><c>exs1024s</c></tag>
+      <item>
+        <p>Xorshift1024*, 64 bits precision and a period of 2^1024-1</p>
+        <p>Jump function: equivalent to 2^512 calls</p>
+      </item>
       <tag><c>exsp</c></tag>
       <item>
         <p>Xorshift116+, 58 bits precision and period of 2^116-1</p>
@@ -103,7 +129,7 @@
     </taglist>
 
     <p>
-      The default algorithm is <c>exrop</c> (Xoroshiro116+).
+      The default algorithm is <c>exsss</c> (Xorshift116**).
       If a specific algorithm is
       required, ensure to always use <seealso marker="#seed-1">
       <c>seed/1</c></seealso> to initialize the state.
@@ -174,19 +200,19 @@ R1 = rand:uniform(),</pre>
     <p>Use a specified algorithm:</p>
 
     <pre>
-_ = rand:seed(exs1024s),
+_ = rand:seed(exs928ss),
 R2 = rand:uniform(),</pre>
 
     <p>Use a specified algorithm with a constant seed:</p>
 
     <pre>
-_ = rand:seed(exs1024s, {123, 123534, 345345}),
+_ = rand:seed(exs928ss, {123, 123534, 345345}),
 R3 = rand:uniform(),</pre>
 
    <p>Use the functional API with a non-constant seed:</p>
 
    <pre>
-S0 = rand:seed_s(exrop),
+S0 = rand:seed_s(exsss),
 {R4, S1} = rand:uniform_s(S0),</pre>
 
    <p>Textbook basic form Box-Muller standard normal deviate</p>
@@ -215,8 +241,9 @@ SND0 = math:sqrt(-2 * math:log(R5)) * math:cos(math:pi() * R6)</pre>
     </note>
 
     <p>
-      For all these generators except <c>exro928ss</c> the lowest bit(s)
-      has got a slightly less random behaviour than all other bits.
+      For all these generators except <c>exro928ss</c> and <c>exsss</c>
+      the lowest bit(s) has got a slightly less
+      random behaviour than all other bits.
       1 bit for <c>exrop</c> (and <c>exsp</c>),
       and 3 bits for <c>exs1024s</c>.
       See for example the explanation in the
@@ -231,7 +258,7 @@ up to (and included) 16TB, with the exception of binary rank tests,
 which fail due to the lowest bit being an LFSR; all other bits pass all
 tests. We suggest to use a sign test to extract a random Boolean value.</pre>
     <p>
-      If this is a problem; to generate a boolean
+      If this is a problem; to generate a boolean with these algorithms
       use something like this:
     </p>
     <pre>(rand:uniform(16) > 8)</pre>
@@ -299,19 +326,19 @@ tests. We suggest to use a sign test to extract a random Boolean value.</pre>
       </desc>
     </datatype>
     <datatype>
-      <name name="exrop_state"/>
+      <name name="exsplus_state"/>
       <desc><p>Algorithm specific internal state</p></desc>
     </datatype>
     <datatype>
-      <name name="exs1024_state"/>
+      <name name="exro928_state"/>
       <desc><p>Algorithm specific internal state</p></desc>
     </datatype>
     <datatype>
-      <name name="exro928_state"/>
+      <name name="exrop_state"/>
       <desc><p>Algorithm specific internal state</p></desc>
     </datatype>
     <datatype>
-      <name name="exsplus_state"/>
+      <name name="exs1024_state"/>
       <desc><p>Algorithm specific internal state</p></desc>
     </datatype>
     <datatype>