1#[cfg(test)]
4mod test;
5#[cfg(test)]
6mod test_data;
7
8use std::num::Wrapping;
9
10pub const M: Wrapping<i64> = Wrapping((1 << 48) - 1);
12
13pub const A: Wrapping<i64> = Wrapping(0x5DEECE66D);
15
16pub const C: Wrapping<i64> = Wrapping(11);
18
19const F32_DIV: f32 = (1u32 << 24) as f32;
20const F64_DIV: f64 = (1u64 << 53) as f64;
21
22#[derive(Debug, Clone)]
23pub struct Random {
24 state: Wrapping<i64>,
25 next_gaussian: Option<f64>
26}
27
28impl Random {
29 pub fn new(seed: u64) -> Self {
30 Random {
31 state: Wrapping((seed as i64) ^ A.0) & M,
32 next_gaussian: None
33 }
34 }
35
36 pub fn set_seed(&mut self, seed: u64) {
38 *self = Random::new(seed);
39 }
40
41 pub fn next(&mut self, bits: u8) -> i32 {
46 if bits > 48 {
47 panic!("Too many bits!")
48 }
49
50 self.state = (self.state * A + C) & M;
51
52 ((self.state.0 as u64) >> (48 - bits)) as i32
53 }
54
55 pub fn next_bytes(&mut self, bytes: &mut [u8]) {
57 for chunk in bytes.chunks_mut(4) {
58 let mut block = self.next_u32();
59
60 for item in chunk {
61 *item = (block & 0xFF) as u8;
62 block >>= 8;
63 }
64 }
65 }
66
67 pub fn next_i32(&mut self) -> i32 {
69 self.next(32) as i32
70 }
71
72 pub fn next_u32(&mut self) -> u32 {
74 self.next(32) as u32
75 }
76
77 pub fn next_i32_bound(&mut self, max: i32) -> i32 {
84 if max <= 0 {
85 panic!("Maximum must be > 0")
86 }
87
88 if (max as u32).is_power_of_two() {
89 let max = max as i64;
90
91 return ((max.wrapping_mul(self.next(31) as i64)) >> 31) as i32;
92 }
93
94 let mut bits = self.next(31);
95 let mut val = bits % max;
96
97 while bits.wrapping_sub(val).wrapping_add(max - 1) < 0 {
98 bits = self.next(31);
99 val = bits % max;
100 }
101
102 val
103 }
104
105 pub fn next_u32_bound(&mut self, max: u32) -> u32 {
113 self.next_i32_bound(max as i32) as u32
114 }
115
116 pub fn next_i64(&mut self) -> i64 {
118 ((self.next(32) as i64) << 32).wrapping_add(self.next(32) as i64)
119 }
120
121 pub fn next_u64(&mut self) -> u64 {
123 self.next_i64() as u64
124 }
125
126 pub fn next_bool(&mut self) -> bool {
128 self.next(1) == 1
129 }
130
131 pub fn next_f32(&mut self) -> f32 {
133 (self.next(24) as f32) / F32_DIV
134 }
135
136 pub fn next_f64(&mut self) -> f64 {
138 let high = (self.next(26) as i64) << 27;
139 let low = self.next(27) as i64;
140
141 (high.wrapping_add(low) as f64) / F64_DIV
142 }
143
144 fn next_gaussian_pair(&mut self) -> (f64, f64) {
146 let mut next_candidate = || {
147 let v = (
148 2.0 * self.next_f64() - 1.0,
149 2.0 * self.next_f64() - 1.0
150 );
151
152 (v, v.0*v.0 + v.1*v.1)
153 };
154
155 let (mut v, mut s) = next_candidate();
156
157 while s >= 1.0 || s == 0.0 {
158 let (vn, sn) = next_candidate();
159 v = vn;
160 s = sn;
161 }
162
163 let multiplier = ((s.log(::std::f64::consts::E) / s) * -2.0).sqrt();
165
166 (v.0 * multiplier, v.1 * multiplier)
167 }
168
169 pub fn next_gaussian(&mut self) -> f64 {
171 match self.next_gaussian.take() {
172 Some(next) => next,
173 None => {
174 let (v0, v1) = self.next_gaussian_pair();
175
176 self.next_gaussian = Some(v1);
177
178 v0
179 }
180 }
181 }
182}