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
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
|
(***********************************************************************)
(* *)
(* 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" 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 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 s_bits s =
s.idx <- (s.idx + 1) mod 55;
let newval = (s.st.((s.idx + 24) mod 55) + s.st.(s.idx)) land 0x3FFFFFFF in
s.st.(s.idx) <- newval;
newval
;;
(* Returns a float 0 <= x < 1 with at most 90 bits of precision. *)
let s_rawfloat s =
let scale = 1073741824.0
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 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 (s_intaux s (0x3FFFFFFF / bound * bound)) mod bound
;;
let s_float s bound = s_rawfloat s *. bound
let s_bool s = (s_bits s land 1 = 0);;
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
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 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. *)
let get_state () = s_copy default;;
let set_state s = assign_state 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, 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 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
;;
********************)
|
