use num_complex::{Complex, Complex64, ComplexFloat};

use crate::filters::{iir::lib::{aberth::{aberth, make_f, make_fp}, kve, newton::{newton, newton_complex}, polyval::{polyval_scalar, polyval_vec}}, BesselCriterion};

fn campos_zeros(n: usize) -> Vec<Complex<f64>> {
    if n == 1 {
        return vec![Complex::new(-1.0, 0.0)];
    }

    // s = polyval(n, [0, 0, 2, 0, -3, 1])
    let s = polyval_scalar(n as f64, &[0.0, 0.0, 2.0, 0.0, -3.0, 1.0]);
    let b3 = polyval_scalar(n as f64, &[16.0, -8.0]) / s;
    let b2 = polyval_scalar(n as f64, &[-24.0, -12.0, 12.0]) / s;
    let b1 = polyval_scalar(n as f64, &[8.0, 24.0, -12.0, -2.0]) / s;
    let b0 = polyval_scalar(n as f64, &[0.0, -6.0, 0.0, 5.0, -1.0]) / s;

    let r = polyval_scalar(n as f64, &[0.0, 0.0, 2.0, 1.0]);
    let a1 = polyval_scalar(n as f64, &[-6.0, -6.0]) / r;
    let a2 = 6.0 / r;

    let k: Vec<f64> = (1..=n).map(|i| i as f64).collect();
    let x_vals = polyval_vec(&k, &[0.0, a1, a2]);
    let y_vals = polyval_vec(&k, &[b0, b1, b2, b3]);

    x_vals
        .into_iter()
        .zip(y_vals.into_iter())
        .map(|(x, y)| Complex::new(x, y))
        .collect()
}

/// 模拟计算 Bessel 多项式的零点
fn bessel_zeros(n: usize) -> Result<Vec<Complex<f64>>, String> {
    if n == 0 {
        return Err("n == 0 !".into());
    }

    let x0 = campos_zeros(n);

    let f = make_f(n);
    let fp = make_fp(n);
    match aberth(&f, &fp, &x0,  1e-15, 50) {
        Ok(mut x) => {
            // Step 2: Newton refine
            x.iter_mut().for_each(|xi: &mut Complex<f64>| {
                let func = |z: Complex<f64>| {
                    let inv = Complex64::new(1.0, 0.0) / z; // real part only for placeholder
                    kve(n as f64 + 0.5, inv, 50)
                };
                let dfunc = |z: Complex<f64>| {
                    let inv = Complex64::new(1.0, 0.0) / z;
                    let x2 = z * z;
                    let term1 = kve(n as f64 - 0.5, inv, 50) / (2.0 * x2);
                    let term2 = kve(n as f64 + 0.5, inv, 50) / x2;
                    let term3 = kve(n as f64 + 1.5, inv, 50) / (2.0 * x2);
                    term1 - term2 + term3
                };
                if let Ok(result) = newton_complex(
                    &func,
                    *xi,
                    Some(&dfunc),
                    1e-8,
                    50,
                    None,
                    None,
                    0.0,
                    true,
                ){
                    *xi = result
                }
            });
            // Step 3: Average conjugate pairs
            let len = x.len();
            for i in 0..len / 2 {
                let avg = (x[i] + x[len - 1 - i].conj()) / 2.0;
                x[i] = avg;
                x[len - 1 - i] = avg.conj();
            }
            // Step 4: Verify sum ≈ -1
            let sum: Complex<f64> = x.iter().copied().sum();
            if (sum + Complex::new(1.0, 0.0)).norm() > 1e-12 {
                panic!("Generated zeros are inaccurate (sum != -1)");
            }
            Ok(x)
        },
        Err(e) => {
            Err(e)
        }
    }
}

fn norm_factor(p: &[Complex<f64>], k: f64) -> Result<f64, String> {
    // 构建 cutoff 闭包
    let g = |w: f64| {
        let jw = Complex::new(0.0, w);
        let prod = p.iter()
            .map(|&pole| jw - pole)
            .fold(Complex::new(1.0, 0.0), |acc, x| acc * x);
        k / prod.norm()
    };

    let cutoff = |w: f64| {
        g(w) - 1.0 / 2.0_f64.sqrt()
    };
    println!("cutoff{:?}", cutoff(1.4));

    // 使用已实现的 newton 方法求解 cutoff(w) = 0，初始值为 1.5
    let result = newton(
        &cutoff,
        1.5,
        None,
        1e-8,
        50,
        None,
        None,
        0.0,
        true,
    )?;
    Ok(result)
}

pub fn besselap(n: usize, norm: BesselCriterion) -> Result<(Vec<Complex<f64>>, Vec<Complex<f64>>, f64), String> {
    if n == 0 {
        return Ok((vec![], vec![], 1.0));
    }

    let mut poles: Vec<Complex<f64>> = bessel_zeros(n)?
        .into_iter()
        .map(|z| z.recip())
        .collect();

    let a_last = falling_factorial((2 * n) as i64, n as i64) / (1 << n);

    let k: f64;
    match norm {
        BesselCriterion::Delay | BesselCriterion::Mag => {
            k = a_last as f64;
            if norm == BesselCriterion::Mag {
                let norm_fac = norm_factor(&poles, k)?;
                poles.iter_mut().for_each(|p| *p = *p / norm_fac);
                let updated_k = norm_fac.powf(-(n as f64)) * a_last as f64;
                return Ok((vec![], poles, updated_k));
            }
        }
        BesselCriterion::Phase => {
            let scale = 10f64.powf(-(a_last as f64).log10() / n as f64);
            poles.iter_mut().for_each(|p| *p *= scale);
            k = 1.0;
        }
    }
    Ok((vec![], poles, k))

}
fn falling_factorial(x: i64, n: i64) -> i64 {
    let mut val = 1;
    for k in (x - n + 1)..=x {
        val *= k;
    }
    val
}