Simple chi-square test for functions from module Random

This test appeared during review of PR #9489, but seems useful
independently of this PR.

testsuite/tests/lib-random/chi2.ml

@@ -0,0 +1,68 @@
+(* TEST
+*)
+
+(* A basic chi-square test to detect simple errors in the Random module. *)
+
+(* Accumulate [n] samples from function [f] and check the chi-square.
+   Only the low 8 bits of the result of [f] are sampled. *)
+
+let chisquare n f =
+  let r = 256 in
+  let freq = Array.make r 0 in
+  for i = 0 to n - 1 do
+    let t = f () land 0xFF in
+    freq.(t) <- freq.(t) + 1
+  done;
+  let expected = float n /. float r in
+  let t =
+    Array.fold_left
+      (fun s x -> let d = float x -. expected in d *. d +. s)
+      0.0 freq in
+  let chi2 = t /. expected in
+  let degfree = float r -. 1.0 in
+  (* The degree of freedom is high, so we approximate as a normal
+     distribution with mean equal to degfree and variance 2 * degfree.
+     Four sigmas correspond to a 99.9968% confidence interval.
+     (Without the approximation, the confidence interval seems to be 99.986%.)
+  *)
+  chi2 <= degfree +. 4.0 *. sqrt (2.0 *. degfree)
+
+let test name f =
+  if not (chisquare 100_000 f)
+  then Printf.printf "%s: suspicious result\n%!" name
+
+let _ =
+  test "Random.bits (bits 0-7)"
+       Random.bits;
+  test "Random.bits (bits 12-19)"
+       (fun () -> Random.bits() lsr 12);
+  test "Random.bits (bits 22-29)"
+       (fun () -> Random.bits() lsr 22);
+  test "Random.int 2^26 (bits 0-7)"
+       (fun () -> Random.int (1 lsl 26));
+  test "Random.int 2^26 (bits 18-25)"
+       (fun () -> Random.int (1 lsl 26) lsr 18);
+  test "Random.int (256 * p) / p"
+       (fun () -> Random.int (256 * 853187) / 853187);
+  test "Random.float 1.0 (first 8 bits)"
+       (fun () -> int_of_float (Random.float 1.0 *. 256.0));
+  test "Random.float 1.0 (next 8 bits)"
+       (fun () -> int_of_float (Random.float 1.0 *. 65536.0));
+  test "Random.int32 2^30 (bits 0-7)"
+       (fun () -> Int32.to_int (Random.int32 0x40000000l));
+  test "Random.int32 2^30 (bits 20-27)"
+       (fun () -> Int32.(to_int (shift_right (Random.int32 0x40000000l) 20)));
+  test "Random.int32 (256 * p) / p"
+       (let p = 7048673l in
+        fun () -> Int32.(to_int (div (Random.int32 (mul 256l p)) p)));
+  test "Random.int64 2^60 (bits 0-7)"
+       (fun () -> Int64.to_int (Random.int64 0x1000000000000000L));
+  test "Random.int64 2^60 (bits 30-37)"
+       (fun () -> Int64.(to_int (shift_right (Random.int64 0x1000000000000000L)
+                                             30)));
+  test "Random.int64 2^60 (bits 52-59)"
+       (fun () -> Int64.(to_int (shift_right (Random.int64 0x1000000000000000L)
+                                             52)));
+  test "Random.int64 (256 * p) / p"
+       (let p = 16430454264262693L in
+        fun () -> Int64.(to_int (div (Random.int64 (mul 256L p)) p)))