changement seeding + ajout fonctions serieuses
git-svn-id: http://caml.inria.fr/svn/ocaml/trunk@5592 f963ae5c-01c2-4b8c-9fe0-0dff7051ff02master
parent
0b2e6f5f7c
commit
60a5460849
1
Changes
1
Changes
|
@ -55,6 +55,7 @@ Libraries:
|
|||
- Module Arg: added Set_*, Symbol.
|
||||
- Module Scanf: %n and %N formats to count characters / items read so far;
|
||||
assorted bug fixes.
|
||||
- Module Random: better seeding; added functions for serious use.
|
||||
|
||||
Runtime system:
|
||||
- output_value/input_value: fixed bug with large blocks (>= 4 Mwords)
|
||||
|
|
204
stdlib/random.ml
204
stdlib/random.ml
|
@ -13,93 +13,124 @@
|
|||
|
||||
(* $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.
|
||||
*)
|
||||
|
||||
(* 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
|
||||
|]
|
||||
type state = { st : int array; mutable idx : int };;
|
||||
|
||||
let index = ref 0
|
||||
(* This is the state you get with [init 27182818] on a 32-bit machine. *)
|
||||
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;
|
||||
;;
|
||||
|
||||
let full_init seed = s_full_init default 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 ();;
|
||||
|
||||
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;
|
||||
;;
|
||||
|
||||
(* Create, initialise, and return a new state value. *)
|
||||
let s_make seed =
|
||||
let result = new_state () in
|
||||
s_full_init result seed;
|
||||
result
|
||||
;;
|
||||
|
||||
let s_copy s =
|
||||
let result = new_state () in
|
||||
assign_state result s;
|
||||
result
|
||||
;;
|
||||
|
||||
(* Manipulating the current state. *)
|
||||
|
||||
type state = { st : int array; idx : int };;
|
||||
|
||||
let get_state () = { st = Array.copy state; idx = !index };;
|
||||
|
||||
let set_state s =
|
||||
Array.blit s.st 0 state 0 55;
|
||||
index := s.idx;
|
||||
;;
|
||||
let get_state () = s_copy default;;
|
||||
let set_state s = assign_state default s;;
|
||||
|
||||
(********************
|
||||
|
||||
|
@ -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;
|
||||
|
|
|
@ -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.
|
||||
*)
|
||||
|
|
Loading…
Reference in New Issue