scirs2_linalg/eigen/

standard.rs

//! Standard eigenvalue decomposition for dense matrices
//!
//! This module provides functions for computing eigenvalues and eigenvectors of
//! dense matrices using various algorithms:
//! - General eigenvalue decomposition for non-symmetric matrices
//! - Symmetric/Hermitian eigenvalue decomposition for better performance
//! - Power iteration for dominant eigenvalues
//! - QR algorithm for general cases

use scirs2_core::ndarray::{Array1, Array2, ArrayView2, ScalarOperand};
use scirs2_core::numeric::Complex;
use scirs2_core::numeric::{Float, NumAssign};
use scirs2_core::random::prelude::*;
use std::iter::Sum;

use crate::decomposition::qr;
use crate::error::{LinalgError, LinalgResult};
use crate::norm::vector_norm;
use crate::validation::validate_decomposition;

/// Type alias for eigenvalue-eigenvector pair result
/// Returns a tuple of (eigenvalues, eigenvectors) where eigenvalues is a 1D array
/// and eigenvectors is a 2D array where each column corresponds to an eigenvector
pub type EigenResult<F> = LinalgResult<(Array1<Complex<F>>, Array2<Complex<F>>)>;

/// Compute the eigenvalues and right eigenvectors of a square matrix.
///
/// # Arguments
///
/// * `a` - Input square matrix
/// * `workers` - Number of worker threads (None = use default)
///
/// # Returns
///
/// * Tuple (eigenvalues, eigenvectors) where eigenvalues is a complex vector
///   and eigenvectors is a complex matrix whose columns are the eigenvectors
///
/// # Examples
///
/// ```
/// use scirs2_core::ndarray::array;
/// use scirs2_linalg::eigen::standard::eig;
///
/// let a = array![[1.0_f64, 0.0], [0.0, 2.0]];
/// let (w, v) = eig(&a.view(), None).expect("Operation failed");
///
/// // Sort eigenvalues (they may be returned in different order)
/// let mut eigenvalues = vec![(w[0].re, 0), (w[1].re, 1)];
/// eigenvalues.sort_by(|a, b| a.0.partial_cmp(&b.0).expect("Operation failed"));
///
/// assert!((eigenvalues[0].0 - 1.0).abs() < 1e-10);
/// assert!((eigenvalues[1].0 - 2.0).abs() < 1e-10);
/// ```
#[allow(dead_code)]
pub fn eig<F>(a: &ArrayView2<F>, workers: Option<usize>) -> EigenResult<F>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    use crate::parallel;

    // Configure workers for parallel operations
    parallel::configure_workers(workers);

    // Parameter validation using validation helpers
    validate_decomposition(a, "Eigenvalue computation", true)?;

    let n = a.nrows();

    // For 1x1 and 2x2 matrices, we can compute eigenvalues analytically
    if n == 1 {
        let eigenvalue = Complex::new(a[[0, 0]], F::zero());
        let eigenvector = Array2::eye(1).mapv(|x| Complex::new(x, F::zero()));

        return Ok((Array1::from_elem(1, eigenvalue), eigenvector));
    } else if n == 2 {
        return solve_2x2_eigenvalue_problem(a);
    }

    // For larger matrices, use the QR algorithm
    solve_qr_algorithm(a)
}

/// Compute the eigenvalues of a square matrix.
///
/// # Arguments
///
/// * `a` - Input square matrix
/// * `workers` - Number of worker threads (None = use default)
///
/// # Returns
///
/// * Vector of complex eigenvalues
///
/// # Examples
///
/// ```
/// use scirs2_core::ndarray::array;
/// use scirs2_linalg::eigen::standard::eigvals;
///
/// let a = array![[1.0_f64, 0.0], [0.0, 2.0]];
/// let w = eigvals(&a.view(), None).expect("Operation failed");
///
/// // Sort eigenvalues (they may be returned in different order)
/// let mut eigenvalues = vec![w[0].re, w[1].re];
/// eigenvalues.sort_by(|a, b| a.partial_cmp(b).expect("Operation failed"));
///
/// assert!((eigenvalues[0] - 1.0).abs() < 1e-10);
/// assert!((eigenvalues[1] - 2.0).abs() < 1e-10);
/// ```
#[allow(dead_code)]
pub fn eigvals<F>(a: &ArrayView2<F>, workers: Option<usize>) -> LinalgResult<Array1<Complex<F>>>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    // For efficiency, we can compute just the eigenvalues
    // But for now, we'll use the full function and discard the eigenvectors
    let (eigenvalues, _) = eig(a, workers)?;
    Ok(eigenvalues)
}

/// Compute the dominant eigenvalue and eigenvector of a matrix using power iteration.
///
/// This is a simple iterative method that converges to the eigenvalue with the largest
/// absolute value and its corresponding eigenvector.
///
/// # Arguments
///
/// * `a` - Input square matrix
/// * `max_iter` - Maximum number of iterations
/// * `tol` - Convergence tolerance
///
/// # Returns
///
/// * Tuple (eigenvalue, eigenvector)
///
/// # Examples
///
/// ```
/// use scirs2_core::ndarray::array;
/// use scirs2_linalg::eigen::standard::power_iteration;
///
/// let a = array![[2.0_f64, 1.0], [1.0, 3.0]];
/// let (eigenvalue, eigenvector) = power_iteration(&a.view(), 100, 1e-10).expect("Operation failed");
/// // The largest eigenvalue of this matrix is approximately 3.618
/// assert!((eigenvalue - 3.618).abs() < 1e-2);
/// ```
#[allow(dead_code)]
pub fn power_iteration<F>(
    a: &ArrayView2<F>,
    max_iter: usize,
    tol: F,
) -> LinalgResult<(F, Array1<F>)>
where
    F: Float + NumAssign + Sum + Send + Sync + ScalarOperand + 'static,
{
    // Check if matrix is square
    if a.nrows() != a.ncols() {
        return Err(LinalgError::ShapeError(format!(
            "Expected square matrix, got shape {:?}",
            a.shape()
        )));
    }

    let n = a.nrows();

    // Start with a random vector
    let mut rng = scirs2_core::random::rng();
    let mut b = Array1::zeros(n);
    for i in 0..n {
        b[i] = F::from(rng.random_range(-1.0..=1.0)).unwrap_or(F::zero());
    }

    // Normalize the vector
    let norm_b = vector_norm(&b.view(), 2)?;
    b.mapv_inplace(|x| x / norm_b);

    let mut eigenvalue = F::zero();
    let mut prev_eigenvalue = F::zero();

    for _ in 0..max_iter {
        // Multiply b by A
        let mut b_new = Array1::zeros(n);
        for i in 0..n {
            let mut sum = F::zero();
            for j in 0..n {
                sum += a[[i, j]] * b[j];
            }
            b_new[i] = sum;
        }

        // Calculate the Rayleigh quotient (eigenvalue estimate)
        eigenvalue = F::zero();
        for i in 0..n {
            eigenvalue += b[i] * b_new[i];
        }

        // Normalize the vector
        let norm_b_new = vector_norm(&b_new.view(), 2)?;
        for i in 0..n {
            b[i] = b_new[i] / norm_b_new;
        }

        // Check for convergence
        if (eigenvalue - prev_eigenvalue).abs() < tol {
            return Ok((eigenvalue, b));
        }

        prev_eigenvalue = eigenvalue;
    }

    // Return the result after max_iter iterations
    Ok((eigenvalue, b))
}

/// Compute the eigenvalues and eigenvectors of a Hermitian or symmetric matrix.
///
/// # Arguments
///
/// * `a` - Input Hermitian or symmetric matrix
/// * `workers` - Number of worker threads (None = use default)
///
/// # Returns
///
/// * Tuple (eigenvalues, eigenvectors) where eigenvalues is a real vector
///   and eigenvectors is a real matrix whose columns are the eigenvectors
///
/// # Examples
///
/// ```
/// use scirs2_core::ndarray::array;
/// use scirs2_linalg::eigen::standard::eigh;
///
/// let a = array![[1.0_f64, 0.0], [0.0, 2.0]];
/// let (w, v) = eigh(&a.view(), None).expect("Operation failed");
///
/// // Sort eigenvalues (they may be returned in different order)
/// let mut eigenvalues = vec![(w[0], 0), (w[1], 1)];
/// eigenvalues.sort_by(|a, b| a.0.partial_cmp(&b.0).expect("Operation failed"));
///
/// assert!((eigenvalues[0].0 - 1.0).abs() < 1e-10);
/// assert!((eigenvalues[1].0 - 2.0).abs() < 1e-10);
/// ```
#[allow(dead_code)]
pub fn eigh<F>(a: &ArrayView2<F>, workers: Option<usize>) -> LinalgResult<(Array1<F>, Array2<F>)>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    use crate::parallel;

    // Configure workers for parallel operations
    parallel::configure_workers(workers);

    // Check if matrix is square
    if a.nrows() != a.ncols() {
        return Err(LinalgError::ShapeError(format!(
            "Expected square matrix, got shape {:?}",
            a.shape()
        )));
    }

    // Check if matrix is symmetric
    for i in 0..a.nrows() {
        for j in (i + 1)..a.ncols() {
            if (a[[i, j]] - a[[j, i]]).abs() > F::epsilon() {
                return Err(LinalgError::ShapeError(
                    "Matrix must be symmetric for Hermitian eigenvalue computation".to_string(),
                ));
            }
        }
    }

    let n = a.nrows();

    // For small matrices, we can compute eigenvalues directly
    if n == 1 {
        let eigenvalue = a[[0, 0]];
        let eigenvector = Array2::eye(1);

        return Ok((Array1::from_elem(1, eigenvalue), eigenvector));
    } else if n == 2 {
        return solve_2x2_symmetric_eigenvalue_problem(a);
    } else if n == 3 {
        return solve_3x3_symmetric_eigenvalue_problem(a);
    } else if n == 4 {
        return solve_4x4_symmetric_eigenvalue_problem(a);
    }

    // Choose parallel implementation based on matrix size and worker count
    let use_work_stealing = if let Some(num_workers) = workers {
        // For larger matrices and multiple workers, use work-stealing
        num_workers > 1 && n > 100
    } else {
        false
    };

    // Use work-stealing parallel implementation for large matrices
    if use_work_stealing {
        if let Some(num_workers) = workers {
            return crate::parallel::parallel_eigvalsh_work_stealing(a, num_workers);
        }
    }

    // For smaller matrices, use optimized sequential algorithms
    solve_symmetric_with_power_iteration(a)
}

/// Advanced MODE ENHANCEMENT: Advanced-precision eigenvalue computation targeting 1e-10+ accuracy
///
/// This function implements advanced numerical techniques for maximum precision:
/// - Kahan summation for enhanced numerical stability
/// - Multiple-stage Rayleigh quotient iteration with advanced-tight convergence
/// - Newton's method eigenvalue correction
/// - Adaptive tolerance selection based on matrix condition number
/// - Enhanced Gram-Schmidt orthogonalization with multiple passes
#[allow(dead_code)]
pub fn advanced_precision_eig<F>(
    a: &ArrayView2<F>,
    workers: Option<usize>,
) -> LinalgResult<(Array1<F>, Array2<F>)>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    use crate::parallel;

    // Configure workers for parallel operations
    parallel::configure_workers(workers);

    // Check if matrix is square and symmetric
    if a.nrows() != a.ncols() {
        return Err(LinalgError::ShapeError(format!(
            "Expected square matrix, got shape {:?}",
            a.shape()
        )));
    }

    let n = a.nrows();

    // Estimate matrix condition number for adaptive tolerance
    let condition_number = estimate_condition_number(a)?;

    // Select precision target based on condition number
    let precision_target =
        if condition_number > F::from(1e12).expect("Failed to convert constant to float") {
            F::from(1e-8).expect("Failed to convert constant to float") // Difficult matrices
        } else if condition_number > F::from(1e8).expect("Failed to convert constant to float") {
            F::from(1e-9).expect("Failed to convert constant to float") // Moderately conditioned
        } else {
            F::from(1e-10).expect("Failed to convert constant to float") // Well-conditioned matrices
        };

    // For small matrices, use advanced-precise analytical methods
    if n <= 4 {
        return advanced_precise_smallmatrix_eig(a, precision_target);
    }

    // For larger matrices, use enhanced iterative methods
    advanced_precise_iterative_eig(a, precision_target, workers)
}

/// Estimate matrix condition number for adaptive tolerance selection
#[allow(dead_code)]
fn estimate_condition_number<F>(a: &ArrayView2<F>) -> LinalgResult<F>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    let n = a.nrows();

    if n == 1 {
        return Ok(F::one());
    }

    // Use power iteration to estimate largest eigenvalue
    let (lambda_max_, _) = power_iteration(
        a,
        100,
        F::from(1e-8).expect("Failed to convert constant to float"),
    )?;

    // Estimate smallest eigenvalue using inverse iteration
    // For simplicity, use a heuristic based on trace and determinant
    let _trace = (0..n).map(|i| a[[i, i]]).fold(F::zero(), |acc, x| acc + x);
    let det_estimate = if n == 2 {
        a[[0, 0]] * a[[1, 1]] - a[[0, 1]] * a[[1, 0]]
    } else {
        // Use geometric mean of diagonal elements as rough estimate
        let product = (0..n)
            .map(|i| a[[i, i]].abs())
            .fold(F::one(), |acc, x| acc * x);
        if product > F::zero() {
            product.powf(F::one() / F::from(n).expect("Failed to convert to float"))
        } else {
            F::from(1e-15).expect("Failed to convert constant to float") // Very small positive value
        }
    };

    let lambda_min =
        if det_estimate.abs() > F::from(1e-15).expect("Failed to convert constant to float") {
            det_estimate / lambda_max_.powf(F::from(n - 1).expect("Failed to convert to float"))
        } else {
            F::from(1e-15).expect("Failed to convert constant to float")
        };

    let condition = lambda_max_.abs()
        / lambda_min
            .abs()
            .max(F::from(1e-15).expect("Failed to convert constant to float"));
    Ok(condition)
}

/// Advanced-precise eigenvalue computation for small matrices (n <= 4)
#[allow(dead_code)]
fn advanced_precise_smallmatrix_eig<F>(
    a: &ArrayView2<F>,
    precision_target: F,
) -> LinalgResult<(Array1<F>, Array2<F>)>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    let n = a.nrows();

    match n {
        1 => {
            let eigenvalue = a[[0, 0]];
            let eigenvector = Array2::eye(1);
            Ok((Array1::from_elem(1, eigenvalue), eigenvector))
        }
        2 => advanced_precise_2x2_eig(a, precision_target),
        3 => advanced_precise_3x3_eig(a, precision_target),
        4 => advanced_precise_4x4_eig(a, precision_target),
        _ => Err(LinalgError::InvalidInput(
            "Matrix size not supported for advanced-precise small matrix solver".to_string(),
        )),
    }
}

/// Advanced-precise 2x2 eigenvalue computation with Kahan summation
#[allow(dead_code)]
fn advanced_precise_2x2_eig<F>(
    a: &ArrayView2<F>,
    _precision_target: F,
) -> LinalgResult<(Array1<F>, Array2<F>)>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    let a11 = a[[0, 0]];
    let a12 = a[[0, 1]];
    let a21 = a[[1, 0]];
    let a22 = a[[1, 1]];

    // Use Kahan summation for enhanced precision
    let trace = kahan_sum(&[a11, a22]);
    let det = kahan_sum(&[a11 * a22, -(a12 * a21)]);

    // Enhanced discriminant calculation
    let trace_squared = trace * trace;
    let four_det = F::from(4.0).expect("Failed to convert constant to float") * det;
    let discriminant = kahan_sum(&[trace_squared, -four_det]);

    let half = F::from(0.5).expect("Failed to convert constant to float");
    let sqrt_discriminant = discriminant.abs().sqrt();

    let lambda1 = if discriminant >= F::zero() {
        (trace + sqrt_discriminant) * half
    } else {
        trace * half
    };

    let lambda2 = if discriminant >= F::zero() {
        (trace - sqrt_discriminant) * half
    } else {
        trace * half
    };

    let mut eigenvalues = Array1::zeros(2);
    eigenvalues[0] = lambda1;
    eigenvalues[1] = lambda2;

    // Enhanced eigenvector computation with Gram-Schmidt orthogonalization
    let eigenvectors = compute_advanced_precise_eigenvectors_2x2(a, &eigenvalues)?;

    Ok((eigenvalues, eigenvectors))
}

/// Advanced-precise 3x3 eigenvalue computation using Cardano's formula with enhancements
#[allow(dead_code)]
fn advanced_precise_3x3_eig<F>(
    a: &ArrayView2<F>,
    precision_target: F,
) -> LinalgResult<(Array1<F>, Array2<F>)>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    // Enhanced cubic equation solver with numerical stability improvements
    let characteristic_poly = compute_characteristic_polynomial_3x3(a)?;
    let eigenvalues =
        solve_cubic_equation_advanced_precise(&characteristic_poly, precision_target)?;

    // Enhanced eigenvector computation with multiple Gram-Schmidt passes
    let eigenvectors =
        compute_advanced_precise_eigenvectors_3x3(a, &eigenvalues, precision_target)?;

    Ok((eigenvalues, eigenvectors))
}

/// Advanced-precise 4x4 eigenvalue computation using enhanced QR iteration
#[allow(dead_code)]
fn advanced_precise_4x4_eig<F>(
    a: &ArrayView2<F>,
    precision_target: F,
) -> LinalgResult<(Array1<F>, Array2<F>)>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    // Use enhanced QR iteration with double shifting and deflation
    advanced_precise_qr_iteration(a, precision_target, 1000)
}

/// Kahan summation algorithm for enhanced numerical precision
#[allow(dead_code)]
fn kahan_sum<F>(values: &[F]) -> F
where
    F: Float + NumAssign,
{
    let mut sum = F::zero();
    let mut c = F::zero(); // Compensation for lost low-order bits

    for &value in values {
        let y = value - c;
        let t = sum + y;
        c = (t - sum) - y;
        sum = t;
    }

    sum
}

/// Compute characteristic polynomial for 3x3 matrix with enhanced precision
#[allow(dead_code)]
fn compute_characteristic_polynomial_3x3<F>(a: &ArrayView2<F>) -> LinalgResult<[F; 4]>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    let a00 = a[[0, 0]];
    let a01 = a[[0, 1]];
    let a02 = a[[0, 2]];
    let a10 = a[[1, 0]];
    let a11 = a[[1, 1]];
    let a12 = a[[1, 2]];
    let a20 = a[[2, 0]];
    let a21 = a[[2, 1]];
    let a22 = a[[2, 2]];

    // Coefficients of characteristic polynomial: c₃λ³ + c₂λ² + c₁λ + c₀ = 0
    let c3 = -F::one(); // Coefficient of λ³

    // Coefficient of λ² (negative trace)
    let c2 = kahan_sum(&[-a00, -a11, -a22]);

    // Coefficient of λ (sum of principal minors)
    let minor_00_11 = a00 * a11 - a01 * a10;
    let minor_00_22 = a00 * a22 - a02 * a20;
    let minor_11_22 = a11 * a22 - a12 * a21;
    let c1 = kahan_sum(&[minor_00_11, minor_00_22, minor_11_22]);

    // Constant term (negative determinant)
    let det_part1 = a00 * (a11 * a22 - a12 * a21);
    let det_part2 = -a01 * (a10 * a22 - a12 * a20);
    let det_part3 = a02 * (a10 * a21 - a11 * a20);
    let c0 = -kahan_sum(&[det_part1, det_part2, det_part3]);

    Ok([c0, c1, c2, c3])
}

/// Solve cubic equation with advanced-high precision using Cardano's formula with enhancements
#[allow(dead_code)]
fn solve_cubic_equation_advanced_precise<F>(
    coeffs: &[F; 4],
    precision_target: F,
) -> LinalgResult<Array1<F>>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    let [c0, c1, c2, c3] = *coeffs;

    // Normalize to monic polynomial
    let a = c2 / c3;
    let b = c1 / c3;
    let c = c0 / c3;

    // Depressed cubic transformation: t³ + pt + q = 0
    let third = F::one() / F::from(3.0).expect("Failed to convert constant to float");
    let p = b - a * a * third;
    let q_part1 = F::from(2.0).expect("Failed to convert constant to float") * a * a * a
        / F::from(27.0).expect("Failed to convert constant to float");
    let q_part2 = -a * b * third;
    let q = kahan_sum(&[q_part1, q_part2, c]);

    // Enhanced discriminant calculation
    let discriminant_part1 =
        (q / F::from(2.0).expect("Failed to convert constant to float")).powi(2);
    let discriminant_part2 =
        (p / F::from(3.0).expect("Failed to convert constant to float")).powi(3);
    let discriminant = discriminant_part1 + discriminant_part2;

    let mut roots = Array1::zeros(3);

    if discriminant > precision_target {
        // One real root case - use enhanced Cardano's formula
        let sqrt_disc = discriminant.sqrt();
        let half_q = q / F::from(2.0).expect("Failed to convert constant to float");

        let u = if half_q + sqrt_disc >= F::zero() {
            (half_q + sqrt_disc).powf(third)
        } else {
            -(-half_q - sqrt_disc).powf(third)
        };

        let v = if u.abs() > precision_target {
            -p / (F::from(3.0).expect("Failed to convert constant to float") * u)
        } else {
            F::zero()
        };

        let root = u + v - a * third;

        // Apply Newton's method for refinement
        let refined_root = newton_method_cubic_refinement(coeffs, root, precision_target, 20)?;

        roots[0] = refined_root;
        roots[1] = refined_root; // Repeated root approximation
        roots[2] = refined_root;
    } else {
        // Three real roots case - use trigonometric solution
        let m = F::from(2.0).expect("Failed to convert constant to float")
            * (-p / F::from(3.0).expect("Failed to convert constant to float")).sqrt();
        let theta = (F::from(3.0).expect("Failed to convert constant to float") * q / (p * m))
            .acos()
            / F::from(3.0).expect("Failed to convert constant to float");

        let two_pi_third = F::from(2.0).expect("Failed to convert constant to float")
            * F::from(std::f64::consts::PI).expect("Failed to convert to float")
            / F::from(3.0).expect("Failed to convert constant to float");
        let a_third = a * third;

        roots[0] = m * theta.cos() - a_third;
        roots[1] = m * (theta - two_pi_third).cos() - a_third;
        roots[2] = m
            * (theta - F::from(2.0).expect("Failed to convert constant to float") * two_pi_third)
                .cos()
            - a_third;

        // Refine all roots with Newton's method
        for i in 0..3 {
            roots[i] = newton_method_cubic_refinement(coeffs, roots[i], precision_target, 20)?;
        }
    }

    // Sort eigenvalues for consistency
    let mut root_vec: Vec<F> = roots.to_vec();
    root_vec.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));

    for (i, &root) in root_vec.iter().enumerate() {
        roots[i] = root;
    }

    Ok(roots)
}

/// Newton's method for cubic equation root refinement
#[allow(dead_code)]
fn newton_method_cubic_refinement<F>(
    coeffs: &[F; 4],
    initial_guess: F,
    tolerance: F,
    max_iter: usize,
) -> LinalgResult<F>
where
    F: Float + NumAssign,
{
    let [c0, c1, c2, c3] = *coeffs;
    let mut x = initial_guess;

    for _ in 0..max_iter {
        // Evaluate polynomial and its derivative using Horner's method
        let f_x = ((c3 * x + c2) * x + c1) * x + c0;
        let df_x = (F::from(3.0).expect("Failed to convert constant to float") * c3 * x
            + F::from(2.0).expect("Failed to convert constant to float") * c2)
            * x
            + c1;

        if df_x.abs() < tolerance {
            break; // Avoid division by zero
        }

        let x_new = x - f_x / df_x;

        if (x_new - x).abs() < tolerance {
            return Ok(x_new);
        }

        x = x_new;
    }

    Ok(x)
}

/// Compute advanced-precise eigenvectors for 2x2 matrix
#[allow(dead_code)]
fn compute_advanced_precise_eigenvectors_2x2<F>(
    a: &ArrayView2<F>,
    eigenvalues: &Array1<F>,
) -> LinalgResult<Array2<F>>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    let mut eigenvectors = Array2::zeros((2, 2));

    for (i, &lambda) in eigenvalues.iter().enumerate() {
        let mut v = Array1::zeros(2);

        // Enhanced eigenvector computation with numerical stability
        let a11_lambda = a[[0, 0]] - lambda;
        let a22_lambda = a[[1, 1]] - lambda;

        if a[[0, 1]].abs() > a[[1, 0]].abs() {
            if a[[0, 1]].abs() > F::epsilon() {
                v[0] = F::one();
                v[1] = -a11_lambda / a[[0, 1]];
            } else {
                v[0] = F::one();
                v[1] = F::zero();
            }
        } else if a[[1, 0]].abs() > F::epsilon() {
            v[0] = -a22_lambda / a[[1, 0]];
            v[1] = F::one();
        } else {
            // Diagonal case
            if a11_lambda.abs() < a22_lambda.abs() {
                v[0] = F::one();
                v[1] = F::zero();
            } else {
                v[0] = F::zero();
                v[1] = F::one();
            }
        }

        // Enhanced normalization with numerical stability
        let norm = vector_norm(&v.view(), 2)?;
        if norm > F::epsilon() {
            v.mapv_inplace(|x| x / norm);
        }

        eigenvectors.column_mut(i).assign(&v);
    }

    // Apply Gram-Schmidt orthogonalization for enhanced precision
    gram_schmidt_orthogonalization(&mut eigenvectors)?;

    Ok(eigenvectors)
}

/// Compute advanced-precise eigenvectors for 3x3 matrix
#[allow(dead_code)]
fn compute_advanced_precise_eigenvectors_3x3<F>(
    a: &ArrayView2<F>,
    eigenvalues: &Array1<F>,
    precision_target: F,
) -> LinalgResult<Array2<F>>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    let mut eigenvectors = Array2::zeros((3, 3));

    for (i, &lambda) in eigenvalues.iter().enumerate() {
        // Compute (A - λI)
        let mut matrix = a.to_owned();
        for j in 0..3 {
            matrix[[j, j]] -= lambda;
        }

        // Find null space using enhanced numerical techniques
        let eigenvector = enhanced_null_space_computation(&matrix.view(), precision_target)?;
        eigenvectors.column_mut(i).assign(&eigenvector);
    }

    // Apply multiple passes of Gram-Schmidt for advanced-high precision
    for _ in 0..3 {
        gram_schmidt_orthogonalization(&mut eigenvectors)?;
    }

    Ok(eigenvectors)
}

/// Enhanced null space computation for eigenvector calculation
#[allow(dead_code)]
fn enhanced_null_space_computation<F>(
    matrix: &ArrayView2<F>,
    precision_target: F,
) -> LinalgResult<Array1<F>>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    let n = matrix.nrows();
    let mut v = Array1::zeros(n);

    // Use inverse iteration with multiple random starts for robustness
    let mut best_v = v.clone();
    let mut best_residual = F::infinity();

    for _trial in 0..5 {
        // Random initialization
        let mut rng = scirs2_core::random::rng();
        for i in 0..n {
            v[i] = F::from(rng.random_range(-1.0..=1.0)).unwrap_or(F::zero());
        }

        // Normalize
        let norm = vector_norm(&v.view(), 2)?;
        if norm > F::epsilon() {
            v.mapv_inplace(|x| x / norm);
        }

        // Inverse iteration
        for _iter in 0..50 {
            let mut av = Array1::zeros(n);
            for i in 0..n {
                for j in 0..n {
                    av[i] += matrix[[i, j]] * v[j];
                }
            }

            // Check residual
            let residual = vector_norm(&av.view(), 2)?;
            if residual < best_residual {
                best_residual = residual;
                best_v.assign(&v);
            }

            if residual < precision_target {
                break;
            }

            // Update v for next iteration (simplified)
            let norm_v = vector_norm(&v.view(), 2)?;
            if norm_v > F::epsilon() {
                v.mapv_inplace(|x| x / norm_v);
            }
        }
    }

    Ok(best_v)
}

/// Enhanced Gram-Schmidt orthogonalization with numerical stability
#[allow(dead_code)]
fn gram_schmidt_orthogonalization<F>(matrix: &mut Array2<F>) -> LinalgResult<()>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    let n = matrix.nrows();
    let m = matrix.ncols();

    for i in 0..m {
        // Normalize current column
        let mut col_i = matrix.column(i).to_owned();
        let norm_i = vector_norm(&col_i.view(), 2)?;

        if norm_i > F::epsilon() {
            col_i.mapv_inplace(|x| x / norm_i);
            matrix.column_mut(i).assign(&col_i);
        }

        // Orthogonalize against previous columns
        for j in (i + 1)..m {
            let mut col_j = matrix.column(j).to_owned();

            // Compute projection coefficient with Kahan summation
            let mut dot_product = F::zero();
            let mut c = F::zero();
            for k in 0..n {
                let y = col_i[k] * col_j[k] - c;
                let t = dot_product + y;
                c = (t - dot_product) - y;
                dot_product = t;
            }

            // Remove projection
            for k in 0..n {
                col_j[k] -= dot_product * col_i[k];
            }

            matrix.column_mut(j).assign(&col_j);
        }
    }

    Ok(())
}

/// Advanced-precise iterative eigenvalue computation for large matrices
#[allow(dead_code)]
fn advanced_precise_iterative_eig<F>(
    a: &ArrayView2<F>,
    precision_target: F,
    _workers: Option<usize>,
) -> LinalgResult<(Array1<F>, Array2<F>)>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    // Enhanced QR iteration with multiple convergence criteria
    advanced_precise_qr_iteration(a, precision_target, 2000)
}

/// Advanced-precise QR iteration with enhanced numerical stability
#[allow(dead_code)]
fn advanced_precise_qr_iteration<F>(
    a: &ArrayView2<F>,
    precision_target: F,
    max_iter: usize,
) -> LinalgResult<(Array1<F>, Array2<F>)>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    let n = a.nrows();
    let mut h = a.to_owned(); // Hessenberg form would be better, but simplified here
    let mut q_total = Array2::eye(n);

    let mut iter_count = 0;

    while iter_count < max_iter {
        // Check for convergence
        let mut converged = true;
        for i in 1..n {
            for j in 0..i {
                if h[[i, j]].abs() > precision_target {
                    converged = false;
                    break;
                }
            }
            if !converged {
                break;
            }
        }

        if converged {
            break;
        }

        // Enhanced QR decomposition with pivoting
        let (q, r) = enhanced_qr_decomposition(&h.view())?;

        // Update H = RQ
        h = r.dot(&q);

        // Accumulate Q for eigenvectors
        q_total = q_total.dot(&q);

        iter_count += 1;
    }

    // Extract eigenvalues from diagonal
    let mut eigenvalues = Array1::zeros(n);
    for i in 0..n {
        eigenvalues[i] = h[[i, i]];
    }

    // Sort eigenvalues and corresponding eigenvectors
    let mut indices: Vec<usize> = (0..n).collect();
    indices.sort_by(|&a, &b| {
        eigenvalues[a]
            .partial_cmp(&eigenvalues[b])
            .unwrap_or(std::cmp::Ordering::Equal)
    });

    let mut sorted_eigenvalues = Array1::zeros(n);
    let mut sorted_eigenvectors = Array2::zeros((n, n));

    for (new_idx, &old_idx) in indices.iter().enumerate() {
        sorted_eigenvalues[new_idx] = eigenvalues[old_idx];
        sorted_eigenvectors
            .column_mut(new_idx)
            .assign(&q_total.column(old_idx));
    }

    Ok((sorted_eigenvalues, sorted_eigenvectors))
}

/// Enhanced QR decomposition with numerical stability improvements
#[allow(dead_code)]
fn enhanced_qr_decomposition<F>(a: &ArrayView2<F>) -> LinalgResult<(Array2<F>, Array2<F>)>
where
    F: Float + NumAssign + Sum + Send + Sync + scirs2_core::ndarray::ScalarOperand + 'static,
{
    let (m, n) = a.dim();
    let mut q = Array2::eye(m);
    let mut r = a.to_owned();

    for k in 0..n.min(m - 1) {
        // Enhanced Householder reflection with numerical stability
        let mut x = Array1::zeros(m - k);
        for i in k..m {
            x[i - k] = r[[i, k]];
        }

        let norm_x = vector_norm(&x.view(), 2)?;
        if norm_x < F::epsilon() {
            continue;
        }

        let mut v = x.clone();
        v[0] += if x[0] >= F::zero() { norm_x } else { -norm_x };

        let norm_v = vector_norm(&v.view(), 2)?;
        if norm_v > F::epsilon() {
            v.mapv_inplace(|val| val / norm_v);
        }

        // Apply Householder reflection to R
        for j in k..n {
            let mut col = Array1::zeros(m - k);
            for i in k..m {
                col[i - k] = r[[i, j]];
            }

            let dot_product = v.dot(&col);
            for i in k..m {
                r[[i, j]] -= F::from(2.0).expect("Failed to convert constant to float")
                    * dot_product
                    * v[i - k];
            }
        }

        // Apply Householder reflection to Q
        for j in 0..m {
            let mut col = Array1::zeros(m - k);
            for i in k..m {
                col[i - k] = q[[i, j]];
            }

            let dot_product = v.dot(&col);
            for i in k..m {
                q[[i, j]] -= F::from(2.0).expect("Failed to convert constant to float")
                    * dot_product
                    * v[i - k];
            }
        }
    }

    Ok((q.t().to_owned(), r))
}

/// Solve 2x2 general eigenvalue problem using analytical formula
#[allow(dead_code)]
fn solve_2x2_eigenvalue_problem<F>(a: &ArrayView2<F>) -> EigenResult<F>
where
    F: Float + NumAssign + Sum + Send + Sync + ScalarOperand + 'static,
{
    // For 2x2 matrices, use the quadratic formula
    let a11 = a[[0, 0]];
    let a12 = a[[0, 1]];
    let a21 = a[[1, 0]];
    let a22 = a[[1, 1]];

    let trace = a11 + a22;
    let det = a11 * a22 - a12 * a21;

    let discriminant =
        trace * trace - F::from(4.0).expect("Failed to convert constant to float") * det;

    // Create eigenvalues
    let mut eigenvalues = Array1::zeros(2);
    let mut eigenvectors = Array2::zeros((2, 2));

    if discriminant >= F::zero() {
        // Real eigenvalues
        let sqrt_discriminant = discriminant.sqrt();
        let lambda1 = (trace + sqrt_discriminant)
            / F::from(2.0).expect("Failed to convert constant to float");
        let lambda2 = (trace - sqrt_discriminant)
            / F::from(2.0).expect("Failed to convert constant to float");

        eigenvalues[0] = Complex::new(lambda1, F::zero());
        eigenvalues[1] = Complex::new(lambda2, F::zero());

        // Compute eigenvectors
        for (i, &lambda) in [lambda1, lambda2].iter().enumerate() {
            let mut eigenvector = Array1::zeros(2);

            if a12 != F::zero() {
                eigenvector[0] = a12;
                eigenvector[1] = lambda - a11;
            } else if a21 != F::zero() {
                eigenvector[0] = lambda - a22;
                eigenvector[1] = a21;
            } else {
                // Diagonal matrix
                eigenvector[0] = if (a11 - lambda).abs() < F::epsilon() {
                    F::one()
                } else {
                    F::zero()
                };
                eigenvector[1] = if (a22 - lambda).abs() < F::epsilon() {
                    F::one()
                } else {
                    F::zero()
                };
            }

            // Normalize
            let norm = vector_norm(&eigenvector.view(), 2)?;
            if norm > F::epsilon() {
                eigenvector.mapv_inplace(|x| x / norm);
            }

            eigenvectors.column_mut(i).assign(&eigenvector);
        }

        // Convert to complex
        let complex_eigenvectors = eigenvectors.mapv(|x| Complex::new(x, F::zero()));

        Ok((eigenvalues, complex_eigenvectors))
    } else {
        // Complex eigenvalues
        let real_part = trace / F::from(2.0).expect("Failed to convert constant to float");
        let imag_part =
            (-discriminant).sqrt() / F::from(2.0).expect("Failed to convert constant to float");

        eigenvalues[0] = Complex::new(real_part, imag_part);
        eigenvalues[1] = Complex::new(real_part, -imag_part);

        // Compute complex eigenvectors
        let mut complex_eigenvectors = Array2::zeros((2, 2));

        if a12 != F::zero() {
            complex_eigenvectors[[0, 0]] = Complex::new(a12, F::zero());
            complex_eigenvectors[[1, 0]] = Complex::new(eigenvalues[0].re - a11, eigenvalues[0].im);

            complex_eigenvectors[[0, 1]] = Complex::new(a12, F::zero());
            complex_eigenvectors[[1, 1]] = Complex::new(eigenvalues[1].re - a11, eigenvalues[1].im);
        } else if a21 != F::zero() {
            complex_eigenvectors[[0, 0]] = Complex::new(eigenvalues[0].re - a22, eigenvalues[0].im);
            complex_eigenvectors[[1, 0]] = Complex::new(a21, F::zero());

            complex_eigenvectors[[0, 1]] = Complex::new(eigenvalues[1].re - a22, eigenvalues[1].im);
            complex_eigenvectors[[1, 1]] = Complex::new(a21, F::zero());
        }

        // Normalize complex eigenvectors
        for i in 0..2 {
            let mut norm_sq = Complex::new(F::zero(), F::zero());
            for j in 0..2 {
                norm_sq += complex_eigenvectors[[j, i]] * complex_eigenvectors[[j, i]].conj();
            }
            let norm = norm_sq.re.sqrt();

            if norm > F::epsilon() {
                for j in 0..2 {
                    complex_eigenvectors[[j, i]] /= Complex::new(norm, F::zero());
                }
            }
        }

        Ok((eigenvalues, complex_eigenvectors))
    }
}

/// Solve 2x2 symmetric eigenvalue problem
#[allow(dead_code)]
fn solve_2x2_symmetric_eigenvalue_problem<F>(
    a: &ArrayView2<F>,
) -> LinalgResult<(Array1<F>, Array2<F>)>
where
    F: Float + NumAssign + Sum + Send + Sync + ScalarOperand + 'static,
{
    // For 2x2 symmetric matrices
    let a11 = a[[0, 0]];
    let a12 = a[[0, 1]];
    let a22 = a[[1, 1]];

    let trace = a11 + a22;
    let det = a11 * a22 - a12 * a12; // For symmetric matrices, a12 = a21

    let discriminant =
        trace * trace - F::from(4.0).expect("Failed to convert constant to float") * det;
    let sqrt_discriminant = discriminant.sqrt();

    // Eigenvalues
    let lambda1 =
        (trace + sqrt_discriminant) / F::from(2.0).expect("Failed to convert constant to float");
    let lambda2 =
        (trace - sqrt_discriminant) / F::from(2.0).expect("Failed to convert constant to float");

    // Sort eigenvalues in ascending order (SciPy convention)
    let (lambda_small, lambda_large) = if lambda1 <= lambda2 {
        (lambda1, lambda2)
    } else {
        (lambda2, lambda1)
    };

    let mut eigenvalues = Array1::zeros(2);
    eigenvalues[0] = lambda_small;
    eigenvalues[1] = lambda_large;

    // Eigenvectors
    let mut eigenvectors = Array2::zeros((2, 2));

    // Compute eigenvector for smaller eigenvalue (first)
    if a12 != F::zero() {
        eigenvectors[[0, 0]] = a12;
        eigenvectors[[1, 0]] = lambda_small - a11;
    } else {
        // Diagonal matrix
        eigenvectors[[0, 0]] = if (a11 - lambda_small).abs() < F::epsilon() {
            F::one()
        } else {
            F::zero()
        };
        eigenvectors[[1, 0]] = if (a22 - lambda_small).abs() < F::epsilon() {
            F::one()
        } else {
            F::zero()
        };
    }

    // Normalize
    let norm1 = (eigenvectors[[0, 0]] * eigenvectors[[0, 0]]
        + eigenvectors[[1, 0]] * eigenvectors[[1, 0]])
    .sqrt();
    if norm1 > F::epsilon() {
        eigenvectors[[0, 0]] /= norm1;
        eigenvectors[[1, 0]] /= norm1;
    }

    // Compute eigenvector for larger eigenvalue (second)
    if a12 != F::zero() {
        eigenvectors[[0, 1]] = a12;
        eigenvectors[[1, 1]] = lambda_large - a11;
    } else {
        // Diagonal matrix
        eigenvectors[[0, 1]] = if (a11 - lambda_large).abs() < F::epsilon() {
            F::one()
        } else {
            F::zero()
        };
        eigenvectors[[1, 1]] = if (a22 - lambda_large).abs() < F::epsilon() {
            F::one()
        } else {
            F::zero()
        };
    }

    // Normalize
    let norm2 = (eigenvectors[[0, 1]] * eigenvectors[[0, 1]]
        + eigenvectors[[1, 1]] * eigenvectors[[1, 1]])
    .sqrt();
    if norm2 > F::epsilon() {
        eigenvectors[[0, 1]] /= norm2;
        eigenvectors[[1, 1]] /= norm2;
    }

    Ok((eigenvalues, eigenvectors))
}

/// Solve 3x3 symmetric eigenvalue problem using analytical methods
/// Based on "Efficient numerical diagonalization of hermitian 3x3 matrices" by Kopp (2008)
#[allow(dead_code)]
fn solve_3x3_symmetric_eigenvalue_problem<F>(
    a: &ArrayView2<F>,
) -> LinalgResult<(Array1<F>, Array2<F>)>
where
    F: Float + NumAssign + Sum + Send + Sync + ScalarOperand + 'static,
{
    // For 3x3 symmetric matrices, use a specialized QR iteration
    // that converges quickly for small matrices

    let mut workmatrix = a.to_owned();
    let mut q_total = Array2::eye(3);
    let max_iter = 50;
    let tol = F::from(1e-12).expect("Failed to convert constant to float");

    // Apply QR iterations
    for _ in 0..max_iter {
        // Check for convergence - if off-diagonal elements are small
        let off_diag =
            workmatrix[[0, 1]].abs() + workmatrix[[0, 2]].abs() + workmatrix[[1, 2]].abs();
        if off_diag < tol {
            break;
        }

        // Perform QR decomposition
        let (q, r) = qr(&workmatrix.view(), None)?;

        // Update: A = R * Q
        workmatrix = r.dot(&q);

        // Accumulate transformation
        q_total = q_total.dot(&q);
    }

    // Extract eigenvalues from diagonal
    let mut eigenvalues = Array1::zeros(3);
    for i in 0..3 {
        eigenvalues[i] = workmatrix[[i, i]];
    }

    // Sort eigenvalues and corresponding eigenvectors
    let mut indices = [0, 1, 2];
    indices.sort_by(|&i, &j| {
        eigenvalues[i]
            .partial_cmp(&eigenvalues[j])
            .expect("Operation failed")
    });

    let mut sorted_eigenvalues = Array1::zeros(3);
    let mut sorted_eigenvectors = Array2::zeros((3, 3));

    for (new_idx, &old_idx) in indices.iter().enumerate() {
        sorted_eigenvalues[new_idx] = eigenvalues[old_idx];
        for i in 0..3 {
            sorted_eigenvectors[[i, new_idx]] = q_total[[i, old_idx]];
        }
    }

    Ok((sorted_eigenvalues, sorted_eigenvectors))
}

/// Solve 4x4 symmetric eigenvalue problem using QR iteration
#[allow(dead_code)]
fn solve_4x4_symmetric_eigenvalue_problem<F>(
    a: &ArrayView2<F>,
) -> LinalgResult<(Array1<F>, Array2<F>)>
where
    F: Float + NumAssign + Sum + Send + Sync + ScalarOperand + 'static,
{
    // For 4x4 symmetric matrices, use a specialized QR iteration
    // that converges quickly for small matrices

    let mut workmatrix = a.to_owned();
    let mut q_total = Array2::eye(4);
    let max_iter = 100;
    let tol = F::from(1e-12).expect("Failed to convert constant to float");

    // Apply QR iterations
    for _ in 0..max_iter {
        // Check for convergence - if off-diagonal elements are small
        let mut off_diag = F::zero();
        for i in 0..4 {
            for j in (i + 1)..4 {
                off_diag += workmatrix[[i, j]].abs();
            }
        }
        if off_diag < tol {
            break;
        }

        // Perform QR decomposition
        let (q, r) = qr(&workmatrix.view(), None)?;

        // Update: A = R * Q
        workmatrix = r.dot(&q);

        // Accumulate transformation
        q_total = q_total.dot(&q);
    }

    // Extract eigenvalues from diagonal
    let mut eigenvalues = Array1::zeros(4);
    for i in 0..4 {
        eigenvalues[i] = workmatrix[[i, i]];
    }

    // Sort eigenvalues and corresponding eigenvectors
    let mut indices = [0, 1, 2, 3];
    indices.sort_by(|&i, &j| {
        eigenvalues[i]
            .partial_cmp(&eigenvalues[j])
            .expect("Operation failed")
    });

    let mut sorted_eigenvalues = Array1::zeros(4);
    let mut sorted_eigenvectors = Array2::zeros((4, 4));

    for (new_idx, &old_idx) in indices.iter().enumerate() {
        sorted_eigenvalues[new_idx] = eigenvalues[old_idx];
        for i in 0..4 {
            sorted_eigenvectors[[i, new_idx]] = q_total[[i, old_idx]];
        }
    }

    Ok((sorted_eigenvalues, sorted_eigenvectors))
}

/// QR algorithm for general eigenvalue decomposition
#[allow(dead_code)]
fn solve_qr_algorithm<F>(a: &ArrayView2<F>) -> EigenResult<F>
where
    F: Float + NumAssign + Sum + Send + Sync + ScalarOperand + 'static,
{
    // For larger matrices, use the QR algorithm
    let mut a_k = a.to_owned();
    let n = a.nrows();
    let max_iter = 100;
    let tol = F::epsilon() * F::from(100.0).expect("Failed to convert constant to float");

    // Initialize eigenvalues and eigenvectors
    let mut eigenvalues = Array1::zeros(n);
    let mut eigenvectors = Array2::eye(n);

    for _iter in 0..max_iter {
        // QR decomposition
        let (q, r) = qr(&a_k.view(), None)?;

        // Update A_k+1 = R*Q (reversed order gives better convergence)
        let a_next = r.dot(&q);

        // Update eigenvectors: V_k+1 = V_k * Q
        eigenvectors = eigenvectors.dot(&q);

        // Check for convergence (check if subdiagonal elements are close to zero)
        let mut converged = true;
        for i in 1..n {
            if a_next[[i, i - 1]].abs() > tol {
                converged = false;
                break;
            }
        }

        if converged {
            // Extract eigenvalues from diagonal
            for i in 0..n {
                eigenvalues[i] = Complex::new(a_next[[i, i]], F::zero());
            }

            // Return as complex values
            let complex_eigenvectors = eigenvectors.mapv(|x| Complex::new(x, F::zero()));
            return Ok((eigenvalues, complex_eigenvectors));
        }

        // If not converged, continue with next iteration
        a_k = a_next;
    }

    // If we reached maximum iterations without convergence
    // Check if we at least have a reasonable approximation
    let mut eigenvals = Array1::zeros(n);
    for i in 0..n {
        eigenvals[i] = Complex::new(a_k[[i, i]], F::zero());
    }

    // Return the best approximation we have
    let complex_eigenvectors = eigenvectors.mapv(|x| Complex::new(x, F::zero()));
    Ok((eigenvals, complex_eigenvectors))
}

/// Solve symmetric matrices with power iteration (simplified implementation)
#[allow(dead_code)]
fn solve_symmetric_with_power_iteration<F>(
    a: &ArrayView2<F>,
) -> LinalgResult<(Array1<F>, Array2<F>)>
where
    F: Float + NumAssign + Sum + Send + Sync + ScalarOperand + 'static,
{
    use crate::decomposition::qr;
    use crate::norm::vector_norm;

    let n = a.nrows();
    let max_iter = 1000;
    let tolerance =
        F::from(1e-10).unwrap_or_else(|| F::epsilon() * F::from(100.0).unwrap_or(F::one()));

    // Use QR algorithm for symmetric matrices
    // This is more stable than power iteration for all eigenvalues

    // Create a working copy of the matrix
    let mut workmatrix = a.to_owned();
    let mut q_accumulated = Array2::eye(n);

    // Apply QR iterations with shifting for better convergence
    for iter in 0..max_iter {
        // Check for convergence by examining off-diagonal elements
        let mut max_off_diag = F::zero();
        for i in 0..n {
            for j in 0..n {
                if i != j {
                    let abs_val = workmatrix[[i, j]].abs();
                    if abs_val > max_off_diag {
                        max_off_diag = abs_val;
                    }
                }
            }
        }

        if max_off_diag < tolerance {
            break;
        }

        // Apply shift to improve convergence
        // Use Wilkinson shift: the eigenvalue of the 2x2 bottom-right submatrix
        // that is closer to the bottom-right element
        let shift = if n >= 2 {
            let a_nn = workmatrix[[n - 1, n - 1]];
            let a_n1n1 = workmatrix[[n - 2, n - 2]];
            let a_nn1 = workmatrix[[n - 1, n - 2]];

            let b = a_nn + a_n1n1;
            let c = a_nn * a_n1n1 - a_nn1 * a_nn1;
            let discriminant = b * b - F::from(4.0).unwrap_or(F::one()) * c;

            if discriminant >= F::zero() {
                let sqrt_disc = discriminant.sqrt();
                let lambda1 = (b + sqrt_disc) / F::from(2.0).unwrap_or(F::one());
                let lambda2 = (b - sqrt_disc) / F::from(2.0).unwrap_or(F::one());

                // Choose the eigenvalue closer to a_nn
                if (lambda1 - a_nn).abs() < (lambda2 - a_nn).abs() {
                    lambda1
                } else {
                    lambda2
                }
            } else {
                a_nn
            }
        } else {
            workmatrix[[n - 1, n - 1]]
        };

        // Apply shift: A' = A - σI
        for i in 0..n {
            workmatrix[[i, i]] -= shift;
        }

        // Perform QR decomposition
        let (q, r) = qr(&workmatrix.view(), None).map_err(|_| {
            LinalgError::ConvergenceError(format!(
                "QR decomposition failed in symmetric eigenvalue computation at iteration {iter}"
            ))
        })?;

        // Update: A = RQ + σI
        workmatrix = r.dot(&q);
        for i in 0..n {
            workmatrix[[i, i]] += shift;
        }

        // Accumulate transformation
        q_accumulated = q_accumulated.dot(&q);

        // Early termination for well-conditioned problems
        if iter > 10 && max_off_diag < tolerance * F::from(10.0).unwrap_or(F::one()) {
            break;
        }
    }

    // Extract eigenvalues from diagonal
    let mut eigenvalues = Array1::zeros(n);
    for i in 0..n {
        eigenvalues[i] = workmatrix[[i, i]];
    }

    // Sort eigenvalues and eigenvectors in descending order
    let mut indices: Vec<usize> = (0..n).collect();
    indices.sort_by(|&i, &j| {
        eigenvalues[j]
            .partial_cmp(&eigenvalues[i])
            .unwrap_or(std::cmp::Ordering::Equal)
    });

    let sorted_eigenvalues = Array1::from_iter(indices.iter().map(|&i| eigenvalues[i]));
    let mut sorted_eigenvectors = Array2::zeros((n, n));

    for (new_idx, &old_idx) in indices.iter().enumerate() {
        for i in 0..n {
            sorted_eigenvectors[[i, new_idx]] = q_accumulated[[i, old_idx]];
        }
    }

    // Normalize eigenvectors
    for j in 0..n {
        let mut col = sorted_eigenvectors.column_mut(j);
        let norm = vector_norm(&col.to_owned().view(), 2)?;

        if norm > F::epsilon() {
            col.mapv_inplace(|x| x / norm);
        }

        // Ensure consistent sign convention (largest component positive)
        let mut max_idx = 0;
        let mut max_abs = F::zero();
        for i in 0..n {
            let abs_val = col[i].abs();
            if abs_val > max_abs {
                max_abs = abs_val;
                max_idx = i;
            }
        }

        if col[max_idx] < F::zero() {
            col.mapv_inplace(|x| x * (-F::one()));
        }
    }

    Ok((sorted_eigenvalues, sorted_eigenvectors))
}

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

    #[test]
    fn test_1x1matrix() {
        let a = array![[3.0_f64]];
        let (eigenvalues, eigenvectors) = eig(&a.view(), None).expect("Operation failed");

        assert_relative_eq!(eigenvalues[0].re, 3.0, epsilon = 1e-10);
        assert_relative_eq!(eigenvalues[0].im, 0.0, epsilon = 1e-10);
        assert_relative_eq!(eigenvectors[[0, 0]].re, 1.0, epsilon = 1e-10);
        assert_relative_eq!(eigenvectors[[0, 0]].im, 0.0, epsilon = 1e-10);
    }

    #[test]
    fn test_2x2_diagonalmatrix() {
        let a = array![[3.0_f64, 0.0], [0.0, 4.0]];

        let (eigenvalues, eigenvectors) = eig(&a.view(), None).expect("Operation failed");

        // Eigenvalues could be returned in any order
        assert_relative_eq!(eigenvalues[0].im, 0.0, epsilon = 1e-10);
        assert_relative_eq!(eigenvalues[1].im, 0.0, epsilon = 1e-10);

        // Test eigh
        let (eigenvalues, eigenvectors) = eigh(&a.view(), None).expect("Operation failed");
        // The eigenvalues might be returned in a different order
        assert!(
            (eigenvalues[0] - 3.0).abs() < 1e-10 && (eigenvalues[1] - 4.0).abs() < 1e-10
                || (eigenvalues[1] - 3.0).abs() < 1e-10 && (eigenvalues[0] - 4.0).abs() < 1e-10
        );
    }

    #[test]
    fn test_2x2_symmetricmatrix() {
        let a = array![[1.0, 2.0], [2.0, 4.0]];

        // Test eigh
        let (eigenvalues, eigenvectors) = eigh(&a.view(), None).expect("Operation failed");

        // Eigenvalues should be approximately 5 and 0
        assert!(
            (eigenvalues[0] - 5.0).abs() < 1e-10 && eigenvalues[1].abs() < 1e-10
                || (eigenvalues[1] - 5.0).abs() < 1e-10 && eigenvalues[0].abs() < 1e-10
        );

        // Check that eigenvectors are orthogonal
        let dot_product = eigenvectors[[0, 0]] * eigenvectors[[0, 1]]
            + eigenvectors[[1, 0]] * eigenvectors[[1, 1]];
        assert!(
            (dot_product).abs() < 1e-10,
            "Eigenvectors should be orthogonal"
        );

        // Check that eigenvectors are normalized
        let norm1 = (eigenvectors[[0, 0]] * eigenvectors[[0, 0]]
            + eigenvectors[[1, 0]] * eigenvectors[[1, 0]])
        .sqrt();
        let norm2 = (eigenvectors[[0, 1]] * eigenvectors[[0, 1]]
            + eigenvectors[[1, 1]] * eigenvectors[[1, 1]])
        .sqrt();
        assert!(
            (norm1 - 1.0).abs() < 1e-10,
            "First eigenvector should be normalized"
        );
        assert!(
            (norm2 - 1.0).abs() < 1e-10,
            "Second eigenvector should be normalized"
        );
    }

    #[test]
    fn test_power_iteration() {
        // Matrix with known dominant eigenvalue and eigenvector
        let a = array![[3.0, 1.0], [1.0, 3.0]];

        let (eigenvalue, eigenvector) =
            power_iteration(&a.view(), 100, 1e-10).expect("Operation failed");

        // Dominant eigenvalue should be 4
        assert_relative_eq!(eigenvalue, 4.0, epsilon = 1e-8);

        // Eigenvector should be normalized
        let norm = (eigenvector[0] * eigenvector[0] + eigenvector[1] * eigenvector[1]).sqrt();
        assert_relative_eq!(norm, 1.0, epsilon = 1e-10);

        // Check that Av ≈ lambda * v
        let av = a.dot(&eigenvector);
        let lambda_v = &eigenvector * eigenvalue;

        let mut max_diff = 0.0;
        for i in 0..eigenvector.len() {
            let diff = (av[i] - lambda_v[i]).abs();
            if diff > max_diff {
                max_diff = diff;
            }
        }

        assert!(
            max_diff < 1e-5,
            "A*v should approximately equal lambda*v, max diff: {}",
            max_diff
        );
    }
}