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