New branch binary-annot by Fabrice Le Fessant from r10789

git-svn-id: http://caml.inria.fr/svn/ocaml/branches/binary-annot@10790 f963ae5c-01c2-4b8c-9fe0-0dff7051ff02

1 files changed, 275 insertions, 0 deletions

diff --git a/ocaml-binary-annot/stdlib/random.ml b/ocaml-binary-annot/stdlib/random.ml
new file mode 100644
index 0000000000..bd3979914a
--- /dev/null
+++ b/ocaml-binary-annot/stdlib/random.ml
@@ -0,0 +1,275 @@
+(***********************************************************************)
+(*                                                                     *)
+(*                           Objective Caml                            *)
+(*                                                                     *)
+(*              Damien Doligez, projet Para, INRIA Rocquencourt        *)
+(*                                                                     *)
+(*  Copyright 1996 Institut National de Recherche en Informatique et   *)
+(*  en Automatique.  All rights reserved.  This file is distributed    *)
+(*  under the terms of the GNU Library General Public License, with    *)
+(*  the special exception on linking described in file ../LICENSE.     *)
+(*                                                                     *)
+(***********************************************************************)
+
+(* $Id$ *)
+
+(* Pseudo-random number generator
+   This is a lagged-Fibonacci F(55, 24, +) with a modified addition
+   function to enhance the mixing of bits.
+   If we use normal addition, the low-order bit fails tests 1 and 7
+   of the Diehard test suite, and bits 1 and 2 also fail test 7.
+   If we use multiplication as suggested by Marsaglia, it doesn't fare
+   much better.
+   By mixing the bits of one of the numbers before addition (XOR the
+   5 high-order bits into the low-order bits), we get a generator that
+   passes all the Diehard tests.
+*)
+
+external random_seed: unit -> int = "caml_sys_random_seed";;
+
+module State = struct
+
+  type t = { st : int array; mutable idx : int };;
+
+  let new_state () = { st = Array.make 55 0; idx = 0 };;
+  let assign st1 st2 =
+    Array.blit st2.st 0 st1.st 0 55;
+    st1.idx <- st2.idx;
+  ;;
+
+  let full_init s seed =
+    let combine accu x = Digest.string (accu ^ string_of_int x) in
+    let extract d =
+      Char.code d.[0] + (Char.code d.[1] lsl 8) + (Char.code d.[2] lsl 16)
+      + (Char.code d.[3] lsl 24)
+    in
+    let seed = if seed = [| |] then [| 0 |] else seed in
+    let l = Array.length seed in
+    for i = 0 to 54 do
+      s.st.(i) <- i;
+    done;
+    let accu = ref "x" in
+    for i = 0 to 54 + max 55 l do
+      let j = i mod 55 in
+      let k = i mod l in
+      accu := combine !accu seed.(k);
+      s.st.(j) <- s.st.(j) lxor extract !accu;
+    done;
+    s.idx <- 0;
+  ;;
+
+  let make seed =
+    let result = new_state () in
+    full_init result seed;
+    result
+  ;;
+
+  let make_self_init () = make [| random_seed () |];;
+
+  let copy s =
+    let result = new_state () in
+    assign result s;
+    result
+  ;;
+
+  (* Returns 30 random bits as an integer 0 <= x < 1073741824 *)
+  let bits s =
+    s.idx <- (s.idx + 1) mod 55;
+    let newval = s.st.((s.idx + 24) mod 55)
+                 + (s.st.(s.idx) lxor ((s.st.(s.idx) lsr 25) land 31)) in
+    s.st.(s.idx) <- newval;
+    newval land 0x3FFFFFFF   (* land is needed for 64-bit arch *)
+  ;;
+
+  let rec intaux s n =
+    let r = bits s in
+    let v = r mod n in
+    if r - v > 0x3FFFFFFF - n + 1 then intaux s n else v
+  ;;
+  let int s bound =
+    if bound > 0x3FFFFFFF || bound <= 0
+    then invalid_arg "Random.int"
+    else intaux s bound
+  ;;
+
+  let rec int32aux s n =
+    let b1 = Int32.of_int (bits s) in
+    let b2 = Int32.shift_left (Int32.of_int (bits s land 1)) 30 in
+    let r = Int32.logor b1 b2 in
+    let v = Int32.rem r n in
+    if Int32.sub r v > Int32.add (Int32.sub Int32.max_int n) 1l
+    then int32aux s n
+    else v
+  ;;
+  let int32 s bound =
+    if bound <= 0l
+    then invalid_arg "Random.int32"
+    else int32aux s bound
+  ;;
+
+  let rec int64aux s n =
+    let b1 = Int64.of_int (bits s) in
+    let b2 = Int64.shift_left (Int64.of_int (bits s)) 30 in
+    let b3 = Int64.shift_left (Int64.of_int (bits s land 7)) 60 in
+    let r = Int64.logor b1 (Int64.logor b2 b3) in
+    let v = Int64.rem r n in
+    if Int64.sub r v > Int64.add (Int64.sub Int64.max_int n) 1L
+    then int64aux s n
+    else v
+  ;;
+  let int64 s bound =
+    if bound <= 0L
+    then invalid_arg "Random.int64"
+    else int64aux s bound
+  ;;
+
+  let nativeint =
+    if Nativeint.size = 32
+    then fun s bound -> Nativeint.of_int32 (int32 s (Nativeint.to_int32 bound))
+    else fun s bound -> Int64.to_nativeint (int64 s (Int64.of_nativeint bound))
+  ;;
+
+  (* Returns a float 0 <= x < 1 with at most 90 bits of precision. *)
+  let rawfloat s =
+    let scale = 1073741824.0
+    and r0 = Pervasives.float (bits s)
+    and r1 = Pervasives.float (bits s)
+    and r2 = Pervasives.float (bits s)
+    in ((r0 /. scale +. r1) /. scale +. r2) /. scale
+  ;;
+
+  let float s bound = rawfloat s *. bound;;
+
+  let bool s = (bits s land 1 = 0);;
+
+end;;
+
+(* This is the state you get with [init 27182818]. *)
+let default = {
+  State.st = [|
+      0x7ae2522b; 0x5d8d4634; 0x15b4fad0; 0x18b14ace; 0x12f8a3c4; 0x7b086c47;
+      0x16d467d6; 0x501d91c7; 0x321df177; 0x4176c193; 0x1ff72bf1; 0x5e889109;
+      0x0b464b18; 0x6b86b97c; 0x4891da48; 0x03137463; 0x485ac5a1; 0x15d61f2f;
+      0x7bced359; 0x69c1c132; 0x7a86766e; 0x366d8c86; 0x1f5b6222; 0x7ce1b59f;
+      0x2ebf78e1; 0x67cd1b86; 0x658f3dc3; 0x789a8194; 0x42e4c44c; 0x58c43f7d;
+      0x0f6e534f; 0x1e7df359; 0x455d0b7e; 0x10e84e7e; 0x126198e4; 0x4e7722cb;
+      0x5cbede28; 0x7391b964; 0x7d40e92a; 0x4c59933d; 0x0b8cd0b7; 0x64efff1c;
+      0x2803fdaa; 0x08ebc72e; 0x4f522e32; 0x45398edc; 0x2144a04c; 0x4aef3cbd;
+      0x41ad4719; 0x75b93cd6; 0x2a559d4f; 0x5e6fd768; 0x66e27f36; 0x186f18c3;
+      0x2fbf967a;
+    |];
+  State.idx = 0;
+};;
+
+let bits () = State.bits default;;
+let int bound = State.int default bound;;
+let int32 bound = State.int32 default bound;;
+let nativeint bound = State.nativeint default bound;;
+let int64 bound = State.int64 default bound;;
+let float scale = State.float default scale;;
+let bool () = State.bool default;;
+
+let full_init seed = State.full_init default seed;;
+let init seed = State.full_init default [| seed |];;
+let self_init () = init (random_seed());;
+
+(* Manipulating the current state. *)
+
+let get_state () = State.copy default;;
+let set_state s = State.assign default s;;
+
+(********************
+
+(* Test functions.  Not included in the library.
+   The [chisquare] function should be called with n > 10r.
+   It returns a triple (low, actual, high).
+   If low <= actual <= high, the [g] function passed the test,
+   otherwise it failed.
+
+  Some results:
+
+init 27182818; chisquare int 100000 1000;;
+init 27182818; chisquare int 100000 100;;
+init 27182818; chisquare int 100000 5000;;
+init 27182818; chisquare int 1000000 1000;;
+init 27182818; chisquare int 100000 1024;;
+init 299792643; chisquare int 100000 1024;;
+init 14142136; chisquare int 100000 1024;;
+init 27182818; init_diff 1024; chisquare diff 100000 1024;;
+init 27182818; init_diff 100; chisquare diff 100000 100;;
+init 27182818; init_diff2 1024; chisquare diff2 100000 1024;;
+init 27182818; init_diff2 100; chisquare diff2 100000 100;;
+init 14142136; init_diff2 100; chisquare diff2 100000 100;;
+init 299792643; init_diff2 100; chisquare diff2 100000 100;;
+- : float * float * float = (936.754446796632465, 997.5, 1063.24555320336754)
+# - : float * float * float = (80., 89.7400000000052387, 120.)
+# - : float * float * float = (4858.57864376269, 5045.5, 5141.42135623731)
+# - : float * float * float =
+(936.754446796632465, 944.805999999982305, 1063.24555320336754)
+# - : float * float * float = (960., 1019.19744000000355, 1088.)
+# - : float * float * float = (960., 1059.31776000000536, 1088.)
+# - : float * float * float = (960., 1039.98463999999512, 1088.)
+# - : float * float * float = (960., 1054.38207999999577, 1088.)
+# - : float * float * float = (80., 90.096000000005, 120.)
+# - : float * float * float = (960., 1076.78720000000612, 1088.)
+# - : float * float * float = (80., 85.1760000000067521, 120.)
+# - : float * float * float = (80., 85.2160000000003492, 120.)
+# - : float * float * float = (80., 80.6220000000030268, 120.)
+
+*)
+
+(* Return the sum of the squares of v[i0,i1[ *)
+let rec sumsq v i0 i1 =
+  if i0 >= i1 then 0.0
+  else if i1 = i0 + 1 then Pervasives.float v.(i0) *. Pervasives.float v.(i0)
+  else sumsq v i0 ((i0+i1)/2) +. sumsq v ((i0+i1)/2) i1
+;;
+
+let chisquare g n r =
+  if n <= 10 * r then invalid_arg "chisquare";
+  let f = Array.make r 0 in
+  for i = 1 to n do
+    let t = g r in
+    f.(t) <- f.(t) + 1
+  done;
+  let t = sumsq f 0 r
+  and r = Pervasives.float r
+  and n = Pervasives.float n in
+  let sr = 2.0 *. sqrt r in
+  (r -. sr,   (r *. t /. n) -. n,   r +. sr)
+;;
+
+(* This is to test for linear dependencies between successive random numbers.
+*)
+let st = ref 0;;
+let init_diff r = st := int r;;
+let diff r =
+  let x1 = !st
+  and x2 = int r
+  in
+  st := x2;
+  if x1 >= x2 then
+    x1 - x2
+  else
+    r + x1 - x2
+;;
+
+let st1 = ref 0
+and st2 = ref 0
+;;
+
+(* This is to test for quadratic dependencies between successive random
+   numbers.
+*)
+let init_diff2 r = st1 := int r; st2 := int r;;
+let diff2 r =
+  let x1 = !st1
+  and x2 = !st2
+  and x3 = int r
+  in
+  st1 := x2;
+  st2 := x3;
+  (x3 - x2 - x2 + x1 + 2*r) mod r
+;;
+
+********************)