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.

3 files changed, 247 insertions, 53 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>
diff --git a/lib/stdlib/src/rand.erl b/lib/stdlib/src/rand.erl
index 9854c778a1..3a9a1e007b 100644
--- a/lib/stdlib/src/rand.erl
+++ b/lib/stdlib/src/rand.erl
@@ -43,13 +43,13 @@
 %% Debug
 -export([make_float/3, float2str/1, bc64/1]).
 
--compile({inline, [exs64_next/1, exsplus_next/1,
+-compile({inline, [exs64_next/1, exsplus_next/1, exsss_next/1,
 		   exs1024_next/1, exs1024_calc/2,
                    exro928_next_state/4,
                    exrop_next/1, exrop_next_s/2,
 		   get_52/1, normal_kiwi/1]}).
 
--define(DEFAULT_ALG_HANDLER, exrop).
+-define(DEFAULT_ALG_HANDLER, exsss).
 -define(SEED_DICT, rand_seed).
 
 %% =====================================================================
@@ -86,7 +86,7 @@
 
 %% This depends on the algorithm handler function
 -type alg_state() ::
-	exrop_state() | exs1024_state() |  exro928_state() |  exsplus_state() |
+	exsplus_state() | exro928_state() |  exrop_state() | exs1024_state() |
 	exs64_state() | term().
 
 %% This is the algorithm handling definition within this module,
@@ -131,7 +131,7 @@
 %% Algorithm state
 -type state() :: {alg_handler(), alg_state()}.
 -type builtin_alg() ::
-	exrop | exs1024s | exro928ss | exsp | exs64 | exsplus | exs1024.
+	exsss | exro928ss | exrop | exs1024s | exsp | exs64 | exsplus | exs1024.
 -type alg() :: builtin_alg() | atom().
 -type export_state() :: {alg(), alg_state()}.
 -type seed() :: [integer()] | integer() | {integer(), integer(), integer()}.
@@ -139,7 +139,7 @@
    [builtin_alg/0, alg/0, alg_handler/0, alg_state/0,
     state/0, export_state/0, seed/0]).
 -export_type(
-   [exrop_state/0, exs1024_state/0, exro928_state/0, exsplus_state/0,
+   [exsplus_state/0, exro928_state/0, exrop_state/0, exs1024_state/0,
     exs64_state/0]).
 
 %% =====================================================================
@@ -618,6 +618,11 @@ mk_alg(exsp) ->
        uniform=>fun exsp_uniform/1, uniform_n=>fun exsp_uniform/2,
        jump=>fun exsplus_jump/1},
      fun exsplus_seed/1};
+mk_alg(exsss) ->
+    {#{type=>exsss, bits=>58, next=>fun exsss_next/1,
+       uniform=>fun exsss_uniform/1, uniform_n=>fun exsss_uniform/2,
+       jump=>fun exsplus_jump/1},
+     fun exsss_seed/1};
 mk_alg(exs1024) ->
     {#{type=>exs1024, max=>?MASK(64), next=>fun exs1024_next/1,
        jump=>fun exs1024_jump/1},
@@ -675,6 +680,36 @@ exs64_next(R) ->
 %% 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.
+%%
+%% Reference C code for Xorshift116+ and Xorshift116**
+%%
+%% #include <stdint.h>
+%%
+%% #define MASK(b, v) (((UINT64_C(1) << (b)) - 1) & (v))
+%% #define BSL(b, v, n) (MASK((b)-(n), (v)) << (n))
+%% #define ROTL(b, v, n) (BSL((b), (v), (n)) | ((v) >> ((b)-(n))))
+%%
+%% uint64_t s[2];
+%%
+%% uint64_t next(void) {
+%%     uint64_t s1 = s[0];
+%%     const uint64_t s0 = s[1];
+%%
+%%     s1 ^= BSL(58, s1, 24); // a
+%%     s1 ^= s0 ^ (s1 >> 11) ^ (s0 >> 41); // b, c
+%%     s[0] = s0;
+%%     s[1] = s1;
+%%
+%%     const uint64_t result_plus = MASK(58, s0 + s1);
+%%     uint64_t result_starstar = s0;
+%%     result_starstar = MASK(58, result_starstar * 5);
+%%     result_starstar = ROTL(58, result_starstar, 7);
+%%     result_starstar = MASK(58, result_starstar * 9);
+%%
+%%     return result_plus;
+%%     return result_starstar;
+%% }
+%%
 %% =====================================================================
 -opaque exsplus_state() :: nonempty_improper_list(uint58(), uint58()).
 
@@ -697,16 +732,62 @@ exsplus_seed({A1, A2, A3}) ->
                  tl(R1)]),
     R2.
 
+-dialyzer({no_improper_lists, exsss_seed/1}).
+
+exsss_seed(L) when is_list(L) ->
+    [S0,S1] = seed58_nz(2, L),
+    [S0|S1];
+exsss_seed(X) when is_integer(X) ->
+    [S0,S1] = seed58(2, ?MASK(64, X)),
+    [S0|S1];
+%%
+%% Seed from traditional integer triple - mix into splitmix
+exsss_seed({A1, A2, A3}) ->
+    {_, X0} = seed58(?MASK(64, A1)),
+    {S0, X1} = seed58(?MASK(64, A2) bxor X0),
+    {S1, _} = seed58(?MASK(64, A3) bxor X1),
+    [S0|S1].
+
+%% Advance Xorshift116 state one step
+-define(
+   exs_next(S0, S1, S1_b),
+   begin
+       S1_b = S1 bxor ?BSL(58, S1, 24),
+       S1_b bxor S0 bxor (S1_b bsr 11) bxor (S0 bsr 41)
+   end).
+
+-define(
+   scramble_starstar(S, V_a, V_b),
+   begin
+       %% The multiply by add shifted trick avoids creating bignums
+       %% which improves performance significantly
+       %%
+       V_a = ?MASK(58, S + ?BSL(58, S, 2)), % V_a = S * 5
+       V_b = ?ROTL(58, V_a, 7),
+       ?MASK(58, V_b + ?BSL(58, V_b, 3)) % V_b * 9
+   end).
+
 -dialyzer({no_improper_lists, exsplus_next/1}).
 
-%% Advance xorshift116+ state for one step and generate 58bit unsigned integer
+%% Advance state 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 ?BSL(58, S1, 24),
-    S12 = S11 bxor S0 bxor (S11 bsr 11) bxor (S0 bsr 41),
-    {?MASK(58, S0 + S12), [S0|S12]}.
+    NewS1 = ?exs_next(S0, S1, S1_1),
+    {?MASK(58, S0 + NewS1), [S0|NewS1]}.
+%%    %% Note: members s0 and s1 are swapped here
+%%    S11 = S1 bxor ?BSL(58, S1, 24),
+%%    S12 = S11 bxor S0 bxor (S11 bsr 11) bxor (S0 bsr 41),
+%%    {?MASK(58, S0 + S12), [S0|S12]}.
+
+-dialyzer({no_improper_lists, exsss_next/1}).
 
+-spec exsss_next(exsplus_state()) -> {uint58(), exsplus_state()}.
+exsss_next([S1|S0]) ->
+    %% Note: members s0 and s1 are swapped here
+    NewS1 = ?exs_next(S0, S1, S1_1),
+    {?scramble_starstar(S0, V_0, V_1), [S0|NewS1]}.
+%%    {?MASK(58, S0 + NewS1), [S0|NewS1]}.
 
 exsp_uniform({Alg, R0}) ->
     {I, R1} = exsplus_next(R0),
@@ -714,18 +795,48 @@ exsp_uniform({Alg, R0}) ->
     %% randomness quality than the others
     {(I bsr (58-53)) * ?TWO_POW_MINUS53, {Alg, R1}}.
 
+exsss_uniform({Alg, R0}) ->
+    {I, R1} = exsss_next(R0),
+    {(I bsr (58-53)) * ?TWO_POW_MINUS53, {Alg, R1}}.
+
 exsp_uniform(Range, {Alg, R}) ->
     {V, R1} = exsplus_next(R),
     MaxMinusRange = ?BIT(58) - Range,
     ?uniform_range(Range, Alg, R1, V, MaxMinusRange, I).
 
+exsss_uniform(Range, {Alg, R}) ->
+    {V, R1} = exsss_next(R),
+    MaxMinusRange = ?BIT(58) - Range,
+    ?uniform_range(Range, Alg, R1, V, MaxMinusRange, I).
+
 
-%% This is the jump function for the exsplus generator, equivalent
+%% This is the jump function for the exs* generators,
+%% i.e the Xorshift116 generators,  equivalent
 %% to 2^64 calls to next/1; it can be used to generate 2^52
 %% non-overlapping subsequences for parallel computations.
 %% Note: the jump function takes 116 times of the execution time of
 %% next/1.
-
+%%
+%% #include <stdint.h>
+%%
+%% void jump(void) {
+%%   static const uint64_t JUMP[] = { 0x02f8ea6bc32c797,
+%% 				   0x345d2a0f85f788c };
+%%   int i, b;
+%%   uint64_t s0 = 0;
+%%   uint64_t s1 = 0;
+%%   for(i = 0; i < sizeof JUMP / sizeof *JUMP; i++)
+%%     for(b = 0; b < 58; b++) {
+%%       if (JUMP[i] & 1ULL << b) {
+%% 	s0 ^= s[0];
+%% 	s1 ^= s[1];
+%%       }
+%%       next();
+%%     }
+%%   s[0] = s0;
+%%   s[1] = s1;
+%% }
+%%
 %% -define(JUMPCONST, 16#000d174a83e17de2302f8ea6bc32c797).
 %% split into 58-bit chunks
 %% and two iterative executions
@@ -973,13 +1084,14 @@ exro928ss_next({[S15,S0|Ss], Rs}) ->
     %% {S, R, T} = {5, 7, 9}
     %% const uint64_t result_starstar = rotl(s0 * S, R) * T;
     %%
-    %% The multiply by add shifted trick avoids creating bignums
-    %% which improves performance significantly
-    %%
-    V0 = ?MASK(58, S0 + ?BSL(58, S0, 2)), % V0 = S0 * 5
-    V1 = ?ROTL(58, V0, 7),
-    V = ?MASK(58, V1 + ?BSL(58, V1, 3)), % V = V1 * 9
-    {V, SR};
+    {?scramble_starstar(S0, V_0, V_1), SR};
+%%    %% The multiply by add shifted trick avoids creating bignums
+%%    %% which improves performance significantly
+%%    %%
+%%    V0 = ?MASK(58, S0 + ?BSL(58, S0, 2)), % V0 = S0 * 5
+%%    V1 = ?ROTL(58, V0, 7),
+%%    V = ?MASK(58, V1 + ?BSL(58, V1, 3)), % V = V1 * 9
+%%    {V, SR};
 exro928ss_next({[S15], Rs}) ->
     exro928ss_next({[S15|lists:reverse(Rs)], []}).
 
diff --git a/lib/stdlib/test/rand_SUITE.erl b/lib/stdlib/test/rand_SUITE.erl
index 6c9fab51dc..7685c17967 100644
--- a/lib/stdlib/test/rand_SUITE.erl
+++ b/lib/stdlib/test/rand_SUITE.erl
@@ -79,7 +79,7 @@ test() ->
       end, Tests).
 
 algs() ->
-    [exrop, exsp, exs1024s, exs64, exsplus, exs1024, exro928ss].
+    [exsss, exrop, exsp, exs1024s, exs64, exsplus, exs1024, exro928ss].
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -276,7 +276,7 @@ gen(Algo) ->
         if
             Algo =:= exs64 -> %% Printed with orig 'C' code and this seed
                 rand:seed_s(exs64, [12345678]);
-            Algo =:= exsplus; Algo =:= exsp; Algo =:= exrop ->
+            Algo =:= exsplus; Algo =:= exsp; Algo =:= exrop; Algo =:= exsss ->
                 %% Printed with orig 'C' code and this seed
                 rand:seed_s(Algo, [12345678,12345678]);
             Algo =:= exs1024; Algo =:= exs1024s; Algo =:= exro928ss ->
@@ -438,7 +438,7 @@ stats_standard_normal_box_muller(Config) when is_list(Config) ->
                     ([S|Z]) ->
                         {Z, [S]}
                 end,
-            State = [rand:seed(exrop)],
+            State = [rand:seed(exsss)],
             stats_standard_normal(NormalS, State, 3)
     catch error:_ ->
             {skip, "math:erfc/1 not supported"}
@@ -463,7 +463,7 @@ stats_standard_normal_box_muller_2(Config) when is_list(Config) ->
                     ([S|Z]) ->
                         {Z, [S]}
                 end,
-            State = [rand:seed(exrop)],
+            State = [rand:seed(exsss)],
             stats_standard_normal(NormalS, State, 3)
     catch error:_ ->
             {skip, "math:erfc/1 not supported"}
@@ -475,7 +475,7 @@ stats_standard_normal(Config) when is_list(Config) ->
     try math:erfc(1.0) of
         _ ->
             stats_standard_normal(
-              fun rand:normal_s/1, rand:seed_s(exrop), Retries)
+              fun rand:normal_s/1, rand:seed_s(exsss), Retries)
     catch error:_ ->
             {skip, "math:erfc/1 not supported"}
     end.
@@ -1071,7 +1071,7 @@ do_measure(_Config) ->
                     end,
                     State)
           end,
-          exrop, TMarkNormalFloat),
+          exsss, TMarkNormalFloat),
     ok.
 
 -define(LOOP_MEASURE, (?LOOP div 5)).
@@ -1170,7 +1170,7 @@ gen_jump_1(Algo) ->
 	    catch
 		error:not_implemented -> [error_not_implemented]
 	    end;
-	_ when Algo =:= exsplus; Algo =:= exsp; Algo =:= exrop ->
+	_ when Algo =:= exsplus; Algo =:= exsp; Algo =:= exrop; Algo =:= exsss ->
 	    %% Printed with orig 'C' code and this seed
 	    gen_jump_2(
 	      rand:seed_s(Algo, [12345678,12345678]));
@@ -1224,7 +1224,7 @@ gen_jump_p1(Algo) ->
 	    catch
 		error:not_implemented -> [error_not_implemented]
 	    end;
-	_ when Algo =:= exsplus; Algo =:= exsp; Algo =:= exrop ->
+	_ when Algo =:= exsplus; Algo =:= exsp; Algo =:= exrop; Algo =:= exsss ->
 	    %% Printed with orig 'C' code and this seed
 	    gen_jump_p2(
 	      rand:seed(Algo, [12345678,12345678]));
@@ -1377,6 +1377,34 @@ reference_val(exsplus) ->
      16#3dd493b8012970f,16#be13bed1e00e5c,16#ceef033b74ae10,16#3da38c6a50abe03,
      16#15cbd1a421c7a8c,16#22794e3ec6ef3b1,16#26154d26e7ea99f,16#3a66681359a6ab6];
 
+reference_val(exsss) ->
+    [16#108e8d5b01,16#33b72092117209a,16#224d4d2961a2d0a,16#2c4c81aac3da48d,
+     16#2f4bc39bfc36f3a,16#41826d4c4d243a,16#19871b8bb4e23ee,16#3e2112cdf9384b1,
+     16#69801943bf91ab,16#2de1a603c31ec45,16#a90ca1991b831e,16#51ca29571a69a7,
+     16#93ce3e511906cf,16#93ebc5768aef75,16#2412f284b902ae7,16#1ac10e758410c52,
+     16#3f32494560368f6,16#39a5e82dcf0de95,16#3f4b14d59cc6a21,16#3174668db0b36ae,
+     16#1449812fb8bd54e,16#eaca1f8ece51e1,16#2564b2545fd23c1,16#3cf3a2d2217e0d7,
+     16#226f4164ba1d054,16#10dac9ae207ceef,16#17f2c4b2d40fcb9,16#1c1b282d386fdcb,
+     16#a264f450ba2912,16#2a0a1dd67e52666,16#2be84eb835cb1e1,16#2a1cd9aa16ccc37,
+     16#7dd5e8c2b3f490,16#254a3db4976c05b,16#2a0a67971ec1e63,16#13a0cbf7c0eed8a,
+     16#3192d7edc0a20bc,16#2705ad756292e84,16#3ec429a18119c81,16#25944b38baa975b,
+     16#291dcc43e3256f4,16#30d10b759237db,16#c1522a652058a,16#8ef1e9378381e6,
+     16#1f442f33c2439f4,16#186087710a73818,16#12887f94b2b8387,16#3e42e8b1f3c9b4b,
+     16#e462859d55f9d8,16#2356ae85be908de,16#15e96a927b3bc52,16#35c6dc52511ce46,
+     16#7bc0624ce66e01,16#33ab7d95b738322,16#26f01effc182aa0,16#1b66ae7eaafea88,
+     16#278f3dc14943b90,16#22178bc8d8faf28,16#396c37d53c11985,16#5e0d79d0b10f18,
+     16#1be3de3b5675ec,16#d4db298f1f4b50,16#2da6cb99bb5c7b1,16#130b2dc17d03be8,
+     16#f1847e7e059e9f,16#2da6591788326e7,16#222e4a18c24211c,16#949213ca49baab,
+     16#b5129fec56f6a2,16#30f25f1e926f43e,16#1ddd8d04445fb4d,16#15995b542514150,
+     16#1595fe879296296,16#e2f237a488453b,16#23e5cd2d6047890,16#3a5dc88fc954666,
+     16#89bca9969b103,16#5e6893cd35dc63,16#1fed534feeeef5a,16#26f40e2147ee558,
+     16#30c131a00625837,16#2618a7e617422e9,16#23630b297e45e7,16#1143b17502f3219,
+     16#15607dac41168da,16#2886bdc314b3fb8,16#465d1cc1536546,16#30b09123e3a02e4,
+     16#245a375f810be52,16#6a1b0792376a03,16#221425f59f2470f,16#867ce16dfac81c,
+     16#9c62d95fae9b58,16#380381db1394426,16#34908dedc01c324,16#1f0ff517089b561,
+     16#1571366dd873d32,16#3ee353dc56e192,16#15a1dee8d889b11,16#41036ad76d9888
+    ];
+
 reference_val(exsp) ->
     reference_val(exsplus);
 reference_val(exs1024s) ->
@@ -1517,6 +1545,33 @@ reference_jump_val(exsplus) ->
      12504080415362731, 45083100453836317, 270968267812126657, 93505647407734103,
      252852934678537969, 258758309277167202, 74250882143432077, 141629095984552833];
 
+reference_jump_val(exsss) ->
+    [16#304ae783d40db2b,16#1dfb196b3a5600a,16#2a24116effc6a0d,16#1f138d68c56725,
+     16#9360a445e2f989,16#32ed8080390e242,16#294ca85a270cff6,16#1418e6296a88bf,
+     16#114fae3dc578ba7,16#479c42c760eb72,16#334a40655df22d6,16#e7a85dd4d37d72,
+     16#181db16c8925c77,16#1b8a5a8afd16cbd,16#329107bf9777a39,16#2fc915c08535e42,
+     16#16696d142c6078,16#2e2a2601c919448,16#2246150d1000568,16#26109007cb3dd44,
+     16#3761360723e3175,16#169abd352db74de,16#1c97d520983684f,16#12455f0adee8c66,
+     16#46719cff00622d,16#1fc92792ed4e437,16#18e2edae21affb5,16#3a67fa9e3e7d46e,
+     16#1313fdc2728aa74,16#1c1a2b577581db8,16#db49357ea196b1,16#10e219a21d93fc7,
+     16#3c43abede083666,16#3eef5055a58bbf9,16#1975056f95d90e3,16#3916c133ab16d87,
+     16#2bc0bea891c26f1,16#391e4b369fc6b36,16#183f83155a359f6,16#1d9f137e9d2e488,
+     16#ef084de5f4cd3c,16#36a9cf7e29e55d3,16#19eca704e0409a7,16#1bdb99902896c69,
+     16#21777e2ad128203,16#5d0369ec0563e4,16#36db40b863bd74a,16#33feb71b7515159,
+     16#208d923ce26f257,16#3841b32891c082d,16#2748f224c2ba226,16#2fcd93b2daf79bb,
+     16#2c8e6cacad58ec4,16#39850131a1a85f,16#134648d6eea624d,16#2e102e197d5725c,
+     16#12ac280fa744758,16#1c18266c7442d16,16#22b5f91b15fe17e,16#316740ca870f7c8,
+     16#720ed4836c426,16#1aac0f738d04f8c,16#34fcd2a647b462c,16#3d430ac755114a3,
+     16#3692e3670fdf2a,16#265279ab0fc0a15,16#10bd883dee80945,16#10e7843413175e4,
+     16#b291deba08cee2,16#3915a8234caf11,16#34b911b96707dbd,16#ae63fcda15fde6,
+     16#b13b9091e82e41,16#29de1b6d70dc04f,16#23fbcbc409617e8,16#1389a0738061066,
+     16#360f39af790f5d1,16#f436da2a7d12f5,16#2d06ba8da21e08,16#3601a6492b887d,
+     16#2b2590b8c6cc186,16#f8d613b6904464,16#e5456786e46b78,16#201b8b1f96ed80c,
+     16#1b75b86d9b843f2,16#2e8bfaa7243a630,16#125ff068a78c3b4,16#3875a28c48bd26e,
+     16#f09a06941fc9d7,16#107c4de8ca77744,16#357c34144bb9ed6,16#3ccc55d3ebb3378,
+     16#28db7cea7d3fdee,16#3197fd0b49f6370,16#11af6fedb708ea6,16#2bde0382e37469e,
+     16#10666171abddb3f,16#1a8876c1f4e78a8,16#169c0efd4422043,16#1501c49abf0440f];
+
 reference_jump_val(exs1024) ->
     [2655961906500790629, 17003395417078685063, 10466831598958356428, 7603399148503548021,
      1650550950190587188, 12294992315080723704, 15743995773860389219, 5492181000145247327,