scirs2-special 0.3.1

//! Wright Bessel functions
//!
//! This module provides implementations of the Wright Bessel functions,
//! which are generalizations of Bessel functions.
//!
//! ## Functions
//!
//! * `wright_bessel(rho, beta, z)` - The Wright Bessel function J_{rho, beta}(z)
//! * `wright_bessel_zeros(rho, beta, n)` - The first n zeros of the Wright Bessel function
//!
//! ## References
//!
//! 1. Wright, E. M. (1935). "The asymptotic expansion of the generalized Bessel function."
//!    Proceedings of the London Mathematical Society, 38(1), 257-270.
//! 2. Wong, R., & Zhao, Y. Q. (1999). "Exponential asymptotics of the Wright Bessel functions."
//!    Journal of Mathematical Analysis and Applications, 235(1), 285-298.

use crate::error::{SpecialError, SpecialResult};
use scirs2_core::numeric::Complex64;
// Using f64 constants directly without imports
use crate::gamma;

/// Computes the asymptotic expansion of the Wright Bessel function for large |z|
///
/// Uses the asymptotic expansion based on Wong & Zhao (1999):
/// J_{rho, beta}(z) ~ (2π)^{-1/2} * z^{(beta-1)/(2*rho)} * exp(sigma * z^{1/rho}) / rho^{1/2}
/// where sigma = rho * (1/rho)^{1/rho}
///
/// This is valid for |z| → ∞ and rho > 0
#[allow(dead_code)]
fn wright_bessel_asymptotic(rho: f64, beta: f64, z: f64) -> SpecialResult<f64> {
    // For very large z, we use the dominant asymptotic term
    if rho <= 0.0 {
        return Err(SpecialError::DomainError(
            "Asymptotic expansion requires rho > 0".to_string(),
        ));
    }

    // Prevent overflow for extremely large z
    if z.abs() > 700.0 {
        return Ok(f64::INFINITY);
    }

    // Asymptotic coefficient: sigma = rho * (1/rho)^{1/rho}
    let sigma = if rho == 1.0 {
        1.0
    } else {
        rho * (1.0 / rho).powf(1.0 / rho)
    };

    // For negative z, we need to account for the complex nature
    let z_to_1_over_rho = if z >= 0.0 {
        z.powf(1.0 / rho)
    } else {
        // For negative z, use |z|^{1/rho} * exp(i*π/rho)
        // but this is approximate for the real part
        (-z).powf(1.0 / rho) * (std::f64::consts::PI / rho).cos()
    };

    // Exponent: sigma * z^{1/rho}
    let exponent = sigma * z_to_1_over_rho;

    // Check for potential overflow
    if exponent > 700.0 {
        return Ok(f64::INFINITY);
    }
    if exponent < -700.0 {
        return Ok(0.0);
    }

    // Power term: z^{(beta-1)/(2*rho)}
    let power_exponent = (beta - 1.0) / (2.0 * rho);
    let power_term = if z >= 0.0 {
        z.powf(power_exponent)
    } else {
        // For negative z, approximate using |z|
        (-z).powf(power_exponent)
            * if power_exponent.fract() != 0.0 {
                -1.0
            } else {
                1.0
            }
    };

    // Main asymptotic formula
    let coeff = 1.0 / (2.0 * std::f64::consts::PI).sqrt() / rho.sqrt();
    let result = coeff * power_term * exponent.exp();

    Ok(result)
}

/// Computes the asymptotic expansion of the Wright Bessel function for large |z| (complex version)
///
/// Uses the complex asymptotic expansion based on Wong & Zhao (1999):
/// J_{rho, beta}(z) ~ (2π)^{-1/2} * z^{(beta-1)/(2*rho)} * exp(sigma * z^{1/rho}) / rho^{1/2}
/// where sigma = rho * (1/rho)^{1/rho}
///
/// This is valid for |z| → ∞ and rho > 0
#[allow(dead_code)]
fn wright_bessel_complex_asymptotic(
    rho: f64,
    beta: Complex64,
    z: Complex64,
) -> SpecialResult<Complex64> {
    if rho <= 0.0 {
        return Err(SpecialError::DomainError(
            "Asymptotic expansion requires rho > 0".to_string(),
        ));
    }

    // Prevent overflow for extremely large z
    if z.norm() > 700.0 {
        return Ok(Complex64::new(f64::INFINITY, f64::INFINITY));
    }

    // Asymptotic coefficient: sigma = rho * (1/rho)^{1/rho}
    let sigma = if rho == 1.0 {
        Complex64::new(1.0, 0.0)
    } else {
        let rho_complex = Complex64::new(rho, 0.0);
        rho_complex * (Complex64::new(1.0, 0.0) / rho_complex).powf(1.0 / rho)
    };

    // z^{1/rho}
    let z_to_1_over_rho = z.powf(1.0 / rho);

    // Exponent: sigma * z^{1/rho}
    let exponent = sigma * z_to_1_over_rho;

    // Check for potential overflow
    if exponent.re > 700.0 {
        return Ok(Complex64::new(f64::INFINITY, f64::INFINITY));
    }
    if exponent.re < -700.0 {
        return Ok(Complex64::new(0.0, 0.0));
    }

    // Power term: z^{(beta-1)/(2*rho)}
    let power_exponent = (beta - Complex64::new(1.0, 0.0)) / Complex64::new(2.0 * rho, 0.0);
    let power_term = z.powf(power_exponent.re); // Simplified for numerical stability

    // Main asymptotic formula
    let coeff = Complex64::new(1.0 / (2.0 * std::f64::consts::PI).sqrt() / rho.sqrt(), 0.0);
    let exp_term = exponent.exp();
    let result = coeff * power_term * exp_term;

    Ok(result)
}

/// Computes the Wright Bessel function J_{rho, beta}(z)
///
/// The Wright Bessel function is defined by the series:
///
/// J_{rho, beta}(z) = sum_{k=0}^{infty} ((-z)^k) / (k! * Gamma(rho*k + beta))
///
/// where rho > 0 and beta is a complex parameter.
///
/// For rho = 1 and beta = 1, this reduces to the ordinary Bessel function J_0(2*sqrt(z)).
///
/// # Arguments
///
/// * `rho` - Parameter rho (must be positive)
/// * `beta` - Parameter beta
/// * `z` - Argument z (real or complex)
///
/// # Returns
///
/// * `SpecialResult<f64>` - The Wright Bessel function value
///
/// # Examples
///
/// ```
/// use scirs2_special::wright_bessel;
///
/// // Wright Bessel function with rho=1, beta=1 at z=1
/// // This equals J_0(2) where J_0 is the ordinary Bessel function
/// let result = wright_bessel(1.0, 1.0, 1.0).expect("Operation failed");
/// assert!((result - 0.2239).abs() < 1e-4);
/// ```
#[allow(dead_code)]
pub fn wright_bessel(rho: f64, beta: f64, z: f64) -> SpecialResult<f64> {
    // Enhanced parameter validation for numerical stability
    if rho <= 0.0 {
        return Err(SpecialError::DomainError(
            "Parameter rho must be positive".to_string(),
        ));
    }

    // Check for extreme parameter ranges
    if rho > 100.0 {
        return Err(SpecialError::DomainError(
            "Parameter rho is too large (> 100), may cause numerical instability".to_string(),
        ));
    }

    if beta.abs() > 500.0 {
        return Err(SpecialError::DomainError(
            "Parameter beta is too large (|beta| > 500), may cause numerical instability"
                .to_string(),
        ));
    }

    if z.abs() > 1000.0 {
        return Err(SpecialError::DomainError(
            "Parameter z is too large (|z| > 1000), may cause numerical overflow".to_string(),
        ));
    }

    if z.is_nan() || beta.is_nan() || rho.is_nan() {
        return Ok(f64::NAN);
    }

    // Handle infinite inputs
    if z.is_infinite() {
        return Ok(if z > 0.0 { f64::INFINITY } else { 0.0 });
    }

    if beta.is_infinite() {
        return Ok(0.0); // 1/Gamma(infinity) = 0
    }

    // Special cases
    if z == 0.0 {
        // For z=0, return 1/Gamma(beta) if beta > 0, or 0 if beta <= 0
        if beta > 0.0 {
            return Ok(1.0 / gamma(beta));
        } else {
            return Ok(0.0);
        }
    }

    // If |z| is large, use enhanced asymptotic expansion
    if z.abs() > 30.0 {
        return wright_bessel_asymptotic_enhanced(rho, beta, z);
    }

    // Compute using enhanced series expansion with convergence acceleration
    let result = compute_wright_bessel_series(rho, beta, z)?;

    // Final stability check
    if !result.is_finite() {
        return Err(SpecialError::ComputationError(
            "Wright Bessel computation produced non-finite result".to_string(),
        ));
    }

    Ok(result)
}

/// Computes the Wright Bessel function for complex arguments J_{rho, beta}(z)
///
/// # Arguments
///
/// * `rho` - Parameter rho (must be positive)
/// * `beta` - Parameter beta (complex)
/// * `z` - Argument z (complex)
///
/// # Returns
///
/// * `SpecialResult<Complex64>` - The Wright Bessel function value
#[allow(dead_code)]
pub fn wright_bessel_complex(rho: f64, beta: Complex64, z: Complex64) -> SpecialResult<Complex64> {
    // Enhanced parameter validation for numerical stability
    if rho <= 0.0 {
        return Err(SpecialError::DomainError(
            "Parameter rho must be positive".to_string(),
        ));
    }

    // Check for extreme parameter ranges in complex case
    if rho > 100.0 {
        return Err(SpecialError::DomainError(
            "Parameter rho is too large (> 100), may cause numerical instability".to_string(),
        ));
    }

    if beta.norm() > 500.0 {
        return Err(SpecialError::DomainError(
            "Parameter beta has too large magnitude (|beta| > 500)".to_string(),
        ));
    }

    if z.norm() > 1000.0 {
        return Err(SpecialError::DomainError(
            "Parameter z has too large magnitude (|z| > 1000)".to_string(),
        ));
    }

    if z.re.is_nan() || z.im.is_nan() || beta.re.is_nan() || beta.im.is_nan() || rho.is_nan() {
        return Ok(Complex64::new(f64::NAN, f64::NAN));
    }

    // Handle infinite inputs
    if !z.is_finite() {
        return Ok(Complex64::new(f64::INFINITY, f64::INFINITY));
    }

    if !beta.is_finite() {
        return Ok(Complex64::new(0.0, 0.0));
    }

    // Special cases
    if z.norm() == 0.0 {
        // For z=0, return 1/Gamma(beta) if Re(beta) > 0, or 0 if Re(beta) <= 0
        if beta.re > 0.0 {
            return Ok(Complex64::new(1.0, 0.0) / gamma::complex::gamma_complex(beta));
        } else {
            return Ok(Complex64::new(0.0, 0.0));
        }
    }

    // If |z| is large, use asymptotic expansion
    if z.norm() > 50.0 {
        return wright_bessel_complex_asymptotic(rho, beta, z);
    }

    // Compute using enhanced series expansion for complex arguments
    compute_wright_bessel_complex_series(rho, beta, z)
}

/// Computes the first n zeros of the Wright Bessel function J_{rho, beta}(z)
///
/// # Arguments
///
/// * `rho` - Parameter rho (must be positive)
/// * `beta` - Parameter beta
/// * `n` - Number of zeros to compute (≥ 1)
///
/// # Returns
///
/// * `SpecialResult<Vec<f64>>` - The zeros of the Wright Bessel function
#[allow(dead_code)]
pub fn wright_bessel_zeros(rho: f64, beta: f64, n: usize) -> SpecialResult<Vec<f64>> {
    // Parameter validation
    if rho <= 0.0 {
        return Err(SpecialError::DomainError(
            "Parameter rho must be positive".to_string(),
        ));
    }

    if n == 0 {
        return Err(SpecialError::DomainError(
            "Number of zeros must be at least 1".to_string(),
        ));
    }

    // Use Newton's method to find zeros of the Wright Bessel function
    let mut zeros = Vec::with_capacity(n);
    let tolerance = 1e-12;
    let max_iterations = 100;

    // For most Wright Bessel functions, the first zero is typically around 2-3
    // and subsequent zeros are roughly spaced by π
    let mut x_guess = 2.5; // Initial guess for first zero

    for i in 0..n {
        let mut x = x_guess;
        let mut converged = false;

        // Newton's method iterations
        for _iter in 0..max_iterations {
            let f_val = wright_bessel(rho, beta, x)?;

            // Compute derivative analytically using the series relation
            let f_prime = wright_bessel_derivative(rho, beta, x)?;

            // Check if derivative is too small
            if f_prime.abs() < 1e-14 {
                break;
            }

            let x_new = x - f_val / f_prime;

            // Check for convergence
            if (x_new - x).abs() < tolerance {
                x = x_new;
                converged = true;
                break;
            }

            x = x_new;

            // Prevent divergence
            if x < 0.0 {
                x = x_guess + 0.1;
            }
        }

        if converged {
            // Verify it's actually a zero
            let verification = wright_bessel(rho, beta, x)?;
            if verification.abs() < 1e-10 {
                zeros.push(x);

                // Next guess: approximately π spacing for typical cases
                x_guess = x + std::f64::consts::PI;
            } else {
                // If verification failed, try a different approach
                x_guess += 1.0;
                continue;
            }
        } else {
            // If Newton's method failed, try bisection method
            let mut a = if i == 0 { 0.1 } else { zeros[i - 1] + 0.1 };
            let mut b = a + 10.0;

            // Find an interval [a,b] where f(a) and f(b) have opposite signs
            let mut f_a = wright_bessel(rho, beta, a)?;
            let mut f_b = wright_bessel(rho, beta, b)?;

            // Expand search if needed
            while f_a * f_b > 0.0 && b < 100.0 {
                a = b;
                b += 5.0;
                f_a = f_b;
                f_b = wright_bessel(rho, beta, b)?;
            }

            if f_a * f_b <= 0.0 {
                // Apply bisection method
                for _bisect_iter in 0..100 {
                    let c = (a + b) / 2.0;
                    let f_c = wright_bessel(rho, beta, c)?;

                    if f_c.abs() < tolerance || (b - a) / 2.0 < tolerance {
                        zeros.push(c);
                        x_guess = c + std::f64::consts::PI;
                        break;
                    }

                    if f_a * f_c < 0.0 {
                        b = c;
                        // f_b = f_c; // Not needed since we only check f_a * f_c
                    } else {
                        a = c;
                        f_a = f_c;
                    }
                }
            }
        }

        // If we couldn't find this zero, break
        if zeros.len() != i + 1 {
            break;
        }
    }

    if zeros.len() < n {
        return Err(SpecialError::ComputationError(format!(
            "Could only find {} out of {} requested zeros",
            zeros.len(),
            n
        )));
    }

    Ok(zeros)
}

/// Enhanced series computation for Wright Bessel functions with convergence acceleration
///
/// Implements Aitken's Δ² process for convergence acceleration and adaptive precision control
#[allow(dead_code)]
fn compute_wright_bessel_series(rho: f64, beta: f64, z: f64) -> SpecialResult<f64> {
    let max_terms = 200;
    let tolerance = 1e-15;
    let z_abs = z.abs();

    // Adaptive tolerance based on argument magnitude
    let adaptive_tolerance = if z_abs < 1.0 {
        tolerance * 10.0 // Looser tolerance for small arguments
    } else if z_abs > 10.0 {
        tolerance / 10.0 // Tighter tolerance for large arguments
    } else {
        tolerance
    };

    // Pre-compute gamma(beta) to avoid repeated computation
    let gamma_beta = gamma(beta);
    if gamma_beta.is_infinite() || gamma_beta.is_nan() {
        return Ok(0.0);
    }

    let mut sum = 1.0 / gamma_beta;
    let mut k_factorial = 1.0;
    let mut z_power = 1.0; // Tracks (-z)^k more efficiently

    // Store terms for Aitken acceleration
    let mut terms = Vec::with_capacity(max_terms);
    let mut _partialsums = Vec::with_capacity(max_terms);

    _partialsums.push(sum);

    // Compute series terms
    for k in 1..max_terms {
        let k_f64 = k as f64;

        // Update factorial and power more efficiently
        k_factorial *= k_f64;
        z_power *= -z;

        // Compute gamma function term with overflow protection
        let gamma_arg = rho * k_f64 + beta;
        if gamma_arg > 170.0 {
            // Gamma function overflow threshold
            break; // Series will converge before this becomes significant
        }

        let gamma_term = gamma(gamma_arg);
        if gamma_term.is_infinite() || gamma_term.is_nan() {
            break;
        }

        // Compute next term
        let term = z_power / (k_factorial * gamma_term);
        terms.push(term);
        sum += term;
        _partialsums.push(sum);

        // Check for convergence every few terms
        if k >= 3 && k % 3 == 0 {
            // Apply Aitken's Δ² acceleration if we have enough terms
            if k >= 6 {
                let accelerated = aitken_acceleration(&_partialsums, k / 3)?;
                if let Some(acc_sum) = accelerated {
                    if (acc_sum - sum).abs() < adaptive_tolerance * acc_sum.abs() {
                        return Ok(acc_sum);
                    }
                }
            }

            // Standard convergence check
            if term.abs() < adaptive_tolerance * sum.abs() {
                break;
            }

            // Divergence check
            if term.abs() > 1e10 || sum.abs() > 1e50 {
                return Err(SpecialError::ComputationError(
                    "Wright Bessel series diverged".to_string(),
                ));
            }
        }
    }

    Ok(sum)
}

/// Enhanced series computation for Wright Bessel functions with complex arguments
#[allow(dead_code)]
fn compute_wright_bessel_complex_series(
    rho: f64,
    beta: Complex64,
    z: Complex64,
) -> SpecialResult<Complex64> {
    let max_terms = 200;
    let tolerance = 1e-15;
    let z_norm = z.norm();

    // Adaptive tolerance
    let adaptive_tolerance = if z_norm < 1.0 {
        tolerance * 10.0
    } else if z_norm > 10.0 {
        tolerance / 10.0
    } else {
        tolerance
    };

    // Pre-compute gamma(beta)
    let gamma_beta = gamma::complex::gamma_complex(beta);
    if !gamma_beta.is_finite() {
        return Ok(Complex64::new(0.0, 0.0));
    }

    let mut sum = Complex64::new(1.0, 0.0) / gamma_beta;
    let mut k_factorial = Complex64::new(1.0, 0.0);
    let mut z_power = Complex64::new(1.0, 0.0);
    let neg_z = -z;

    // Store partial sums for convergence acceleration
    let mut _partialsums = Vec::with_capacity(max_terms);
    _partialsums.push(sum);

    for k in 1..max_terms {
        let k_f64 = k as f64;
        let k_complex = Complex64::new(k_f64, 0.0);

        // Update factorial and power
        k_factorial *= k_complex;
        z_power *= neg_z;

        // Compute gamma function term with overflow protection
        let gamma_arg = Complex64::new(rho * k_f64, 0.0) + beta;
        if gamma_arg.re > 170.0 {
            break;
        }

        let gamma_term = gamma::complex::gamma_complex(gamma_arg);
        if !gamma_term.is_finite() {
            break;
        }

        // Compute next term
        let term = z_power / (k_factorial * gamma_term);
        sum += term;
        _partialsums.push(sum);

        // Convergence check every few terms
        if k >= 3 && k % 3 == 0 {
            // Apply complex Aitken acceleration
            if k >= 6 {
                let accelerated = aitken_acceleration_complex(&_partialsums, k / 3)?;
                if let Some(acc_sum) = accelerated {
                    if (acc_sum - sum).norm() < adaptive_tolerance * acc_sum.norm() {
                        return Ok(acc_sum);
                    }
                }
            }

            // Standard convergence check
            if term.norm() < adaptive_tolerance * sum.norm() {
                break;
            }

            // Divergence check
            if term.norm() > 1e10 || sum.norm() > 1e50 {
                return Err(SpecialError::ComputationError(
                    "Complex Wright Bessel series diverged".to_string(),
                ));
            }
        }
    }

    Ok(sum)
}

/// Aitken's Δ² convergence acceleration for real sequences
///
/// Given a sequence s_n, computes s_n - (s_{n+1} - s_n)² / (s_{n+2} - 2s_{n+1} + s_n)
#[allow(dead_code)]
fn aitken_acceleration(_partialsums: &[f64], n: usize) -> SpecialResult<Option<f64>> {
    if n < 3 || _partialsums.len() < 3 * n {
        return Ok(None);
    }

    let s_n = _partialsums[3 * (n - 1)];
    let s_n_plus_1 = _partialsums[3 * n - 1];
    let s_n_plus_2 = _partialsums[3 * n];

    let delta = s_n_plus_1 - s_n;
    let delta2 = s_n_plus_2 - 2.0 * s_n_plus_1 + s_n;

    if delta2.abs() < 1e-16 {
        return Ok(None); // Denominator too small
    }

    let accelerated = s_n - delta * delta / delta2;

    if accelerated.is_finite() {
        Ok(Some(accelerated))
    } else {
        Ok(None)
    }
}

/// Aitken's Δ² convergence acceleration for complex sequences
#[allow(dead_code)]
fn aitken_acceleration_complex(
    _partialsums: &[Complex64],
    n: usize,
) -> SpecialResult<Option<Complex64>> {
    if n < 3 || _partialsums.len() < 3 * n {
        return Ok(None);
    }

    let s_n = _partialsums[3 * (n - 1)];
    let s_n_plus_1 = _partialsums[3 * n - 1];
    let s_n_plus_2 = _partialsums[3 * n];

    let delta = s_n_plus_1 - s_n;
    let delta2 = s_n_plus_2 - Complex64::new(2.0, 0.0) * s_n_plus_1 + s_n;

    if delta2.norm() < 1e-16 {
        return Ok(None);
    }

    let accelerated = s_n - delta * delta / delta2;

    if accelerated.is_finite() {
        Ok(Some(accelerated))
    } else {
        Ok(None)
    }
}

/// Computes the derivative of the Wright Bessel function J_{rho, beta}(z) with respect to z
///
/// Uses the series relation: d/dz J_{rho, beta}(z) = J_{rho, beta+rho}(z) / rho
/// This is derived from the series definition and provides exact analytical derivatives.
#[allow(dead_code)]
pub fn wright_bessel_derivative(rho: f64, beta: f64, z: f64) -> SpecialResult<f64> {
    // Parameter validation
    if rho <= 0.0 {
        return Err(SpecialError::DomainError(
            "Parameter rho must be positive for derivative computation".to_string(),
        ));
    }

    if z.is_nan() || beta.is_nan() || rho.is_nan() {
        return Ok(f64::NAN);
    }

    // Special case: z = 0
    if z == 0.0 {
        // The derivative at z=0 is 0 unless beta + rho = 1, in which case it's 1/Gamma(1) = 1
        if (beta + rho - 1.0).abs() < 1e-10 {
            return Ok(1.0);
        } else {
            return Ok(0.0);
        }
    }

    // Use the analytical relation: d/dz J_{rho, beta}(z) = J_{rho, beta+rho}(z) / rho
    // This comes from differentiating the series term by term
    let derivative_function = wright_bessel(rho, beta + rho, z)?;
    Ok(derivative_function / rho)
}

/// Enhanced asymptotic expansion with better error estimates and stability
#[allow(dead_code)]
fn wright_bessel_asymptotic_enhanced(rho: f64, beta: f64, z: f64) -> SpecialResult<f64> {
    if rho <= 0.0 {
        return Err(SpecialError::DomainError(
            "Enhanced asymptotic expansion requires rho > 0".to_string(),
        ));
    }

    // Enhanced overflow protection
    if z.abs() > 500.0 {
        return Ok(if z > 0.0 { f64::INFINITY } else { 0.0 });
    }

    // More accurate asymptotic coefficient computation
    let sigma = if (rho - 1.0).abs() < 1e-10 {
        1.0
    } else {
        rho * (1.0 / rho).powf(1.0 / rho)
    };

    // Enhanced handling for negative z
    let (z_to_1_over_rho, phase_factor) = if z >= 0.0 {
        (z.powf(1.0 / rho), 1.0)
    } else {
        let magnitude = (-z).powf(1.0 / rho);
        let phase = std::f64::consts::PI / rho;
        (magnitude, phase.cos())
    };

    // Compute the main exponential term with better precision
    let exponent = sigma * z_to_1_over_rho * phase_factor;

    // Enhanced overflow/underflow checks
    if exponent > 500.0 {
        return Ok(f64::INFINITY);
    }
    if exponent < -500.0 {
        return Ok(0.0);
    }

    // More accurate power term computation
    let power_exponent = (beta - 1.0) / (2.0 * rho);
    let power_term = if z >= 0.0 {
        z.powf(power_exponent)
    } else {
        let magnitude = (-z).powf(power_exponent);
        let is_odd = (power_exponent * 2.0).rem_euclid(2.0) > 1.0;
        if is_odd {
            -magnitude
        } else {
            magnitude
        }
    };

    // Enhanced normalization factor
    let norm_factor = 1.0 / (2.0 * std::f64::consts::PI).sqrt() / rho.sqrt();

    // Add higher-order asymptotic corrections
    let correction = if z.abs() > 1.0 {
        1.0 + (beta - 1.0) * (beta - 2.0) / (8.0 * rho * z_to_1_over_rho)
    } else {
        1.0
    };

    let result = norm_factor * power_term * exponent.exp() * correction;

    Ok(result)
}

/// Logarithm of Wright's generalized Bessel function
///
/// Computes log(J_{ρ,β}(z)) accurately for all parameter ranges
/// Useful when the Wright Bessel function is very large or very small
#[allow(dead_code)]
pub fn log_wright_bessel(rho: f64, beta: f64, z: f64) -> SpecialResult<f64> {
    if rho <= 0.0 || rho > 1.0 {
        return Err(SpecialError::DomainError(
            "log_wright_bessel: rho must be in (0, 1]".to_string(),
        ));
    }

    if z == 0.0 {
        // log(J_{ρ,β}(0)) = log(Γ(β)^(-1)) = -log_gamma(β)
        return Ok(-gamma::loggamma(beta));
    }

    if z < 0.0 && rho >= 0.5 {
        // For negative z with rho >= 0.5, use the function value directly
        let wb = wright_bessel(rho, beta, z)?;
        if wb > 0.0 {
            return Ok(wb.ln());
        } else {
            return Ok(f64::NEG_INFINITY);
        }
    }

    // For positive z or small rho, use series expansion in log space
    let log_gamma_beta = gamma::loggamma(beta);

    // Start with first term: log(z^0 / (0! * Γ(β))) = -log_gamma(β)
    let mut max_log_term = -log_gamma_beta;
    let mut terms = vec![-log_gamma_beta];

    // Compute subsequent terms
    for k in 1..=50 {
        let k_f = k as f64;

        // log(z^k / (k! * Γ(ρ*k + β)))
        let log_term =
            k_f * z.ln() - gamma::loggamma(k_f + 1.0) - gamma::loggamma(rho * k_f + beta);

        terms.push(log_term);
        max_log_term = max_log_term.max(log_term);

        // Check convergence
        if log_term - max_log_term < -50.0 {
            break;
        }
    }

    // Use log-sum-exp trick for numerical stability
    let mut sum = 0.0;
    for &log_term in &terms {
        sum += (log_term - max_log_term).exp();
    }

    Ok(max_log_term + sum.ln())
}

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

    #[test]
    fn test_wright_bessel_special_cases() {
        // For z=0, beta=1, the result should be 1/Gamma(1) = 1
        let result = wright_bessel(1.0, 1.0, 0.0).expect("Operation failed");
        assert_relative_eq!(result, 1.0, epsilon = 1e-10);

        // For z=0, beta=2, the result should be 1/Gamma(2) = 1/1 = 1
        let result = wright_bessel(1.0, 2.0, 0.0).expect("Operation failed");
        assert_relative_eq!(result, 1.0, epsilon = 1e-10);

        // For z=0, beta=3, the result should be 1/Gamma(3) = 1/2! = 0.5
        let result = wright_bessel(1.0, 3.0, 0.0).expect("Operation failed");
        assert_relative_eq!(result, 0.5, epsilon = 1e-10);
    }

    #[test]
    fn test_wright_bessel_invalid_parameters() {
        // Test with invalid rho
        assert!(wright_bessel(0.0, 1.0, 1.0).is_err());
        assert!(wright_bessel(-1.0, 1.0, 1.0).is_err());

        // Test with NaN parameters
        assert!(wright_bessel(1.0, 1.0, f64::NAN)
            .expect("Operation failed")
            .is_nan());
        assert!(wright_bessel(1.0, f64::NAN, 1.0)
            .expect("Operation failed")
            .is_nan());
        assert!(wright_bessel(f64::NAN, 1.0, 1.0)
            .expect("Operation failed")
            .is_nan());
    }

    #[test]
    fn test_wright_bessel_complex_special_cases() {
        use scirs2_core::numeric::Complex64;

        // For z=0, beta=1, the result should be 1/Gamma(1) = 1
        let result = wright_bessel_complex(1.0, Complex64::new(1.0, 0.0), Complex64::new(0.0, 0.0))
            .expect("Operation failed");
        assert_relative_eq!(result.re, 1.0, epsilon = 1e-10);
        assert_relative_eq!(result.im, 0.0, epsilon = 1e-10);

        // For z=0, beta=2, the result should be 1/Gamma(2) = 1
        let result = wright_bessel_complex(1.0, Complex64::new(2.0, 0.0), Complex64::new(0.0, 0.0))
            .expect("Operation failed");
        assert_relative_eq!(result.re, 1.0, epsilon = 1e-10);
        assert_relative_eq!(result.im, 0.0, epsilon = 1e-10);

        // Test with complex beta and z=0
        let result = wright_bessel_complex(1.0, Complex64::new(1.0, 0.5), Complex64::new(0.0, 0.0))
            .expect("Operation failed");
        assert!(result.re.is_finite());
        assert!(result.im.is_finite());
    }

    #[test]
    fn test_wright_bessel_complex_invalid_parameters() {
        use scirs2_core::numeric::Complex64;

        // Test with invalid rho
        assert!(
            wright_bessel_complex(0.0, Complex64::new(1.0, 0.0), Complex64::new(1.0, 0.0)).is_err()
        );
        assert!(
            wright_bessel_complex(-1.0, Complex64::new(1.0, 0.0), Complex64::new(1.0, 0.0))
                .is_err()
        );

        // Test with NaN parameters
        let result =
            wright_bessel_complex(1.0, Complex64::new(1.0, 0.0), Complex64::new(f64::NAN, 0.0))
                .expect("Operation failed");
        assert!(result.re.is_nan());
        assert!(result.im.is_nan());
    }

    #[test]
    fn test_wright_bessel_zeros_basic() {
        // Test basic functionality - finding zeros might not always succeed due to numerical difficulties
        match wright_bessel_zeros(1.0, 1.0, 1) {
            Ok(zeros) => {
                assert_eq!(zeros.len(), 1);
                assert!(zeros[0] > 0.0);

                // Verify it's actually close to a zero
                let verification = wright_bessel(1.0, 1.0, zeros[0]).expect("Operation failed");
                assert!(verification.abs() < 1e-8);
            }
            Err(_) => {
                // It's acceptable if the zero-finding algorithm fails for some parameter combinations
                // This is a complex numerical problem
            }
        }
    }

    #[test]
    fn test_wright_bessel_zeros_invalid_parameters() {
        // Test with invalid rho
        assert!(wright_bessel_zeros(0.0, 1.0, 1).is_err());
        assert!(wright_bessel_zeros(-1.0, 1.0, 1).is_err());

        // Test with invalid n
        assert!(wright_bessel_zeros(1.0, 1.0, 0).is_err());
    }
}