nd_300/speedtest/
stat_primitives.rs

1//! Pinned statistical primitives for SpeedQX Methodology v4.
2//!
3//! These primitives are the portability contract between the TypeScript
4//! (website + app WebView) and Rust (CLI) implementations. They MUST behave
5//! identically across languages so the committed `golden-vectors.json` fixture
6//! asserts parity:
7//!
8//!   1. [`quantile`]            — type-7 linear interpolation (BIT-EXACT)
9//!   2. [`Pcg32`] + [`lemire_bounded`] — deterministic PRNG index streams (BIT-EXACT)
10//!   3. [`inv_normal`] / [`phi`] — Acklam / West rational approximations (1e-9 rel)
11//!   4. [`t975`]                — hardcoded Student-t 0.975 quantiles (BIT-EXACT)
12//!   5. [`sum`] / [`sample_mean`] / [`sample_variance`] — fixed-order summation
13//!
14//! FP discipline (governing rule): fixed left-to-right summation order, no FMA /
15//! fast-math (never `mul_add`), `f64` everywhere, and `u64` wrapping arithmetic
16//! for the PCG32 state. This mirrors `stat-primitives.ts` in the SpeedQX web
17//! repo byte-for-byte in behavior. See METHODOLOGY.md §11.
18
19// ── type-7 quantile ─────────────────────────────────────────────────────
20
21/// Type-7 (linear-interpolation) quantile on a pre-sorted ascending slice.
22///
23/// `h = (n - 1) * p`, interpolate between the two bracketing order statistics.
24/// BIT-EXACT primitive — identical formula in TS and Rust. `p` must be in
25/// `[0, 1]`. Returns `NaN` for an empty slice (mirrors the TS reference).
26pub fn quantile(sorted: &[f64], p: f64) -> f64 {
27    let n = sorted.len();
28    if n == 0 {
29        return f64::NAN;
30    }
31    if n == 1 {
32        return sorted[0];
33    }
34    let h = (n - 1) as f64 * p;
35    let lo = h.floor() as usize;
36    let hi = h.ceil() as usize;
37    sorted[lo] + (h - lo as f64) * (sorted[hi] - sorted[lo])
38}
39
40// ── PCG32 + Lemire bounded index ─────────────────────────────────────────
41
42/// PCG32 LCG multiplier (Melissa O'Neill's reference constant).
43pub const PCG32_MULT: u64 = 6364136223846793005;
44/// Default initial state (PCG reference `PCG32_INITIALIZER`). Pinned for golden vectors.
45pub const PCG32_DEFAULT_STATE: u64 = 0x853c49e6748fea9b;
46/// Default stream increment (odd). Pinned for golden vectors.
47pub const PCG32_DEFAULT_INC: u64 = 0xda3e39cb94b95bdb;
48
49/// Lemire multiply-shift bounded index: `floor(u32 * n / 2^32)`.
50///
51/// Uses the full 64-bit product then shifts right 32. Deterministic across TS
52/// and Rust. Returns a value in `[0, n)`.
53pub fn lemire_bounded(u32_val: u32, n: usize) -> usize {
54    ((u32_val as u64 * n as u64) >> 32) as usize
55}
56
57/// PCG32 (permuted congruential generator, XSH-RR 64/32 variant).
58///
59/// The 64-bit state advances with wrapping `u64` arithmetic; the output
60/// permutation (xorshift + rotate) is truncated to 32 bits, exactly mirroring
61/// the C reference and the TS BigInt implementation. `inc` must be odd.
62#[derive(Debug, Clone)]
63pub struct Pcg32 {
64    state: u64,
65    inc: u64,
66}
67
68impl Pcg32 {
69    /// Construct with the pinned default state + increment (golden-vector stream).
70    pub fn new() -> Self {
71        Self {
72            state: PCG32_DEFAULT_STATE,
73            inc: PCG32_DEFAULT_INC,
74        }
75    }
76
77    /// Construct with an explicit state + increment. `inc` is used verbatim
78    /// (the caller is responsible for it being odd, as in the reference).
79    pub fn with_seed(state: u64, inc: u64) -> Self {
80        Self { state, inc }
81    }
82
83    /// Advance the state and return the next uniformly-distributed `u32`.
84    pub fn next_u32(&mut self) -> u32 {
85        let old = self.state;
86        self.state = old.wrapping_mul(PCG32_MULT).wrapping_add(self.inc);
87        let xorshifted = ((((old >> 18) ^ old) >> 27) & 0xffff_ffff) as u32;
88        let rot = (old >> 59) as u32; // 0..=31
89        xorshifted.rotate_right(rot)
90    }
91
92    /// Lemire bounded index into `[0, n)` drawn from the next `u32`.
93    pub fn bounded_index(&mut self, n: usize) -> usize {
94        lemire_bounded(self.next_u32(), n)
95    }
96
97    /// Next double in `[0, 1)`: `next_u32() / 2^32`. Convenience for input synthesis.
98    pub fn next_f64(&mut self) -> f64 {
99        self.next_u32() as f64 / 4_294_967_296.0
100    }
101}
102
103impl Default for Pcg32 {
104    fn default() -> Self {
105        Self::new()
106    }
107}
108
109// ── invNormal (Acklam) + phi (West) ──────────────────────────────────────
110
111/// Inverse standard-normal CDF via Peter Acklam's rational approximation.
112///
113/// Max relative error ≈ 1.15e-9 vs the true quantile (no Halley refinement, so
114/// TS and Rust stay identical from the same constants). Returns ±∞ at the open
115/// boundaries. 1e-9-relative primitive.
116#[allow(clippy::excessive_precision)]
117pub fn inv_normal(p: f64) -> f64 {
118    if p <= 0.0 {
119        return f64::NEG_INFINITY;
120    }
121    if p >= 1.0 {
122        return f64::INFINITY;
123    }
124
125    let a1 = -3.969683028665376e1;
126    let a2 = 2.209460984245205e2;
127    let a3 = -2.759285104469687e2;
128    let a4 = 1.38357751867269e2;
129    let a5 = -3.066479806614716e1;
130    let a6 = 2.506628277459239e0;
131    let b1 = -5.447609879822406e1;
132    let b2 = 1.615858368580409e2;
133    let b3 = -1.556989798598866e2;
134    let b4 = 6.680131188771972e1;
135    let b5 = -1.328068155288572e1;
136    let c1 = -7.784894002430293e-3;
137    let c2 = -3.223964580411365e-1;
138    let c3 = -2.400758277161838e0;
139    let c4 = -2.549732539343734e0;
140    let c5 = 4.374664141464968e0;
141    let c6 = 2.938163982698783e0;
142    let d1 = 7.784695709041462e-3;
143    let d2 = 3.224671290700398e-1;
144    let d3 = 2.445134137142996e0;
145    let d4 = 3.754408661907416e0;
146
147    let p_low = 0.02425;
148    let p_high = 1.0 - p_low;
149
150    if p < p_low {
151        let q = (-2.0 * p.ln()).sqrt();
152        (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6)
153            / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1.0)
154    } else if p <= p_high {
155        let q = p - 0.5;
156        let r = q * q;
157        (((((a1 * r + a2) * r + a3) * r + a4) * r + a5) * r + a6) * q
158            / (((((b1 * r + b2) * r + b3) * r + b4) * r + b5) * r + 1.0)
159    } else {
160        let q = (-2.0 * (1.0 - p).ln()).sqrt();
161        -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6)
162            / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1.0)
163    }
164}
165
166/// Standard-normal CDF Φ(x) via Graeme West's `cumnorm` (Hart-style rational
167/// approximation), accurate to ~1e-15 using only `exp` + arithmetic.
168/// Φ(0) = 0.5 exactly. 1e-9-relative primitive; the companion of [`inv_normal`]
169/// for BCa.
170#[allow(clippy::excessive_precision)]
171pub fn phi(x: f64) -> f64 {
172    let xa = x.abs();
173    if xa > 37.0 {
174        return if x > 0.0 { 1.0 } else { 0.0 };
175    }
176    let e = (-xa * xa / 2.0).exp();
177    let c: f64;
178    if xa < 7.07106781186547 {
179        let mut b = 3.52624965998911e-2 * xa + 0.700383064443688;
180        b = b * xa + 6.37396220353165;
181        b = b * xa + 33.912866078383;
182        b = b * xa + 112.079291497871;
183        b = b * xa + 221.213596169931;
184        b = b * xa + 220.206867912376;
185        let num = e * b;
186        b = 8.83883476483184e-2 * xa + 1.75566716318264;
187        b = b * xa + 16.064177579207;
188        b = b * xa + 86.7807322029461;
189        b = b * xa + 296.564248779674;
190        b = b * xa + 637.333633378831;
191        b = b * xa + 793.826512519948;
192        b = b * xa + 440.413735824752;
193        c = num / b;
194    } else {
195        let mut b = xa + 0.65;
196        b = xa + 4.0 / b;
197        b = xa + 3.0 / b;
198        b = xa + 2.0 / b;
199        b = xa + 1.0 / b;
200        c = e / b / 2.506628274631;
201    }
202    if x > 0.0 {
203        1.0 - c
204    } else {
205        c
206    }
207}
208
209// ── Student-t 0.975 table ────────────────────────────────────────────────
210
211/// Student-t 0.975 quantile for `df` degrees of freedom.
212///
213/// The canonical table (METHODOLOGY.md §6) is pinned for `df ∈ [1, 7]`, which
214/// covers every reachable provider count. `df ≤ 1` returns the `df = 1` value;
215/// `df ≥ 7` clamps to `df = 7` (2.365) — conservative (wider CI) and
216/// unreachable with the current provider registry. BIT-EXACT.
217pub fn t975(df: usize) -> f64 {
218    match df {
219        0 | 1 => 12.706,
220        2 => 4.303,
221        3 => 3.182,
222        4 => 2.776,
223        5 => 2.571,
224        6 => 2.447,
225        _ => 2.365,
226    }
227}
228
229// ── fixed-order sum helpers ──────────────────────────────────────────────
230
231/// Deterministic left-to-right sum. No pairwise / Kahan reassociation — the
232/// summation order is part of the golden-vector contract.
233pub fn sum(xs: &[f64]) -> f64 {
234    let mut s = 0.0;
235    for &x in xs {
236        s += x;
237    }
238    s
239}
240
241/// Arithmetic mean via fixed-order [`sum`]. Returns 0 for the empty slice.
242pub fn sample_mean(xs: &[f64]) -> f64 {
243    if xs.is_empty() {
244        return 0.0;
245    }
246    sum(xs) / xs.len() as f64
247}
248
249/// Unbiased sample variance (Bessel's `n-1` correction), fixed summation order.
250pub fn sample_variance(xs: &[f64]) -> f64 {
251    let n = xs.len();
252    if n < 2 {
253        return 0.0;
254    }
255    let m = sample_mean(xs);
256    let mut s = 0.0;
257    for &x in xs {
258        let d = x - m;
259        s += d * d;
260    }
261    s / (n - 1) as f64
262}
263
264#[cfg(test)]
265mod tests {
266    use super::*;
267
268    #[test]
269    fn quantile_hand_computed() {
270        assert_eq!(quantile(&[1.0, 2.0, 3.0, 4.0], 0.5), 2.5);
271        assert_eq!(quantile(&[1.0, 2.0, 3.0, 4.0], 0.25), 1.75);
272        assert_eq!(quantile(&[1.0, 2.0, 3.0, 4.0], 0.0), 1.0);
273        assert_eq!(quantile(&[1.0, 2.0, 3.0, 4.0], 1.0), 4.0);
274        assert_eq!(quantile(&[7.0], 0.5), 7.0);
275        assert!(quantile(&[], 0.5).is_nan());
276    }
277
278    #[test]
279    fn pcg32_first_outputs() {
280        let mut r = Pcg32::new();
281        assert_eq!(
282            [r.next_u32(), r.next_u32(), r.next_u32(), r.next_u32()],
283            [355248013, 41705475, 3406281715, 4186697710]
284        );
285    }
286
287    #[test]
288    fn lemire_bounds() {
289        assert_eq!(lemire_bounded(0, 100), 0);
290        assert_eq!(lemire_bounded(0xffff_ffff, 100), 99);
291        assert_eq!(lemire_bounded(0x8000_0000, 10), 5);
292    }
293
294    #[test]
295    fn inv_normal_and_phi_anchors() {
296        assert_eq!(inv_normal(0.5), 0.0);
297        assert_eq!(phi(0.0), 0.5);
298        assert_eq!(inv_normal(0.0), f64::NEG_INFINITY);
299        assert_eq!(inv_normal(1.0), f64::INFINITY);
300        assert!((inv_normal(0.975) - 1.959963984540054).abs() < 1e-8);
301        assert!((phi(1.959963984540054) - 0.975).abs() < 1e-9);
302    }
303
304    #[test]
305    fn t975_table_and_clamps() {
306        assert_eq!(t975(1), 12.706);
307        assert_eq!(t975(0), 12.706);
308        assert_eq!(t975(2), 4.303);
309        assert_eq!(t975(7), 2.365);
310        assert_eq!(t975(99), 2.365);
311    }
312
313    #[test]
314    fn fixed_order_sums() {
315        let xs = [2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0];
316        assert_eq!(sum(&xs), 40.0);
317        assert_eq!(sample_mean(&xs), 5.0);
318        assert_eq!(sample_variance(&xs), 32.0 / 7.0);
319        // The famous non-associative residue must be preserved.
320        assert_eq!(sum(&[0.1, 0.1, 0.1]), 0.30000000000000004);
321    }
322}
nd_300/speedtest/stat_primitives.rs

nd_300/speedtest/
stat_primitives.rs
