rand_float/pekkizen.rs
1// Derived from the uniFloats wiki https://github.com/pekkizen/prng (functions
2// Float64_64, Float64_117 and Float64full), distributed under the following
3// license:
4//
5// Copyright (c) 2020, Pekka Pulkkinen
6//
7// Copying and distribution of this article, ideas and source code are permitted
8// worldwide, without royalty, in any medium, provided the copyright notice is
9// preserved and a reference to this article is given.
10
11//! Pekka Pulkkinen’s leading-zeros technique.
12//!
13//! A Rust port of `Float64_64`, `Float64_117` and `Float64full` from [Pekka
14//! Pulkkinen’s uniFloats
15//! wiki](https://github.com/pekkizen/prng/wiki/uniFloats). One 64-bit word is
16//! interpreted as a uniform fixed-point real in [0 . . 1): the count of leading
17//! zeros picks the binade (a geometric distribution), and the remaining bits
18//! are shifted into the mantissa.
19//!
20//! The three variants differ in how far down the fixed point extends:
21//!
22//! - [`f64_64`] stops at the first word: a uniform real rounded down to the
23//! 2⁻⁶⁴ grid, complete (every float reachable) in [2⁻¹² . . 1);
24//! - [`f64_117`] draws a second word when the first has 12 or more leading
25//! zeros (probability 2⁻¹²), extending completeness to [2⁻⁶⁵ . . 1)
26//! with a 2⁻¹¹⁷ grid below;
27//! - [`f64_full`] additionally keeps drawing words while they are zero,
28//! reaching every float in [0 . . 1), subnormals included.
29//!
30//! All variants consume one word per call with probability 1 − 2⁻¹² and, in
31//! the typical case, cost only a couple of operations more than
32//! [`standard`](crate::standard) scaling.
33
34/// Returns a random `f64` distributed as a uniform 64-bit fixed-point real
35/// in [0 . . 1) rounded down to the nearest representable value: every float
36/// in [2⁻¹² . . 1), and the 2⁵² multiples of 2⁻⁶⁴ below 2⁻¹².
37#[inline]
38pub fn f64_64(mut bits: impl FnMut() -> u64) -> f64 {
39 let u = bits();
40 if u == 0 {
41 return 0.0;
42 }
43 let z = u.leading_zeros() as u64 + 1;
44 // The Go original computes `u << z` with z possibly 64, which Go defines
45 // as 0; Rust declares 64-bit shifts overflow, so the shift is split as
46 // `(u << (z - 1)) << 1` (z ≥ 1 always).
47 f64::from_bits((1023 - z) << 52 | ((u << (z - 1)) << 1) >> 12)
48}
49
50/// Returns 2⁻ⁿ as an `f64`; `n` must be at most 1022.
51#[inline]
52const fn two_to_minus(n: u64) -> f64 {
53 debug_assert!(n <= 1022);
54 f64::from_bits((1023 - n) << 52)
55}
56
57/// Returns a random `f64` distributed as a uniform 117-bit fixed-point real
58/// in [0 . . 1) rounded down to the nearest representable value: every
59/// float in [2⁻⁶⁵ . . 1), and the multiples of 2⁻¹¹⁷ below 2⁻⁶⁵.
60///
61/// A port of `Float64_117`: as [`f64_64`], except that when the first word
62/// has 12 or more leading zeros (probability 2⁻¹²) a second word extends the
63/// fixed point to 117 bits.
64#[inline]
65pub fn f64_117(mut bits: impl FnMut() -> u64) -> f64 {
66 let u = bits();
67 let z = u.leading_zeros() as u64 + 1;
68 if z <= 12 {
69 // 99.975% of cases.
70 return f64::from_bits((1023 - z) << 52 | ((u << (z - 1)) << 1) >> 12);
71 }
72 // Kluge; see [`crate::cold::cold_barrier`].
73 crate::cold::cold_barrier();
74 let z = z - 1;
75 // The Go original computes `u << z` with z possibly 64 (first word 0),
76 // which Go defines as 0; `checked_shl` + `unwrap_or` mirrors that
77 // (as would `unbounded_shl`, which however requires Rust 1.87).
78 let u = u.checked_shl(z as u32).unwrap_or(0) | bits() >> (64 - z);
79 (u >> 11) as f64 * two_to_minus(53 + z)
80}
81
82/// Returns a random `f64` distributed as a uniform real in [0 . . 1)
83/// rounded down to the nearest representable value: every float in
84/// [0 . . 1) is reachable, subnormals and 0 included, with probability
85/// equal to the measure of the reals that round down to it.
86///
87/// A port of `Float64full`: as [`f64_117`], except that zero words keep the
88/// zoom going, 64 binades at a time, down to the bottom of the subnormals.
89#[inline]
90pub fn f64_full(mut bits: impl FnMut() -> u64) -> f64 {
91 let mut u = bits();
92 let mut z = u.leading_zeros() as u64 + 1;
93 if z <= 12 {
94 // 99.975% of cases.
95 return f64::from_bits((1023 - z) << 52 | ((u << (z - 1)) << 1) >> 12);
96 }
97 // Kluge; see [`crate::cold::cold_barrier`].
98 crate::cold::cold_barrier();
99 z -= 1;
100 let mut exp = 0u64;
101 while u == 0 {
102 u = bits();
103 z = u.leading_zeros() as u64;
104 exp += 64;
105 if exp + z >= 1074 {
106 return 0.0;
107 }
108 }
109 // The Go original computes `bits() >> (64 - z)` with z possibly 0 after
110 // the loop, which Go defines as 0; `checked_shr` + `unwrap_or` mirrors
111 // that (as would `unbounded_shr`, which however requires Rust 1.87).
112 let u = u << z | bits().checked_shr(64 - z as u32).unwrap_or(0);
113 exp += z;
114 if exp < 1022 {
115 return f64::from_bits((1022 - exp) << 52 | (u << 1) >> 12);
116 }
117 // The 2⁵² subnormal floats.
118 f64::from_bits(u >> (exp - 1022) >> 12)
119}
120
121#[cfg(test)]
122mod tests {
123 use super::*;
124 use crate::sources::Weyl;
125
126 /// The bit-building form must agree with the wiki’s division form,
127 /// `float64(u << z >> 11) / 2^53 / 2^z` with z = leadingZeros(u).
128 #[test]
129 fn test_matches_reference_division_form() {
130 let division_form = |u: u64| {
131 let z = u.leading_zeros();
132 let m = ((u << z) >> 11) as f64;
133 m / (1u64 << 53) as f64 / (1u64 << z) as f64
134 };
135 let mut rng = Weyl(7);
136 for _ in 0..100_000 {
137 let u = rng.0.wrapping_add(0x9E3779B97F4A7C15);
138 assert_eq!(f64_64(|| rng.next_u64()), division_form(u));
139 }
140 }
141
142 #[test]
143 fn test_range() {
144 let mut rng = Weyl(42);
145 for _ in 0..100_000 {
146 let x = f64_64(|| rng.next_u64());
147 assert!((0.0..1.0).contains(&x), "f64_64: {x}");
148 }
149 }
150
151 /// Extremes of the leading-zeros count.
152 #[test]
153 fn test_edge_cases() {
154 assert_eq!(f64_64(|| 1 << 63), 0.5);
155 assert_eq!(f64_64(|| u64::MAX), f64::from_bits(1.0f64.to_bits() - 1));
156 assert_eq!(f64_64(|| 1), f64::from_bits((1023 - 64) << 52)); // 2^-64
157 assert_eq!(f64_64(|| 0), 0.0);
158 }
159 /// A source that replays a fixed sequence of words, then panics.
160 fn replay(words: &[u64]) -> impl FnMut() -> u64 + '_ {
161 let mut iter = words.iter();
162 move || *iter.next().expect("source exhausted")
163 }
164
165 /// Golden values produced by the verbatim Go code of the wiki, one per
166 /// branch (fast path, slow path, extreme leading-zero counts, zero
167 /// words).
168 #[test]
169 fn test_known_values_117() {
170 for (words, expected) in [
171 (&[0x8000000000000000, 0][..], 0x3fe0000000000000u64),
172 (&[0xDEADBEEFDEADBEEF, 0][..], 0x3febd5b7ddfbd5b7),
173 (&[1 << 52, 0][..], 0x3f30000000000000),
174 (&[1 << 51, 0xC0FFEE0DDF00D5ED][..], 0x3f20000000000001),
175 (&[1, 0xC0FFEE0DDF00D5ED][..], 0x3bfc0ffee0ddf00d),
176 (&[3, 0xFFFFFFFFFFFFFFFF][..], 0x3c0fffffffffffff),
177 (&[0, 0xC0FFEE0DDF00D5ED][..], 0x3be81ffdc1bbe01a),
178 (&[0, 0][..], 0),
179 ] {
180 assert_eq!(
181 f64_117(replay(words)).to_bits(),
182 expected,
183 "words {words:x?}"
184 );
185 }
186 }
187
188 /// Golden values produced by the verbatim Go code of the wiki; the deep
189 /// sequences exercise the zero-word zoom, the subnormal construction and
190 /// the all-zeros cutoff (which fires before the final exponent add).
191 #[test]
192 fn test_known_values_full() {
193 for (words, expected) in [
194 (vec![0x8000000000000000], 0x3fe0000000000000u64),
195 (vec![0xDEADBEEFDEADBEEF], 0x3febd5b7ddfbd5b7),
196 (vec![1 << 52], 0x3f30000000000000),
197 (vec![1 << 51, 0xC0FFEE0DDF00D5ED], 0x3f20000000000001),
198 (vec![1, 0xC0FFEE0DDF00D5ED], 0x3bfc0ffee0ddf00d),
199 (vec![3, 0xFFFFFFFFFFFFFFFF], 0x3c0fffffffffffff),
200 (vec![0, 0xC0FFEE0DDF00D5ED, 0xABCD], 0x3be81ffdc1bbe01a),
201 (vec![0, 0, 0xC0FFEE0DDF00D5ED, 0xAB], 0x37e81ffdc1bbe01a),
202 (
203 [vec![0; 15], vec![2, u64::MAX]].concat(),
204 0x000bffffffffffff,
205 ),
206 (
207 [vec![0; 15], vec![1, u64::MAX]].concat(),
208 0x0007ffffffffffff,
209 ),
210 (
211 [vec![0; 16], vec![1 << 50, 0xFFFFFFFFFFFFF]].concat(),
212 0x0000001000000000,
213 ),
214 ([vec![0; 16], vec![3, 0]].concat(), 0),
215 (vec![0; 18], 0),
216 ] {
217 assert_eq!(
218 f64_full(replay(&words)).to_bits(),
219 expected,
220 "words {words:x?}"
221 );
222 }
223 }
224
225 /// XOR-rotate hash over 10⁶ draws, cross-checked against the verbatim Go
226 /// code of the wiki driven by the same Weyl source. The two variants
227 /// agree on any zero-free source stream.
228 #[test]
229 fn test_matches_go_reference_hash() {
230 let mut rng = Weyl(42);
231 let mut h = 0u64;
232 for _ in 0..1_000_000 {
233 h = h.rotate_left(1) ^ f64_117(|| rng.next_u64()).to_bits();
234 }
235 assert_eq!(h, 0xf9b6db6017240be7);
236
237 let mut rng = Weyl(42);
238 let mut h = 0u64;
239 for _ in 0..1_000_000 {
240 h = h.rotate_left(1) ^ f64_full(|| rng.next_u64()).to_bits();
241 }
242 assert_eq!(h, 0xf9b6db6017240be7);
243 }
244
245 #[test]
246 fn test_range_117_full() {
247 let mut rng = Weyl(42);
248 for _ in 0..100_000 {
249 let x = f64_117(|| rng.next_u64());
250 assert!((0.0..1.0).contains(&x), "f64_117: {x}");
251 let x = f64_full(|| rng.next_u64());
252 assert!((0.0..1.0).contains(&x), "f64_full: {x}");
253 }
254 }
255}