Module Random: revised initialization of PRNG from array of integers

Sometimes we are given a single integer (as in Random.init) and
sometimes we are given an array of 12 bytes (as in Random.self_init
with the /dev/urandom implementation).

In the first case, from a single integer we need to come up with 4
values for the 4 components of the PRNG state, avoiding bad values
like 0, 0 for the x component.

In the second case, we need to collect the 96 bits of entropy spread
among these 12 bytes, then distribute them on the 4 components of the
PRNG state.

This commit treats the array as a string of 64-bit characters and applies
a hash function to this string, producing a 256-bit hash, which is then
used as the initial PRNG state.

The hash function used in FNV1a, because it supports 256-bit outputs
and it is relatively easy to implement.

runtime/prng.c

@@ -13,6 +13,7 @@
 /*                                                                        */
 /**************************************************************************/
 
+#include <string.h>
 #include "caml/alloc.h"
 #include "caml/bigarray.h"
 #include "caml/mlvalues.h"
@@ -34,7 +35,7 @@ struct LXM_state {
 
 #define LXM_val(v) ((struct LXM_state *) Caml_ba_data_val(v))
 
-static inline uint64_t rotl(const uint64_t x, int k) {
+Caml_inline uint64_t rotl(const uint64_t x, int k) {
   return (x << k) | (x >> (64 - k));
 }
 
@@ -77,17 +78,58 @@ static void caml_lxm_set(value v, uint64_t i1, uint64_t i2,
   st->s = i4;
 }
 
+static void add256(uint64_t x[4], uint64_t y[4])
+{
+  int i;
+  unsigned int carry = 0;
+  for (i = 0; i < 4; i++) {
+    uint64_t t1 = x[i];
+    uint64_t t2 = t1 + y[i];
+    uint64_t t3 = t2 + carry;
+    x[i] = t3;
+    carry = (t2 < t1) + (t3 < t2);
+  }
+}
+
+static void shl256(uint64_t x[4], int amount)
+{
+  while (amount >= 64) {
+    x[3] = x[2]; x[2] = x[1]; x[1] = x[0]; x[0] = 0;
+    amount -= 64;
+  }
+  if (amount == 0) return;
+  x[3] = x[3] << amount | x[2] >> (64 - amount);
+  x[2] = x[2] << amount | x[1] >> (64 - amount);
+  x[1] = x[1] << amount | x[0] >> (64 - amount);
+  x[0] = x[0] << amount;
+}
+
 CAMLprim value caml_lxm_init(value v, value a)
 {
-  const uint64_t mix = 6364136223846793005;
-  /* Multiplier taken from the MMIX LCG, Knoth TAOCP vol 2, 1998 edition */
-  uint64_t d[4] = {0, 0, 0, 0};
+  /* The FNV1-256 offset basis,
+     1000292579580525809070709686206257048370927960
+     14241193945225284501741471925557,
+     as four 64-bit digits, little endian */
+  uint64_t h[4] = { 0x1023b4c8caee0535,
+                    0xc8b1536847b6bbb3,
+                    0x2d98c384c4e576cc,
+                    0xdd268dbcaac55036 };
+  uint64_t t[4];
   mlsize_t i, len;
 
   for (i = 0, len = Wosize_val(a); i < len; i++) {
-    d[i % 4] = d[i % 4] * mix + Long_val(Field(a, i));
+    /* On 32-bit hosts, force sign-extension to 64 bits, so that
+       the results are the same on 32-bit and 64-bit platforms */
+    h[0] ^= (int64_t) Long_val(Field(a, i));
+    /* Multiply by the FNV1-256 prime, 2^168 + 2^8 + 2^6 + 2^5 + 2^1 + 2^0 */
+    memcpy(t, h, sizeof(t));
+    shl256(t, 1-0); add256(h, t);
+    shl256(t, 5-1); add256(h, t);
+    shl256(t, 6-5); add256(h, t);
+    shl256(t, 8-6); add256(h, t);
+    shl256(t, 168-8); add256(h, t);
   }
-  caml_lxm_set(v, d[0], d[1], d[2], d[3]);
+  caml_lxm_set(v, h[0], h[1], h[2], h[3]);
   return Val_unit;
 }
 

testsuite/tests/lib-dynlink-native/main.reference

@@ -6,7 +6,7 @@ Registering module Plugin2
 1
 2
 6
-6
+1
 XXX
 Loading plugin_thread.so
 Registering module Plugin_thread

testsuite/tests/lib-random/defaultinit.reference

@@ -1,4 +1,4 @@
- 175 168 780 564 484 867 123 178 768 467 576 542 797 603 992 635 874 323 778 676 422
+ 192 796 468 294 544 121 972 623 517 515 878 696 2 285 379 415 680 7 351 147 607
 
- 669.156953125 242.875937515 704.511016963 148.823985021 483.099782622 777.74287 771.413854095 894.790411 761.141039066 921.111045515 984.657815283 468.838002847 891.446431758 87.8734596194 708.49973828 110.822253298 484.719704342 157.513223036 714.773419366 59.3660539904 701.591141363
+ 511.211801957 844.958249829 499.607291677 491.690183172 363.251949921 511.447587009 742.101717446 159.235602873 217.702834638 93.2460661129 791.464061354 58.9841280761 839.865622773 681.503997651 738.538364284 898.529623312 190.903176011 431.537064919 270.998675184 430.84362776 424.00611019
 All tests succeeded.

testsuite/tests/lib-random/testvectors.ml

@@ -4,12 +4,14 @@
 (* Check the numbers drawn from a known state against the numbers
    obtained from the reference Java implementation. *)
 
-open Random
+open Bigarray
 
 let _ =
-  full_init [| 1; 2; 3; 4 |];
+  let a = Array1.of_array Int64 C_layout [| 1L; 2L; 3L; 4L |] in
+  (* Violate abstraction of type Random.State.t to manipulate state directly *)
+  let r = (Obj.magic a : Random.State.t) in
   for i = 0 to 49 do
-    Printf.printf "%Ld\n" (bits64 ())
+    Printf.printf "%Ld\n" (Random.State.bits64 r)
   done
 
 let _ = exit 0

testsuite/tests/lib-random/testvectors.reference

@@ -1,50 +1,50 @@
-1125592859273023139
--8555116859351185690
-1740974640465038003
--114103598616990594
--4537459191886299567
-6291660581237232798
-4399537237262424059
-1444832433108046408
--885086767605962159
-5035407203018967627
-2454027010810206102
--3730687192620176031
--9121837924364225174
-5158779247817520871
-346180386384953991
--5994913815963856552
-6869762371190108983
--351309903469929100
--4064370530550403097
--1705171716212891727
--2597460413620866941
-3726722190751020339
--2161674832373023869
-1104736721682310151
--8637978141014867568
-563925799384532606
--8148729055587940695
--2600730985292293862
-8948433074316439393
-227565377938431895
--8692872471483886118
-4158043793302624760
--1494134816437625629
-3067733291085251934
--2202687045387381141
-5220625884742653315
--3801895474121595209
-8988626187385988242
--25423021318571012
--1519542409018449719
-8318221973461969899
-8866230809488490189
-8260650303702392390
--1241814289009151420
-4356707657081988900
--855677221024118431
--5700738296854965036
--2490506808380844303
--7083650095061909347
--6592815571846300453
+3860816457867857678
+21223322560256856
+3500884966496404595
+1500794501795755166
+-4410813772481188553
+-4572900151748500348
+4265851130944421541
+7572018898959715966
+-2937174041541593693
+-1802830867415637455
+8366519798672692173
+-6590954168183217171
+-4921483866549552021
+-3689136988381952376
+1051817605112975897
+-7328239262545604981
+5123201494011052613
+-2341724972498217376
+9115420906531662800
+-2254527559866817705
+6561621740404805009
+-3304584699295016645
+-6306474873117248843
+-8254191423720223116
+7918980772536573383
+3812821519051744912
+3169704758133872700
+1303179780519235243
+2791474158728480712
+5710006355063646940
+-4119479146064124870
+4742660959001555540
+-1339660087226824925
+-3233387961267533196
+-5600632561554148951
+6742170938733010945
+6703547594041469408
+6019528866072981890
+1653175532392381808
+1251959135604991018
+-3454440012186344425
+-6919427552847598775
+-9047964058899192553
+-4707136377285073135
+-5588355652467780140
+1107195072897197378
+-3891759467946419528
+97248838865565210
+-5756171654209134086
+1043655509106157291