diff options
Diffstat (limited to 'stdlib')
| -rw-r--r-- | stdlib/random.ml | 202 | ||||
| -rw-r--r-- | stdlib/random.mli | 38 |
2 files changed, 149 insertions, 91 deletions
diff --git a/stdlib/random.ml b/stdlib/random.ml index 35d19ff0c..2da5e98f0 100644 --- a/stdlib/random.ml +++ b/stdlib/random.ml @@ -13,94 +13,125 @@ (* $Id$ *) -(* "Linear feedback shift register" random number generator. *) +(* "Linear feedback shift register" pseudo-random number generator. *) (* References: Robert Sedgewick, "Algorithms", Addison-Wesley *) (* The PRNG is a linear feedback shift register. It is seeded by a MD5-based PRNG. *) +type state = { st : int array; mutable idx : int };; + (* This is the state you get with [init 27182818] on a 32-bit machine. *) -let state = [| - 561073064; 1051173471; 764306064; 9858203; 1023641486; 615350359; - 552627506; 486882977; 147054819; 951240904; 869261341; 71648846; 848741663; - 337696531; 66770770; 473370118; 998499212; 477485839; 814302728; 281896889; - 206134737; 796925167; 762624501; 971004788; 878960411; 233350272; - 965168955; 933858406; 572927557; 708896334; 32881167; 462134267; 868098973; - 768795410; 567327260; 4136554; 268309077; 804670393; 854580894; 781847598; - 310632349; 22990936; 187230644; 714526560; 146577263; 979459837; 514922558; - 414383108; 21528564; 896816596; 33747835; 180326017; 414576093; 124177607; - 440266690 -|] - -let index = ref 0 +let default = { + st = [| + 561073064; 1051173471; 764306064; 9858203; 1023641486; 615350359; + 552627506; 486882977; 147054819; 951240904; 869261341; 71648846; + 848741663; 337696531; 66770770; 473370118; 998499212; 477485839; + 814302728; 281896889; 206134737; 796925167; 762624501; 971004788; + 878960411; 233350272; 965168955; 933858406; 572927557; 708896334; + 32881167; 462134267; 868098973; 768795410; 567327260; 4136554; + 268309077; 804670393; 854580894; 781847598; 310632349; 22990936; + 187230644; 714526560; 146577263; 979459837; 514922558; 414383108; + 21528564; 896816596; 33747835; 180326017; 414576093; 124177607; + 440266690; + |]; + idx = 0; +};; (* Returns 30 random bits as an integer 0 <= x < 1073741824 *) -let bits () = - index := (!index + 1) mod 55; - let newval = - state.((!index + 24) mod 55) + state.(!index) in - state.(!index) <- newval; +let s_bits s = + s.idx <- (s.idx + 1) mod 55; + let newval = s.st.((s.idx + 24) mod 55) + s.st.(s.idx) in + s.st.(s.idx) <- newval; newval land 0x3FFFFFFF +;; (* Returns a float 0 <= x < 1 with at most 90 bits of precision. *) -let rawfloat () = +let s_rawfloat s = let scale = 1073741824.0 - and r0 = float (bits ()) - and r1 = float (bits ()) - and r2 = float (bits ()) + and r0 = Pervasives.float (s_bits s) + and r1 = Pervasives.float (s_bits s) + and r2 = Pervasives.float (s_bits s) in ((r0 /. scale +. r1) /. scale +. r2) /. scale +;; -let rec intaux n = - let r = bits () in - if r >= n then intaux n else r -let int bound = +let rec s_intaux s n = + let r = s_bits s in + if r >= n then s_intaux s n else r +;; +let s_int s bound = if bound > 0x3FFFFFFF || bound <= 0 then invalid_arg "Random.int" - else (intaux (0x3FFFFFFF / bound * bound)) mod bound + else (s_intaux s (0x3FFFFFFF / bound * bound)) mod bound +;; -let float bound = rawfloat () *. bound +let s_float s bound = s_rawfloat s *. bound -let bool () = (bits () land 1 = 0);; +let s_bool s = (s_bits s land 1 = 0);; -(* Simple initialisation. The seed is an integer. *) -let init seed = - let mdg i = - let d = Digest.string (string_of_int i ^ string_of_int seed) in +let bits () = s_bits default;; +let int bound = s_int default bound;; +let float scale = s_float default scale;; +let bool () = s_bool default;; + +(* Full initialisation. The seed is an array of integers. *) +let s_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)) lxor (Char.code d.[3] lsl 22) in + let l = Array.length seed in for i = 0 to 54 do - state.(i) <- (mdg i) + s.st.(i) <- i; done; - index := 0 - -(* Full initialisation. The seed is an array of integers. *) -let full_init seed = - init 27182818; - for i = 0 to Array.length (seed) - 1 do + let accu = ref "x" in + for i = 0 to 54 + max 55 l do let j = i mod 55 in - state.(j) <- state.(j) + seed.(i) - done + let k = i mod l in + accu := combine !accu seed.(k); + s.st.(j) <- s.st.(j) lxor extract !accu; + done; + s.idx <- 0; +;; -(* Low-entropy system-dependent initialisation. *) +let full_init seed = s_full_init default seed;; -external random_seed: unit -> int = "sys_random_seed";; +(* Simple initialisation. The seed is an integer. *) +let init seed = s_full_init default [| seed |];; +(* Low-entropy system-dependent initialisation. *) +external random_seed: unit -> int = "sys_random_seed";; let self_init () = init (random_seed());; +(* The default PRNG is initialised with self_init. *) +self_init ();; -(* Manipulating the current state. *) - -type state = { st : int array; idx : int };; +let new_state () = { st = Array.make 55 0; idx = 0 };; +let assign_state st1 st2 = + Array.blit st2.st 0 st1.st 0 55; + st1.idx <- st2.idx; +;; -let get_state () = { st = Array.copy state; idx = !index };; +(* Create, initialise, and return a new state value. *) +let s_make seed = + let result = new_state () in + s_full_init result seed; + result +;; -let set_state s = - Array.blit s.st 0 state 0 55; - index := s.idx; +let s_copy s = + let result = new_state () in + assign_state result s; + result ;; +(* Manipulating the current state. *) + +let get_state () = s_copy default;; +let set_state s = assign_state default s;; + (******************** (* Test functions. Not included in the library. @@ -111,39 +142,40 @@ let set_state s = Some results: -Random.init 27182818; chisquare Random.int 100000 1000;; -Random.init 27182818; chisquare Random.int 100000 100;; -Random.init 27182818; chisquare Random.int 100000 5000;; -Random.init 27182818; chisquare Random.int 1000000 1000;; -Random.init 27182818; chisquare Random.int 100000 1024;; -Random.init 299792643; chisquare Random.int 100000 1024;; -Random.init 14142136; chisquare Random.int 100000 1024;; -Random.init 27182818; init_diff 1024; chisquare diff 100000 1024;; -Random.init 27182818; init_diff 100; chisquare diff 100000 100;; -Random.init 27182818; init_diff2 1024; chisquare diff2 100000 1024;; -Random.init 27182818; init_diff2 100; chisquare diff2 100000 100;; -Random.init 14142136; init_diff2 100; chisquare diff2 100000 100;; -Random.init 299792643; init_diff2 100; chisquare diff2 100000 100;; -- : float * float * float = 936.754446797, 948.8, 1063.2455532 -#- : float * float * float = 80, 80.076, 120 -#- : float * float * float = 4858.57864376, 4767.5, 5141.42135624 ********* -#- : float * float * float = 936.754446797, 951.2, 1063.2455532 -#- : float * float * float = 960, 1028.31104, 1088 -#- : float * float * float = 960, 1012.64384, 1088 -#- : float * float * float = 960, 970.25024, 1088 -#- : float * float * float = 960, 982.29248, 1088 -#- : float * float * float = 80, 110.418, 120 -#- : float * float * float = 960, 1022.76096, 1088 -#- : float * float * float = 80, 96.894, 120 -#- : float * float * float = 80, 83.864, 120 -#- : float * float * float = 80, 89.956, 120 +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, 1032., 1063.24555320336754) +# - : float * float * float = (80., 91.3699999999953434, 120.) +# - : float * float * float = (4858.57864376269026, 4982., 5141.42135623730974) +# - : float * float * float = +(936.754446796632465, 1017.99399999994785, 1063.24555320336754) +# - : float * float * float = (960., 984.565759999997681, 1088.) +# - : float * float * float = (960., 1003.40735999999742, 1088.) +# - : float * float * float = (960., 1035.23328000000038, 1088.) +# - : float * float * float = (960., 1026.79551999999967, 1088.) +# - : float * float * float = (80., 110.194000000003143, 120.) +# - : float * float * float = (960., 1067.98080000000482, 1088.) +# - : float * float * float = (80., 107.292000000001281, 120.) +# - : float * float * float = (80., 85.1180000000022119, 120.) +# - : float * float * float = (80., 86.614000000001397, 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 float v.(i0) *. float v.(i0) + 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 ;; @@ -155,8 +187,8 @@ let chisquare g n r = f.(t) <- f.(t) + 1 done; let t = sumsq f 0 r - and r = float r - and n = float n in + 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) ;; @@ -164,10 +196,10 @@ let chisquare g n r = (* This is to test for linear dependencies between successive random numbers. *) let st = ref 0;; -let init_diff r = st := Random.int r;; +let init_diff r = st := int r;; let diff r = let x1 = !st - and x2 = Random.int r + and x2 = int r in st := x2; if x1 >= x2 then @@ -183,11 +215,11 @@ and st2 = ref 0 (* This is to test for quadratic dependencies between successive random numbers. *) -let init_diff2 r = st1 := Random.int r; st2 := Random.int r;; +let init_diff2 r = st1 := int r; st2 := int r;; let diff2 r = let x1 = !st1 and x2 = !st2 - and x3 = Random.int r + and x3 = int r in st1 := x2; st2 := x3; diff --git a/stdlib/random.mli b/stdlib/random.mli index 32854a659..0f1ccc459 100644 --- a/stdlib/random.mli +++ b/stdlib/random.mli @@ -13,7 +13,9 @@ (* $Id$ *) -(** Pseudo-random number generator (PRNG). *) +(** Pseudo-random number generators (PRNG). *) + +(** {6 functions for casual users} *) val init : int -> unit (** Initialize the generator, using the argument as a seed. @@ -24,7 +26,8 @@ val full_init : int array -> unit val self_init : unit -> unit (** Initialize the generator with a more-or-less random seed chosen - in a system-dependent way. *) + in a system-dependent way. The generator is initialised with this + function at the start of the program. *) val bits : unit -> int (** Return 30 random bits in a nonnegative integer. *) @@ -37,8 +40,8 @@ val int : int -> int val float : float -> float (** [Random.float bound] returns a random floating-point number between 0 (inclusive) and [bound] (exclusive). If [bound] is - negative, the result is negative. If [bound] is 0, the result - is 0. *) + negative, the result is negative or zero. If [bound] is 0, + the result is 0. *) val bool : unit -> bool (** [Random.bool ()] returns [true] or [false] with probability 0.5 each. *) @@ -48,10 +51,33 @@ type state generator. *) val get_state : unit -> state -(** Returns the current state of the generator. This is useful for +(** Return the current state of the generator. This is useful for checkpointing computations that use the PRNG. *) val set_state : state -> unit -(** Resets the state of the generator to some previous state returned by +(** Reset the state of the generator to some previous state returned by {!Random.get_state}. *) + +(** {6 functions for serious users} *) + +(** These function manipulate the current state explicitely. + This allows you to use one or several deterministic PRNGs, + even in a multi-threaded program, without interference from + other parts of the program (for example, the Filename module + and some object-oriented primitives use the default PRNG). +*) + +val s_make : int array -> state;; +(** Create a new state and initialize it with the given seed. *) + +val s_copy : state -> state;; +(** Make a copy of the given state. *) + +val s_bits : state -> int;; +val s_int : state -> int -> int;; +val s_float : state -> float -> float;; +val s_bool : state -> bool;; +(** These functions are the same as the above versions, except that they + use (and update) the given PRNG state instead of the default one. +*) |