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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
//! Random real module - port of Modules/random_real.c
//!
//! Provides high-quality floating-point random numbers.
use crate::random;
/// Generate a random double in [0, 1] uniformly distributed.
/// Direct port of `random_real` from src/zsh/Src/Modules/random_real.c
/// lines 145-212 (Campbell's algorithm, see
/// <http://mumble.net/~campbell/2014/04/28/uniform-random-float>).
///
/// Algorithm:
/// 1. Read 64-bit chunks until we see a non-zero one. Each
/// all-zero chunk shifts the exponent down by 64.
/// 2. Count leading zeros in the first non-zero chunk; shift
/// them into the exponent and refill the low bits of the
/// significand from another 64-bit draw.
/// 3. Set the sticky bit (significand |= 1) so the round-to-
/// nearest-even doesn't bias toward even when the trailing
/// bits would have decided ties.
/// 4. ldexp(significand, exponent) — convert to double.
///
/// This is the only correct way to generate a uniform double
/// from a uniform bit source. The naïve "53-bit fraction" approach
/// (random_real_53 below) skews ~3% of the [0, 1) interval toward
/// values just below 0.5.
pub fn random_real() -> f64 {
let mut exponent: i32 = 0;
let mut significand: u64 = 0;
// random_real.c:158-175 — read zeros into exponent until
// we hit a non-zero chunk.
while significand == 0 {
exponent -= 64;
significand = random::get_random_u64();
// random_real.c:172-174 — exp below -1074 means it would
// round to zero anyway (smallest subnormal exponent).
if exponent < -1074 {
return 0.0;
}
}
// random_real.c:185-196 — leading-zero shift.
let shift = significand.leading_zeros() as i32;
if shift != 0 {
let r = random::get_random_u64();
exponent -= shift;
significand <<= shift;
significand |= r >> (64 - shift);
}
// random_real.c:205 — sticky bit so round-to-nearest doesn't
// false-tie.
significand |= 1;
// random_real.c:211 — ldexp(significand, exponent).
(significand as f64) * (exponent as f64).exp2()
}
/// Generate a random double in [0, max)
pub fn random_real_max(max: f64) -> f64 {
random_real() * max
}
/// Generate a random double in [min, max)
pub fn random_real_range(min: f64, max: f64) -> f64 {
min + random_real() * (max - min)
}
/// Generate high-precision random in [0, 1) using 53 bits.
/// Slightly biased compared to `random_real`; kept for callers
/// that need deterministic 53-bit precision.
pub fn random_real_53() -> f64 {
let a = random::get_random_u32() >> 5;
let b = random::get_random_u32() >> 6;
(a as f64 * 67108864.0 + b as f64) * (1.0 / 9007199254740992.0)
}
/// Math function for random real
pub fn math_random_real(args: &[f64]) -> Result<f64, String> {
match args.len() {
0 => Ok(random_real_53()),
1 => {
let max = args[0];
if max <= 0.0 {
return Err("random: max must be positive".to_string());
}
Ok(random_real_max(max))
}
2 => {
let min = args[0];
let max = args[1];
if max <= min {
return Err("random: max must be greater than min".to_string());
}
Ok(random_real_range(min, max))
}
_ => Err("random: too many arguments".to_string()),
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_random_real_range() {
for _ in 0..100 {
let r = random_real();
assert!(r >= 0.0 && r < 1.0);
}
}
#[test]
fn test_random_real_max() {
for _ in 0..100 {
let r = random_real_max(10.0);
assert!(r >= 0.0 && r < 10.0);
}
}
#[test]
fn test_random_real_min_max() {
for _ in 0..100 {
let r = random_real_range(5.0, 10.0);
assert!(r >= 5.0 && r < 10.0);
}
}
#[test]
fn test_random_real_53() {
for _ in 0..100 {
let r = random_real_53();
assert!(r >= 0.0 && r < 1.0);
}
}
#[test]
fn test_math_random_real_no_args() {
let result = math_random_real(&[]);
assert!(result.is_ok());
let r = result.unwrap();
assert!(r >= 0.0 && r < 1.0);
}
#[test]
fn test_math_random_real_one_arg() {
let result = math_random_real(&[100.0]);
assert!(result.is_ok());
let r = result.unwrap();
assert!(r >= 0.0 && r < 100.0);
}
#[test]
fn test_math_random_real_two_args() {
let result = math_random_real(&[10.0, 20.0]);
assert!(result.is_ok());
let r = result.unwrap();
assert!(r >= 10.0 && r < 20.0);
}
#[test]
fn test_math_random_real_invalid() {
assert!(math_random_real(&[-1.0]).is_err());
assert!(math_random_real(&[10.0, 5.0]).is_err());
assert!(math_random_real(&[1.0, 2.0, 3.0]).is_err());
}
}