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
|
(***********************************************************************)
(* *)
(* Objective Caml *)
(* *)
(* Damien Doligez, projet Para, INRIA Rocquencourt *)
(* *)
(* Copyright 1996 Institut National de Recherche en Informatique et *)
(* Automatique. Distributed only by permission. *)
(* *)
(***********************************************************************)
(* $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 j = ref 0
(* Returns 30 random bits as an integer 0 <= x < 1073741824 *)
let bits () =
j := (!j + 1) mod 55;
let newval =
Array.unsafe_get state ((!j+24) mod 55) + Array.unsafe_get state !j in
Array.unsafe_set state !j 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
(* Simple initialisation. The seed is an integer.
Two seeds that are close enough will not produce uncorrelated
pseudo-random sequences.
*)
let init seed =
let st = ref seed in
let mdg () =
st := !st + 1;
let d = Digest.string (string_of_int !st) 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
Array.unsafe_set state i (mdg ())
done;
j := 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
Array.unsafe_set state j
(Array.unsafe_get state j + Array.unsafe_get seed i)
done
(* Low-entropy system-dependent initialisation. *)
let self_init () = init (int_of_float (Sys.date () -. 2e9));;
(********************
(* 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
;;
********************)
|
