22528
|
1 |
(* ID: $Id$
|
|
2 |
Author: Florian Haftmann, TU Muenchen
|
|
3 |
*)
|
|
4 |
|
26265
|
5 |
header {* A HOL random engine *}
|
22528
|
6 |
|
|
7 |
theory Random
|
26265
|
8 |
imports State_Monad Code_Index
|
22528
|
9 |
begin
|
|
10 |
|
26265
|
11 |
subsection {* Auxiliary functions *}
|
|
12 |
|
|
13 |
definition
|
|
14 |
inc_shift :: "index \<Rightarrow> index \<Rightarrow> index"
|
|
15 |
where
|
|
16 |
"inc_shift v k = (if v = k then 1 else k + 1)"
|
|
17 |
|
|
18 |
definition
|
|
19 |
minus_shift :: "index \<Rightarrow> index \<Rightarrow> index \<Rightarrow> index"
|
|
20 |
where
|
|
21 |
"minus_shift r k l = (if k < l then r + k - l else k - l)"
|
|
22 |
|
|
23 |
function
|
|
24 |
log :: "index \<Rightarrow> index \<Rightarrow> index"
|
|
25 |
where
|
|
26 |
"log b i = (if b \<le> 1 \<or> i < b then 1 else 1 + log b (i div b))"
|
|
27 |
by pat_completeness auto
|
|
28 |
termination
|
|
29 |
by (relation "measure (nat_of_index o snd)")
|
|
30 |
(auto simp add: index)
|
|
31 |
|
|
32 |
|
|
33 |
subsection {* Random seeds *}
|
26038
|
34 |
|
|
35 |
types seed = "index \<times> index"
|
22528
|
36 |
|
26265
|
37 |
primrec
|
26038
|
38 |
"next" :: "seed \<Rightarrow> index \<times> seed"
|
|
39 |
where
|
26265
|
40 |
"next (v, w) = (let
|
|
41 |
k = v div 53668;
|
|
42 |
v' = minus_shift 2147483563 (40014 * (v mod 53668)) (k * 12211);
|
|
43 |
l = w div 52774;
|
|
44 |
w' = minus_shift 2147483399 (40692 * (w mod 52774)) (l * 3791);
|
|
45 |
z = minus_shift 2147483562 v' (w' + 1) + 1
|
|
46 |
in (z, (v', w')))"
|
|
47 |
|
|
48 |
lemma next_not_0:
|
|
49 |
"fst (next s) \<noteq> 0"
|
|
50 |
apply (cases s)
|
|
51 |
apply (auto simp add: minus_shift_def Let_def)
|
|
52 |
done
|
|
53 |
|
|
54 |
primrec
|
|
55 |
seed_invariant :: "seed \<Rightarrow> bool"
|
|
56 |
where
|
|
57 |
"seed_invariant (v, w) \<longleftrightarrow> 0 < v \<and> v < 9438322952 \<and> 0 < w \<and> True"
|
|
58 |
|
|
59 |
lemma if_same:
|
|
60 |
"(if b then f x else f y) = f (if b then x else y)"
|
|
61 |
by (cases b) simp_all
|
|
62 |
|
|
63 |
(*lemma seed_invariant:
|
|
64 |
assumes "seed_invariant (index_of_nat v, index_of_nat w)"
|
|
65 |
and "(index_of_nat z, (index_of_nat v', index_of_nat w')) = next (index_of_nat v, index_of_nat w)"
|
|
66 |
shows "seed_invariant (index_of_nat v', index_of_nat w')"
|
|
67 |
using assms
|
|
68 |
apply (auto simp add: seed_invariant_def)
|
|
69 |
apply (auto simp add: minus_shift_def Let_def)
|
|
70 |
apply (simp_all add: if_same cong del: if_cong)
|
|
71 |
apply safe
|
|
72 |
unfolding not_less
|
|
73 |
oops*)
|
22528
|
74 |
|
|
75 |
definition
|
26038
|
76 |
split_seed :: "seed \<Rightarrow> seed \<times> seed"
|
|
77 |
where
|
|
78 |
"split_seed s = (let
|
|
79 |
(v, w) = s;
|
|
80 |
(v', w') = snd (next s);
|
26265
|
81 |
v'' = inc_shift 2147483562 v;
|
26038
|
82 |
s'' = (v'', w');
|
26265
|
83 |
w'' = inc_shift 2147483398 w;
|
26038
|
84 |
s''' = (v', w'')
|
|
85 |
in (s'', s'''))"
|
|
86 |
|
|
87 |
|
26265
|
88 |
subsection {* Base selectors *}
|
22528
|
89 |
|
26038
|
90 |
function
|
|
91 |
range_aux :: "index \<Rightarrow> index \<Rightarrow> seed \<Rightarrow> index \<times> seed"
|
|
92 |
where
|
|
93 |
"range_aux k l s = (if k = 0 then (l, s) else
|
|
94 |
let (v, s') = next s
|
|
95 |
in range_aux (k - 1) (v + l * 2147483561) s')"
|
|
96 |
by pat_completeness auto
|
|
97 |
termination
|
|
98 |
by (relation "measure (nat_of_index o fst)")
|
|
99 |
(auto simp add: index)
|
22528
|
100 |
|
|
101 |
definition
|
26038
|
102 |
range :: "index \<Rightarrow> seed \<Rightarrow> index \<times> seed"
|
|
103 |
where
|
26265
|
104 |
"range k = (do
|
|
105 |
v \<leftarrow> range_aux (log 2147483561 k) 1;
|
|
106 |
return (v mod k)
|
|
107 |
done)"
|
|
108 |
|
|
109 |
lemma range:
|
|
110 |
assumes "k > 0"
|
|
111 |
shows "fst (range k s) < k"
|
|
112 |
proof -
|
|
113 |
obtain v w where range_aux:
|
|
114 |
"range_aux (log 2147483561 k) 1 s = (v, w)"
|
|
115 |
by (cases "range_aux (log 2147483561 k) 1 s")
|
|
116 |
with assms show ?thesis
|
26589
|
117 |
by (simp add: range_def run_def scomp_def split_def del: range_aux.simps log.simps)
|
26265
|
118 |
qed
|
26038
|
119 |
|
22528
|
120 |
definition
|
26038
|
121 |
select :: "'a list \<Rightarrow> seed \<Rightarrow> 'a \<times> seed"
|
|
122 |
where
|
26265
|
123 |
"select xs = (do
|
|
124 |
k \<leftarrow> range (index_of_nat (length xs));
|
|
125 |
return (nth xs (nat_of_index k))
|
|
126 |
done)"
|
|
127 |
|
|
128 |
lemma select:
|
|
129 |
assumes "xs \<noteq> []"
|
|
130 |
shows "fst (select xs s) \<in> set xs"
|
|
131 |
proof -
|
|
132 |
from assms have "index_of_nat (length xs) > 0" by simp
|
|
133 |
with range have
|
|
134 |
"fst (range (index_of_nat (length xs)) s) < index_of_nat (length xs)" by best
|
|
135 |
then have
|
|
136 |
"nat_of_index (fst (range (index_of_nat (length xs)) s)) < length xs" by simp
|
|
137 |
then show ?thesis
|
26589
|
138 |
by (auto simp add: select_def run_def scomp_def split_def)
|
26265
|
139 |
qed
|
22528
|
140 |
|
26038
|
141 |
definition
|
26265
|
142 |
select_default :: "index \<Rightarrow> 'a \<Rightarrow> 'a \<Rightarrow> seed \<Rightarrow> 'a \<times> seed"
|
26038
|
143 |
where
|
26265
|
144 |
[code func del]: "select_default k x y = (do
|
|
145 |
l \<leftarrow> range k;
|
|
146 |
return (if l + 1 < k then x else y)
|
|
147 |
done)"
|
|
148 |
|
|
149 |
lemma select_default_zero:
|
|
150 |
"fst (select_default 0 x y s) = y"
|
26589
|
151 |
by (simp add: run_def scomp_def split_def select_default_def)
|
26038
|
152 |
|
26265
|
153 |
lemma select_default_code [code]:
|
|
154 |
"select_default k x y = (if k = 0 then do
|
|
155 |
_ \<leftarrow> range 1;
|
|
156 |
return y
|
|
157 |
done else do
|
|
158 |
l \<leftarrow> range k;
|
|
159 |
return (if l + 1 < k then x else y)
|
|
160 |
done)"
|
|
161 |
proof (cases "k = 0")
|
|
162 |
case False then show ?thesis by (simp add: select_default_def)
|
22528
|
163 |
next
|
26265
|
164 |
case True then show ?thesis
|
26589
|
165 |
by (simp add: run_def scomp_def split_def select_default_def expand_fun_eq range_def)
|
26265
|
166 |
qed
|
22528
|
167 |
|
26265
|
168 |
|
|
169 |
subsection {* @{text ML} interface *}
|
22528
|
170 |
|
|
171 |
ML {*
|
26265
|
172 |
structure Random_Engine =
|
22528
|
173 |
struct
|
|
174 |
|
26038
|
175 |
type seed = int * int;
|
22528
|
176 |
|
|
177 |
local
|
26038
|
178 |
|
26265
|
179 |
val seed = ref
|
|
180 |
(let
|
|
181 |
val now = Time.toMilliseconds (Time.now ());
|
26038
|
182 |
val (q, s1) = IntInf.divMod (now, 2147483562);
|
|
183 |
val s2 = q mod 2147483398;
|
26265
|
184 |
in (s1 + 1, s2 + 1) end);
|
|
185 |
|
22528
|
186 |
in
|
26038
|
187 |
|
|
188 |
fun run f =
|
|
189 |
let
|
26265
|
190 |
val (x, seed') = f (! seed);
|
26038
|
191 |
val _ = seed := seed'
|
|
192 |
in x end;
|
|
193 |
|
22528
|
194 |
end;
|
|
195 |
|
|
196 |
end;
|
|
197 |
*}
|
|
198 |
|
26038
|
199 |
end
|