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 z = x[i] + y[i] + carry;
+    carry = carry ? (z <= x[i]) : (z < x[i]);
+    x[i] = z;
+  }
+}
+
+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 4 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; 4L; 2L; 3L |] 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