wowlab_engine_rng/stochastic/
prd.rs1const PRD_MAX_ATTEMPTS: u32 = 100_000;
2const PRD_NEWTON_ITERATIONS: u8 = 20;
3const PRD_PRECISION: f64 = 1e-12;
4const PRD_INITIAL_GUESS_SCALE: f64 = 1.25;
5const PRD_INITIAL_GUESS_CAP: f64 = 0.99;
6
7#[must_use]
9pub fn prd_constant(average_proc_rate: f64) -> f64 {
10 prd_constant_capped(average_proc_rate, 0)
11}
12
13#[must_use]
15pub fn prd_constant_capped(average_proc_rate: f64, cap: u32) -> f64 {
16 let rate = average_proc_rate.clamp(0.0, 1.0);
17 let cap = if cap == 0 {
18 PRD_MAX_ATTEMPTS
19 } else {
20 cap.min(PRD_MAX_ATTEMPTS)
21 };
22
23 if rate <= 0.0 {
24 return 0.0;
25 }
26
27 if rate <= 1.0 / f64::from(cap) {
28 return f64::MIN_POSITIVE;
29 }
30
31 if rate >= 1.0 {
32 return 1.0;
33 }
34
35 let target_attempts = 1.0 / rate;
36 let mut constant = (PRD_INITIAL_GUESS_SCALE * rate * rate).min(PRD_INITIAL_GUESS_CAP);
37
38 for _ in 0..PRD_NEWTON_ITERATIONS {
39 let (expected, derivative) = prd_expected_attempts(constant, cap);
40 let error = expected - target_attempts;
41
42 if error.abs() < PRD_PRECISION {
43 break;
44 }
45
46 constant -= error / derivative;
47 }
48
49 constant
50}
51
52fn prd_expected_attempts(constant: f64, cap: u32) -> (f64, f64) {
53 let mut chain = 1.0;
54 let mut chain_derivative = 0.0;
55 let mut expected = 0.0;
56 let mut expected_derivative = 0.0;
57
58 for attempt in 1..=cap {
59 expected += chain;
60 expected_derivative += chain_derivative;
61
62 let attempt = f64::from(attempt);
63 let chance = attempt * constant;
64
65 if chance >= 1.0 {
66 break;
67 }
68
69 chain_derivative = chain_derivative * (1.0 - chance) - chain * attempt;
70 chain *= 1.0 - chance;
71 }
72
73 (expected, expected_derivative)
74}
75
76#[cfg(test)]
77mod tests {
78 use googletest::prelude::*;
79 use rstest::rstest;
80
81 use super::*;
82
83 #[gtest]
84 fn constant_preserves_average_proc_rate() -> Result<()> {
85 let constant = prd_constant(0.15);
86
87 verify_that!(constant, near(0.032_220_914, 1e-9))?;
88 let (expected_attempts, _) = prd_expected_attempts(constant, PRD_MAX_ATTEMPTS);
89
90 verify_that!(expected_attempts.recip(), near(0.15, 1e-12))
91 }
92
93 #[gtest]
94 fn capped_constant_preserves_average_proc_rate() -> Result<()> {
95 let constant = prd_constant_capped(0.12, 13);
96 let (expected_attempts, _) = prd_expected_attempts(constant, 13);
97
98 verify_that!(expected_attempts.recip(), near(0.12, 1e-12))?;
99
100 verify_that!(constant, lt(prd_constant(0.12)))
101 }
102
103 #[gtest]
104 #[rstest]
105 #[case::never(0.0, 0.0)]
106 #[case::always(1.0, 1.0)]
107 fn boundary_contract(#[case] average_proc_rate: f64, #[case] expected: f64) -> Result<()> {
108 verify_that!(prd_constant(average_proc_rate), eq(expected))
109 }
110
111 #[gtest]
112 fn capped_unrepresentable_rate_returns_smallest_positive_value() -> Result<()> {
113 verify_that!(
114 prd_constant_capped(0.01, 10).to_bits(),
115 eq(f64::MIN_POSITIVE.to_bits())
116 )
117 }
118}