-
Notifications
You must be signed in to change notification settings - Fork 80
/
random.lean
114 lines (101 loc) · 4.02 KB
/
random.lean
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
/-
Copyright (c) 2018 Microsoft Corporation. All rights reserved.
Released under Apache 2.0 license as described in the file LICENSE.
Authors: Leonardo de Moura
-/
universes u
/-!
# Basic random number generator support
based on the one available on the Haskell library
-/
/-- Interface for random number generators. -/
class random_gen (g : Type u) :=
/- `range` returns the range of values returned by
the generator. -/
(range : g → nat × nat)
/- `next` operation returns a natural number that is uniformly distributed
the range returned by `range` (including both end points),
and a new generator. -/
(next : g → nat × g)
/-
The 'split' operation allows one to obtain two distinct random number
generators. This is very useful in functional programs (for example, when
passing a random number generator down to recursive calls). -/
(split : g → g × g)
/-- "Standard" random number generator. -/
structure std_gen :=
(s1 : nat) (s2 : nat)
def std_range := (1, 2147483562)
instance : has_repr std_gen :=
{ repr := λ ⟨s1, s2⟩, "⟨" ++ to_string s1 ++ ", " ++ to_string s2 ++ "⟩" }
def std_next : std_gen → nat × std_gen
| ⟨s1, s2⟩ :=
let k : int := s1 / 53668,
s1' : int := 40014 * ((s1 : int) - k * 53668) - k * 12211,
s1'' : int := if s1' < 0 then s1' + 2147483563 else s1',
k' : int := s2 / 52774,
s2' : int := 40692 * ((s2 : int) - k' * 52774) - k' * 3791,
s2'' : int := if s2' < 0 then s2' + 2147483399 else s2',
z : int := s1'' - s2'',
z' : int := if z < 1 then z + 2147483562 else z % 2147483562
in (z'.to_nat, ⟨s1''.to_nat, s2''.to_nat⟩)
def std_split : std_gen → std_gen × std_gen
| g@⟨s1, s2⟩ :=
let new_s1 := if s1 = 2147483562 then 1 else s1 + 1,
new_s2 := if s2 = 1 then 2147483398 else s2 - 1,
new_g := (std_next g).2,
left_g := std_gen.mk new_s1 new_g.2,
right_g := std_gen.mk new_g.1 new_s2
in (left_g, right_g)
instance : random_gen std_gen :=
{range := λ _, std_range,
next := std_next,
split := std_split}
/-- Return a standard number generator. -/
def mk_std_gen (s : nat := 0) : std_gen :=
let q := s / 2147483562,
s1 := s % 2147483562,
s2 := q % 2147483398 in
⟨s1 + 1, s2 + 1⟩
/--
Auxiliary function for random_nat_val.
Generate random values until we exceed the target magnitude.
`gen_lo` and `gen_mag` are the generator lower bound and magnitude.
The parameter `r` is the "remaining" magnitude.
-/
private def rand_nat_aux {gen : Type u} [random_gen gen] (gen_lo gen_mag : nat) (h : gen_mag > 0) : nat → nat → gen → nat × gen
| 0 v g := (v, g)
| r'@(r+1) v g :=
let (x, g') := random_gen.next g,
v' := v*gen_mag + (x - gen_lo)
in have r' / gen_mag - 1 < r',
begin
by_cases h : (r + 1) / gen_mag = 0,
{ rw [h], simp, apply nat.zero_lt_succ },
{ have : (r + 1) / gen_mag > 0, from nat.pos_of_ne_zero h,
have h₁ : (r + 1) / gen_mag - 1 < (r + 1) / gen_mag, { apply nat.sub_lt, assumption, tactic.comp_val },
have h₂ : (r + 1) / gen_mag ≤ r + 1, { apply nat.div_le_self },
exact lt_of_lt_of_le h₁ h₂ }
end,
rand_nat_aux (r' / gen_mag - 1) v' g'
/-- Generate a random natural number in the interval [lo, hi]. -/
def rand_nat {gen : Type u} [random_gen gen] (g : gen) (lo hi : nat) : nat × gen :=
let lo' := if lo > hi then hi else lo,
hi' := if lo > hi then lo else hi,
(gen_lo, gen_hi) := random_gen.range g,
gen_mag := gen_hi - gen_lo + 1,
/-
Probabilities of the most likely and least likely result
will differ at most by a factor of (1 +- 1/q). Assuming the RandomGen
is uniform, of course
-/
q := 1000,
k := hi' - lo' + 1,
tgt_mag := k * q,
(v, g') := rand_nat_aux gen_lo gen_mag (nat.zero_lt_succ _) tgt_mag 0 g,
v' := lo' + (v % k)
in (v', g')
/-- Generate a random Boolean. -/
def rand_bool {gen : Type u} [random_gen gen] (g : gen) : bool × gen :=
let (v, g') := rand_nat g 0 1
in (v = 1, g')