stdlib/random.ml

   1 (***********************************************************************)
   2 (*                                                                     *)
   3 (*                           Objective Caml                            *)
   4 (*                                                                     *)
   5 (*              Damien Doligez, projet Para, INRIA Rocquencourt        *)
   6 (*                                                                     *)
   7 (*  Copyright 1996 Institut National de Recherche en Informatique et   *)
   8 (*  en Automatique.  All rights reserved.  This file is distributed    *)
   9 (*  under the terms of the GNU Library General Public License, with    *)
  10 (*  the special exception on linking described in file ../LICENSE.     *)
  11 (*                                                                     *)
  12 (***********************************************************************)
  13
  14 (* $Id$ *)
  15
  16 (* "Linear feedback shift register" pseudo-random number generator. *)
  17 (* References: Robert Sedgewick, "Algorithms", Addison-Wesley *)
  18
  19 (* The PRNG is a linear feedback shift register.
  20    It is seeded by a MD5-based PRNG.
  21 *)
  22
  23 external random_seed: unit -> int = "caml_sys_random_seed";;
  24
  25 module State = struct
  26
  27   type t = { st : int array; mutable idx : int };;
  28
  29   let new_state () = { st = Array.make 55 0; idx = 0 };;
  30   let assign st1 st2 =
  31     Array.blit st2.st 0 st1.st 0 55;
  32     st1.idx <- st2.idx;
  33   ;;
  34
  35   let full_init s seed =
  36     let combine accu x = Digest.string (accu ^ string_of_int x) in
  37     let extract d =
  38       (Char.code d.[0] + (Char.code d.[1] lsl 8) + (Char.code d.[2] lsl 16))
  39       lxor (Char.code d.[3] lsl 22)
  40     in
  41     let l = Array.length seed in
  42     for i = 0 to 54 do
  43       s.st.(i) <- i;
  44     done;
  45     let accu = ref "x" in
  46     for i = 0 to 54 + max 55 l do
  47       let j = i mod 55 in
  48       let k = i mod l in
  49       accu := combine !accu seed.(k);
  50       s.st.(j) <- s.st.(j) lxor extract !accu;
  51     done;
  52     s.idx <- 0;
  53   ;;
  54
  55   let make seed =
  56     let result = new_state () in
  57     full_init result seed;
  58     result
  59   ;;
  60
  61   let make_self_init () = make [| random_seed () |];;
  62
  63   let copy s =
  64     let result = new_state () in
  65     assign result s;
  66     result
  67   ;;
  68
  69   (* Returns 30 random bits as an integer 0 <= x < 1073741824 *)
  70   let bits s =
  71     s.idx <- (s.idx + 1) mod 55;
  72     let newval = (s.st.((s.idx + 24) mod 55) + s.st.(s.idx)) land 0x3FFFFFFF in
  73     s.st.(s.idx) <- newval;
  74     newval
  75   ;;
  76
  77   let rec intaux s n =
  78     let r = bits s in
  79     let v = r mod n in
  80     if r - v > 0x3FFFFFFF - n + 1 then intaux s n else v
  81   ;;
  82   let int s bound =
  83     if bound > 0x3FFFFFFF || bound <= 0
  84     then invalid_arg "Random.int"
  85     else intaux s bound
  86   ;;
  87
  88   let rec int32aux s n =
  89     let b1 = Int32.of_int (bits s) in
  90     let b2 = Int32.shift_left (Int32.of_int (bits s land 1)) 30 in
  91     let r = Int32.logor b1 b2 in
  92     let v = Int32.rem r n in
  93     if Int32.sub r v > Int32.add (Int32.sub Int32.max_int n) 1l
  94     then int32aux s n
  95     else v
  96   ;;
  97   let int32 s bound =
  98     if bound <= 0l
  99     then invalid_arg "Random.int32"
 100     else int32aux s bound
 101   ;;
 102
 103   let rec int64aux s n =
 104     let b1 = Int64.of_int (bits s) in
 105     let b2 = Int64.shift_left (Int64.of_int (bits s)) 30 in
 106     let b3 = Int64.shift_left (Int64.of_int (bits s land 7)) 60 in
 107     let r = Int64.logor b1 (Int64.logor b2 b3) in
 108     let v = Int64.rem r n in
 109     if Int64.sub r v > Int64.add (Int64.sub Int64.max_int n) 1L
 110     then int64aux s n
 111     else v
 112   ;;
 113   let int64 s bound =
 114     if bound <= 0L
 115     then invalid_arg "Random.int64"
 116     else int64aux s bound
 117   ;;
 118
 119   let nativeint =
 120     if Nativeint.size = 32
 121     then fun s bound -> Nativeint.of_int32 (int32 s (Nativeint.to_int32 bound))
 122     else fun s bound -> Int64.to_nativeint (int64 s (Int64.of_nativeint bound))
 123   ;;
 124
 125   (* Returns a float 0 <= x < 1 with at most 90 bits of precision. *)
 126   let rawfloat s =
 127     let scale = 1073741824.0
 128     and r0 = Pervasives.float (bits s)
 129     and r1 = Pervasives.float (bits s)
 130     and r2 = Pervasives.float (bits s)
 131     in ((r0 /. scale +. r1) /. scale +. r2) /. scale
 132   ;;
 133
 134   let float s bound = rawfloat s *. bound;;
 135
 136   let bool s = (bits s land 1 = 0);;
 137
 138 end;;
 139
 140 (* This is the state you get with [init 27182818] on a 32-bit machine. *)
 141 let default = {
 142   State.st = [|
 143       509760043; 399328820; 99941072; 112282318; 611886020; 516451399;
 144       626288598; 337482183; 748548471; 808894867; 657927153; 386437385;
 145       42355480; 977713532; 311548488; 13857891; 307938721; 93724463;
 146       1041159001; 444711218; 1040610926; 233671814; 664494626; 1071756703;
 147       188709089; 420289414; 969883075; 513442196; 275039308; 918830973;
 148       598627151; 134083417; 823987070; 619204222; 81893604; 871834315;
 149       398384680; 475117924; 520153386; 324637501; 38588599; 435158812;
 150       168033706; 585877294; 328347186; 293179100; 671391820; 846150845;
 151       283985689; 502873302; 718642511; 938465128; 962756406; 107944131;
 152       192910970;
 153     |];
 154   State.idx = 0;
 155 };;
 156
 157 let bits () = State.bits default;;
 158 let int bound = State.int default bound;;
 159 let int32 bound = State.int32 default bound;;
 160 let nativeint bound = State.nativeint default bound;;
 161 let int64 bound = State.int64 default bound;;
 162 let float scale = State.float default scale;;
 163 let bool () = State.bool default;;
 164
 165 let full_init seed = State.full_init default seed;;
 166 let init seed = State.full_init default [| seed |];;
 167 let self_init () = init (random_seed());;
 168
 169 (* Manipulating the current state. *)
 170
 171 let get_state () = State.copy default;;
 172 let set_state s = State.assign default s;;
 173
 174 (********************
 175
 176 (* Test functions.  Not included in the library.
 177    The [chisquare] function should be called with n > 10r.
 178    It returns a triple (low, actual, high).
 179    If low <= actual <= high, the [g] function passed the test,
 180    otherwise it failed.
 181
 182   Some results:
 183
 184 init 27182818; chisquare int 100000 1000;;
 185 init 27182818; chisquare int 100000 100;;
 186 init 27182818; chisquare int 100000 5000;;
 187 init 27182818; chisquare int 1000000 1000;;
 188 init 27182818; chisquare int 100000 1024;;
 189 init 299792643; chisquare int 100000 1024;;
 190 init 14142136; chisquare int 100000 1024;;
 191 init 27182818; init_diff 1024; chisquare diff 100000 1024;;
 192 init 27182818; init_diff 100; chisquare diff 100000 100;;
 193 init 27182818; init_diff2 1024; chisquare diff2 100000 1024;;
 194 init 27182818; init_diff2 100; chisquare diff2 100000 100;;
 195 init 14142136; init_diff2 100; chisquare diff2 100000 100;;
 196 init 299792643; init_diff2 100; chisquare diff2 100000 100;;
 197 - : float * float * float = (936.754446796632465, 1032., 1063.24555320336754)
 198 # - : float * float * float = (80., 91.3699999999953434, 120.)
 199 # - : float * float * float = (4858.57864376269026, 4982., 5141.42135623730974)
 200 # - : float * float * float =
 201 (936.754446796632465, 1017.99399999994785, 1063.24555320336754)
 202 # - : float * float * float = (960., 984.565759999997681, 1088.)
 203 # - : float * float * float = (960., 1003.40735999999742, 1088.)
 204 # - : float * float * float = (960., 1035.23328000000038, 1088.)
 205 # - : float * float * float = (960., 1026.79551999999967, 1088.)
 206 # - : float * float * float = (80., 110.194000000003143, 120.)
 207 # - : float * float * float = (960., 1067.98080000000482, 1088.)
 208 # - : float * float * float = (80., 107.292000000001281, 120.)
 209 # - : float * float * float = (80., 85.1180000000022119, 120.)
 210 # - : float * float * float = (80., 86.614000000001397, 120.)
 211
 212 *)
 213
 214 (* Return the sum of the squares of v[i0,i1[ *)
 215 let rec sumsq v i0 i1 =
 216   if i0 >= i1 then 0.0
 217   else if i1 = i0 + 1 then Pervasives.float v.(i0) *. Pervasives.float v.(i0)
 218   else sumsq v i0 ((i0+i1)/2) +. sumsq v ((i0+i1)/2) i1
 219 ;;
 220
 221 let chisquare g n r =
 222   if n <= 10 * r then invalid_arg "chisquare";
 223   let f = Array.make r 0 in
 224   for i = 1 to n do
 225     let t = g r in
 226     f.(t) <- f.(t) + 1
 227   done;
 228   let t = sumsq f 0 r
 229   and r = Pervasives.float r
 230   and n = Pervasives.float n in
 231   let sr = 2.0 *. sqrt r in
 232   (r -. sr,   (r *. t /. n) -. n,   r +. sr)
 233 ;;
 234
 235 (* This is to test for linear dependencies between successive random numbers.
 236 *)
 237 let st = ref 0;;
 238 let init_diff r = st := int r;;
 239 let diff r =
 240   let x1 = !st
 241   and x2 = int r
 242   in
 243   st := x2;
 244   if x1 >= x2 then
 245     x1 - x2
 246   else
 247     r + x1 - x2
 248 ;;
 249
 250 let st1 = ref 0
 251 and st2 = ref 0
 252 ;;
 253
 254 (* This is to test for quadratic dependencies between successive random
 255    numbers.
 256 *)
 257 let init_diff2 r = st1 := int r; st2 := int r;;
 258 let diff2 r =
 259   let x1 = !st1
 260   and x2 = !st2
 261   and x3 = int r
 262   in
 263   st1 := x2;
 264   st2 := x3;
 265   (x3 - x2 - x2 + x1 + 2*r) mod r
 266 ;;
 267
 268 ********************)