//! Interactive Proof Explorer for Special Functions
//!
//! This module provides a comprehensive system for exploring mathematical proofs

#![allow(clippy::all)]
//! interactively, with step-by-step guidance, visualization, and verification.
//!
//! Features:
//! - Step-by-step proof navigation with detailed explanations
//! - Interactive mathematical derivations with user input validation
//! - Visual aids and geometric interpretations
//! - Cross-references to related theorems and applications
//! - Computational verification of theoretical results
//! - Adaptive difficulty based on user understanding
//!
//! Run with: cargo run --example interactive_proof_explorer

use scirs2_core::ndarray::{Array1, Array2};
use scirs2_core::Complex64;
use scirs2_special::*;
use std::collections::HashMap;
use std::f64::consts::PI;
use std::io::{self, Write};
use std::time::{Duration, Instant};

#[allow(dead_code)]
#[derive(Debug, Clone)]
struct ProofExplorer {
    available_proofs: Vec<ProofModule>,
    user_progress: UserProgress,
    current_proof: Option<ProofSession>,
    difficulty_settings: DifficultySettings,
}

#[allow(dead_code)]
#[derive(Debug, Clone)]
struct ProofModule {
    id: String,
    title: String,
    description: String,
    mathematical_level: MathematicalLevel,
    prerequisites: Vec<String>,
    main_theorem: Theorem,
    proof_steps: Vec<ProofStep>,
    applications: Vec<Application>,
    related_proofs: Vec<String>,
    computational_examples: Vec<ComputationalExample>,
}

#[allow(dead_code)]
#[derive(Debug, Clone)]
enum MathematicalLevel {
    Undergraduate,
    Graduate,
    Advanced,
    Research,
}

#[allow(dead_code)]
#[derive(Debug, Clone)]
struct Theorem {
    statement: String,
    mathematical_formulation: String,
    conditions: Vec<String>,
    significance: String,
    historical_context: String,
}

#[allow(dead_code)]
#[derive(Debug, Clone)]
struct ProofStep {
    id: String,
    title: String,
    explanation: String,
    mathematical_content: String,
    justification: String,
    interactive_elements: Vec<InteractiveElement>,
    verification_method: Option<VerificationMethod>,
    common_confusions: Vec<String>,
    hints: Vec<String>,
}

#[allow(dead_code)]
#[derive(Debug, Clone)]
enum InteractiveElement {
    UserInput {
        prompt: String,
        expected_answer: String,
        validation_tolerance: f64,
        hints: Vec<String>,
    },
    Visualization {
        title: String,
        description: String,
        plot_type: PlotType,
    },
    Computation {
        description: String,
        code_example: String,
        expected_result: String,
    },
    ConceptualQuestion {
        question: String,
        options: Vec<String>,
        correct_answer: usize,
        explanation: String,
    },
}

#[allow(dead_code)]
#[derive(Debug, Clone)]
enum PlotType {
    FunctionGraph {
        domain: (f64, f64),
        functions: Vec<String>,
    },
    ComplexPlane {
        radius: f64,
        function: String,
    },
    Convergence {
        series: String,
        terms: usize,
    },
    AsymptoticComparison {
        exact: String,
        approximation: String,
        domain: (f64, f64),
    },
}

#[allow(dead_code)]
#[derive(Debug, Clone)]
enum VerificationMethod {
    NumericalCheck {
        description: String,
        test_cases: Vec<TestCase>,
    },
    SymbolicVerification {
        description: String,
        identity: String,
    },
    AsymptoticVerification {
        description: String,
        limit_expression: String,
    },
}

#[allow(dead_code)]
#[derive(Debug, Clone)]
struct TestCase {
    input: f64,
    expected: f64,
    tolerance: f64,
    description: String,
}

#[allow(dead_code)]
#[derive(Debug, Clone)]
struct Application {
    title: String,
    description: String,
    mathematical_context: String,
    practical_example: String,
}

#[allow(dead_code)]
#[derive(Debug, Clone)]
struct ComputationalExample {
    title: String,
    description: String,
    code: String,
    expected_output: String,
    learning_objective: String,
}

#[allow(dead_code)]
#[derive(Debug, Clone)]
struct UserProgress {
    completed_proofs: Vec<String>,
    difficulty_level: MathematicalLevel,
    understanding_scores: HashMap<String, f64>,
    time_spent: HashMap<String, Duration>,
    favorite_topics: Vec<String>,
    achievements: Vec<String>,
}

#[allow(dead_code)]
#[derive(Debug, Clone)]
struct ProofSession {
    proof_id: String,
    current_step: usize,
    start_time: Instant,
    step_completion_times: Vec<Duration>,
    user_responses: Vec<UserResponse>,
    understanding_score: f64,
}

#[allow(dead_code)]
#[derive(Debug, Clone)]
struct UserResponse {
    step_id: String,
    response: String,
    correct: bool,
    attempts: u32,
    time_taken: Duration,
}

#[allow(dead_code)]
#[derive(Debug, Clone)]
struct DifficultySettings {
    show_hints_automatically: bool,
    require_step_verification: bool,
    include_computational_checks: bool,
    mathematical_rigor_level: MathematicalLevel,
}

impl ProofExplorer {
    fn new() -> Self {
        let mut explorer = Self {
            available_proofs: Vec::new(),
            user_progress: UserProgress {
                completed_proofs: Vec::new(),
                difficulty_level: MathematicalLevel::Undergraduate,
                understanding_scores: HashMap::new(),
                time_spent: HashMap::new(),
                favorite_topics: Vec::new(),
                achievements: Vec::new(),
            },
            current_proof: None,
            difficulty_settings: DifficultySettings {
                show_hints_automatically: true,
                require_step_verification: true,
                include_computational_checks: true,
                mathematical_rigor_level: MathematicalLevel::Undergraduate,
            },
        };

        explorer.initialize_proof_library();
        explorer
    }

    fn initialize_proof_library(&mut self) {
        // Gamma Function Reflection Formula
        self.available_proofs.push(ProofModule {
            id: "gamma_reflection".to_string(),
            title: "Gamma Function Reflection Formula".to_string(),
            description: "A beautiful proof connecting the gamma function to trigonometry through complex analysis".to_string(),
            mathematical_level: MathematicalLevel::Graduate,
            prerequisites: vec![
                "complex_analysis".to_string(),
                "gamma_function_basics".to_string(),
                "residue_calculus".to_string(),
            ],
            main_theorem: Theorem {
                statement: "For all z not a non-positive integer, Γ(z)Γ(1-z) = π/sin(πz)".to_string(),
                mathematical_formulation: "\\Gamma(z)\\Gamma(1-z) = \\frac{\\pi}{\\sin(\\pi z)}".to_string(),
                conditions: vec![
                    "z ∉ {..., -2, -1, 0}".to_string(),
                    "The formula extends by analytic continuation".to_string(),
                ],
                significance: "This formula connects discrete factorials to continuous trigonometric functions, demonstrating the deep unity of mathematics".to_string(),
                historical_context: "Discovered by Euler as part of his investigation into extending factorials to non-integer values".to_string(),
            },
            proof_steps: create_reflection_formula_proof_steps(),
            applications: vec![
                Application {
                    title: "Beta Function Simplification".to_string(),
                    description: "Simplifies beta function evaluations for special arguments".to_string(),
                    mathematical_context: "B(z, 1-z) = Γ(z)Γ(1-z) = π/sin(πz)".to_string(),
                    practical_example: "Computing B(1/3, 2/3) = 2π/√3".to_string(),
                },
                Application {
                    title: "Special Values of Gamma".to_string(),
                    description: "Provides a method to compute Γ(1/2), Γ(1/3), etc.".to_string(),
                    mathematical_context: "For z = 1/2: Γ(1/2)² = π".to_string(),
                    practical_example: "Γ(1/2) = √π ≈ 1.77245385".to_string(),
                },
            ],
            related_proofs: vec![
                "gamma_duplication".to_string(),
                "beta_function_properties".to_string(),
            ],
            computational_examples: create_reflection_formula_examples(),
        });

        // Stirling's Approximation with Complete Asymptotic Series
        self.available_proofs.push(ProofModule {
            id: "stirling_asymptotic".to_string(),
            title: "Stirling's Asymptotic Series".to_string(),
            description: "Derivation of the complete asymptotic expansion of the gamma function using the Euler-Maclaurin formula".to_string(),
            mathematical_level: MathematicalLevel::Advanced,
            prerequisites: vec![
                "euler_maclaurin_formula".to_string(),
                "asymptotic_series".to_string(),
                "bernoulli_numbers".to_string(),
            ],
            main_theorem: Theorem {
                statement: "For large |z| with |arg z| < π: ln Γ(z) ~ (z-1/2)ln z - z + (1/2)ln(2π) + asymptotic series".to_string(),
                mathematical_formulation: "\\ln \\Gamma(z) \\sim \\left(z - \\frac{1}{2}\\right) \\ln z - z + \\frac{1}{2}\\ln(2\\pi) + \\sum_{k=1}^{\\infty} \\frac{B_{2k}}{2k(2k-1)z^{2k-1}}".to_string(),
                conditions: vec![
                    "|z| → ∞".to_string(),
                    "|arg z| < π".to_string(),
                    "Series is asymptotic, not convergent".to_string(),
                ],
                significance: "Provides highly accurate approximations for large arguments and connects gamma function to Bernoulli numbers".to_string(),
                historical_context: "Stirling's original approximation was refined by many mathematicians to give the complete asymptotic series".to_string(),
            },
            proof_steps: create_stirling_proof_steps(),
            applications: vec![
                Application {
                    title: "High-Precision Gamma Computation".to_string(),
                    description: "Enables accurate computation of gamma function for large arguments".to_string(),
                    mathematical_context: "Each term in the series provides additional precision".to_string(),
                    practical_example: "Computing Γ(100) with 15-digit accuracy".to_string(),
                },
            ],
            related_proofs: vec!["euler_maclaurin_derivation".to_string()],
            computational_examples: create_stirling_examples(),
        });

        // Add more proofs...
        self.add_additional_proofs();
    }

    fn add_additional_proofs(&mut self) {
        // Bessel Function Generating Function
        self.available_proofs.push(ProofModule {
            id: "bessel_generating".to_string(),
            title: "Bessel Function Generating Function".to_string(),
            description: "Proof that exp(x(t-1/t)/2) generates the Bessel functions".to_string(),
            mathematical_level: MathematicalLevel::Graduate,
            prerequisites: vec![
                "power_series".to_string(),
                "bessel_functions".to_string(),
                "complex_analysis".to_string(),
            ],
            main_theorem: Theorem {
                statement: "The exponential exp(x(t-1/t)/2) is the generating function for Bessel functions of the first kind".to_string(),
                mathematical_formulation: "e^{\\frac{x}{2}(t-\\frac{1}{t})} = \\sum_{n=-\\infty}^{\\infty} J_n(x) t^n".to_string(),
                conditions: vec![
                    "x ∈ ℂ".to_string(),
                    "t ≠ 0".to_string(),
                    "Series converges for all finite x and t ≠ 0".to_string(),
                ],
                significance: "Provides a unified approach to deriving properties of all Bessel functions simultaneously".to_string(),
                historical_context: "This generating function was crucial in Bessel's original work on planetary motion".to_string(),
            },
            proof_steps: create_bessel_generating_proof_steps(),
            applications: vec![
                Application {
                    title: "Bessel Function Identities".to_string(),
                    description: "Generates recurrence relations and addition formulas".to_string(),
                    mathematical_context: "Differentiating the generating function gives identities".to_string(),
                    practical_example: "Deriving J_{n-1}(x) + J_{n+1}(x) = (2n/x)J_n(x)".to_string(),
                },
            ],
            related_proofs: vec!["bessel_recurrence_relations".to_string()],
            computational_examples: create_bessel_generating_examples(),
        });

        // Error Function Series Expansion
        self.available_proofs.push(ProofModule {
            id: "error_function_series".to_string(),
            title: "Error Function Series Expansion".to_string(),
            description: "Derivation of the power series for the error function through term-by-term integration".to_string(),
            mathematical_level: MathematicalLevel::Undergraduate,
            prerequisites: vec![
                "power_series".to_string(),
                "integration".to_string(),
                "gaussian_integral".to_string(),
            ],
            main_theorem: Theorem {
                statement: "The error function has the series expansion erf(x) = (2/√π) Σ (-1)^n x^(2n+1) / (n!(2n+1))".to_string(),
                mathematical_formulation: "\\text{erf}(x) = \\frac{2}{\\sqrt{\\pi}} \\sum_{n=0}^{\\infty} \\frac{(-1)^n x^{2n+1}}{n!(2n+1)}".to_string(),
                conditions: vec![
                    "x ∈ ℂ".to_string(),
                    "Series converges for all finite x".to_string(),
                    "Radius of convergence is infinite".to_string(),
                ],
                significance: "Enables accurate computation of error function values and provides insight into its analytical structure".to_string(),
                historical_context: "This series was developed as part of the study of the normal distribution in probability theory".to_string(),
            },
            proof_steps: create_error_function_proof_steps(),
            applications: vec![
                Application {
                    title: "Probability Calculations".to_string(),
                    description: "Computing probabilities for normal distributions".to_string(),
                    mathematical_context: "P(X ≤ x) = (1/2)[1 + erf(x/√2)]".to_string(),
                    practical_example: "Computing P(|X| ≤ 1) for standard normal X".to_string(),
                },
            ],
            related_proofs: vec!["gaussian_integral_derivation".to_string()],
            computational_examples: create_error_function_examples(),
        });
    }

    fn run(&mut self) -> Result<(), Box<dyn std::error::Error>> {
        println!("🧮 Interactive Proof Explorer for Special Functions");
        println!("=====================================================\n");

        println!("Welcome to the advanced mathematical proof exploration system!");
        println!("Here you can learn and verify complex mathematical proofs step by step.\n");

        self.display_main_menu()?;

        loop {
            println!("\n🔍 Main Menu:");
            println!("1. 📚 Browse available proofs");
            println!("2. 🎯 Continue a proof in progress");
            println!("3. 📊 View your progress and achievements");
            println!("4. ⚙️ Adjust difficulty settings");
            println!("5. 🔍 Search proofs by topic");
            println!("6. 💡 Get proof recommendations");
            println!("7. 📖 Mathematical reference library");
            println!("8. 💾 Save progress and exit");

            let choice = get_user_input("Enter your choice (1-8): ")?;

            match choice.parse::<u32>() {
                Ok(1) => self.browse_proofs()?,
                Ok(2) => self.continue_proof()?,
                Ok(3) => self.display_progress()?,
                Ok(4) => self.adjust_difficulty()?,
                Ok(5) => self.search_proofs()?,
                Ok(6) => self.get_recommendations()?,
                Ok(7) => self.reference_library()?,
                Ok(8) => {
                    self.save_progress()?;
                    println!("👋 Thank you for exploring mathematical proofs!");
                    break;
                }
                _ => println!("❌ Invalid choice. Please try again."),
            }
        }

        Ok(())
    }

    fn display_main_menu(&self) -> Result<(), Box<dyn std::error::Error>> {
        println!("📈 Your Current Status:");
        println!(
            "• Mathematical Level: {:?}",
            self.user_progress.difficulty_level
        );
        println!(
            "• Completed Proofs: {}",
            self.user_progress.completed_proofs.len()
        );
        println!("• Available Proofs: {}", self.available_proofs.len());

        if let Some(current) = &self.current_proof {
            let proof = self
                .get_proof_by_id(&current.proof_id)
                .expect("Operation failed");
            println!(
                "• Current Proof: {} (Step {}/{})",
                proof.title,
                current.current_step + 1,
                proof.proof_steps.len()
            );
        }

        Ok(())
    }

    fn browse_proofs(&self) -> Result<(), Box<dyn std::error::Error>> {
        println!("\n📚 Available Mathematical Proofs");
        println!("==================================\n");

        for (i, proof) in self.available_proofs.iter().enumerate() {
            let completed = if self.user_progress.completed_proofs.contains(&proof.id) {
                "✅"
            } else {
                "🔒"
            };

            let level_indicator = match proof.mathematical_level {
                MathematicalLevel::Undergraduate => "🟢",
                MathematicalLevel::Graduate => "🟡",
                MathematicalLevel::Advanced => "🟠",
                MathematicalLevel::Research => "🔴",
            };

            println!(
                "{}. {} {} {}",
                i + 1,
                completed,
                level_indicator,
                proof.title
            );
            println!("   Level: {:?}", proof.mathematical_level);
            println!("   {}", proof.description);
            println!("   Prerequisites: {}", proof.prerequisites.join(", "));
            println!();
        }

        println!("Legend: ✅ Completed, 🔒 Available, 🟢 Undergraduate, 🟡 Graduate, 🟠 Advanced, 🔴 Research\n");

        let choice = get_user_input("Enter proof number to start, or 'back' to return: ")?;

        if choice.to_lowercase() == "back" {
            return Ok(());
        }

        if let Ok(index) = choice.parse::<usize>() {
            if index > 0 && index <= self.available_proofs.len() {
                let proof = &self.available_proofs[index - 1];
                self.start_proof_exploration(&proof.id)?;
            } else {
                println!("❌ Invalid proof number.");
            }
        }

        Ok(())
    }

    fn start_proof_exploration(&self, proofid: &str) -> Result<(), Box<dyn std::error::Error>> {
        let proof = self.get_proof_by_id(proofid).expect("Operation failed");

        println!("\n🎯 Starting Proof: {}", proof.title);
        println!("=================={}", "=".repeat(proof.title.len()));
        println!();

        // Display theorem
        println!("📜 Theorem Statement:");
        println!("{}", proof.main_theorem.statement);
        println!();

        println!("🧮 Mathematical Formulation:");
        println!("{}", proof.main_theorem.mathematical_formulation);
        println!();

        println!("📋 Conditions:");
        for condition in &proof.main_theorem.conditions {
            println!("  • {}", condition);
        }
        println!();

        println!("🌟 Significance:");
        println!("{}", proof.main_theorem.significance);
        println!();

        println!("📚 Historical Context:");
        println!("{}", proof.main_theorem.historical_context);
        println!();

        let start = get_user_input("Ready to begin the proof? (y/n): ")?;
        if start.to_lowercase() == "y" {
            self.run_proof_session(proof)?;
        }

        Ok(())
    }

    fn run_proof_session(&self, proof: &ProofModule) -> Result<(), Box<dyn std::error::Error>> {
        println!("\n🚀 Beginning Proof Session");
        println!("===========================\n");

        let session_start = Instant::now();
        let mut current_step = 0;
        let mut step_times = Vec::new();
        let mut user_responses = Vec::new();

        while current_step < proof.proof_steps.len() {
            let step = &proof.proof_steps[current_step];
            let step_start = Instant::now();

            println!(
                "📝 Step {} of {}: {}",
                current_step + 1,
                proof.proof_steps.len(),
                step.title
            );
            println!("{}", "=".repeat(50));
            println!();

            println!("💡 Explanation:");
            println!("{}", step.explanation);
            println!();

            println!("🧮 Mathematical Content:");
            println!("{}", step.mathematical_content);
            println!();

            println!("📖 Justification:");
            println!("{}", step.justification);
            println!();

            // Handle interactive elements
            for (i, element) in step.interactive_elements.iter().enumerate() {
                match self.handle_interactive_element(element, current_step, i)? {
                    ElementResult::Correct(response) => {
                        user_responses.push(UserResponse {
                            step_id: step.id.clone(),
                            response,
                            correct: true,
                            attempts: 1,
                            time_taken: step_start.elapsed(),
                        });
                    }
                    ElementResult::Incorrect(response, attempts) => {
                        user_responses.push(UserResponse {
                            step_id: step.id.clone(),
                            response,
                            correct: false,
                            attempts,
                            time_taken: step_start.elapsed(),
                        });
                    }
                    ElementResult::Skipped => {}
                }
            }

            // Verification if enabled
            if self.difficulty_settings.require_step_verification {
                if let Some(verification) = &step.verification_method {
                    self.run_verification(verification)?;
                }
            }

            step_times.push(step_start.elapsed());

            let continue_choice =
                get_user_input("\nPress Enter to continue to next step, or 'quit' to exit: ")?;
            if continue_choice.to_lowercase() == "quit" {
                break;
            }

            current_step += 1;
            println!("\n{}\n", "=".repeat(80));
        }

        if current_step == proof.proof_steps.len() {
            self.celebrate_completion(proof, session_start.elapsed())?;
        }

        Ok(())
    }

    fn handle_interactive_element(
        &self,
        element: &InteractiveElement,
        step: usize,
        element_index: usize,
    ) -> Result<ElementResult, Box<dyn std::error::Error>> {
        match element {
            InteractiveElement::UserInput {
                prompt,
                expected_answer,
                validation_tolerance,
                hints,
            } => {
                println!("❓ {}", prompt);

                let mut attempts = 0;
                loop {
                    attempts += 1;
                    let response = get_user_input("Your answer: ")?;

                    if self.validate_numerical_response(
                        &response,
                        expected_answer,
                        *validation_tolerance,
                    ) {
                        println!("✅ Correct!");
                        return Ok(ElementResult::Correct(response));
                    } else {
                        println!("❌ Not quite right.");

                        if attempts <= hints.len() && !hints.is_empty() {
                            println!("💡 Hint: {}", hints[attempts - 1]);
                        }

                        if attempts >= 3 {
                            println!("🔍 Expected answer: {}", expected_answer);
                            return Ok(ElementResult::Incorrect(response, attempts as u32));
                        }
                    }
                }
            }

            InteractiveElement::Visualization {
                title,
                description,
                plot_type,
            } => {
                println!("📊 Visualization: {}", title);
                println!("{}", description);
                self.create_visualization(plot_type)?;
                Ok(ElementResult::Skipped)
            }

            InteractiveElement::Computation {
                description,
                code_example,
                expected_result,
            } => {
                println!("💻 Computational Verification:");
                println!("{}", description);
                println!("Code example:");
                println!("```rust");
                println!("{}", code_example);
                println!("```");
                println!("Expected result: {}", expected_result);

                let run_code = get_user_input("Run computation? (y/n): ")?;
                if run_code.to_lowercase() == "y" {
                    self.execute_computational_example(code_example)?;
                }
                Ok(ElementResult::Skipped)
            }

            InteractiveElement::ConceptualQuestion {
                question,
                options,
                correct_answer,
                explanation,
            } => {
                println!("🤔 Conceptual Question:");
                println!("{}", question);
                println!();

                for (i, option) in options.iter().enumerate() {
                    println!("  {}. {}", (b'A' + i as u8) as char, option);
                }

                let response = get_user_input("\nYour answer (A, B, C, etc.): ")?;
                let answer_index = response.to_uppercase().chars().next().and_then(|c: char| {
                    if c.is_ascii_uppercase() {
                        Some((c as u8 - b'A') as usize)
                    } else {
                        None
                    }
                });

                if let Some(idx) = answer_index {
                    if idx == *correct_answer {
                        println!("✅ Correct!");
                        println!("💡 Explanation: {}", explanation);
                        return Ok(ElementResult::Correct(response));
                    } else {
                        println!("❌ Incorrect.");
                        println!("💡 Explanation: {}", explanation);
                        return Ok(ElementResult::Incorrect(response, 1));
                    }
                }

                Ok(ElementResult::Incorrect(response, 1))
            }
        }
    }

    fn validate_numerical_response(&self, response: &str, expected: &str, tolerance: f64) -> bool {
        if let (Ok(user_val), Ok(expected_val)) = (response.parse::<f64>(), expected.parse::<f64>())
        {
            (user_val - expected_val).abs() <= tolerance
        } else {
            response.trim() == expected.trim()
        }
    }

    fn create_visualization(&self, plottype: &PlotType) -> Result<(), Box<dyn std::error::Error>> {
        match plottype {
            PlotType::FunctionGraph { domain, functions } => {
                println!("📈 Function Graph:");
                println!("Domain: [{}, {}]", domain.0, domain.1);
                println!("Functions: {}", functions.join(", "));

                // Simple ASCII visualization
                self.create_ascii_plot(domain, functions)?;
            }

            PlotType::ComplexPlane { radius, function } => {
                println!("🌀 Complex Plane Visualization:");
                println!("Function: {}", function);
                println!("Radius: {}", radius);

                // Simple complex plane representation
                self.create_complex_plane_plot(*radius, function)?;
            }

            PlotType::Convergence { series, terms } => {
                println!("📊 Convergence Visualization:");
                println!("Series: {}", series);
                println!("Terms shown: {}", terms);

                self.demonstrate_convergence(series, *terms)?;
            }

            PlotType::AsymptoticComparison {
                exact,
                approximation,
                domain,
            } => {
                println!("📏 Asymptotic Comparison:");
                println!("Exact: {}", exact);
                println!("Approximation: {}", approximation);
                println!("Domain: [{}, {}]", domain.0, domain.1);

                self.compare_asymptotic_behavior(exact, approximation, domain)?;
            }
        }

        wait_for_user_input()?;
        Ok(())
    }

    fn create_ascii_plot(
        &self,
        domain: &(f64, f64),
        functions: &[String],
    ) -> Result<(), Box<dyn std::error::Error>> {
        println!("\nASCII Plot:");

        let num_points = 40;
        let step = (domain.1 - domain.0) / (num_points - 1) as f64;

        println!("x     f(x)     Visualization");
        println!("────────────────────────────────────");

        for i in 0..num_points {
            let x = domain.0 + i as f64 * step;

            // For demonstration, use gamma function
            let y = if functions.contains(&"gamma".to_string()) {
                gamma(x)
            } else if functions.contains(&"sin".to_string()) {
                x.sin()
            } else {
                x // Default to identity
            };

            let normalized_y = ((y + 2.0) * 10.0) as usize;
            let display_pos = normalized_y.min(20);

            let mut line = [' '; 21];
            line[10] = '|'; // Zero line
            if display_pos < line.len() {
                line[display_pos] = '●';
            }

            let display: String = line.iter().collect();
            println!("{:5.1} {:8.3}   {}", x, y, display);
        }

        Ok(())
    }

    fn create_complex_plane_plot(
        &self,
        radius: f64,
        function: &str,
    ) -> Result<(), Box<dyn std::error::Error>> {
        println!("\nComplex Plane (|z| ≤ {}):", radius);

        for row in 0..15 {
            for col in 0..30 {
                let re = (col as f64 - 15.0) * radius / 15.0;
                let im = (7.5 - row as f64) * radius / 7.5;
                let z = Complex64::new(re, im);

                if z.norm() <= radius {
                    // Evaluate function at z
                    let result = if function == "gamma" {
                        gamma_complex(z)
                    } else {
                        z // Default
                    };

                    let magnitude = result.norm();
                    let char = if magnitude < 1.0 {
                        '·'
                    } else if magnitude < 2.0 {
                        '▒'
                    } else if magnitude < 5.0 {
                        '▓'
                    } else {
                        '█'
                    };
                    print!("{}", char);
                } else {
                    print!(" ");
                }
            }
            println!();
        }

        Ok(())
    }

    fn demonstrate_convergence(
        &self,
        series: &str,
        terms: usize,
    ) -> Result<(), Box<dyn std::error::Error>> {
        println!("\nConvergence Demonstration for {}:", series);
        println!("Term    Partial Sum    Error");
        println!("──────────────────────────────");

        if series.contains("gamma") {
            // Demonstrate Stirling series convergence
            let x: f64 = 10.0;
            let exact = gamma(x);
            let mut sum = (x - 0.5) * x.ln() - x + 0.5 * (2.0 * PI).ln();

            println!(
                "{:4} {:12.6} {:10.2e}",
                0,
                sum.exp(),
                (sum.exp() - exact).abs()
            );

            // Add first few terms of Stirling series
            for k in 1..=terms.min(5) {
                let bernoulli_coeff = match k {
                    1 => 1.0 / 12.0,
                    2 => -1.0 / 360.0,
                    3 => 1.0 / 1260.0,
                    4 => -1.0 / 1680.0,
                    _ => 0.0,
                };

                sum += bernoulli_coeff / x.powi(2 * k as i32 - 1);
                println!(
                    "{:4} {:12.6} {:10.2e}",
                    k,
                    sum.exp(),
                    (sum.exp() - exact).abs()
                );
            }
        }

        Ok(())
    }

    fn compare_asymptotic_behavior(
        &self,
        exact: &str,
        approximation: &str,
        domain: &(f64, f64),
    ) -> Result<(), Box<dyn std::error::Error>> {
        println!("\nAsymptotic Comparison:");
        println!("x        Exact      Approx     Relative Error");
        println!("────────────────────────────────────────────");

        let num_points = 10;
        let step = (domain.1 - domain.0) / (num_points - 1) as f64;

        for i in 0..num_points {
            let x = domain.0 + i as f64 * step;

            let exact_val = if exact.contains("gamma") { gamma(x) } else { x };

            let approx_val = if approximation.contains("stirling") {
                // Stirling approximation
                ((2.0 * PI / x).sqrt() * (x / std::f64::consts::E).powf(x)).exp()
            } else {
                exact_val
            };

            let rel_error = if exact_val != 0.0 {
                ((exact_val - approx_val) / exact_val).abs()
            } else {
                0.0
            };

            println!(
                "{:8.1} {:10.4} {:10.4} {:12.2e}",
                x, exact_val, approx_val, rel_error
            );
        }

        Ok(())
    }

    fn execute_computational_example(&self, code: &str) -> Result<(), Box<dyn std::error::Error>> {
        println!("\n💻 Executing computational example...");

        // Parse and execute simple mathematical expressions
        if code.contains("gamma(") {
            // Extract gamma function calls and evaluate
            println!("Result: (executing gamma function computations)");

            // Example evaluations
            if code.contains("1/2") {
                let result = gamma(0.5);
                println!("gamma(1/2) = {:.10}", result);
                println!("√π = {:.10}", PI.sqrt());
                println!("Difference: {:.2e}", (result - PI.sqrt()).abs());
            }
        }

        Ok(())
    }

    fn run_verification(
        &self,
        verification: &VerificationMethod,
    ) -> Result<(), Box<dyn std::error::Error>> {
        println!("\n🔍 Verification Step:");

        match verification {
            VerificationMethod::NumericalCheck {
                description,
                test_cases,
            } => {
                println!("{}", description);
                println!("\nRunning numerical tests:");

                for (i, test_case) in test_cases.iter().enumerate() {
                    println!("Test {}: {}", i + 1, test_case.description);

                    // Evaluate based on test case
                    let actual = gamma(test_case.input); // Example
                    let error = (actual - test_case.expected).abs();
                    let passed = error <= test_case.tolerance;

                    println!("  Input: {}", test_case.input);
                    println!("  Expected: {:.6}", test_case.expected);
                    println!("  Actual: {:.6}", actual);
                    println!("  Error: {:.2e}", error);
                    println!("  Result: {}", if passed { "✅ PASS" } else { "❌ FAIL" });
                    println!();
                }
            }

            VerificationMethod::SymbolicVerification {
                description,
                identity,
            } => {
                println!("{}", description);
                println!("Identity to verify: {}", identity);
                println!("(Symbolic verification would be performed here)");
            }

            VerificationMethod::AsymptoticVerification {
                description,
                limit_expression,
            } => {
                println!("{}", description);
                println!("Limit expression: {}", limit_expression);
                println!("(Asymptotic verification would be performed here)");
            }
        }

        wait_for_user_input()?;
        Ok(())
    }

    fn celebrate_completion(
        &self,
        proof: &ProofModule,
        total_time: Duration,
    ) -> Result<(), Box<dyn std::error::Error>> {
        println!("\n🎉 Proof Completed Successfully!");
        println!("===================================\n");

        println!("🏆 Congratulations! You have successfully completed:");
        println!("📜 {}", proof.title);
        println!(
            "⏱️ Total _time: {:.1} minutes",
            total_time.as_secs_f64() / 60.0
        );
        println!("🎯 Mathematical level: {:?}", proof.mathematical_level);
        println!();

        println!("🌟 What you've learned:");
        println!("  • Understanding of: {}", proof.main_theorem.statement);
        println!("  • Complex analysis techniques");
        println!("  • Connection between discrete and continuous mathematics");
        println!();

        println!("🔗 Related topics you might explore next:");
        for related in &proof.related_proofs {
            println!("  • {}", related);
        }
        println!();

        println!("🏅 Achievement unlocked: Advanced Proof Explorer!");

        wait_for_user_input()?;
        Ok(())
    }

    fn continue_proof(&self) -> Result<(), Box<dyn std::error::Error>> {
        if let Some(_current) = &self.current_proof {
            println!("Continuing proof session...");
            // Implementation would restore the session state
        } else {
            println!("No proof session in progress.");
        }
        Ok(())
    }

    fn display_progress(&self) -> Result<(), Box<dyn std::error::Error>> {
        println!("\n📊 Your Learning Progress");
        println!("==========================\n");

        println!("🎯 Overall Statistics:");
        println!(
            "• Completed Proofs: {}",
            self.user_progress.completed_proofs.len()
        );
        println!(
            "• Mathematical Level: {:?}",
            self.user_progress.difficulty_level
        );
        println!(
            "• Total Study Time: {:.1} hours",
            self.user_progress
                .time_spent
                .values()
                .map(|d| d.as_secs_f64())
                .sum::<f64>()
                / 3600.0
        );

        println!("\n🏆 Achievements:");
        for achievement in &self.user_progress.achievements {
            println!("• {}", achievement);
        }

        if self.user_progress.achievements.is_empty() {
            println!("• Complete your first proof to earn achievements!");
        }

        wait_for_user_input()?;
        Ok(())
    }

    fn adjust_difficulty(&mut self) -> Result<(), Box<dyn std::error::Error>> {
        println!("\n⚙️ Difficulty Settings");
        println!("=======================\n");

        println!("Current settings:");
        println!(
            "• Mathematical rigor: {:?}",
            self.difficulty_settings.mathematical_rigor_level
        );
        println!(
            "• Show hints automatically: {}",
            self.difficulty_settings.show_hints_automatically
        );
        println!(
            "• Require step verification: {}",
            self.difficulty_settings.require_step_verification
        );
        println!(
            "• Include computational checks: {}",
            self.difficulty_settings.include_computational_checks
        );

        // Implementation for adjusting settings would go here

        wait_for_user_input()?;
        Ok(())
    }

    fn search_proofs(&self) -> Result<(), Box<dyn std::error::Error>> {
        println!("\n🔍 Search Proofs");
        println!("=================\n");

        let query = get_user_input("Enter search term (theorem name, topic, or keyword): ")?;

        println!("Search results for '{}':", query);

        for proof in &self.available_proofs {
            if proof.title.to_lowercase().contains(&query.to_lowercase())
                || proof
                    .description
                    .to_lowercase()
                    .contains(&query.to_lowercase())
                || proof
                    .prerequisites
                    .iter()
                    .any(|p| p.to_lowercase().contains(&query.to_lowercase()))
            {
                println!("• {} ({:?})", proof.title, proof.mathematical_level);
                println!("  {}", proof.description);
            }
        }

        wait_for_user_input()?;
        Ok(())
    }

    fn get_recommendations(&self) -> Result<(), Box<dyn std::error::Error>> {
        println!("\n💡 Personalized Recommendations");
        println!("=================================\n");

        println!("Based on your current level and progress, we recommend:");

        // Simple recommendation logic
        match self.user_progress.difficulty_level {
            MathematicalLevel::Undergraduate => {
                println!("• Start with the Error Function Series proof");
                println!("• Build foundation with basic special function properties");
            }
            MathematicalLevel::Graduate => {
                println!("• Explore the Gamma Function Reflection Formula");
                println!("• Learn complex analysis techniques through Bessel functions");
            }
            _ => {
                println!("• Challenge yourself with Stirling's Asymptotic Series");
                println!("• Explore advanced topics in special function theory");
            }
        }

        wait_for_user_input()?;
        Ok(())
    }

    fn reference_library(&self) -> Result<(), Box<dyn std::error::Error>> {
        println!("\n📖 Mathematical Reference Library");
        println!("===================================\n");

        println!("Available references:");
        println!("1. 📚 Special Functions Handbook");
        println!("2. 🧮 Complex Analysis Review");
        println!("3. 📊 Asymptotic Methods");
        println!("4. 🔢 Number Theory Connections");
        println!("5. 🌐 Historical Development");

        let choice = get_user_input("Select reference (1-5) or 'back': ")?;

        match choice.as_str() {
            "1" => self.show_special_functions_reference()?,
            "2" => self.show_complex_analysis_reference()?,
            "3" => self.show_asymptotic_methods_reference()?,
            "4" => self.show_number_theory_reference()?,
            "5" => self.show_historical_reference()?,
            _ => return Ok(()),
        }

        Ok(())
    }

    fn show_special_functions_reference(&self) -> Result<(), Box<dyn std::error::Error>> {
        println!("\n📚 Special Functions Handbook");
        println!("==============================\n");

        println!("🎲 Gamma Function:");
        println!("• Definition: Γ(z) = ∫₀^∞ t^(z-1) e^(-t) dt");
        println!("• Recurrence: Γ(z+1) = z·Γ(z)");
        println!("• Reflection: Γ(z)Γ(1-z) = π/sin(πz)");
        println!("• Special values: Γ(1/2) = √π, Γ(n) = (n-1)!");
        println!();

        println!("🌊 Bessel Functions:");
        println!("• Differential equation: x²y'' + xy' + (x²-ν²)y = 0");
        println!("• Generating function: exp(x(t-1/t)/2) = Σ J_n(x) t^n");
        println!("• Orthogonality: ∫₀¹ J_n(αx)J_n(βx)x dx = 0 if α ≠ β");
        println!();

        println!("📊 Error Function:");
        println!("• Definition: erf(x) = (2/√π) ∫₀^x e^(-t²) dt");
        println!("• Complement: erfc(x) = 1 - erf(x)");
        println!("• Series: erf(x) = (2/√π) Σ (-1)^n x^(2n+1) / (n!(2n+1))");

        wait_for_user_input()?;
        Ok(())
    }

    fn show_complex_analysis_reference(&self) -> Result<(), Box<dyn std::error::Error>> {
        println!("\n🧮 Complex Analysis Review");
        println!("============================\n");

        println!("🔍 Key Concepts for Special Functions:");
        println!("• Analytic continuation extends functions beyond original domains");
        println!("• Residue calculus enables evaluation of complex integrals");
        println!("• Branch cuts and Riemann surfaces handle multivalued functions");
        println!("• Steepest descent method provides asymptotic expansions");
        println!();

        println!("🎯 Applications:");
        println!("• Gamma function poles and residues");
        println!("• Bessel function integral representations");
        println!("• Hypergeometric function transformation formulas");

        wait_for_user_input()?;
        Ok(())
    }

    fn show_asymptotic_methods_reference(&self) -> Result<(), Box<dyn std::error::Error>> {
        println!("\n📊 Asymptotic Methods");
        println!("======================\n");

        println!("🎯 Key Methods:");
        println!("• Stirling's method for factorial-like functions");
        println!("• Saddle point method for exponential integrals");
        println!("• WKB approximation for differential equations");
        println!("• Euler-Maclaurin formula for sum-integral relationships");

        wait_for_user_input()?;
        Ok(())
    }

    fn show_number_theory_reference(&self) -> Result<(), Box<dyn std::error::Error>> {
        println!("\n🔢 Number Theory Connections");
        println!("==============================\n");

        println!("🌟 Special Functions in Number Theory:");
        println!("• Riemann zeta function: ζ(s) = Σ 1/n^s");
        println!("• Gamma function and factorials");
        println!("• Bernoulli numbers in asymptotic series");
        println!("• Euler's constant γ in logarithmic integrals");

        wait_for_user_input()?;
        Ok(())
    }

    fn show_historical_reference(&self) -> Result<(), Box<dyn std::error::Error>> {
        println!("\n🌐 Historical Development");
        println!("==========================\n");

        println!("📅 Timeline:");
        println!("• 1720s: Euler begins work on factorial extension");
        println!("• 1812: Legendre introduces Γ notation");
        println!("• 1824: Bessel develops his functions for astronomy");
        println!("• 1856: Weierstrass proves infinite product formulas");
        println!("• 1900s: Modern asymptotic theory develops");

        wait_for_user_input()?;
        Ok(())
    }

    fn save_progress(&self) -> Result<(), Box<dyn std::error::Error>> {
        println!("💾 Saving progress...");
        // In a real implementation, this would save to a file or database
        println!("✅ Progress saved successfully!");
        Ok(())
    }

    fn get_proof_by_id(&self, id: &str) -> Option<&ProofModule> {
        self.available_proofs.iter().find(|p| p.id == id)
    }
}

#[allow(dead_code)]
#[derive(Debug)]
enum ElementResult {
    Correct(String),
    Incorrect(String, u32),
    Skipped,
}

// Helper functions for creating proof steps
#[allow(dead_code)]
fn create_reflection_formula_proof_steps() -> Vec<ProofStep> {
    vec![
        ProofStep {
            id: "step1".to_string(),
            title: "Beta Function Connection".to_string(),
            explanation: "We start by connecting the gamma function to the beta function, which will provide a pathway to the trigonometric identity.".to_string(),
            mathematical_content: "B(z, 1-z) = ∫₀¹ t^(z-1)(1-t)^(-z) dt = Γ(z)Γ(1-z)/Γ(1) = Γ(z)Γ(1-z)".to_string(),
            justification: "This follows from the integral representation of the beta function and the fact that Γ(1) = 1.".to_string(),
            interactive_elements: vec![
                InteractiveElement::UserInput {
                    prompt: "What is the value of Γ(1)?".to_string(),
                    expected_answer: "1".to_string(),
                    validation_tolerance: 0.01,
                    hints: vec!["Remember that Γ(n) = (n-1)! for positive integers".to_string()],
                },
            ],
            verification_method: Some(VerificationMethod::NumericalCheck {
                description: "Verify the beta function identity numerically".to_string(),
                test_cases: vec![
                    TestCase {
                        input: 0.5,
                        expected: PI,
                        tolerance: 1e-10,
                        description: "B(1/2, 1/2) should equal π".to_string(),
                    },
                ],
            }),
            common_confusions: vec![
                "Don't confuse B(z,w) with Γ(z)Γ(w)".to_string(),
            ],
            hints: vec![
                "The beta function is symmetric: B(z,w) = B(w,z)".to_string(),
            ],
        },
        ProofStep {
            id: "step2".to_string(),
            title: "Integral Transformation".to_string(),
            explanation: "We transform the beta integral using substitution to make it amenable to complex analysis techniques.".to_string(),
            mathematical_content: "Using t = u/(1+u), we get B(z,1-z) = ∫₀^∞ u^(z-1)/(1+u) du".to_string(),
            justification: "This substitution transforms the finite interval [0,1] to the infinite interval [0,∞), which is better suited for complex analysis.".to_string(),
            interactive_elements: vec![
                InteractiveElement::Computation {
                    description: "Verify the substitution t = u/(1+u) transforms the integral correctly".to_string(),
                    code_example: "// When t = u/(1+u), then dt = du/(1+u)²\n// and (1-t) = 1/(1+u)".to_string(),
                    expected_result: "The Jacobian gives the correct transformation".to_string(),
                },
            ],
            verification_method: None,
            common_confusions: vec![
                "Be careful with the Jacobian of the transformation".to_string(),
            ],
            hints: vec![
                "Work out dt in terms of du carefully".to_string(),
            ],
        },
        // Additional steps would be added here...
    ]
}

#[allow(dead_code)]
fn create_stirling_proof_steps() -> Vec<ProofStep> {
    vec![
        ProofStep {
            id: "stirling1".to_string(),
            title: "Euler-Maclaurin Setup".to_string(),
            explanation: "We begin with the Euler-Maclaurin formula to relate sums to integrals, which is key to understanding the asymptotic behavior of the gamma function.".to_string(),
            mathematical_content: "∑_{k=1}^n f(k) = ∫₁ⁿ f(x)dx + ½[f(1)+f(n)] + ∑_{k=1}^p B₂ₖ/(2k)! [f^(2k-1)(n)-f^(2k-1)(1)] + Rp".to_string(),
            justification: "The Euler-Maclaurin formula provides the connection between discrete sums and continuous integrals through Bernoulli numbers.".to_string(),
            interactive_elements: vec![
                InteractiveElement::ConceptualQuestion {
                    question: "What is the key insight behind the Euler-Maclaurin formula?".to_string(),
                    options: vec![
                        "It approximates integrals with sums".to_string(),
                        "It relates sums to integrals plus correction terms".to_string(),
                        "It computes derivatives numerically".to_string(),
                        "It solves differential equations".to_string(),
                    ],
                    correct_answer: 1,
                    explanation: "The formula expresses sums as integrals plus correction terms involving Bernoulli numbers".to_string(),
                },
            ],
            verification_method: None,
            common_confusions: vec![
                "The Bernoulli numbers appear with specific powers".to_string(),
            ],
            hints: vec![
                "Think about approximating discrete with continuous".to_string(),
            ],
        },
    ]
}

#[allow(dead_code)]
fn create_bessel_generating_proof_steps() -> Vec<ProofStep> {
    vec![
        ProofStep {
            id: "bessel1".to_string(),
            title: "Exponential Series Expansion".to_string(),
            explanation: "We expand the exponential function exp(x(t-1/t)/2) in powers of t to identify the coefficients as Bessel functions.".to_string(),
            mathematical_content: "exp(x(t-1/t)/2) = exp(xt/2)exp(-x/(2t)) = ∑_{m=0}^∞ (xt/2)^m/m! ∑_{n=0}^∞ (-x/(2t))^n/n!".to_string(),
            justification: "We use the product of two exponential series and collect terms by powers of t.".to_string(),
            interactive_elements: vec![
                InteractiveElement::Visualization {
                    title: "Bessel Function Generating Function".to_string(),
                    description: "Visualize how the generating function creates oscillating Bessel functions".to_string(),
                    plot_type: PlotType::FunctionGraph {
                        domain: (0.0, 20.0),
                        functions: vec!["J_0(x)".to_string(), "J_1(x)".to_string()],
                    },
                },
            ],
            verification_method: Some(VerificationMethod::NumericalCheck {
                description: "Verify the generating function for small values".to_string(),
                test_cases: vec![
                    TestCase {
                        input: 1.0,
                        expected: 0.765198,
                        tolerance: 1e-5,
                        description: "J_0(1) from generating function".to_string(),
                    },
                ],
            }),
            common_confusions: vec![
                "Keep track of positive and negative powers of t".to_string(),
            ],
            hints: vec![
                "Use the Cauchy product formula for the series".to_string(),
            ],
        },
    ]
}

#[allow(dead_code)]
fn create_error_function_proof_steps() -> Vec<ProofStep> {
    vec![
        ProofStep {
            id: "erf1".to_string(),
            title: "Series Expansion of Exponential".to_string(),
            explanation: "We start with the series expansion of e^(-t²) and integrate term by term to obtain the error function series.".to_string(),
            mathematical_content: "e^(-t²) = ∑_{n=0}^∞ (-1)^n t^(2n)/n!".to_string(),
            justification: "This is the standard power series expansion of the exponential function with argument -t².".to_string(),
            interactive_elements: vec![
                InteractiveElement::UserInput {
                    prompt: "What is the coefficient of t⁴ in the expansion of e^(-t²)?".to_string(),
                    expected_answer: "1/2".to_string(),
                    validation_tolerance: 0.01,
                    hints: vec!["Use the formula (-1)^n/n! with n=2".to_string()],
                },
                InteractiveElement::Visualization {
                    title: "Error Function Convergence".to_string(),
                    description: "Watch how the series converges to the error function".to_string(),
                    plot_type: PlotType::Convergence {
                        series: "erf(x) series".to_string(),
                        terms: 10,
                    },
                },
            ],
            verification_method: Some(VerificationMethod::NumericalCheck {
                description: "Verify series convergence for erf(1)".to_string(),
                test_cases: vec![
                    TestCase {
                        input: 1.0,
                        expected: 0.842701,
                        tolerance: 1e-5,
                        description: "erf(1) from series".to_string(),
                    },
                ],
            }),
            common_confusions: vec![
                "Don't forget the 2/√π factor in the error function definition".to_string(),
            ],
            hints: vec![
                "Integrate each term of the series separately".to_string(),
            ],
        },
    ]
}

#[allow(dead_code)]
fn create_reflection_formula_examples() -> Vec<ComputationalExample> {
    vec![
        ComputationalExample {
            title: "Verifying the Reflection Formula".to_string(),
            description: "Numerical verification of Γ(z)Γ(1-z) = π/sin(πz)".to_string(),
            code: "let z = 0.3;\nlet left_side = gamma(z) * gamma(1.0 - z);\nlet right_side = PI / (PI * z).sin();\nprintln!(\"Left: {}, Right: {}\", left_side, right_side);".to_string(),
            expected_output: "Both sides should be approximately equal".to_string(),
            learning_objective: "Understand how the reflection formula connects gamma and trigonometric functions".to_string(),
        },
    ]
}

#[allow(dead_code)]
fn create_stirling_examples() -> Vec<ComputationalExample> {
    vec![
        ComputationalExample {
            title: "Stirling Approximation Accuracy".to_string(),
            description: "Compare Stirling approximation with exact gamma values".to_string(),
            code: "for n in 5..20 {\n  let exact = gamma(n as f64);\n  let stirling = stirling_approx(n as f64);\n  println!(\"n={}: exact={:.2e}, stirling={:.2e}\", n, exact, stirling);\n}".to_string(),
            expected_output: "Stirling approximation becomes increasingly accurate for larger n".to_string(),
            learning_objective: "See how asymptotic approximations work in practice".to_string(),
        },
    ]
}

#[allow(dead_code)]
fn create_bessel_generating_examples() -> Vec<ComputationalExample> {
    vec![
        ComputationalExample {
            title: "Generating Function Verification".to_string(),
            description: "Verify that the generating function produces correct Bessel function values".to_string(),
            code: "let x = 2.0;\nlet t = 1.0;\nlet gen_func = (x * (t - 1.0/t) / 2.0).exp();\nlet series_sum = j0(x) + j1(x) * t + jminus1(x) / t;\nprintln!(\"Gen func: {}, Series: {}\", gen_func, series_sum);".to_string(),
            expected_output: "The generating function equals the series sum".to_string(),
            learning_objective: "Understand how generating functions encode infinite families of functions".to_string(),
        },
    ]
}

#[allow(dead_code)]
fn create_error_function_examples() -> Vec<ComputationalExample> {
    vec![
        ComputationalExample {
            title: "Error Function Series Convergence".to_string(),
            description: "Demonstrate convergence of the error function power series".to_string(),
            code: "let x = 1.0;\nlet mut sum = 0.0;\nfor n in 0..20 {\n  let term = (-1.0_f64).powi(n) * x.powi(2*n+1) / (factorial(n) * (2*n+1));\n  sum += term;\n  println!(\"n={}: term={:.6e}, sum={:.6}\", n, term, 2.0/PI.sqrt() * sum);\n}".to_string(),
            expected_output: "Series converges to erf(1) ≈ 0.8427".to_string(),
            learning_objective: "See how infinite series represent transcendental functions".to_string(),
        },
    ]
}

#[allow(dead_code)]
fn wait_for_user_input() -> Result<(), Box<dyn std::error::Error>> {
    get_user_input("Press Enter to continue...")?;
    Ok(())
}

#[allow(dead_code)]
fn get_user_input(prompt: &str) -> Result<String, Box<dyn std::error::Error>> {
    print!("{}", prompt);
    io::stdout().flush()?;
    let mut input = String::new();
    io::stdin().read_line(&mut input)?;
    Ok(input.trim().to_string())
}

#[allow(dead_code)]
fn gamma_complex(z: Complex64) -> Complex64 {
    // Complex gamma function implementation using Lanczos approximation
    if z.re > 0.0 {
        // Use Lanczos approximation for positive real part
        let g = 7.0;
        let coef = [
            0.999_999_999_999_809_9,
            676.5203681218851,
            -1259.1392167224028,
            771.323_428_777_653_1,
            -176.615_029_162_140_6,
            12.507343278686905,
            -0.13857109526572012,
            9.984_369_578_019_572e-6,
            1.5056327351493116e-7,
        ];

        let z_shifted = z - Complex64::new(1.0, 0.0);
        let mut x = Complex64::new(coef[0], 0.0);

        for i in 1..coef.len() {
            x += Complex64::new(coef[i], 0.0) / (z_shifted + Complex64::new(i as f64, 0.0));
        }

        let t = z_shifted + Complex64::new(g + 0.5, 0.0);
        let sqrt_2pi = Complex64::new((2.0 * PI).sqrt(), 0.0);

        sqrt_2pi * t.powc(z_shifted + Complex64::new(0.5, 0.0)) * (-t).exp() * x
    } else {
        // Use reflection formula for negative real part
        let pi_z = Complex64::new(PI, 0.0) * z;
        let sin_pi_z = pi_z.sin();

        if sin_pi_z.norm() < 1e-15 {
            // Pole at negative integer
            Complex64::new(f64::INFINITY, 0.0)
        } else {
            Complex64::new(PI, 0.0) / (sin_pi_z * gamma_complex(Complex64::new(1.0, 0.0) - z))
        }
    }
}

#[allow(dead_code)]
fn main() -> Result<(), Box<dyn std::error::Error>> {
    let mut explorer = ProofExplorer::new();
    explorer.run()
}
