doryen_extra/random/
algorithms.rs

1/* BSD 3-Clause License
2 *
3 * Copyright © 2019, Alexander Krivács Schrøder <alexschrod@gmail.com>.
4 * Copyright © 2008-2019, Jice and the libtcod contributors.
5 * All rights reserved.
6 *
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions are met:
9 *
10 * 1. Redistributions of source code must retain the above copyright notice,
11 *    this list of conditions and the following disclaimer.
12 *
13 * 2. Redistributions in binary form must reproduce the above copyright notice,
14 *    this list of conditions and the following disclaimer in the documentation
15 *    and/or other materials provided with the distribution.
16 *
17 * 3. Neither the name of the copyright holder nor the names of its
18 *    contributors may be used to endorse or promote products derived from
19 *    this software without specific prior written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
25 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 * POSSIBILITY OF SUCH DAMAGE.
32 */
33
34//! Random number generator algorithms.
35
36use std::mem::{transmute, MaybeUninit};
37
38const RAND_DIV: f32 = 1.0 / 0xffff_ffff_u32 as f32; // u32::MAX
39#[allow(clippy::unnecessary_cast)]
40const RAND_DIV_DOUBLE: f64 = 1.0 / 0xffff_ffff_u32 as f64; // u32::MAX
41
42/// Random number generator algorithm trait.
43pub trait Algorithm {
44    /// Generate a 32-bit integer.
45    fn get_int(&mut self) -> u32;
46
47    /// Generate a 32-bit floating point number.
48    fn get_float(&mut self) -> f32 {
49        if cfg!(feature = "libtcod-compat") {
50            self.get_int() as f32 * RAND_DIV
51        } else {
52            // Generate a random 32-bit floating point number between 0 and 1
53            // using the Allen Downey algorithm described here:
54            // <https://allendowney.com/research/rand/downey07randfloat.pdf>
55            //
56            // The desirable properties of this algorithm is that it can produce
57            // every representable floating-point value in a given range
58            // (as opposed to about 7% of them with the "libtcod-compat" method
59            // above) and it produces uniformly distributed values.
60            //
61            // Thanks to <https://github.com/orlp> for the link to the article.
62
63            let low_exp = 0;
64            let high_exp = 127;
65
66            let mut bits = Bits::new(self);
67            let mut exp = high_exp - 1;
68            while exp > low_exp {
69                if bits.get_bit() != 0 {
70                    break;
71                }
72                exp -= 1;
73            }
74
75            let mantissa = bits.algorithm.get_int() & 0x7FFFFF;
76            if mantissa == 0 && bits.get_bit() != 0 {
77                exp += 1;
78            }
79
80            let ans = (exp << 23) | mantissa;
81
82            f32::from_bits(ans)
83        }
84    }
85
86    /// Generate a 64-bit floating point number.
87    fn get_double(&mut self) -> f64 {
88        if cfg!(feature = "libtcod-compat") {
89            f64::from(self.get_int()) * RAND_DIV_DOUBLE
90        } else {
91            // Generate a random 64-bit floating point number between 0 and 1
92            // using the Allen Downey algorithm described here:
93            // <https://allendowney.com/research/rand/downey07randfloat.pdf>
94            //
95            // The desirable properties of this algorithm is that it can produce
96            // every representable floating-point value in a given range
97            // (as opposed to a minuscule amount of them with the
98            // "libtcod-compat" method above) and it produces uniformly
99            // distributed values.
100            //
101            // Thanks to <https://github.com/orlp> for the link to the article.
102
103            let low_exp = 0;
104            let high_exp = 1023;
105
106            let mut bits = Bits::new(self);
107            let mut exp = high_exp - 1;
108            while exp > low_exp {
109                if bits.get_bit() != 0 {
110                    break;
111                }
112                exp -= 1;
113            }
114
115            let mantissa = (u64::from(bits.algorithm.get_int()) << 32
116                | u64::from(bits.algorithm.get_int()))
117                & 0xFFFFFFFFFFFFF;
118            if mantissa == 0 && bits.get_bit() != 0 {
119                exp += 1;
120            }
121
122            let ans = (exp << 52) | mantissa;
123
124            f64::from_bits(ans)
125        }
126    }
127}
128
129/// Mersenne Twister algorithm.
130#[derive(Clone, Copy)]
131pub struct MersenneTwister {
132    mt: [u32; Self::MT19937_RECURRENCE_DEGREE],
133    cur_mt: usize,
134}
135
136impl MersenneTwister {
137    const MT19937: u32 = 1_812_433_253;
138    const MT19937_WORD_SIZE: usize = 32;
139    const MT19937_RECURRENCE_DEGREE: usize = 624;
140    const MT19937_SEPARATION_POINT: usize = 31;
141    const MT19937_MIDDLE_WORD: usize = 397;
142    const MT19937_RATIONAL_NORMAL_FORM_TWIST_MATRIX_COEFFICIENTS: u32 = 0x9908_B0DF;
143    const MT19937_TGFSR_R_TEMPERING_BIT_MASKS: (u32, u32) = (0x9D2C_5680, 0xEFC6_0000);
144    const MT19937_TGFSR_R_TEMPERING_BIT_SHIFTS: (u32, u32) = (7, 15);
145    const MT19937_ADDITIONAL_TEMPERING: (u32, u32, u32) = (11, 0xFFFF_FFFF, 18);
146    const MT19937_LOWER_MASK: u32 = (1 << (Self::MT19937_SEPARATION_POINT)) as u32;
147    const MT19937_UPPER_MASK: u32 = !Self::MT19937_LOWER_MASK;
148
149    /// Create a new Mersenne Twister algorithm instance.
150    pub fn new(seed: u32) -> Self {
151        Self {
152            cur_mt: 624,
153            mt: Self::mt_init(seed),
154        }
155    }
156
157    /* initialize the mersenne twister array */
158    #[allow(unsafe_code)]
159    fn mt_init(seed: u32) -> [u32; Self::MT19937_RECURRENCE_DEGREE] {
160        let mut mt: [MaybeUninit<u32>; Self::MT19937_RECURRENCE_DEGREE] =
161            unsafe { MaybeUninit::uninit().assume_init() };
162        mt[0] = MaybeUninit::new(seed);
163        for i in 1..mt.len() {
164            mt[i] = MaybeUninit::new(unsafe {
165                Self::MT19937.wrapping_mul(
166                    (*mt[i - 1].as_ptr()
167                        ^ (*mt[i - 1].as_ptr() >> (Self::MT19937_WORD_SIZE as u32 - 2)))
168                        .wrapping_add(i as u32),
169                )
170            });
171        }
172
173        unsafe { transmute(mt) }
174    }
175
176    /* get the next random value from the mersenne twister array */
177    fn mt_rand(mt: &mut [u32; Self::MT19937_RECURRENCE_DEGREE], cur_mt: &mut usize) -> u32 {
178        if *cur_mt == Self::MT19937_RECURRENCE_DEGREE {
179            /* our 624 sequence is finished. generate a new one */
180            for i in 0..Self::MT19937_RECURRENCE_DEGREE - 1 {
181                let y = (mt[i] & Self::MT19937_LOWER_MASK) | (mt[i + 1] & Self::MT19937_UPPER_MASK);
182                if y & 1 == 0 {
183                    /* even y */
184                    mt[i] = mt[(i + Self::MT19937_MIDDLE_WORD) % Self::MT19937_RECURRENCE_DEGREE]
185                        ^ (y >> 1);
186                } else {
187                    /* odd y */
188                    mt[i] = mt[(i + Self::MT19937_MIDDLE_WORD) % Self::MT19937_RECURRENCE_DEGREE]
189                        ^ (y >> 1)
190                        ^ Self::MT19937_RATIONAL_NORMAL_FORM_TWIST_MATRIX_COEFFICIENTS;
191                }
192            }
193
194            let y = (mt[Self::MT19937_RECURRENCE_DEGREE - 1] & Self::MT19937_LOWER_MASK)
195                | (mt[0] & Self::MT19937_UPPER_MASK);
196            if y & 1 == 0 {
197                /* even y */
198                mt[Self::MT19937_RECURRENCE_DEGREE - 1] =
199                    mt[Self::MT19937_MIDDLE_WORD - 1] ^ (y >> 1);
200            } else {
201                /* odd y */
202                mt[Self::MT19937_RECURRENCE_DEGREE - 1] = mt[Self::MT19937_MIDDLE_WORD - 1]
203                    ^ (y >> 1)
204                    ^ Self::MT19937_RATIONAL_NORMAL_FORM_TWIST_MATRIX_COEFFICIENTS;
205            }
206
207            *cur_mt = 0;
208        }
209
210        let mut y = mt[*cur_mt];
211        *cur_mt += 1;
212
213        y ^= y >> Self::MT19937_ADDITIONAL_TEMPERING.0;
214        y ^= (y << Self::MT19937_TGFSR_R_TEMPERING_BIT_SHIFTS.0)
215            & Self::MT19937_TGFSR_R_TEMPERING_BIT_MASKS.0;
216        y ^= (y << Self::MT19937_TGFSR_R_TEMPERING_BIT_SHIFTS.1)
217            & Self::MT19937_TGFSR_R_TEMPERING_BIT_MASKS.1;
218        y ^= y >> Self::MT19937_ADDITIONAL_TEMPERING.2;
219
220        y
221    }
222}
223
224impl std::fmt::Debug for MersenneTwister {
225    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> Result<(), std::fmt::Error> {
226        write!(f, "MersenneTwister {{ cur_mt: {} }}", self.cur_mt)
227    }
228}
229
230impl Algorithm for MersenneTwister {
231    fn get_int(&mut self) -> u32 {
232        Self::mt_rand(&mut self.mt, &mut self.cur_mt)
233    }
234}
235
236/// Complementary-Multiply-With-Carry algorithm.
237#[derive(Clone, Copy)]
238pub struct ComplementaryMultiplyWithCarry {
239    q: [u32; 4096],
240    c: u32,
241    cur: usize,
242}
243
244impl ComplementaryMultiplyWithCarry {
245    /// Create a new Complementary-Multiply-With-Carry algorithm instance.
246    #[allow(unsafe_code)]
247    pub fn new(seed: u32) -> Self {
248        let mut s = seed;
249        let mut q: [MaybeUninit<u32>; 4096] = unsafe { MaybeUninit::uninit().assume_init() };
250        for qe in &mut q[..] {
251            s = s.wrapping_mul(1_103_515_245).wrapping_add(12345); /* glibc LCG */
252            unsafe {
253                qe.as_mut_ptr().write(s);
254            }
255        }
256        let c = s.wrapping_mul(1_103_515_245).wrapping_add(12345) % 809_430_660; /* this max value is recommended by George Marsaglia */
257        let cur = 0;
258
259        Self {
260            q: unsafe { transmute(q) },
261            c,
262            cur,
263        }
264    }
265
266    fn get_number(&mut self) -> u32 {
267        self.cur = (self.cur + 1) & 4095;
268        let t = 18782_u64 * u64::from(self.q[self.cur]) + u64::from(self.c);
269        self.c = (t >> 32) as u32;
270        let mut x = (t + u64::from(self.c)) as u32;
271        if x < self.c {
272            x += 1;
273            self.c += 1;
274        }
275        if x.wrapping_add(1) == 0 {
276            self.c += 1;
277            x = 0;
278        }
279        self.q[self.cur] = 0xffff_fffe - x;
280
281        self.q[self.cur]
282    }
283}
284
285impl std::fmt::Debug for ComplementaryMultiplyWithCarry {
286    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> Result<(), std::fmt::Error> {
287        write!(
288            f,
289            "ComplementaryMultiplyWithCarry {{ c: {}, cur: {} }}",
290            self.c, self.cur
291        )
292    }
293}
294
295impl Algorithm for ComplementaryMultiplyWithCarry {
296    fn get_int(&mut self) -> u32 {
297        self.get_number()
298    }
299}
300
301struct Bits<'a, A: Algorithm + ?Sized> {
302    algorithm: &'a mut A,
303    bits: u32,
304    bits_left: u32,
305}
306
307impl<'a, A: Algorithm + ?Sized> Bits<'a, A> {
308    fn new(algorithm: &'a mut A) -> Self {
309        Self {
310            algorithm,
311            bits: 0,
312            bits_left: 0,
313        }
314    }
315
316    fn get_bit(&mut self) -> u32 {
317        if self.bits_left == 0 {
318            self.bits = self.algorithm.get_int();
319            self.bits_left = 32;
320        }
321
322        let bit = self.bits & 1;
323        self.bits >>= 1;
324        self.bits_left -= 1;
325        bit
326    }
327}
doryen_extra/random/algorithms.rs

doryen_extra/random/
algorithms.rs
