src/HOL/ex/Random.thy
author haftmann
Sat Sep 06 14:02:36 2008 +0200 (2008-09-06)
changeset 28145 af3923ed4786
parent 28042 1471f2974eb1
child 28562 4e74209f113e
permissions -rw-r--r--
dropped "run" marker in monad syntax
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
krauss@28042
    23
fun
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
haftmann@26265
    28
subsection {* Random seeds *}
haftmann@26038
    29
haftmann@26038
    30
types seed = "index \<times> index"
haftmann@22528
    31
haftmann@26265
    32
primrec
haftmann@26038
    33
  "next" :: "seed \<Rightarrow> index \<times> seed"
haftmann@26038
    34
where
haftmann@26265
    35
  "next (v, w) = (let
haftmann@26265
    36
     k =  v div 53668;
haftmann@26265
    37
     v' = minus_shift 2147483563 (40014 * (v mod 53668)) (k * 12211);
haftmann@26265
    38
     l =  w div 52774;
haftmann@26265
    39
     w' = minus_shift 2147483399 (40692 * (w mod 52774)) (l * 3791);
haftmann@26265
    40
     z =  minus_shift 2147483562 v' (w' + 1) + 1
haftmann@26265
    41
   in (z, (v', w')))"
haftmann@26265
    42
haftmann@26265
    43
lemma next_not_0:
haftmann@26265
    44
  "fst (next s) \<noteq> 0"
haftmann@26265
    45
apply (cases s)
haftmann@26265
    46
apply (auto simp add: minus_shift_def Let_def)
haftmann@26265
    47
done
haftmann@26265
    48
haftmann@26265
    49
primrec
haftmann@26265
    50
  seed_invariant :: "seed \<Rightarrow> bool"
haftmann@26265
    51
where
haftmann@26265
    52
  "seed_invariant (v, w) \<longleftrightarrow> 0 < v \<and> v < 9438322952 \<and> 0 < w \<and> True"
haftmann@26265
    53
haftmann@26265
    54
lemma if_same:
haftmann@26265
    55
  "(if b then f x else f y) = f (if b then x else y)"
haftmann@26265
    56
  by (cases b) simp_all
haftmann@26265
    57
haftmann@22528
    58
definition
haftmann@26038
    59
  split_seed :: "seed \<Rightarrow> seed \<times> seed"
haftmann@26038
    60
where
haftmann@26038
    61
  "split_seed s = (let
haftmann@26038
    62
     (v, w) = s;
haftmann@26038
    63
     (v', w') = snd (next s);
haftmann@26265
    64
     v'' = inc_shift 2147483562 v;
haftmann@26038
    65
     s'' = (v'', w');
haftmann@26265
    66
     w'' = inc_shift 2147483398 w;
haftmann@26038
    67
     s''' = (v', w'')
haftmann@26038
    68
   in (s'', s'''))"
haftmann@26038
    69
haftmann@26038
    70
haftmann@26265
    71
subsection {* Base selectors *}
haftmann@22528
    72
haftmann@26038
    73
function
haftmann@26038
    74
  range_aux :: "index \<Rightarrow> index \<Rightarrow> seed \<Rightarrow> index \<times> seed"
haftmann@26038
    75
where
haftmann@26038
    76
  "range_aux k l s = (if k = 0 then (l, s) else
haftmann@26038
    77
    let (v, s') = next s
haftmann@26038
    78
  in range_aux (k - 1) (v + l * 2147483561) s')"
haftmann@26038
    79
by pat_completeness auto
haftmann@26038
    80
termination
haftmann@26038
    81
  by (relation "measure (nat_of_index o fst)")
haftmann@26038
    82
    (auto simp add: index)
haftmann@22528
    83
haftmann@22528
    84
definition
haftmann@26038
    85
  range :: "index \<Rightarrow> seed \<Rightarrow> index \<times> seed"
haftmann@26038
    86
where
haftmann@26265
    87
  "range k = (do
haftmann@26265
    88
     v \<leftarrow> range_aux (log 2147483561 k) 1;
haftmann@26265
    89
     return (v mod k)
haftmann@26265
    90
   done)"
haftmann@26265
    91
haftmann@26265
    92
lemma range:
haftmann@26265
    93
  assumes "k > 0"
haftmann@26265
    94
  shows "fst (range k s) < k"
haftmann@26265
    95
proof -
haftmann@26265
    96
  obtain v w where range_aux:
haftmann@26265
    97
    "range_aux (log 2147483561 k) 1 s = (v, w)"
haftmann@26265
    98
    by (cases "range_aux (log 2147483561 k) 1 s")
haftmann@26265
    99
  with assms show ?thesis
haftmann@28145
   100
    by (simp add: monad_collapse range_def del: range_aux.simps log.simps)
haftmann@26265
   101
qed
haftmann@26038
   102
haftmann@22528
   103
definition
haftmann@26038
   104
  select :: "'a list \<Rightarrow> seed \<Rightarrow> 'a \<times> seed"
haftmann@26038
   105
where
haftmann@26265
   106
  "select xs = (do
haftmann@26265
   107
     k \<leftarrow> range (index_of_nat (length xs));
haftmann@26265
   108
     return (nth xs (nat_of_index k))
haftmann@26265
   109
   done)"
haftmann@26265
   110
haftmann@26265
   111
lemma select:
haftmann@26265
   112
  assumes "xs \<noteq> []"
haftmann@26265
   113
  shows "fst (select xs s) \<in> set xs"
haftmann@26265
   114
proof -
haftmann@26265
   115
  from assms have "index_of_nat (length xs) > 0" by simp
haftmann@26265
   116
  with range have
haftmann@26265
   117
    "fst (range (index_of_nat (length xs)) s) < index_of_nat (length xs)" by best
haftmann@26265
   118
  then have
haftmann@26265
   119
    "nat_of_index (fst (range (index_of_nat (length xs)) s)) < length xs" by simp
haftmann@26265
   120
  then show ?thesis
haftmann@28145
   121
    by (auto simp add: monad_collapse select_def)
haftmann@26265
   122
qed
haftmann@22528
   123
haftmann@26038
   124
definition
haftmann@26265
   125
  select_default :: "index \<Rightarrow> 'a \<Rightarrow> 'a \<Rightarrow> seed \<Rightarrow> 'a \<times> seed"
haftmann@26038
   126
where
haftmann@26265
   127
  [code func del]: "select_default k x y = (do
haftmann@26265
   128
     l \<leftarrow> range k;
haftmann@26265
   129
     return (if l + 1 < k then x else y)
haftmann@26265
   130
   done)"
haftmann@26265
   131
haftmann@26265
   132
lemma select_default_zero:
haftmann@26265
   133
  "fst (select_default 0 x y s) = y"
haftmann@28145
   134
  by (simp add: monad_collapse select_default_def)
haftmann@26038
   135
haftmann@26265
   136
lemma select_default_code [code]:
haftmann@26265
   137
  "select_default k x y = (if k = 0 then do
haftmann@26265
   138
     _ \<leftarrow> range 1;
haftmann@26265
   139
     return y
haftmann@26265
   140
   done else do
haftmann@26265
   141
     l \<leftarrow> range k;
haftmann@26265
   142
     return (if l + 1 < k then x else y)
haftmann@26265
   143
   done)"
haftmann@26265
   144
proof (cases "k = 0")
haftmann@26265
   145
  case False then show ?thesis by (simp add: select_default_def)
haftmann@22528
   146
next
haftmann@26265
   147
  case True then show ?thesis
haftmann@28145
   148
    by (simp add: monad_collapse select_default_def range_def)
haftmann@26265
   149
qed
haftmann@22528
   150
haftmann@26265
   151
haftmann@26265
   152
subsection {* @{text ML} interface *}
haftmann@22528
   153
haftmann@22528
   154
ML {*
haftmann@26265
   155
structure Random_Engine =
haftmann@22528
   156
struct
haftmann@22528
   157
haftmann@26038
   158
type seed = int * int;
haftmann@22528
   159
haftmann@22528
   160
local
haftmann@26038
   161
haftmann@26265
   162
val seed = ref 
haftmann@26265
   163
  (let
haftmann@26265
   164
    val now = Time.toMilliseconds (Time.now ());
haftmann@26038
   165
    val (q, s1) = IntInf.divMod (now, 2147483562);
haftmann@26038
   166
    val s2 = q mod 2147483398;
haftmann@26265
   167
  in (s1 + 1, s2 + 1) end);
haftmann@26265
   168
haftmann@22528
   169
in
haftmann@26038
   170
haftmann@26038
   171
fun run f =
haftmann@26038
   172
  let
haftmann@26265
   173
    val (x, seed') = f (! seed);
haftmann@26038
   174
    val _ = seed := seed'
haftmann@26038
   175
  in x end;
haftmann@26038
   176
haftmann@22528
   177
end;
haftmann@22528
   178
haftmann@22528
   179
end;
haftmann@22528
   180
*}
haftmann@22528
   181
haftmann@26038
   182
end
haftmann@28145
   183