1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
|
(***********************************************************************)
(* *)
(* 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$ *)
(* "Linear feedback shift register" 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
|]
let index = ref 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;
newval land 0x3FFFFFFF
(* Returns a float 0 <= x < 1 with at most 90 bits of precision. *)
let rawfloat () =
let scale = 1073741824.0
and r0 = float (bits ())
and r1 = float (bits ())
and r2 = float (bits ())
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 =
if bound > 0x3FFFFFFF || bound <= 0
then invalid_arg "Random.int"
else (intaux (0x3FFFFFFF / bound * bound)) mod bound
let float bound = rawfloat () *. bound
let bool () = (bits () 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
(Char.code d.[0] + (Char.code d.[1] lsl 8) + (Char.code d.[2] lsl 16))
lxor (Char.code d.[3] lsl 22)
in
for i = 0 to 54 do
state.(i) <- (mdg 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 j = i mod 55 in
state.(j) <- state.(j) + seed.(i)
done
(* Low-entropy system-dependent initialisation. *)
external random_seed: unit -> int = "sys_random_seed";;
let self_init () = init (random_seed());;
(* 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;
;;
(********************
(* 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:
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
*)
(* 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 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 = float r
and n = 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 := Random.int r;;
let diff r =
let x1 = !st
and x2 = Random.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 := Random.int r; st2 := Random.int r;;
let diff2 r =
let x1 = !st1
and x2 = !st2
and x3 = Random.int r
in
st1 := x2;
st2 := x3;
(x3 - x2 - x2 + x1 + 2*r) mod r
;;
********************)
|
