alice_analytics/

privacy.rs

//! Local Differential Privacy (LDP) Mechanisms
//!
//! Privacy-preserving data collection where noise is added at the source.
//! Individual data points are deniable, but aggregate statistics emerge.

#[cfg(not(feature = "std"))]
use crate::math::FloatExt;
use crate::math::{f64_i64, hash_index, i64_f64, u64_f64};
use alice_det_math::ln64;
// the only caller (`RandomizedResponse::new`) needs system entropy
#[cfg(feature = "std")]
use alice_det_math::exp64;

// ============================================================================
// Random Number Generation (ChaCha20-based for determinism)
// ============================================================================

/// Simple xorshift64 PRNG for fast random numbers
///
/// Not cryptographically secure, but fast and sufficient for noise injection.
#[derive(Clone, Debug)]
pub struct XorShift64 {
    state: u64,
}

impl XorShift64 {
    /// Create a new PRNG with given seed
    #[inline]
    #[must_use]
    pub const fn new(seed: u64) -> Self {
        // Ensure non-zero state
        Self {
            state: if seed == 0 {
                0x853c_49e6_748f_ea9b
            } else {
                seed
            },
        }
    }

    /// Create from system entropy (`std` only: `SystemTime`)
    ///
    /// `no_std` には entropy source が無いので存在しない (2026-09-15 まで固定 seed を返す
    /// stub があり、差分プライバシーの noise が決定論になっていた) — `no_std` では
    /// `XorShift64::new(seed)` / `*::with_seed` / `with_probability` / `with_seed_params` に
    /// caller が entropy を渡す
    #[cfg(feature = "std")]
    #[must_use]
    pub fn from_entropy() -> Self {
        use std::time::{SystemTime, UNIX_EPOCH};
        let seed = SystemTime::now()
            .duration_since(UNIX_EPOCH)
            // 下位 64 bit で十分 (seed 用途、u128 → u64 の切り捨ては意図)
            .map_or(0x853c_49e6_748f_ea9b, |d| {
                #[allow(clippy::cast_possible_truncation)]
                let low = d.as_nanos() as u64;
                low
            });
        Self::new(seed)
    }

    /// Generate next u64
    #[inline]
    pub const fn next_u64(&mut self) -> u64 {
        let mut x = self.state;
        x ^= x << 13;
        x ^= x >> 7;
        x ^= x << 17;
        self.state = x;
        x
    }

    /// Generate uniform f64 in [0, 1)
    #[inline(always)]
    pub fn next_f64(&mut self) -> f64 {
        // Multiply by 1/2^53 instead of dividing
        const RCP_POW2_53: f64 = 1.0 / u64_f64(1u64 << 53);
        u64_f64(self.next_u64() >> 11) * RCP_POW2_53
    }

    /// Generate uniform f64 in [low, high)
    #[inline]
    pub fn next_f64_range(&mut self, low: f64, high: f64) -> f64 {
        (high - low).mul_add(self.next_f64(), low)
    }

    /// Generate a random boolean with given probability of true
    #[inline]
    pub fn next_bool(&mut self, p: f64) -> bool {
        self.next_f64() < p
    }
}

#[cfg(feature = "std")]
impl Default for XorShift64 {
    fn default() -> Self {
        Self::from_entropy()
    }
}

// ============================================================================
// Laplace Distribution
// ============================================================================

/// Generate Laplace-distributed random variable
///
/// Laplace(μ, b) has PDF: f(x) = (1/2b) * exp(-|x-μ|/b)
///
/// Used for ε-differential privacy with sensitivity Δf:
/// noise scale b = Δf / ε
#[derive(Clone, Debug)]
pub struct LaplaceNoise {
    /// Scale parameter b
    scale: f64,
    /// PRNG
    rng: XorShift64,
}

impl LaplaceNoise {
    /// Create a new Laplace noise generator
    ///
    /// # Arguments
    /// * `sensitivity` - Maximum change in output for one input change (Δf)
    /// * `epsilon` - Privacy parameter ε (smaller = more privacy)
    ///
    /// `std` only (system entropy); `no_std` は [`Self::with_seed`]
    #[cfg(feature = "std")]
    #[must_use]
    pub fn new(sensitivity: f64, epsilon: f64) -> Self {
        let scale = sensitivity / epsilon;
        Self {
            scale,
            rng: XorShift64::from_entropy(),
        }
    }

    /// Create with explicit seed
    #[must_use]
    pub fn with_seed(sensitivity: f64, epsilon: f64, seed: u64) -> Self {
        let scale = sensitivity / epsilon;
        Self {
            scale,
            rng: XorShift64::new(seed),
        }
    }

    /// Generate Laplace noise
    #[inline]
    pub fn sample(&mut self) -> f64 {
        // Inverse transform sampling: X = μ - b * sign(U - 0.5) * ln(1 - 2|U - 0.5|)
        let u = self.rng.next_f64() - 0.5;
        let sign = if u < 0.0 { -1.0 } else { 1.0 };
        -sign * self.scale * ln64(2.0f64.mul_add(-u.abs(), 1.0))
    }

    /// Add noise to a value
    #[inline]
    pub fn privatize(&mut self, value: f64) -> f64 {
        value + self.sample()
    }

    /// Add noise to an integer value (rounds result)
    #[inline]
    pub fn privatize_int(&mut self, value: i64) -> i64 {
        f64_i64((i64_f64(value) + self.sample()).round())
    }

    /// Get the scale parameter
    #[inline]
    #[must_use]
    pub const fn scale(&self) -> f64 {
        self.scale
    }
}

// ============================================================================
// Randomized Response (RAPPOR-style)
// ============================================================================

/// Randomized Response for binary data
///
/// Classic technique: to report a sensitive bit b:
/// - With probability p, report b truthfully
/// - With probability 1-p, report a random bit
///
/// This provides ε-differential privacy where ε = ln((p + 0.5(1-p)) / (0.5(1-p)))
#[derive(Clone, Debug)]
pub struct RandomizedResponse {
    /// Probability of truthful response
    p_true: f64,
    /// PRNG
    rng: XorShift64,
}

impl RandomizedResponse {
    /// Create from privacy parameter epsilon
    ///
    /// Higher epsilon = more accuracy, less privacy
    ///
    /// `std` only (system entropy); `no_std` は [`Self::with_probability`]
    ///
    /// # Panics
    ///
    /// `epsilon` が有限でない、または負の場合 検査するのは入力の `epsilon` で、
    /// `p_true` はそこからの導出値 負の ε は `p_true < 0.5` を作り、
    /// [`Self::estimate_proportion`] の推定値の符号が反転する
    /// (実測: `epsilon = -5.0` で `p_true = 6.69e-3`) ⚠️ 2026-10-07 までは
    /// [`Self::with_probability`] だけが定義域を検査していたので、`new` 経由だと
    /// 同じ不変条件を破れる状態だった
    ///
    /// `epsilon = 0.0` は許す (`p_true = 0.5` = 常に無作為に答える、
    /// `with_probability(0.5)` と同じ意味のある端)
    ///
    /// 上端は切らない: `p = 1/(1 + e^(−ε))` (logistic の数値安定形) で計算するので
    /// `ε → ∞` の極限 `1.0` をそのまま返す ⚠️ 旧式の `e^ε/(1 + e^ε)` は
    /// `e^ε` が `ε ≳ 709.79` で `inf` になり `inf/inf = NaN` を返していた (実測)
    /// NaN は**代数的な書き方の副産物で定義域の問題ではない**ので、ε の上限を
    /// assert で切るのでなく式を直した (ε = 0 / 0.5 / 1 / 50 / 709 では旧式と
    /// bit 一致、ε = 3 で 1 ulp 異なる)
    #[cfg(feature = "std")]
    #[must_use]
    pub fn new(epsilon: f64) -> Self {
        assert!(
            epsilon.is_finite() && epsilon >= 0.0,
            "epsilon must be finite and non-negative"
        );
        // p = e^ε / (1 + e^ε) = 1 / (1 + e^(−ε)): for ε ≥ 0 the exponential
        // cannot overflow in this form, so the whole domain is representable
        let p_true = 1.0 / (1.0 + exp64(-epsilon));
        Self {
            p_true,
            rng: XorShift64::from_entropy(),
        }
    }

    /// Create with explicit probability and seed
    ///
    /// # Panics
    ///
    /// `p_true` が `[0.5, 1.0]` の外 (`NaN` / 無限大を含む、どちらも比較が
    /// 偽になるので同じ assert で弾かれる) の場合
    ///
    /// 定義域の扱いは同 crate の [`ExponentialMovingAverage::new`] に合わせた
    /// (範囲外は黙って丸めず panic する) 2026-10-07 までは
    /// `clamp(0.5, 1.0)` だったので、`0.0` を渡すと黙って `0.5` に、`2.0` は
    /// `1.0` になり、**`NaN` は `clamp` が `NaN` を返すので丸められずに残り**、
    /// `p_true()` が `NaN` を報告する機構ができていた 丸めは呼び出し側が
    /// 頼んだのと違う ε の機構を作るので、推定値の意味が変わる
    ///
    /// `0.5` は「常に無作為に答える」という意味のある下端なので許す
    /// (`new(0.0)` がちょうど `0.5` を作る)
    ///
    /// [`ExponentialMovingAverage::new`]: crate::streaming_ops::ExponentialMovingAverage::new
    #[must_use]
    // `RangeInclusive::contains` is not a `const fn`, so the range cannot be
    // spelled the way clippy prefers without dropping `const` from a public
    // constructor.
    #[allow(clippy::manual_range_contains)]
    pub const fn with_probability(p_true: f64, seed: u64) -> Self {
        assert!(
            p_true >= 0.5 && p_true <= 1.0,
            "p_true must be in [0.5, 1.0]"
        );
        Self {
            p_true,
            rng: XorShift64::new(seed),
        }
    }

    /// Privatize a boolean value
    #[inline]
    pub fn privatize(&mut self, value: bool) -> bool {
        if self.rng.next_bool(self.p_true) {
            // Report truthfully
            value
        } else {
            // Report random
            self.rng.next_bool(0.5)
        }
    }

    /// Privatize a bit (0 or 1)
    #[inline]
    pub fn privatize_bit(&mut self, bit: u8) -> u8 {
        u8::from(self.privatize(bit != 0))
    }

    /// Get the probability of truthful response
    #[inline]
    #[must_use]
    pub const fn p_true(&self) -> f64 {
        self.p_true
    }

    /// Estimate true proportion from noisy counts
    ///
    /// Given N total responses with K positive responses,
    /// estimate the true proportion of positive values.
    #[must_use]
    pub fn estimate_proportion(p_true: f64, n: u64, k: u64) -> f64 {
        if n == 0 {
            return 0.0;
        }
        let inv_n = 1.0 / u64_f64(n);
        let observed_rate = u64_f64(k) * inv_n;
        // Debiasing: true_rate = (observed_rate - 0.5*(1-p)) / (p - 0.5*(1-p))
        // Simplifies to: true_rate = (observed_rate - 0.5 + 0.5*p) / (p - 0.5 + 0.5*p)
        //              = (2*observed_rate - 1 + p) / (2p - 1 + p)
        //              = (2*observed_rate - 1 + p) / (3p - 1)
        // Wait, let me recalculate...
        // P(report=1) = p*true_rate + (1-p)*0.5
        // observed_rate = p*true_rate + 0.5 - 0.5*p
        // true_rate = (observed_rate - 0.5 + 0.5*p) / p
        let inv_p_true = 1.0 / p_true;
        let true_rate = 0.5f64.mul_add(p_true, observed_rate - 0.5) * inv_p_true;
        true_rate.clamp(0.0, 1.0)
    }
}

// ============================================================================
// RAPPOR (Randomized Aggregatable Privacy-Preserving Ordinal Response)
// ============================================================================

/// Default RAPPOR bit size
pub const RAPPOR_BITS: usize = 64;

/// RAPPOR for categorical data with multiple bits
///
/// Encodes categorical values into a Bloom filter, then applies
/// randomized response to each bit.
#[derive(Clone, Debug)]
pub struct Rappor {
    /// Permanent randomized response (for longitudinal studies)
    f: f64,
    /// Instantaneous randomized response parameters
    p: f64,
    q: f64,
    /// PRNG
    rng: XorShift64,
}

impl Rappor {
    /// Number of bits in the Bloom filter
    pub const BITS: usize = RAPPOR_BITS;

    /// Create a new RAPPOR encoder
    ///
    /// # Arguments
    /// * `f` - Probability of flipping a bit in permanent response (0.0 to 0.5)
    /// * `p` - Probability of setting a 1 bit to 1 in instantaneous response
    /// * `q` - Probability of setting a 0 bit to 1 in instantaneous response
    ///
    /// `std` only (system entropy); `no_std` は [`Self::with_seed_params`]
    #[cfg(feature = "std")]
    #[must_use]
    pub fn new(f: f64, p: f64, q: f64) -> Self {
        Self::with_seed_params(f, p, q, XorShift64::from_entropy().next_u64())
    }

    /// Create with explicit parameters and RNG seed (`no_std` 用、entropy は caller が供給)
    #[must_use]
    pub fn with_seed_params(f: f64, p: f64, q: f64, seed: u64) -> Self {
        Self {
            f: f.clamp(0.0, 0.5),
            p: p.clamp(0.0, 1.0),
            q: q.clamp(0.0, 1.0),
            rng: XorShift64::new(seed),
        }
    }

    /// Create with typical parameters for ε-differential privacy
    ///
    /// Uses f=0.5, p=0.75, q=0.25 for approximately ε=2 privacy
    #[cfg(feature = "std")]
    #[must_use]
    pub fn default_params() -> Self {
        Self::new(0.5, 0.75, 0.25)
    }

    /// Encode a value into a Bloom filter (simple hash-based)
    #[allow(clippy::unused_self)]
    fn encode_bloom(&self, value: u64) -> [u8; RAPPOR_BITS] {
        use crate::sketch::FnvHasher;

        let mut bloom = [0u8; RAPPOR_BITS];
        // Use multiple hash functions
        for i in 0..3u128 {
            let h = FnvHasher::hash_u128(u128::from(value) | (i << 64));
            let idx = hash_index(h) % RAPPOR_BITS;
            bloom[idx] = 1;
        }
        bloom
    }

    /// Apply permanent randomized response
    fn permanent_response(&mut self, bloom: &[u8; RAPPOR_BITS]) -> [u8; RAPPOR_BITS] {
        let mut result = [0u8; RAPPOR_BITS];
        for i in 0..RAPPOR_BITS {
            if self.rng.next_bool(self.f) {
                // Flip with probability f
                result[i] = u8::from(self.rng.next_bool(0.5));
            } else {
                // Keep original
                result[i] = bloom[i];
            }
        }
        result
    }

    /// Apply instantaneous randomized response
    fn instantaneous_response(&mut self, permanent: &[u8; RAPPOR_BITS]) -> [u8; RAPPOR_BITS] {
        let mut result = [0u8; RAPPOR_BITS];
        for i in 0..RAPPOR_BITS {
            if permanent[i] == 1 {
                result[i] = u8::from(self.rng.next_bool(self.p));
            } else {
                result[i] = u8::from(self.rng.next_bool(self.q));
            }
        }
        result
    }

    /// Encode and privatize a value
    ///
    /// Returns the privatized bit array that can be sent to the aggregator.
    pub fn privatize(&mut self, value: u64) -> [u8; RAPPOR_BITS] {
        let bloom = self.encode_bloom(value);
        let permanent = self.permanent_response(&bloom);
        self.instantaneous_response(&permanent)
    }

    /// Get privacy parameters
    #[must_use]
    pub const fn params(&self) -> (f64, f64, f64) {
        (self.f, self.p, self.q)
    }
}

// ============================================================================
// Privacy Budget Tracker
// ============================================================================

/// Track cumulative privacy budget (composition theorem)
///
/// Under sequential composition, ε values add up.
/// Under parallel composition over disjoint data, use max.
#[derive(Clone, Debug)]
pub struct PrivacyBudget {
    /// Total epsilon spent
    total_epsilon: f64,
    /// Maximum allowed epsilon
    max_epsilon: f64,
    /// Number of queries made
    query_count: u64,
}

impl PrivacyBudget {
    /// Create a new privacy budget tracker
    #[must_use]
    pub const fn new(max_epsilon: f64) -> Self {
        Self {
            total_epsilon: 0.0,
            max_epsilon,
            query_count: 0,
        }
    }

    /// Try to spend epsilon from budget
    ///
    /// Returns `true` and charges the ledger if the request is well formed and
    /// fits in the remaining budget; returns `false` and **leaves every field
    /// untouched** otherwise.
    ///
    /// A request is well formed when `epsilon` is finite and not negative.
    /// `0.0` (and `-0.0`) is well formed: it charges nothing and counts a
    /// query.
    ///
    /// Rejecting a negative `epsilon` is the whole point of the type. Until
    /// 2026-10-07 the guard was only the budget comparison
    /// `total + epsilon <= max`, which a negative value passes trivially, so
    /// `try_spend(-10.0)` returned `true` and *raised* the remaining budget
    /// from 1 to 11; `-inf` raised it to infinity. A caller that spends a
    /// negative epsilon between two honest queries escapes the limit the
    /// budget exists to enforce. `NaN` and `+inf` were already refused,
    /// though only as a side effect of the comparison being false — they are
    /// now refused by the same explicit check.
    pub fn try_spend(&mut self, epsilon: f64) -> bool {
        // `!(epsilon >= 0.0)` would also catch NaN, but spelling both
        // conditions out keeps the intent readable: a non-finite or negative
        // request is not a spend at all, and never reaches the ledger.
        if !epsilon.is_finite() || epsilon < 0.0 {
            return false;
        }
        if self.total_epsilon + epsilon <= self.max_epsilon {
            self.total_epsilon += epsilon;
            self.query_count += 1;
            true
        } else {
            false
        }
    }

    /// Get remaining budget
    #[inline]
    #[must_use]
    pub fn remaining(&self) -> f64 {
        (self.max_epsilon - self.total_epsilon).max(0.0)
    }

    /// Get total spent
    #[inline]
    #[must_use]
    pub const fn spent(&self) -> f64 {
        self.total_epsilon
    }

    /// Get query count
    #[inline]
    #[must_use]
    pub const fn query_count(&self) -> u64 {
        self.query_count
    }

    /// Check if budget is exhausted
    #[inline]
    #[must_use]
    pub fn is_exhausted(&self) -> bool {
        self.total_epsilon >= self.max_epsilon
    }

    /// Reset the budget
    pub const fn reset(&mut self) {
        self.total_epsilon = 0.0;
        self.query_count = 0;
    }
}

// ============================================================================
// Differentially Private Aggregator
// ============================================================================

/// Aggregates noisy reports and estimates true statistics
#[derive(Clone, Debug)]
pub struct PrivateAggregator {
    /// Sum of noisy values
    noisy_sum: f64,
    /// Count of reports
    count: u64,
    /// Laplace noise scale used
    noise_scale: f64,
}

impl PrivateAggregator {
    /// Create a new aggregator
    #[must_use]
    pub const fn new(noise_scale: f64) -> Self {
        Self {
            noisy_sum: 0.0,
            count: 0,
            noise_scale,
        }
    }

    /// Add a noisy report
    #[inline]
    pub fn add(&mut self, noisy_value: f64) {
        self.noisy_sum += noisy_value;
        self.count += 1;
    }

    /// Estimate the true mean
    ///
    /// As count increases, noise averages out to zero.
    #[must_use]
    pub fn estimate_mean(&self) -> f64 {
        if self.count == 0 {
            0.0
        } else {
            let inv_count = 1.0 / u64_f64(self.count);
            self.noisy_sum * inv_count
        }
    }

    /// Estimate the true sum
    #[must_use]
    pub const fn estimate_sum(&self) -> f64 {
        self.noisy_sum
    }

    /// Get the standard error of the mean estimate
    ///
    /// SE = scale * sqrt(2) / sqrt(n)
    #[must_use]
    pub fn standard_error(&self) -> f64 {
        if self.count == 0 {
            f64::INFINITY
        } else {
            let inv_sqrt_n = 1.0 / (u64_f64(self.count)).sqrt();
            self.noise_scale * core::f64::consts::SQRT_2 * inv_sqrt_n
        }
    }

    /// Get the count
    #[inline]
    #[must_use]
    pub const fn count(&self) -> u64 {
        self.count
    }

    /// Reset the aggregator
    pub const fn reset(&mut self) {
        self.noisy_sum = 0.0;
        self.count = 0;
    }
}

// ============================================================================
// Tests
// ============================================================================

#[cfg(test)]
#[allow(clippy::float_cmp)]
mod tests {
    use super::*;

    #[test]
    fn test_xorshift() {
        let mut rng = XorShift64::new(12345);
        let v1 = rng.next_u64();
        let v2 = rng.next_u64();
        assert_ne!(v1, v2);

        // Test reproducibility
        let mut rng2 = XorShift64::new(12345);
        assert_eq!(v1, rng2.next_u64());
    }

    #[test]
    fn test_laplace_noise() {
        let mut noise = LaplaceNoise::with_seed(1.0, 1.0, 42);

        // Generate many samples and check mean is close to 0
        let mut sum = 0.0;
        let n = 10000;
        for _ in 0..n {
            sum += noise.sample();
        }
        let mean = sum / f64::from(n);
        assert!(mean.abs() < 0.1, "mean = {mean}");
    }

    #[test]
    fn test_randomized_response() {
        let mut rr = RandomizedResponse::with_probability(0.75, 42);

        // With high p_true, most responses should match truth
        let mut correct = 0;
        let n = 1000;
        for i in 0..n {
            let truth = i % 2 == 0;
            let response = rr.privatize(truth);
            if response == truth {
                correct += 1;
            }
        }

        // Should be correct more than 50% of the time
        assert!(correct > n / 2, "correct = {correct}");
    }

    #[test]
    fn test_rr_proportion_estimation() {
        // Simulate: true rate = 0.3, n = 10000
        let p_true = 0.8;
        let true_rate = 0.3;
        let n = 10000u64;

        let mut rr = RandomizedResponse::with_probability(p_true, 42);

        let mut positive_reports = 0u64;
        for i in 0..n {
            // Simulate true value with 30% positive rate
            let truth = (u64_f64(i) / u64_f64(n)) < true_rate;
            if rr.privatize(truth) {
                positive_reports += 1;
            }
        }

        let estimated = RandomizedResponse::estimate_proportion(p_true, n, positive_reports);
        // Should be within 0.1 of true rate
        assert!(
            (estimated - true_rate).abs() < 0.1,
            "estimated = {estimated}, true = {true_rate}"
        );
    }

    #[test]
    fn test_privacy_budget() {
        let mut budget = PrivacyBudget::new(1.0);

        assert!(budget.try_spend(0.3));
        assert!(budget.try_spend(0.3));
        assert!(budget.try_spend(0.3));
        assert!(!budget.try_spend(0.3)); // Would exceed

        assert_eq!(budget.query_count(), 3);
        assert!((budget.spent() - 0.9).abs() < 0.001);
    }

    #[test]
    fn test_private_aggregator() {
        let noise_scale = 1.0;
        let mut aggregator = PrivateAggregator::new(noise_scale);
        let mut noise = LaplaceNoise::with_seed(1.0, 1.0, 42);

        // True values are all 100
        let true_value = 100.0;
        let n = 10000;

        for _ in 0..n {
            let privatized = noise.privatize(true_value);
            aggregator.add(privatized);
        }

        let estimated_mean = aggregator.estimate_mean();
        // Should be close to 100
        assert!(
            (estimated_mean - true_value).abs() < 1.0,
            "estimated = {estimated_mean}"
        );
    }

    // `format!` / `Vec` / the entropy-seeded constructors are `std`-only,
    // so this case cannot run on the `no_std` build
    #[cfg(feature = "std")]
    #[test]
    fn test_rappor() {
        let mut rappor = Rappor::default_params();

        let encoded1 = rappor.privatize(12345);
        let encoded2 = rappor.privatize(12345);

        // Different instances should produce different outputs (randomized)
        // but with same structure (64 bits)
        assert_eq!(encoded1.len(), 64);
        assert_eq!(encoded2.len(), 64);
    }

    #[test]
    fn test_xorshift_deterministic() {
        let mut rng1 = XorShift64::new(42);
        let mut rng2 = XorShift64::new(42);
        assert_eq!(rng1.next_u64(), rng2.next_u64());
        assert_eq!(rng1.next_u64(), rng2.next_u64());
    }

    #[test]
    fn test_xorshift_zero_seed() {
        let mut rng = XorShift64::new(0);
        // Should use fallback seed, not produce zeros
        assert_ne!(rng.next_u64(), 0);
    }

    #[test]
    fn test_xorshift_f64_range() {
        let mut rng = XorShift64::new(42);
        for _ in 0..100 {
            let v = rng.next_f64();
            assert!((0.0..1.0).contains(&v), "v = {v}");
        }
    }

    #[test]
    fn test_xorshift_f64_range_custom() {
        let mut rng = XorShift64::new(42);
        for _ in 0..100 {
            let v = rng.next_f64_range(10.0, 20.0);
            assert!((10.0..20.0).contains(&v), "v = {v}");
        }
    }

    #[test]
    fn test_xorshift_bool() {
        let mut rng = XorShift64::new(42);
        let mut trues = 0;
        let n = 10000;
        for _ in 0..n {
            if rng.next_bool(0.5) {
                trues += 1;
            }
        }
        // Should be roughly 50%
        assert!(trues > 4000 && trues < 6000, "trues = {trues}");
    }

    // `format!` / `Vec` / the entropy-seeded constructors are `std`-only,
    // so this case cannot run on the `no_std` build
    #[cfg(feature = "std")]
    #[test]
    fn test_xorshift_default() {
        let mut rng = XorShift64::default();
        let v = rng.next_u64();
        assert_ne!(v, 0);
    }

    #[test]
    fn test_laplace_scale() {
        let noise = LaplaceNoise::with_seed(2.0, 0.5, 42);
        // scale = sensitivity / epsilon = 2.0 / 0.5 = 4.0
        assert!((noise.scale() - 4.0).abs() < f64::EPSILON);
    }

    #[test]
    fn test_laplace_privatize_int() {
        let mut noise = LaplaceNoise::with_seed(1.0, 1.0, 42);
        let result = noise.privatize_int(100);
        // Should be close to 100 but not exactly
        assert!((result - 100).abs() < 20, "result = {result}");
    }

    #[test]
    fn test_rr_privatize_bit() {
        let mut rr = RandomizedResponse::with_probability(0.9, 42);
        let bit = rr.privatize_bit(1);
        assert!(bit == 0 || bit == 1);
    }

    #[test]
    fn test_rr_p_true() {
        let rr = RandomizedResponse::with_probability(0.8, 42);
        assert!((rr.p_true() - 0.8).abs() < f64::EPSILON);
    }

    #[test]
    fn test_rr_estimate_zero_n() {
        let result = RandomizedResponse::estimate_proportion(0.8, 0, 0);
        assert_eq!(result, 0.0);
    }

    #[test]
    fn test_privacy_budget_exhausted() {
        let mut budget = PrivacyBudget::new(1.0);
        budget.try_spend(0.5);
        budget.try_spend(0.5);
        assert!(budget.is_exhausted());
    }

    #[test]
    fn test_privacy_budget_remaining() {
        let mut budget = PrivacyBudget::new(1.0);
        budget.try_spend(0.3);
        assert!((budget.remaining() - 0.7).abs() < 0.001);
    }

    #[test]
    fn test_privacy_budget_reset() {
        let mut budget = PrivacyBudget::new(1.0);
        budget.try_spend(0.5);
        budget.reset();
        assert_eq!(budget.query_count(), 0);
        assert!((budget.remaining() - 1.0).abs() < f64::EPSILON);
    }

    #[test]
    fn test_private_aggregator_empty() {
        let agg = PrivateAggregator::new(1.0);
        assert_eq!(agg.estimate_mean(), 0.0);
        assert_eq!(agg.count(), 0);
        assert!(agg.standard_error().is_infinite());
    }

    #[test]
    fn test_private_aggregator_reset() {
        let mut agg = PrivateAggregator::new(1.0);
        agg.add(10.0);
        agg.add(20.0);
        agg.reset();
        assert_eq!(agg.count(), 0);
        assert_eq!(agg.estimate_mean(), 0.0);
    }

    #[test]
    fn test_private_aggregator_sum() {
        let mut agg = PrivateAggregator::new(1.0);
        agg.add(10.0);
        agg.add(20.0);
        assert!((agg.estimate_sum() - 30.0).abs() < f64::EPSILON);
    }

    // `format!` / `Vec` / the entropy-seeded constructors are `std`-only,
    // so this case cannot run on the `no_std` build
    #[cfg(feature = "std")]
    #[test]
    fn test_rappor_params() {
        let rappor = Rappor::new(0.3, 0.7, 0.2);
        let (f, p, q) = rappor.params();
        assert!((f - 0.3).abs() < f64::EPSILON);
        assert!((p - 0.7).abs() < f64::EPSILON);
        assert!((q - 0.2).abs() < f64::EPSILON);
    }
}