1 (* Author: Florian Haftmann, TU Muenchen *) |
|
2 |
|
3 header {* A HOL random engine *} |
|
4 |
|
5 theory Random |
|
6 imports Code_Index |
|
7 begin |
|
8 |
|
9 notation fcomp (infixl "o>" 60) |
|
10 notation scomp (infixl "o\<rightarrow>" 60) |
|
11 |
|
12 |
|
13 subsection {* Auxiliary functions *} |
|
14 |
|
15 definition inc_shift :: "index \<Rightarrow> index \<Rightarrow> index" where |
|
16 "inc_shift v k = (if v = k then 1 else k + 1)" |
|
17 |
|
18 definition minus_shift :: "index \<Rightarrow> index \<Rightarrow> index \<Rightarrow> index" where |
|
19 "minus_shift r k l = (if k < l then r + k - l else k - l)" |
|
20 |
|
21 fun log :: "index \<Rightarrow> index \<Rightarrow> index" where |
|
22 "log b i = (if b \<le> 1 \<or> i < b then 1 else 1 + log b (i div b))" |
|
23 |
|
24 |
|
25 subsection {* Random seeds *} |
|
26 |
|
27 types seed = "index \<times> index" |
|
28 |
|
29 primrec "next" :: "seed \<Rightarrow> index \<times> seed" where |
|
30 "next (v, w) = (let |
|
31 k = v div 53668; |
|
32 v' = minus_shift 2147483563 (40014 * (v mod 53668)) (k * 12211); |
|
33 l = w div 52774; |
|
34 w' = minus_shift 2147483399 (40692 * (w mod 52774)) (l * 3791); |
|
35 z = minus_shift 2147483562 v' (w' + 1) + 1 |
|
36 in (z, (v', w')))" |
|
37 |
|
38 lemma next_not_0: |
|
39 "fst (next s) \<noteq> 0" |
|
40 by (cases s) (auto simp add: minus_shift_def Let_def) |
|
41 |
|
42 primrec seed_invariant :: "seed \<Rightarrow> bool" where |
|
43 "seed_invariant (v, w) \<longleftrightarrow> 0 < v \<and> v < 9438322952 \<and> 0 < w \<and> True" |
|
44 |
|
45 lemma if_same: "(if b then f x else f y) = f (if b then x else y)" |
|
46 by (cases b) simp_all |
|
47 |
|
48 definition split_seed :: "seed \<Rightarrow> seed \<times> seed" where |
|
49 "split_seed s = (let |
|
50 (v, w) = s; |
|
51 (v', w') = snd (next s); |
|
52 v'' = inc_shift 2147483562 v; |
|
53 s'' = (v'', w'); |
|
54 w'' = inc_shift 2147483398 w; |
|
55 s''' = (v', w'') |
|
56 in (s'', s'''))" |
|
57 |
|
58 |
|
59 subsection {* Base selectors *} |
|
60 |
|
61 fun iterate :: "index \<Rightarrow> ('b \<Rightarrow> 'a \<Rightarrow> 'b \<times> 'a) \<Rightarrow> 'b \<Rightarrow> 'a \<Rightarrow> 'b \<times> 'a" where |
|
62 "iterate k f x = (if k = 0 then Pair x else f x o\<rightarrow> iterate (k - 1) f)" |
|
63 |
|
64 definition range :: "index \<Rightarrow> seed \<Rightarrow> index \<times> seed" where |
|
65 "range k = iterate (log 2147483561 k) |
|
66 (\<lambda>l. next o\<rightarrow> (\<lambda>v. Pair (v + l * 2147483561))) 1 |
|
67 o\<rightarrow> (\<lambda>v. Pair (v mod k))" |
|
68 |
|
69 lemma range: |
|
70 "k > 0 \<Longrightarrow> fst (range k s) < k" |
|
71 by (simp add: range_def scomp_apply split_def del: log.simps iterate.simps) |
|
72 |
|
73 definition select :: "'a list \<Rightarrow> seed \<Rightarrow> 'a \<times> seed" where |
|
74 "select xs = range (Code_Index.of_nat (length xs)) |
|
75 o\<rightarrow> (\<lambda>k. Pair (nth xs (Code_Index.nat_of k)))" |
|
76 |
|
77 lemma select: |
|
78 assumes "xs \<noteq> []" |
|
79 shows "fst (select xs s) \<in> set xs" |
|
80 proof - |
|
81 from assms have "Code_Index.of_nat (length xs) > 0" by simp |
|
82 with range have |
|
83 "fst (range (Code_Index.of_nat (length xs)) s) < Code_Index.of_nat (length xs)" by best |
|
84 then have |
|
85 "Code_Index.nat_of (fst (range (Code_Index.of_nat (length xs)) s)) < length xs" by simp |
|
86 then show ?thesis |
|
87 by (simp add: scomp_apply split_beta select_def) |
|
88 qed |
|
89 |
|
90 definition select_default :: "index \<Rightarrow> 'a \<Rightarrow> 'a \<Rightarrow> seed \<Rightarrow> 'a \<times> seed" where |
|
91 [code del]: "select_default k x y = range k |
|
92 o\<rightarrow> (\<lambda>l. Pair (if l + 1 < k then x else y))" |
|
93 |
|
94 lemma select_default_zero: |
|
95 "fst (select_default 0 x y s) = y" |
|
96 by (simp add: scomp_apply split_beta select_default_def) |
|
97 |
|
98 lemma select_default_code [code]: |
|
99 "select_default k x y = (if k = 0 |
|
100 then range 1 o\<rightarrow> (\<lambda>_. Pair y) |
|
101 else range k o\<rightarrow> (\<lambda>l. Pair (if l + 1 < k then x else y)))" |
|
102 proof |
|
103 fix s |
|
104 have "snd (range (Code_Index.of_nat 0) s) = snd (range (Code_Index.of_nat 1) s)" |
|
105 by (simp add: range_def scomp_Pair scomp_apply split_beta) |
|
106 then show "select_default k x y s = (if k = 0 |
|
107 then range 1 o\<rightarrow> (\<lambda>_. Pair y) |
|
108 else range k o\<rightarrow> (\<lambda>l. Pair (if l + 1 < k then x else y))) s" |
|
109 by (cases "k = 0") (simp_all add: select_default_def scomp_apply split_beta) |
|
110 qed |
|
111 |
|
112 |
|
113 subsection {* @{text ML} interface *} |
|
114 |
|
115 ML {* |
|
116 structure Random_Engine = |
|
117 struct |
|
118 |
|
119 type seed = int * int; |
|
120 |
|
121 local |
|
122 |
|
123 val seed = ref |
|
124 (let |
|
125 val now = Time.toMilliseconds (Time.now ()); |
|
126 val (q, s1) = IntInf.divMod (now, 2147483562); |
|
127 val s2 = q mod 2147483398; |
|
128 in (s1 + 1, s2 + 1) end); |
|
129 |
|
130 in |
|
131 |
|
132 fun run f = |
|
133 let |
|
134 val (x, seed') = f (! seed); |
|
135 val _ = seed := seed' |
|
136 in x end; |
|
137 |
|
138 end; |
|
139 |
|
140 end; |
|
141 *} |
|
142 |
|
143 no_notation fcomp (infixl "o>" 60) |
|
144 no_notation scomp (infixl "o\<rightarrow>" 60) |
|
145 |
|
146 end |
|
147 |
|