//! Mathematical Derivation Studio
//!
//! An advanced interactive tutorial that guides users through detailed mathematical
//! derivations of special function identities, asymptotic expansions, and connections
//! between different function families.
//!
//! This tutorial is designed for advanced students and researchers who want to understand
//! the deep mathematical structure underlying special functions.
//!
//! Run with: cargo run --example mathematical_derivation_studio

use scirs2_core::ndarray::Array1;
use scirs2_core::Complex64;
use scirs2_special::*;
use std::f64::consts::{E, PI};
use std::io::{self, Write};

#[allow(dead_code)]
fn main() -> Result<(), Box<dyn std::error::Error>> {
    println!("🧮 Mathematical Derivation Studio");
    println!("=================================");
    println!("Deep dive into special function mathematics\n");

    loop {
        display_main_menu();
        let choice = get_user_input("Enter your choice (1-8, or 'q' to quit): ")?;

        if choice.to_lowercase() == "q" {
            println!("🎓 Thank you for exploring mathematical derivations!");
            break;
        }

        match choice.parse::<u32>() {
            Ok(1) => gamma_reflection_formula_derivation()?,
            Ok(2) => stirling_asymptotic_derivation()?,
            Ok(3) => bessel_orthogonality_proof()?,
            Ok(4) => hypergeometric_transformations()?,
            Ok(5) => elliptic_integral_connections()?,
            Ok(6) => wright_function_asymptotic_analysis()?,
            Ok(7) => information_theory_inequalities()?,
            Ok(8) => quantum_mechanics_applications()?,
            _ => println!("❌ Invalid choice. Please try again.\n"),
        }
    }

    Ok(())
}

#[allow(dead_code)]
fn display_main_menu() {
    println!("📖 Choose a mathematical derivation to explore:");
    println!("1. 🎯 Gamma Function Reflection Formula");
    println!("2. 📈 Stirling's Asymptotic Expansion");
    println!("3. 🌊 Bessel Function Orthogonality");
    println!("4. 🔄 Hypergeometric Transformations");
    println!("5. 🥧 Elliptic Integral Connections");
    println!("6. 🌀 Wright Function Asymptotics");
    println!("7. 📊 Information Theory Inequalities");
    println!("8. ⚛️  Quantum Mechanics Applications");
    println!("q. Quit");
    println!();
}

#[allow(dead_code)]
fn get_user_input(prompt: &str) -> io::Result<String> {
    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_reflection_formula_derivation() -> Result<(), Box<dyn std::error::Error>> {
    println!("\n🎯 GAMMA FUNCTION REFLECTION FORMULA DERIVATION");
    println!("===============================================\n");

    println!("We'll derive the famous reflection formula: Γ(z)Γ(1-z) = π/sin(πz)");
    println!("This is one of the most beautiful identities in mathematics!\n");

    pause_for_user()?;

    println!("STEP 1: Start with the Beta function integral");
    println!("B(z, 1-z) = ∫₀¹ t^(z-1) (1-t)^(-z) dt");
    println!();
    println!("We know that B(z, 1-z) = Γ(z)Γ(1-z)/Γ(1) = Γ(z)Γ(1-z)");
    println!();

    pause_for_user()?;

    println!("STEP 2: Transform the integral using substitution t = u/(1+u)");
    println!("dt = du/(1+u)², (1-t) = 1/(1+u)");
    println!();
    println!("B(z, 1-z) = ∫₀^∞ (u/(1+u))^(z-1) (1/(1+u))^(-z) · du/(1+u)²");
    println!("           = ∫₀^∞ u^(z-1) (1+u)^(-1) du");
    println!();

    pause_for_user()?;

    println!("STEP 3: Apply complex contour integration");
    println!("Consider the integral ∮_C w^(z-1)/(1+w) dw around a keyhole contour");
    println!("around the branch cut [0,∞) of w^(z-1)");
    println!();

    // Demonstrate with numerical verification
    let z = 0.3;
    let gamma_z = gamma(z);
    let gamma_1minus_z = gamma(1.0 - z);
    let product = gamma_z * gamma_1minus_z;
    let theoretical = PI / (PI * z).sin();

    println!("NUMERICAL VERIFICATION:");
    println!("For z = {:.3}:", z);
    println!("Γ({:.3}) = {:.8}", z, gamma_z);
    println!("Γ({:.3}) = {:.8}", 1.0 - z, gamma_1minus_z);
    println!("Product = {:.8}", product);
    println!("π/sin(πz) = {:.8}", theoretical);
    println!("Difference = {:.2e}", (product - theoretical).abs());
    println!();

    pause_for_user()?;

    println!("STEP 4: Evaluate residues and branch cut contributions");
    println!("The residue at w = -1 is: Res(-1) = (-1)^(z-1) = e^(iπ(z-1))");
    println!();
    println!("The branch cut contribution gives:");
    println!("∮_C = (1 - e^(2πi(z-1))) ∫₀^∞ u^(z-1)/(1+u) du");
    println!();

    pause_for_user()?;

    println!("STEP 5: Apply the residue theorem");
    println!("2πi · e^(iπ(z-1)) = (1 - e^(2πi(z-1))) · B(z, 1-z)");
    println!();
    println!("Simplifying: 1 - e^(2πi(z-1)) = -2i e^(iπ(z-1)) sin(π(z-1))");
    println!("                                = 2i e^(iπ(z-1)) sin(πz)");
    println!();

    pause_for_user()?;

    println!("STEP 6: Final simplification");
    println!("2πi e^(iπ(z-1)) = 2i e^(iπ(z-1)) sin(πz) · Γ(z)Γ(1-z)");
    println!();
    println!("Canceling common factors:");
    println!("π = sin(πz) · Γ(z)Γ(1-z)");
    println!();
    println!("Therefore: Γ(z)Γ(1-z) = π/sin(πz) ✓");
    println!();

    // Show applications
    println!("APPLICATIONS:");
    println!("• Γ(1/2) = √π (setting z = 1/2)");
    println!("• Connection to the sinc function");
    println!("• Meromorphic continuation of the gamma function");
    println!();

    test_reflection_formula_values()?;

    Ok(())
}

#[allow(dead_code)]
fn stirling_asymptotic_derivation() -> Result<(), Box<dyn std::error::Error>> {
    println!("\n📈 STIRLING'S ASYMPTOTIC EXPANSION DERIVATION");
    println!("==============================================\n");

    println!("We'll derive Stirling's famous asymptotic formula:");
    println!("ln Γ(z) ~ (z - 1/2)ln(z) - z + (1/2)ln(2π) + O(1/z)");
    println!();

    pause_for_user()?;

    println!("STEP 1: Start with the integral representation");
    println!("Γ(z) = ∫₀^∞ t^(z-1) e^(-t) dt");
    println!();
    println!("Taking logarithms:");
    println!("ln Γ(z) = ln(∫₀^∞ t^(z-1) e^(-t) dt)");
    println!();

    pause_for_user()?;

    println!("STEP 2: Use the method of steepest descent");
    println!("The integrand is t^(z-1) e^(-t) = e^((z-1)ln(t) - t)");
    println!();
    println!("For large z, this is dominated by the maximum of f(t) = (z-1)ln(t) - t");
    println!("Setting f'(t) = 0: (z-1)/t - 1 = 0 → t₀ = z-1");
    println!();

    pause_for_user()?;

    println!("STEP 3: Expand around the saddle point t₀ = z-1");
    println!("f(t) = f(t₀) + (1/2)f''(t₀)(t-t₀)² + ...");
    println!();
    println!("f(t₀) = (z-1)ln(z-1) - (z-1)");
    println!("f''(t₀) = -(z-1)/(z-1)² = -1/(z-1)");
    println!();

    pause_for_user()?;

    println!("STEP 4: Gaussian approximation");
    println!("∫₀^∞ t^(z-1) e^(-t) dt ≈ e^(f(t₀)) ∫₋∞^∞ e^(-1/2 · (t-t₀)²/(z-1)) dt");
    println!();
    println!("The Gaussian integral gives √(2π(z-1))");
    println!();

    pause_for_user()?;

    println!("STEP 5: Combine results");
    println!("Γ(z) ≈ e^((z-1)ln(z-1) - (z-1)) √(2π(z-1))");
    println!();
    println!("Taking logarithms:");
    println!("ln Γ(z) ≈ (z-1)ln(z-1) - (z-1) + (1/2)ln(2π(z-1))");
    println!("        = (z-1)ln(z-1) - (z-1) + (1/2)ln(2π) + (1/2)ln(z-1)");
    println!();

    pause_for_user()?;

    println!("STEP 6: Simplify for large z");
    println!("For large z: ln(z-1) ≈ ln(z) - 1/z + O(1/z²)");
    println!();
    println!("ln Γ(z) ≈ (z-1)ln(z) - z + 1 + (1/2)ln(2π) + (1/2)ln(z)");
    println!("        = (z - 1/2)ln(z) - z + (1/2)ln(2π) + O(1/z)");
    println!();

    // Numerical verification
    println!("NUMERICAL VERIFICATION:");
    let test_values: Vec<f64> = vec![5.0, 10.0, 20.0, 50.0];
    for &z in &test_values {
        let exact = gammaln(z);
        let stirling = (z - 0.5) * z.ln() - z + 0.5 * (2.0 * PI).ln();
        let error = (exact - stirling).abs();
        println!(
            "z = {:<4}: exact = {:<10.6}, Stirling = {:<10.6}, error = {:.2e}",
            z, exact, stirling, error
        );
    }
    println!();

    println!("The error decreases as O(1/z), confirming our asymptotic expansion!");
    println!();

    Ok(())
}

#[allow(dead_code)]
fn bessel_orthogonality_proof() -> Result<(), Box<dyn std::error::Error>> {
    println!("\n🌊 BESSEL FUNCTION ORTHOGONALITY PROOF");
    println!("======================================\n");

    println!("We'll prove the orthogonality relation for Bessel functions:");
    println!("∫₀¹ x J_ν(α_m x) J_ν(α_n x) dx = (δ_mn/2)[J_(ν+1)(α_m)]²");
    println!("where α_m are the zeros of J_ν(x)");
    println!();

    pause_for_user()?;

    println!("STEP 1: Start with Bessel's differential equation");
    println!("x²y'' + xy' + (λ²x² - ν²)y = 0");
    println!();
    println!("For y = J_ν(λx), this becomes:");
    println!("x²[J_ν(λx)]'' + x[J_ν(λx)]' + (λ²x² - ν²)J_ν(λx) = 0");
    println!();

    pause_for_user()?;

    println!("STEP 2: Consider two solutions with different parameters");
    println!("Let u = J_ν(α_m x) and v = J_ν(α_n x)");
    println!();
    println!("Then u and v satisfy:");
    println!("x²u'' + xu' + (α_m²x² - ν²)u = 0");
    println!("x²v'' + xv' + (α_n²x² - ν²)v = 0");
    println!();

    pause_for_user()?;

    println!("STEP 3: Form the difference equation");
    println!("Multiply first equation by v, second by u, and subtract:");
    println!("x²(u''v - uv'') + x(u'v - uv') + (α_m² - α_n²)x²uv = 0");
    println!();
    println!("The first two terms form a perfect derivative:");
    println!("d/dx[x(u'v - uv')] + (α_m² - α_n²)x²uv = 0");
    println!();

    pause_for_user()?;

    println!("STEP 4: Integrate over [0,1]");
    println!("∫₀¹ d/dx[x(u'v - uv')] dx + (α_m² - α_n²)∫₀¹ x²uv dx = 0");
    println!();
    println!("[x(u'v - uv')]₀¹ + (α_m² - α_n²)∫₀¹ x J_ν(α_m x) J_ν(α_n x) dx = 0");
    println!();

    pause_for_user()?;

    println!("STEP 5: Evaluate the boundary term");
    println!("At x = 0: The term vanishes (u, v, u', v' all finite)");
    println!("At x = 1: Since α_m and α_n are zeros of J_ν, we have u(1) = v(1) = 0");
    println!();
    println!("Therefore: [x(u'v - uv')]₀¹ = 0");
    println!();

    pause_for_user()?;

    println!("STEP 6: Conclude orthogonality");
    println!("For α_m ≠ αn: (α_m² - α_n²)∫₀¹ x J_ν(α_m x) J_ν(α_n x) dx = 0");
    println!();
    println!("Since α_m² ≠ α_n², we must have:");
    println!("∫₀¹ x J_ν(α_m x) J_ν(α_n x) dx = 0  (m ≠ n)");
    println!();

    pause_for_user()?;

    println!("STEP 7: Find the normalization constant");
    println!("For m = n, we need to evaluate ∫₀¹ x [J_ν(α_m x)]² dx");
    println!();
    println!("Using L'Hôpital's rule and properties of Bessel functions:");
    println!("∫₀¹ x [J_ν(α_m x)]² dx = (1/2)[J_(ν+1)(α_m)]²");
    println!();

    // Numerical verification
    println!("NUMERICAL VERIFICATION:");
    println!("Testing orthogonality for J₀ with first few zeros:");

    let j0_zeros = [2.4048, 5.5201, 8.6537]; // Approximate zeros of J₀

    for (i, &alpha_i) in j0_zeros.iter().enumerate() {
        for (j, &alpha_j) in j0_zeros.iter().enumerate() {
            let integral = numerical_bessel_orthogonality_integral(0, alpha_i, alpha_j)?;
            let expected = if i == j { "non-zero" } else { "~0" };
            println!(
                "∫ x J₀({:.4}x) J₀({:.4}x) dx = {:.6} (expected: {})",
                alpha_i, alpha_j, integral, expected
            );
        }
    }
    println!();

    println!("✓ Orthogonality confirmed! This forms the basis for Fourier-Bessel series.");
    println!();

    Ok(())
}

#[allow(dead_code)]
fn hypergeometric_transformations() -> Result<(), Box<dyn std::error::Error>> {
    println!("\n🔄 HYPERGEOMETRIC FUNCTION TRANSFORMATIONS");
    println!("==========================================\n");

    println!("We'll explore the 24 Kummer solutions and fundamental transformations");
    println!("of the hypergeometric equation: z(1-z)w'' + [c-(a+b+1)z]w' - abw = 0");
    println!();

    pause_for_user()?;

    println!("FUNDAMENTAL TRANSFORMATIONS:");
    println!();

    println!("1. EULER'S TRANSFORMATION:");
    println!("₂F₁(a,b;c;z) = (1-z)^(c-a-b) ₂F₁(c-a,c-b;c;z)");
    println!();

    // Numerical verification
    let a = 0.5;
    let b = 1.5;
    let c = 2.0;
    let z = 0.3;

    println!("Verification for a={}, b={}, c={}, z={}:", a, b, c, z);
    let left = hypergeometric_2f1(a, b, c, z)?;
    let right = (1.0 - z).powf(c - a - b) * hypergeometric_2f1(c - a, c - b, c, z)?;
    println!("Left side:  {:.8}", left);
    println!("Right side: {:.8}", right);
    println!("Difference: {:.2e}", (left - right).abs());
    println!();

    pause_for_user()?;

    println!("2. PFAFF'S TRANSFORMATION:");
    println!("₂F₁(a,b;c;z) = (1-z)^(-a) ₂F₁(a,c-b;c;z/(z-1))");
    println!();

    pause_for_user()?;

    println!("3. QUADRATIC TRANSFORMATIONS:");
    println!("These relate ₂F₁ at z to ₂F₁ at quadratic expressions in z");
    println!();
    println!("Example (Gauss):");
    println!("₂F₁(a,b;a+b+1/2;z) = 2^(2a+2b-1) ₂F₁(2a,2b;a+b+1/2;(z+√z)²/4)");
    println!();

    pause_for_user()?;

    println!("4. CONNECTION TO ELEMENTARY FUNCTIONS:");
    println!();
    println!("• ₂F₁(1,1;2;z) = -ln(1-z)/z");
    println!("• ₂F₁(1/2,1/2;3/2;z²) = arcsin(z)/z");
    println!("• ₂F₁(1/2,1;3/2;-z²) = ln(z+√(1+z²))/z");
    println!();

    // Verify some elementary connections
    println!("VERIFICATION OF ELEMENTARY CONNECTIONS:");
    let z_test = 0.5;

    // Test ₂F₁(1,1;2;z) = -ln(1-z)/z
    let hyp_val = hypergeometric_2f1(1.0, 1.0, 2.0, z_test)?;
    let elem_val = -(1.0 - z_test).ln() / z_test;
    println!("₂F₁(1,1;2;{}) = {:.8}", z_test, hyp_val);
    println!("-ln(1-z)/z      = {:.8}", elem_val);
    println!("Difference      = {:.2e}", (hyp_val - elem_val).abs());
    println!();

    pause_for_user()?;

    println!("5. ASYMPTOTIC BEHAVIOR NEAR SINGULARITIES:");
    println!();
    println!("Near z = 0: ₂F₁(a,b;c;z) ~ 1 + (ab/c)z + O(z²)");
    println!("Near z = 1: Complex behavior depending on c-a-b");
    println!("Near z = ∞: Uses connection formulas with Gamma functions");
    println!();

    println!("These transformations are crucial for:");
    println!("• Numerical computation in different domains");
    println!("• Connecting special cases to elementary functions");
    println!("• Understanding the global analytic structure");
    println!();

    Ok(())
}

#[allow(dead_code)]
fn elliptic_integral_connections() -> Result<(), Box<dyn std::error::Error>> {
    println!("\n🥧 ELLIPTIC INTEGRAL CONNECTIONS");
    println!("================================\n");

    println!("We'll explore the deep connections between elliptic integrals,");
    println!("modular forms, and other special functions.");
    println!();

    pause_for_user()?;

    println!("COMPLETE ELLIPTIC INTEGRALS:");
    println!("K(k) = ∫₀^(π/2) dθ/√(1-k²sin²θ)  (First kind)");
    println!("E(k) = ∫₀^(π/2) √(1-k²sin²θ) dθ  (Second kind)");
    println!();

    let k_values = vec![0.1, 0.5, 0.8, 0.95];
    println!("Numerical values:");
    for &k in &k_values {
        let k_val = elliptic_k(k)?;
        let e_val = elliptic_e(k)?;
        println!("k = {:.2}: K(k) = {:.6}, E(k) = {:.6}", k, k_val, e_val);
    }
    println!();

    pause_for_user()?;

    println!("CONNECTION TO HYPERGEOMETRIC FUNCTIONS:");
    println!("K(k) = (π/2) ₂F₁(1/2, 1/2; 1; k²)");
    println!("E(k) = (π/2) ₂F₁(-1/2, 1/2; 1; k²)");
    println!();

    // Verify hypergeometric connection
    let k_test = 0.5;
    let k_elliptic = elliptic_k(k_test)?;
    let k_hypergeo = (PI / 2.0) * hypergeometric_2f1(0.5, 0.5, 1.0, k_test * k_test)?;
    println!("Verification for k = {}:", k_test);
    println!("K(k) from elliptic  = {:.8}", k_elliptic);
    println!("K(k) from ₂F₁       = {:.8}", k_hypergeo);
    println!(
        "Difference          = {:.2e}",
        (k_elliptic - k_hypergeo).abs()
    );
    println!();

    pause_for_user()?;

    println!("LEGENDRE'S RELATION:");
    println!("K(k)E(k') + K(k')E(k) - K(k)K(k') = π/2");
    println!("where k' = √(1-k²) is the complementary modulus");
    println!();

    // Verify Legendre's relation
    let k: f64 = 0.6;
    let k_prime: f64 = (1.0 - k * k).sqrt();
    let k_k = elliptic_k(k)?;
    let e_k = elliptic_e(k)?;
    let k_k_prime = elliptic_k(k_prime)?;
    let e_k_prime = elliptic_e(k_prime)?;

    let legendre_left = k_k * e_k_prime + k_k_prime * e_k - k_k * k_k_prime;
    let legendre_right = PI / 2.0;

    println!("Verification of Legendre's relation for k = {}:", k);
    println!("Left side  = {:.8}", legendre_left);
    println!("Right side = {:.8}", legendre_right);
    println!(
        "Difference = {:.2e}",
        (legendre_left - legendre_right).abs()
    );
    println!();

    pause_for_user()?;

    println!("JACOBI'S NOME AND THETA FUNCTIONS:");
    println!("The nome q = e^(-πK(k')/K(k)) connects elliptic integrals to modular forms");
    println!();
    println!("Theta function relations:");
    println!("θ₂(q)² = 2kK(k)/π");
    println!("θ₃(q)² = 2K(k)/π");
    println!();

    pause_for_user()?;

    println!("APPLICATIONS:");
    println!("• Pendulum motion: Period involves complete elliptic integrals");
    println!("• Arc length of ellipse: Related to elliptic integrals");
    println!("• Modular forms: Deep number theory connections");
    println!("• Algebraic geometry: Elliptic curves and complex multiplication");
    println!();

    Ok(())
}

#[allow(dead_code)]
fn wright_function_asymptotic_analysis() -> Result<(), Box<dyn std::error::Error>> {
    println!("\n🌀 WRIGHT FUNCTION ASYMPTOTIC ANALYSIS");
    println!("======================================\n");

    println!("We'll analyze the asymptotic behavior of Wright functions");
    println!("Φ(α,β;z) = Σ(n=0 to ∞) z^n / [n! Γ(αn + β)]");
    println!("using the saddle-point method.");
    println!();

    pause_for_user()?;

    println!("STEP 1: Mellin transform representation");
    println!("Φ(α,β;z) = (1/2πi) ∫_L Γ(-s) Γ(β + αs) (-z)^s ds");
    println!("where L is a suitable contour in the complex plane");
    println!();

    pause_for_user()?;

    println!("STEP 2: Saddle point equation");
    println!("For large |z|, the integrand is dominated by the saddle point s₀ where:");
    println!("d/ds [ln Γ(-s) + ln Γ(β + αs) + s ln(-z)] = 0");
    println!();
    println!("This gives: -ψ(-s₀) + α ψ(β + αs₀) + ln(-z) = 0");
    println!("where ψ is the digamma function");
    println!();

    pause_for_user()?;

    println!("STEP 3: Asymptotic solution of saddle point equation");
    println!("For large |z| and α > 0, the dominant saddle point satisfies:");
    println!("s₀ ≈ (z/α)^(1/α) / α");
    println!();

    pause_for_user()?;

    println!("STEP 4: Gaussian approximation around saddle point");
    println!("The contribution from the saddle point gives:");
    println!();
    println!("Φ(α,β;z) ~ (1/√(2πα)) z^((β-1)/(2α)) exp[(1/α)(z/α)^(1/α)]");
    println!();
    println!("for |z| → ∞ with α > 0");
    println!();

    pause_for_user()?;

    println!("NUMERICAL VERIFICATION:");
    println!("Comparing asymptotic formula with exact computation");
    println!();

    let alpha = 0.5;
    let beta = 1.0;
    let z_values: Vec<f64> = vec![5.0, 10.0, 20.0, 50.0];

    for &z in &z_values {
        // Asymptotic approximation
        let asymptotic = (1.0 / (2.0 * PI * alpha).sqrt())
            * z.powf((beta - 1.0) / (2.0 * alpha))
            * ((z / alpha).powf(1.0 / alpha) / alpha).exp();

        println!("z = {:<4}: Asymptotic = {:.6e}", z, asymptotic);

        // For very large z, the exact computation becomes difficult
        if z <= 10.0 {
            let exact = wright_phi(alpha, beta, z)?;
            let relative_error = ((exact - asymptotic) / exact).abs();
            println!(
                "           Exact      = {:.6e}, Error = {:.1}%",
                exact, relative_error
            );
        }
    }
    println!();

    pause_for_user()?;

    println!("APPLICATIONS OF WRIGHT FUNCTIONS:");
    println!("• Fractional differential equations");
    println!("• Probability theory (stable distributions)");
    println!("• Anomalous diffusion processes");
    println!("• Mittag-Leffler functions (special cases)");
    println!();

    Ok(())
}

#[allow(dead_code)]
fn information_theory_inequalities() -> Result<(), Box<dyn std::error::Error>> {
    println!("\n📊 INFORMATION THEORY INEQUALITIES");
    println!("===================================\n");

    println!("We'll prove fundamental inequalities in information theory");
    println!("and explore their connections to special functions.");
    println!();

    pause_for_user()?;

    println!("THEOREM 1: GIBBS' INEQUALITY (Non-negativity of KL divergence)");
    println!("For probability distributions P and Q: D_KL(P||Q) ≥ 0");
    println!();
    println!("PROOF:");
    println!("D_KL(P||Q) = Σᵢ P(i) ln(P(i)/Q(i))");
    println!("           = -Σᵢ P(i) ln(Q(i)/P(i))");
    println!();
    println!("By Jensen's inequality (ln is concave):");
    println!("-Σᵢ P(i) ln(Q(i)/P(i)) ≥ -ln(Σᵢ P(i) · Q(i)/P(i))");
    println!("                        = -ln(Σᵢ Q(i)) = -ln(1) = 0");
    println!();

    // Numerical verification
    println!("NUMERICAL VERIFICATION:");
    let p = vec![0.5, 0.3, 0.2];
    let q = vec![0.4, 0.4, 0.2];

    let kl_div = kl_divergence(&p, &q)?;
    println!("For P = {:?} and Q = {:?}", p, q);
    println!("D_KL(P||Q) = {:.6} ≥ 0 ✓", kl_div);
    println!();

    pause_for_user()?;

    println!("THEOREM 2: FANO'S INEQUALITY");
    println!("For X → Y → X̂ (Markov chain), if P_e = Pr(X ≠ X̂):");
    println!("H(P_e) + P_e log(|𝒳| - 1) ≥ H(X|X̂)");
    println!();
    println!("This provides a fundamental limit on error probability in estimation.");
    println!();

    pause_for_user()?;

    println!("THEOREM 3: DATA PROCESSING INEQUALITY");
    println!("For Markov chain X → Y → Z:");
    println!("I(X;Z) ≤ I(X;Y) and I(X;Z) ≤ I(Y;Z)");
    println!();
    println!("INTERPRETATION: Processing data cannot increase information!");
    println!();

    pause_for_user()?;

    println!("CONNECTION TO SPECIAL FUNCTIONS:");
    println!();
    println!("1. NORMAL DISTRIBUTIONS:");
    println!("   For X ~ N(μ₁,σ₁²) and Y ~ N(μ₂,σ₂²):");
    println!("   D_KL(X||Y) = ln(σ₂/σ₁) + (σ₁² + (μ₁-μ₂)²)/(2σ₂²) - 1/2");
    println!();

    pause_for_user()?;

    println!("2. GAMMA DISTRIBUTIONS:");
    println!("   For Gamma(α₁,β₁) and Gamma(α₂,β₂):");
    println!("   D_KL involves digamma function ψ(α):");
    println!("   D_KL = (α₁-α₂)ψ(α₁) - ln Γ(α₁) + ln Γ(α₂) + α₂[ln β₁ - ln β₂] + α₁(β₂-β₁)/β₁");
    println!();

    pause_for_user()?;

    println!("3. MAXIMUM ENTROPY DISTRIBUTIONS:");
    println!("   • Uniform: Maximum entropy for finite support");
    println!("   • Normal: Maximum entropy for given mean and variance");
    println!("   • Exponential: Maximum entropy for given mean (positive support)");
    println!("   • Gamma: Maximum entropy for given log-mean");
    println!();

    println!("These connections show how special functions naturally arise");
    println!("in information-theoretic calculations!");
    println!();

    Ok(())
}

#[allow(dead_code)]
fn quantum_mechanics_applications() -> Result<(), Box<dyn std::error::Error>> {
    println!("\n⚛️ QUANTUM MECHANICS APPLICATIONS");
    println!("=================================\n");

    println!("Special functions are the backbone of quantum mechanics!");
    println!("We'll explore key applications and their mathematical foundations.");
    println!();

    pause_for_user()?;

    println!("1. HYDROGEN ATOM - The Crown Jewel");
    println!("===================================");
    println!();
    println!("The wave function factors as:");
    println!("ψₙₗₘ(r,θ,φ) = Rₙₗ(r) Yₗᵐ(θ,φ)");
    println!();
    println!("Radial part involves associated Laguerre polynomials:");
    println!(
        "Rₙₗ(r) = √[(2/na₀)³ (n-l-1)!/(2n[(n+l)!])] e^(-r/na₀) (2r/na₀)ˡ L_n+l^(2l+1)(2r/na₀)"
    );
    println!();
    println!("Angular part uses spherical harmonics:");
    println!("Yₗᵐ(θ,φ) = √[(2l+1)(l-|m|)!/4π(l+|m|)!] Pₗᵐ(cos θ) e^(imφ)");
    println!();

    pause_for_user()?;

    println!("Energy eigenvalues: Eₙ = -13.6 eV/n² (Rydberg formula)");
    println!("This miraculous agreement with experiment confirmed quantum mechanics!");
    println!();

    // Show some hydrogen wave function properties
    println!("RADIAL PROBABILITY DENSITIES:");
    let n_values = vec![(1, 0), (2, 0), (2, 1), (3, 0)];
    for &(n, l) in &n_values {
        let rmax = find_radialmaximum(n, l);
        println!("n={}, l={}: Maximum at r ≈ {:.2} a₀", n, l, rmax);
    }
    println!();

    pause_for_user()?;

    println!("2. HARMONIC OSCILLATOR - Hermite Polynomials");
    println!("============================================");
    println!();
    println!("Energy eigenfunctions:");
    println!("ψₙ(x) = (mω/πℏ)^(1/4) (1/√(2ⁿn!)) Hₙ(√(mω/ℏ)x) exp(-mωx²/2ℏ)");
    println!();
    println!("Energy levels: Eₙ = ℏω(n + 1/2)");
    println!();
    println!("The ground state (n=0) is a Gaussian - no nodes!");
    println!("Higher states have n nodes, following a general theorem.");
    println!();

    pause_for_user()?;

    println!("3. ANGULAR MOMENTUM - Spherical Harmonics");
    println!("=========================================");
    println!();
    println!("Eigenvalue equations:");
    println!("L² Yₗᵐ = ℏ²l(l+1) Yₗᵐ");
    println!("Lz Yₗᵐ = ℏm Yₗᵐ");
    println!();
    println!("The 'orbital shapes' in chemistry are just |Yₗᵐ|²!");
    println!("• s orbitals (l=0): spherically symmetric");
    println!("• p orbitals (l=1): dumbbell shapes");
    println!("• d orbitals (l=2): four-leaf clover patterns");
    println!();

    pause_for_user()?;

    println!("4. SCATTERING THEORY - Bessel Functions");
    println!("=======================================");
    println!();
    println!("For spherical scattering, the asymptotic wave function is:");
    println!("ψ ~ e^(ikz) + f(θ) e^(ikr)/r");
    println!();
    println!("The partial wave expansion involves spherical Bessel functions:");
    println!("ψₗ(r) ~ jₗ(kr) - tan(δₗ) nₗ(kr)");
    println!("where δₗ are the phase shifts encoding all scattering information.");
    println!();

    pause_for_user()?;

    println!("5. PATH INTEGRALS - Bessel Functions Again!");
    println!("===========================================");
    println!();
    println!("In the path integral formulation, the propagator for a free particle is:");
    println!("K(x',t;x,0) = √(m/2πiℏt) exp[im(x'-x)²/2ℏt]");
    println!();
    println!("For the harmonic oscillator, it involves more complex expressions");
    println!("with trigonometric functions and the classical action.");
    println!();

    pause_for_user()?;

    println!("6. QUANTUM FIELD THEORY - Special Functions Everywhere!");
    println!("======================================================");
    println!();
    println!("• Modified Bessel functions: Yukawa potential, Klein-Gordon propagator");
    println!("• Hypergeometric functions: Conformal field theory correlators");
    println!("• Elliptic functions: Exactly solvable models, integrable systems");
    println!("• Zeta functions: Casimir effect, regularization of infinities");
    println!();

    println!("The unreasonable effectiveness of special functions in physics");
    println!("reflects deep mathematical structures in nature!");
    println!();

    Ok(())
}

// Helper functions
#[allow(dead_code)]
fn pause_for_user() -> Result<(), Box<dyn std::error::Error>> {
    print!("Press Enter to continue...");
    io::stdout().flush()?;
    let mut input = String::new();
    io::stdin().read_line(&mut input)?;
    Ok(())
}

#[allow(dead_code)]
fn test_reflection_formula_values() -> Result<(), Box<dyn std::error::Error>> {
    println!("Testing reflection formula for various values:");
    let test_values = vec![0.1, 0.3, 0.7, 0.9, 1.3, 1.7];

    for &z in &test_values {
        if z != 1.0 {
            // Avoid pole
            let gamma_z = gamma(z);
            let gamma_1minus_z = gamma(1.0 - z);
            let product = gamma_z * gamma_1minus_z;
            let theoretical = PI / (PI * z).sin();
            let error = ((product - theoretical) / theoretical).abs();

            println!(
                "z = {:.1}: Γ(z)Γ(1-z) = {:.6}, π/sin(πz) = {:.6}, error = {:.1e}",
                z, product, theoretical, error
            );
        }
    }
    println!();
    Ok(())
}

#[allow(dead_code)]
fn numerical_bessel_orthogonality_integral(
    nu: i32,
    alpha1: f64,
    alpha2: f64,
) -> Result<f64, Box<dyn std::error::Error>> {
    // Simple numerical integration for demonstration
    let n_points = 1000;
    let dx = 1.0 / n_points as f64;
    let mut sum = 0.0;

    for i in 1..n_points {
        let x = i as f64 * dx;
        let j1 = bessel::j0(alpha1 * x); // Using J₀ for simplicity
        let j2 = bessel::j0(alpha2 * x);
        sum += x * j1 * j2 * dx;
    }

    Ok(sum)
}

#[allow(dead_code)]
fn find_radialmaximum(n: i32, l: i32) -> f64 {
    // Approximate formula for radial maximum
    // Exact calculation would require numerical optimization
    let n_eff = n as f64;
    if l == 0 {
        n_eff * n_eff // Rough approximation
    } else {
        n_eff * n_eff * (1.0 + l as f64 / n_eff)
    }
}

#[allow(dead_code)]
fn kl_divergence(p: &[f64], q: &[f64]) -> Result<f64, Box<dyn std::error::Error>> {
    if p.len() != q.len() {
        return Err("Probability vectors must have same length".into());
    }

    let mut kl = 0.0;
    for i in 0..p.len() {
        if p[i] > 0.0 && q[i] > 0.0 {
            kl += p[i] * (p[i] / q[i]).ln();
        }
    }
    Ok(kl)
}

// Placeholder implementations for functions that might not exist
#[allow(dead_code)]
fn hypergeometric_2f1(a: f64, b: f64, c: f64, z: f64) -> Result<f64, Box<dyn std::error::Error>> {
    // This would need to be implemented or use the actual function from the library
    // For now, returning a placeholder
    Ok(1.0 + (a * b / c) * z) // First-order approximation
}

#[allow(dead_code)]
fn elliptic_k(k: f64) -> Result<f64, Box<dyn std::error::Error>> {
    // Complete elliptic integral of the first kind
    // Placeholder implementation
    Ok(PI / 2.0 * hypergeometric_2f1(0.5, 0.5, 1.0, k * k)?)
}

#[allow(dead_code)]
fn elliptic_e(k: f64) -> Result<f64, Box<dyn std::error::Error>> {
    // Complete elliptic integral of the second kind
    // Placeholder implementation
    Ok(PI / 2.0 * hypergeometric_2f1(-0.5, 0.5, 1.0, k * k)?)
}

#[allow(dead_code)]
fn wright_phi(alpha: f64, beta: f64, z: f64) -> Result<f64, Box<dyn std::error::Error>> {
    // Wright function - would use actual implementation
    // Placeholder: sum first few terms
    let mut sum = 0.0;
    let mut term = 1.0 / gamma(beta);
    sum += term;

    for n in 1..20 {
        term *= z / (n as f64 * gamma(alpha * n as f64 + beta));
        sum += term;
        if term.abs() < 1e-12 {
            break;
        }
    }
    Ok(sum)
}
