src/HOL/Library/Random.thy
author chaieb
Mon Feb 09 17:21:46 2009 +0000 (2009-02-09)
changeset 29847 af32126ee729
parent 29823 0ab754d13ccd
child 30495 a5f1e4f46d14
permissions -rw-r--r--
added Determinants to Library
     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 subsection {* Random seeds *}
    25 
    26 types seed = "index \<times> index"
    27 
    28 primrec "next" :: "seed \<Rightarrow> index \<times> seed" where
    29   "next (v, w) = (let
    30      k =  v div 53668;
    31      v' = minus_shift 2147483563 (40014 * (v mod 53668)) (k * 12211);
    32      l =  w div 52774;
    33      w' = minus_shift 2147483399 (40692 * (w mod 52774)) (l * 3791);
    34      z =  minus_shift 2147483562 v' (w' + 1) + 1
    35    in (z, (v', w')))"
    36 
    37 lemma next_not_0:
    38   "fst (next s) \<noteq> 0"
    39   by (cases s) (auto simp add: minus_shift_def Let_def)
    40 
    41 primrec seed_invariant :: "seed \<Rightarrow> bool" where
    42   "seed_invariant (v, w) \<longleftrightarrow> 0 < v \<and> v < 9438322952 \<and> 0 < w \<and> True"
    43 
    44 lemma if_same: "(if b then f x else f y) = f (if b then x else y)"
    45   by (cases b) simp_all
    46 
    47 definition split_seed :: "seed \<Rightarrow> seed \<times> seed" where
    48   "split_seed s = (let
    49      (v, w) = s;
    50      (v', w') = snd (next s);
    51      v'' = inc_shift 2147483562 v;
    52      s'' = (v'', w');
    53      w'' = inc_shift 2147483398 w;
    54      s''' = (v', w'')
    55    in (s'', s'''))"
    56 
    57 
    58 subsection {* Base selectors *}
    59 
    60 function range_aux :: "index \<Rightarrow> index \<Rightarrow> seed \<Rightarrow> index \<times> seed" where
    61   "range_aux k l s = (if k = 0 then (l, s) else
    62     let (v, s') = next s
    63   in range_aux (k - 1) (v + l * 2147483561) s')"
    64 by pat_completeness auto
    65 termination
    66   by (relation "measure (Code_Index.nat_of o fst)")
    67     (auto simp add: index)
    68 
    69 definition range :: "index \<Rightarrow> seed \<Rightarrow> index \<times> seed" where
    70   "range k = range_aux (log 2147483561 k) 1
    71     o\<rightarrow> (\<lambda>v. Pair (v mod k))"
    72 
    73 lemma range:
    74   assumes "k > 0"
    75   shows "fst (range k s) < k"
    76 proof -
    77   obtain v w where range_aux:
    78     "range_aux (log 2147483561 k) 1 s = (v, w)"
    79     by (cases "range_aux (log 2147483561 k) 1 s")
    80   with assms show ?thesis
    81     by (simp add: scomp_apply range_def del: range_aux.simps log.simps)
    82 qed
    83 
    84 definition select :: "'a list \<Rightarrow> seed \<Rightarrow> 'a \<times> seed" where
    85   "select xs = range (Code_Index.of_nat (length xs))
    86     o\<rightarrow> (\<lambda>k. Pair (nth xs (Code_Index.nat_of k)))"
    87      
    88 lemma select:
    89   assumes "xs \<noteq> []"
    90   shows "fst (select xs s) \<in> set xs"
    91 proof -
    92   from assms have "Code_Index.of_nat (length xs) > 0" by simp
    93   with range have
    94     "fst (range (Code_Index.of_nat (length xs)) s) < Code_Index.of_nat (length xs)" by best
    95   then have
    96     "Code_Index.nat_of (fst (range (Code_Index.of_nat (length xs)) s)) < length xs" by simp
    97   then show ?thesis
    98     by (simp add: scomp_apply split_beta select_def)
    99 qed
   100 
   101 definition select_default :: "index \<Rightarrow> 'a \<Rightarrow> 'a \<Rightarrow> seed \<Rightarrow> 'a \<times> seed" where
   102   [code del]: "select_default k x y = range k
   103      o\<rightarrow> (\<lambda>l. Pair (if l + 1 < k then x else y))"
   104 
   105 lemma select_default_zero:
   106   "fst (select_default 0 x y s) = y"
   107   by (simp add: scomp_apply split_beta select_default_def)
   108 
   109 lemma select_default_code [code]:
   110   "select_default k x y = (if k = 0
   111     then range 1 o\<rightarrow> (\<lambda>_. Pair y)
   112     else range k o\<rightarrow> (\<lambda>l. Pair (if l + 1 < k then x else y)))"
   113 proof
   114   fix s
   115   have "snd (range (Code_Index.of_nat 0) s) = snd (range (Code_Index.of_nat 1) s)"
   116     by (simp add: range_def scomp_Pair scomp_apply split_beta)
   117   then show "select_default k x y s = (if k = 0
   118     then range 1 o\<rightarrow> (\<lambda>_. Pair y)
   119     else range k o\<rightarrow> (\<lambda>l. Pair (if l + 1 < k then x else y))) s"
   120     by (cases "k = 0") (simp_all add: select_default_def scomp_apply split_beta)
   121 qed
   122 
   123 
   124 subsection {* @{text ML} interface *}
   125 
   126 ML {*
   127 structure Random_Engine =
   128 struct
   129 
   130 type seed = int * int;
   131 
   132 local
   133 
   134 val seed = ref 
   135   (let
   136     val now = Time.toMilliseconds (Time.now ());
   137     val (q, s1) = IntInf.divMod (now, 2147483562);
   138     val s2 = q mod 2147483398;
   139   in (s1 + 1, s2 + 1) end);
   140 
   141 in
   142 
   143 fun run f =
   144   let
   145     val (x, seed') = f (! seed);
   146     val _ = seed := seed'
   147   in x end;
   148 
   149 end;
   150 
   151 end;
   152 *}
   153 
   154 no_notation fcomp (infixl "o>" 60)
   155 no_notation scomp (infixl "o\<rightarrow>" 60)
   156 
   157 end
   158