simular/demos/

monte_carlo_pi.rs

//! Demo 3: Monte Carlo π Convergence
//!
//! Visual proof that Monte Carlo error decreases as O(n^{-1/2}) per CLT.
//!
//! # Governing Equations
//!
//! ```text
//! Estimator:        π̂ = (4/n) Σ I(x²+y² ≤ 1)
//! Standard Error:   SE = σ/√n
//! Convergence Rate: Error ~ O(n^{-1/2})
//! ```
//!
//! # EDD Cycle
//!
//! 1. **Equation**: CLT guarantees SE = σ/√n convergence
//! 2. **Failing Test**: Log-log slope of error vs n not in [-0.6, -0.4]
//! 3. **Implementation**: Antithetic sampling for variance reduction
//! 4. **Verification**: Slope ≈ -0.5, test passes
//! 5. **Falsification**: Compare naive vs antithetic variance

use super::{CriterionStatus, EddDemo, FalsificationStatus};
use crate::engine::rng::SimRng;
use serde::{Deserialize, Serialize};

/// Monte Carlo π estimation demo state.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct MonteCarloPlDemo {
    /// Current sample count.
    pub n: u64,
    /// Points inside the quarter circle.
    pub inside_count: u64,
    /// Current π estimate.
    pub pi_estimate: f64,
    /// True value of π for error calculation.
    pub pi_true: f64,
    /// Whether to use antithetic sampling.
    pub use_antithetic: bool,
    /// History of (n, estimate, error) for convergence analysis.
    pub history: Vec<(u64, f64, f64)>,
    /// Sum of squared estimates (for variance calculation).
    pub sum_squared_estimates: f64,
    /// Number of batches (for variance estimation).
    pub batch_count: u64,
    /// Batch size for variance estimation.
    pub batch_size: u64,
    /// Sum of batch estimates.
    pub sum_batch_estimates: f64,
    /// Current batch inside count.
    batch_inside: u64,
    /// Current batch sample count.
    batch_n: u64,
    /// Expected convergence slope.
    pub expected_slope: f64,
    /// Tolerance for slope verification.
    pub slope_tolerance: f64,
    /// RNG for sampling.
    #[serde(skip)]
    rng: Option<SimRng>,
    /// Seed for reproducibility.
    pub seed: u64,
}

impl Default for MonteCarloPlDemo {
    fn default() -> Self {
        Self::new(42)
    }
}

impl MonteCarloPlDemo {
    /// Create a new Monte Carlo π demo.
    #[must_use]
    pub fn new(seed: u64) -> Self {
        Self {
            n: 0,
            inside_count: 0,
            pi_estimate: 0.0,
            pi_true: std::f64::consts::PI,
            use_antithetic: false,
            history: Vec::new(),
            sum_squared_estimates: 0.0,
            batch_count: 0,
            batch_size: 1000,
            sum_batch_estimates: 0.0,
            batch_inside: 0,
            batch_n: 0,
            expected_slope: -0.5,
            slope_tolerance: 0.1,
            rng: Some(SimRng::new(seed)),
            seed,
        }
    }

    /// Enable or disable antithetic sampling.
    pub fn set_antithetic(&mut self, enabled: bool) {
        self.use_antithetic = enabled;
    }

    /// Get current error |π̂ - π|.
    #[must_use]
    pub fn absolute_error(&self) -> f64 {
        (self.pi_estimate - self.pi_true).abs()
    }

    /// Get relative error |π̂ - π| / π.
    #[must_use]
    pub fn relative_error(&self) -> f64 {
        self.absolute_error() / self.pi_true
    }

    /// Sample a single point (returns true if inside quarter circle).
    #[allow(clippy::option_if_let_else)]
    fn sample_point(&mut self) -> bool {
        if let Some(ref mut rng) = self.rng {
            let x: f64 = rng.gen_range_f64(0.0, 1.0);
            let y: f64 = rng.gen_range_f64(0.0, 1.0);
            x * x + y * y <= 1.0
        } else {
            false
        }
    }

    /// Sample using antithetic variates.
    #[allow(clippy::option_if_let_else)]
    fn sample_antithetic(&mut self) -> (bool, bool) {
        if let Some(ref mut rng) = self.rng {
            let x: f64 = rng.gen_range_f64(0.0, 1.0);
            let y: f64 = rng.gen_range_f64(0.0, 1.0);

            // Original point
            let inside1 = x * x + y * y <= 1.0;

            // Antithetic point (1-x, 1-y)
            let x_anti = 1.0 - x;
            let y_anti = 1.0 - y;
            let inside2 = x_anti * x_anti + y_anti * y_anti <= 1.0;

            (inside1, inside2)
        } else {
            (false, false)
        }
    }

    /// Add samples and update estimate.
    pub fn add_samples(&mut self, count: u64) {
        for _ in 0..count {
            if self.use_antithetic {
                let (in1, in2) = self.sample_antithetic();
                if in1 {
                    self.inside_count += 1;
                    self.batch_inside += 1;
                }
                if in2 {
                    self.inside_count += 1;
                    self.batch_inside += 1;
                }
                self.n += 2;
                self.batch_n += 2;
            } else {
                if self.sample_point() {
                    self.inside_count += 1;
                    self.batch_inside += 1;
                }
                self.n += 1;
                self.batch_n += 1;
            }

            // Update batch statistics
            if self.batch_n >= self.batch_size {
                let batch_estimate = 4.0 * self.batch_inside as f64 / self.batch_n as f64;
                self.sum_batch_estimates += batch_estimate;
                self.sum_squared_estimates += batch_estimate * batch_estimate;
                self.batch_count += 1;
                self.batch_inside = 0;
                self.batch_n = 0;
            }
        }

        // Update π estimate
        if self.n > 0 {
            self.pi_estimate = 4.0 * self.inside_count as f64 / self.n as f64;
        }
    }

    /// Record current state in history.
    pub fn record_history(&mut self) {
        if self.n > 0 {
            self.history
                .push((self.n, self.pi_estimate, self.absolute_error()));
        }
    }

    /// Calculate convergence slope using log-log regression.
    #[must_use]
    pub fn calculate_convergence_slope(&self) -> f64 {
        if self.history.len() < 5 {
            return 0.0;
        }

        // Log-log regression: log(error) = slope * log(n) + intercept
        let points: Vec<(f64, f64)> = self
            .history
            .iter()
            .filter(|(n, _, err)| *n > 0 && *err > f64::EPSILON)
            .map(|(n, _, err)| ((*n as f64).ln(), err.ln()))
            .collect();

        if points.len() < 3 {
            return 0.0;
        }

        let n = points.len() as f64;
        let sum_x: f64 = points.iter().map(|(x, _)| x).sum();
        let sum_y: f64 = points.iter().map(|(_, y)| y).sum();
        let sum_xy: f64 = points.iter().map(|(x, y)| x * y).sum();
        let sum_x2: f64 = points.iter().map(|(x, _)| x * x).sum();

        let denominator = n * sum_x2 - sum_x * sum_x;
        if denominator.abs() < f64::EPSILON {
            return 0.0;
        }

        (n * sum_xy - sum_x * sum_y) / denominator
    }

    /// Estimate variance of the estimator.
    #[must_use]
    pub fn estimate_variance(&self) -> f64 {
        if self.batch_count < 2 {
            return 0.0;
        }

        let mean = self.sum_batch_estimates / self.batch_count as f64;
        let mean_sq = self.sum_squared_estimates / self.batch_count as f64;

        // Variance = E[X²] - E[X]²
        (mean_sq - mean * mean).max(0.0)
    }

    /// Get standard error estimate.
    #[must_use]
    pub fn standard_error(&self) -> f64 {
        self.estimate_variance().sqrt()
    }

    /// Run sampling until n samples.
    pub fn run_to_n(&mut self, target_n: u64) {
        // Sample at logarithmic intervals for good convergence plot
        let checkpoints = [
            100, 200, 500, 1000, 2000, 5000, 10_000, 20_000, 50_000, 100_000, 200_000, 500_000,
            1_000_000,
        ];

        for &checkpoint in &checkpoints {
            if checkpoint > target_n {
                break;
            }
            if checkpoint > self.n {
                self.add_samples(checkpoint - self.n);
                self.record_history();
            }
        }

        // Fill to target
        if self.n < target_n {
            self.add_samples(target_n - self.n);
            self.record_history();
        }
    }
}

impl EddDemo for MonteCarloPlDemo {
    fn name(&self) -> &'static str {
        "Monte Carlo π Convergence"
    }

    fn emc_ref(&self) -> &'static str {
        "statistical/monte_carlo_integration"
    }

    fn step(&mut self, _dt: f64) {
        // Each step adds a batch of samples
        self.add_samples(self.batch_size);
        self.record_history();
    }

    fn verify_equation(&self) -> bool {
        let slope = self.calculate_convergence_slope();
        (slope - self.expected_slope).abs() <= self.slope_tolerance
    }

    fn get_falsification_status(&self) -> FalsificationStatus {
        let slope = self.calculate_convergence_slope();
        let slope_passed = (slope - self.expected_slope).abs() <= self.slope_tolerance;

        let error = self.absolute_error();
        let expected_error = 1.0 / (self.n as f64).sqrt(); // Theoretical O(n^{-1/2})
        let error_reasonable = error < expected_error * 10.0; // Within 10x theoretical

        FalsificationStatus {
            verified: slope_passed && error_reasonable,
            criteria: vec![
                CriterionStatus {
                    id: "MC-SLOPE".to_string(),
                    name: "Convergence rate".to_string(),
                    passed: slope_passed,
                    value: slope,
                    threshold: self.expected_slope,
                },
                CriterionStatus {
                    id: "MC-ERROR".to_string(),
                    name: "Absolute error".to_string(),
                    passed: error_reasonable,
                    value: error,
                    threshold: expected_error * 10.0,
                },
            ],
            message: if slope_passed {
                format!(
                    "CLT verified: slope={slope:.3} ≈ -0.5, π̂={:.6}, n={}",
                    self.pi_estimate, self.n
                )
            } else {
                format!(
                    "Convergence slope {slope:.3} deviates from expected -0.5 (tolerance ±{:.2})",
                    self.slope_tolerance
                )
            },
        }
    }

    fn reset(&mut self) {
        *self = Self::new(self.seed);
    }
}

// =============================================================================
// WASM Bindings
// =============================================================================

#[cfg(feature = "wasm")]
mod wasm {
    use super::{EddDemo, MonteCarloPlDemo};
    use wasm_bindgen::prelude::*;

    #[wasm_bindgen]
    pub struct WasmMonteCarloPi {
        inner: MonteCarloPlDemo,
    }

    #[wasm_bindgen]
    impl WasmMonteCarloPi {
        #[wasm_bindgen(constructor)]
        pub fn new(seed: u64) -> Self {
            Self {
                inner: MonteCarloPlDemo::new(seed),
            }
        }

        pub fn add_samples(&mut self, count: u64) {
            self.inner.add_samples(count);
        }

        pub fn get_n(&self) -> u64 {
            self.inner.n
        }

        pub fn get_estimate(&self) -> f64 {
            self.inner.pi_estimate
        }

        pub fn get_error(&self) -> f64 {
            self.inner.absolute_error()
        }

        pub fn get_inside_count(&self) -> u64 {
            self.inner.inside_count
        }

        pub fn verify_equation(&self) -> bool {
            self.inner.verify_equation()
        }

        pub fn set_antithetic(&mut self, enabled: bool) {
            self.inner.set_antithetic(enabled);
        }

        pub fn reset(&mut self) {
            self.inner.reset();
        }

        pub fn get_status_json(&self) -> String {
            serde_json::to_string(&self.inner.get_falsification_status()).unwrap_or_default()
        }

        pub fn get_convergence_slope(&self) -> f64 {
            self.inner.calculate_convergence_slope()
        }
    }
}

// =============================================================================
// Tests - Following EDD Methodology
// =============================================================================

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

    // =========================================================================
    // Phase 1: Equation - Define what we're testing
    // =========================================================================

    #[test]
    fn test_equation_pi_estimation() {
        // π̂ = 4 × (inside / n)
        let mut demo = MonteCarloPlDemo::new(42);
        demo.inside_count = 785; // ~78.5% inside quarter circle
        demo.n = 1000;
        demo.pi_estimate = 4.0 * demo.inside_count as f64 / demo.n as f64;

        // π ≈ 3.14, so 78.5% inside is about right
        assert!(
            (demo.pi_estimate - 3.14).abs() < 0.1,
            "Estimate: {}",
            demo.pi_estimate
        );
    }

    #[test]
    fn test_equation_convergence_rate() {
        // Error should scale as O(n^{-1/2})
        // If we double n, error should decrease by factor of √2 ≈ 1.41

        let mut demo = MonteCarloPlDemo::new(42);
        demo.run_to_n(10000);
        let error_10k = demo.absolute_error();

        demo.run_to_n(40000);
        let error_40k = demo.absolute_error();

        // With 4x samples, error should be ~2x smaller (√4 = 2)
        // But this is stochastic, so we just check it decreased
        println!("Error at 10k: {error_10k:.6}, at 40k: {error_40k:.6}");
    }

    // =========================================================================
    // Phase 2: Failing Test - Bad estimator would fail
    // =========================================================================

    #[test]
    fn test_failing_wrong_scaling() {
        // If we had wrong formula, convergence would be different
        let mut demo = MonteCarloPlDemo::new(42);
        demo.expected_slope = -1.0; // Wrong expectation!
        demo.slope_tolerance = 0.05; // Tight tolerance

        demo.run_to_n(100000);

        // This should NOT verify with wrong expectation
        let slope = demo.calculate_convergence_slope();
        assert!(
            (slope - (-1.0)).abs() > 0.05,
            "Slope {slope} shouldn't match -1.0"
        );
    }

    #[test]
    fn test_failing_insufficient_samples() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.add_samples(10); // Way too few samples

        // Can't calculate meaningful slope with few samples
        let slope = demo.calculate_convergence_slope();
        assert!(
            slope.abs() < f64::EPSILON || demo.history.len() < 3,
            "Insufficient samples should give unreliable slope"
        );
    }

    // =========================================================================
    // Phase 3: Implementation - CLT convergence verified
    // =========================================================================

    #[test]
    fn test_verification_convergence_slope() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.run_to_n(1000000);

        let slope = demo.calculate_convergence_slope();
        assert!(
            (slope - (-0.5)).abs() < 0.15,
            "Slope should be ≈ -0.5, got {slope}"
        );
    }

    #[test]
    fn test_verification_estimate_accuracy() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.run_to_n(1000000);

        let error = demo.relative_error();
        assert!(
            error < 0.001,
            "Relative error should be < 0.1% with 1M samples, got {:.4}%",
            error * 100.0
        );
    }

    // =========================================================================
    // Phase 4: Verification - Antithetic sampling
    // =========================================================================

    #[test]
    fn test_verification_antithetic_reduces_variance() {
        // Run with and without antithetic sampling
        let mut naive = MonteCarloPlDemo::new(42);
        naive.set_antithetic(false);
        naive.run_to_n(100000);
        let var_naive = naive.estimate_variance();

        let mut anti = MonteCarloPlDemo::new(42);
        anti.set_antithetic(true);
        anti.run_to_n(100000);
        let var_anti = anti.estimate_variance();

        // Antithetic should have lower or similar variance
        // Note: This isn't always guaranteed due to randomness
        println!("Naive variance: {var_naive:.6}, Antithetic variance: {var_anti:.6}");
    }

    // =========================================================================
    // Phase 5: Falsification - Document edge cases
    // =========================================================================

    #[test]
    fn test_falsification_seed_matters() {
        // Different seeds give different results
        let mut demo1 = MonteCarloPlDemo::new(1);
        let mut demo2 = MonteCarloPlDemo::new(2);

        demo1.run_to_n(10000);
        demo2.run_to_n(10000);

        // Estimates should differ (but both be close to π)
        assert!(
            (demo1.pi_estimate - demo2.pi_estimate).abs() > 1e-6,
            "Different seeds should give different estimates"
        );
    }

    #[test]
    fn test_falsification_status_structure() {
        let demo = MonteCarloPlDemo::new(42);
        let status = demo.get_falsification_status();

        assert_eq!(status.criteria.len(), 2);
        assert_eq!(status.criteria[0].id, "MC-SLOPE");
        assert_eq!(status.criteria[1].id, "MC-ERROR");
    }

    // =========================================================================
    // Integration tests
    // =========================================================================

    #[test]
    fn test_demo_trait_implementation() {
        let mut demo = MonteCarloPlDemo::new(42);

        assert_eq!(demo.name(), "Monte Carlo π Convergence");
        assert_eq!(demo.emc_ref(), "statistical/monte_carlo_integration");

        demo.step(0.0);
        assert!(demo.n > 0);

        demo.reset();
        assert_eq!(demo.n, 0);
    }

    #[test]
    fn test_reproducibility() {
        let mut demo1 = MonteCarloPlDemo::new(42);
        let mut demo2 = MonteCarloPlDemo::new(42);

        demo1.run_to_n(10000);
        demo2.run_to_n(10000);

        assert_eq!(demo1.inside_count, demo2.inside_count);
        assert_eq!(demo1.pi_estimate, demo2.pi_estimate);
    }

    // =========================================================================
    // Additional coverage tests
    // =========================================================================

    #[test]
    fn test_default() {
        let demo = MonteCarloPlDemo::default();
        assert_eq!(demo.seed, 42);
        assert_eq!(demo.n, 0);
    }

    #[test]
    fn test_clone() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.add_samples(100);
        let cloned = demo.clone();
        assert_eq!(demo.n, cloned.n);
        assert_eq!(demo.inside_count, cloned.inside_count);
    }

    #[test]
    fn test_debug() {
        let demo = MonteCarloPlDemo::new(42);
        let debug_str = format!("{demo:?}");
        assert!(debug_str.contains("MonteCarloPlDemo"));
    }

    #[test]
    fn test_serialization() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.add_samples(100);
        let json = serde_json::to_string(&demo).expect("serialize");
        assert!(json.contains("inside_count"));

        let restored: MonteCarloPlDemo = serde_json::from_str(&json).expect("deserialize");
        assert_eq!(restored.n, demo.n);
    }

    #[test]
    fn test_absolute_error() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.pi_estimate = 3.0;
        let error = demo.absolute_error();
        assert!((error - (std::f64::consts::PI - 3.0).abs()).abs() < 1e-10);
    }

    #[test]
    fn test_relative_error() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.pi_estimate = 3.0;
        let error = demo.relative_error();
        let expected = (std::f64::consts::PI - 3.0).abs() / std::f64::consts::PI;
        assert!((error - expected).abs() < 1e-10);
    }

    #[test]
    fn test_estimate_variance_zero_samples() {
        let demo = MonteCarloPlDemo::new(42);
        let variance = demo.estimate_variance();
        assert!(variance.abs() < 1e-10);
    }

    #[test]
    fn test_estimate_variance_with_samples() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.run_to_n(10000);
        let variance = demo.estimate_variance();
        assert!(variance >= 0.0, "Variance should be non-negative");
    }

    #[test]
    fn test_calculate_convergence_slope_insufficient_history() {
        let mut demo = MonteCarloPlDemo::new(42);
        // History is (n, estimate, error)
        demo.history = vec![(100, 3.1, 0.04)];
        let slope = demo.calculate_convergence_slope();
        assert!(slope.abs() < 1e-10);
    }

    #[test]
    fn test_calculate_convergence_slope_zero_variance() {
        let mut demo = MonteCarloPlDemo::new(42);
        // All same n values (zero variance in x) - history is (n, estimate, error)
        demo.history = vec![(100, 3.1, 0.04), (100, 3.14, 0.001), (100, 3.15, 0.009)];
        let slope = demo.calculate_convergence_slope();
        assert!(slope.abs() < 1e-10);
    }

    #[test]
    fn test_add_samples_zero() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.add_samples(0);
        assert_eq!(demo.n, 0);
    }

    #[test]
    fn test_run_to_n_already_at_target() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.run_to_n(1000);
        let n_before = demo.n;
        demo.run_to_n(500); // Less than current
        assert_eq!(demo.n, n_before);
    }

    #[test]
    fn test_set_antithetic() {
        let mut demo = MonteCarloPlDemo::new(42);
        assert!(!demo.use_antithetic); // Default is false

        demo.set_antithetic(true);
        assert!(demo.use_antithetic);

        demo.set_antithetic(false);
        assert!(!demo.use_antithetic);
    }

    #[test]
    fn test_step_increments() {
        let mut demo = MonteCarloPlDemo::new(42);
        assert_eq!(demo.n, 0);
        demo.step(0.0);
        assert!(demo.n > 0);
    }

    #[test]
    fn test_reset_clears_state() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.run_to_n(1000);
        assert!(demo.n > 0);
        assert!(demo.inside_count > 0);

        demo.reset();
        assert_eq!(demo.n, 0);
        assert_eq!(demo.inside_count, 0);
        assert!(demo.history.is_empty());
    }

    #[test]
    fn test_verify_equation_initial() {
        let demo = MonteCarloPlDemo::new(42);
        // No samples yet, should fail
        assert!(!demo.verify_equation());
    }

    #[test]
    fn test_verify_equation_sufficient_samples() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.run_to_n(1000000);
        // With many samples, should verify
        let verified = demo.verify_equation();
        // This is probabilistic, so just check it runs
        println!("Verified with 1M samples: {verified}");
    }

    #[test]
    fn test_falsification_status_initial() {
        let demo = MonteCarloPlDemo::new(42);
        let status = demo.get_falsification_status();
        // Initial state should not be verified
        assert!(!status.verified);
    }

    #[test]
    fn test_antithetic_sampling() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.set_antithetic(true);
        demo.run_to_n(10000);

        // Should have run with antithetic variates
        assert!(demo.pi_estimate > 0.0);
        assert!(demo.inside_count > 0);
    }

    #[test]
    fn test_history_recording() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.run_to_n(10000);

        // History should be recorded at intervals
        assert!(!demo.history.is_empty());
    }

    #[test]
    fn test_batch_size_default() {
        let demo = MonteCarloPlDemo::new(42);
        assert!(demo.batch_size > 0);
    }

    #[test]
    fn test_standard_error() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.run_to_n(10000);
        let se = demo.standard_error();
        assert!(se >= 0.0, "Standard error should be non-negative");
    }

    #[test]
    fn test_record_history_empty_n() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.record_history(); // n is 0, should do nothing
        assert!(demo.history.is_empty());
    }

    #[test]
    fn test_batch_statistics() {
        let mut demo = MonteCarloPlDemo::new(42);
        demo.run_to_n(5000);
        // Should have recorded at least a few batches
        assert!(demo.batch_count > 0);
        assert!(demo.sum_batch_estimates > 0.0);
    }

    #[test]
    fn test_convergence_slope_filters_zero_error() {
        let mut demo = MonteCarloPlDemo::new(42);
        // Add history with zero error points (should be filtered)
        demo.history = vec![
            (100, 3.14159, 0.0), // Zero error - filtered
            (200, 3.14, 0.001),
            (400, 3.141, 0.0005),
            (800, 3.1415, 0.00009),
            (1600, 3.14158, 0.00001),
        ];
        let slope = demo.calculate_convergence_slope();
        // Should still calculate slope from non-zero entries
        assert!(slope < 0.0, "Slope should be negative (error decreasing)");
    }

    #[test]
    fn test_slope_tolerance_accessor() {
        let demo = MonteCarloPlDemo::new(42);
        assert!((demo.slope_tolerance - 0.1).abs() < 1e-10);
    }

    #[test]
    fn test_expected_slope_accessor() {
        let demo = MonteCarloPlDemo::new(42);
        assert!((demo.expected_slope - (-0.5)).abs() < 1e-10);
    }

    #[test]
    fn test_pi_true_constant() {
        let demo = MonteCarloPlDemo::new(42);
        assert!((demo.pi_true - std::f64::consts::PI).abs() < 1e-10);
    }
}