use std::{ops::Deref, rc::Rc};

use crate::integrater::IntegrationError;
use fynsul_core::EvalError;
// consider not allowing results for integration

const SIMPSON_1_PRECISION: u64 = 1000000;
const SIMPSON_2_PRECISION: u64 = 3 * 100000;

fn simpson_0(
    f: &impl Fn(f64) -> Result<f64, EvalError>,
    a: f64,
    fa: f64,
    b: f64,
    fb: f64,
) -> Result<(f64, f64, f64), IntegrationError> {
    let m = (a + b) / 2.0;
    match f(m) {
        Ok(fm) => Ok((m, fm, (b - a) / 6.0 * (fa + 4.0 * fm + fb))),
        Err(_e) => Err(IntegrationError::DivergentIntegral),
    }
}

pub fn simpson_1(
    f: &impl Fn(f64) -> Result<f64, EvalError>,
    a: f64,
    b: f64,
) -> Result<f64, IntegrationError> {
    let n = ((b - a).abs().round() as u64).max(2) * SIMPSON_1_PRECISION;
    let h = (b - a) / n as f64;

    let mut integral = 0.0;
    if let Ok(left_end) = f(a) {
        integral += left_end;
    } else {
        return Err(IntegrationError::DivergentIntegral);
    }

    if let Ok(right_end) = f(b) {
        integral += right_end;
    } else {
        return Err(IntegrationError::DivergentIntegral);
    }

    for i in 0..n {
        let x = a + i as f64 * h;
        if let Ok(res) = f(x) {
            if i % 2 == 0 {
                integral += 2.0 * res;
            } else {
                integral += 4.0 * res;
            }
        } else {
            return Err(IntegrationError::DivergentIntegral);
        }
    }
    integral = integral * h / 3.0;
    Ok(integral)
}

// consider not allowing results for integration
pub fn simpson_2(
    f: &impl Fn(f64) -> Result<f64, EvalError>,
    a: f64,
    b: f64,
) -> Result<f64, IntegrationError> {
    let n = SIMPSON_2_PRECISION * ((b - a).abs().round() as u64).max(3);
    let h = (b - a) / n as f64;

    let mut integral = 0.0;
    if let Ok(left_end) = f(a) {
        integral += left_end;
    } else {
        return Err(IntegrationError::DivergentIntegral);
    }

    if let Ok(right_end) = f(b) {
        integral += right_end;
    } else {
        return Err(IntegrationError::DivergentIntegral);
    }

    for i in 0..n {
        let x = a + i as f64 * h;
        if let Ok(res) = f(x) {
            if i % 3 == 0 {
                integral += 2.0 * res;
            } else {
                integral += 3.0 * res;
            }
        } else {
            return Err(IntegrationError::DivergentIntegral);
        }
    }
    integral = integral * 3.0 * h / 8.0;
    Ok(integral)
}

#[allow(clippy::too_many_arguments)]
fn recursive_adap(
    f: Rc<impl Fn(f64) -> Result<f64, EvalError>>,
    a: f64,
    fa: f64,
    b: f64,
    fb: f64,
    eps: f64,
    whole: f64,
    m: f64,
    fm: f64,
) -> Result<f64, IntegrationError> {
    let (lm, flm, left) = simpson_0(f.deref(), a, fa, m, fm)?;
    let (rm, frm, right) = simpson_0(f.deref(), m, fm, b, fb)?;
    let delta = left + right - whole;
    if delta.abs() <= 15.0 * eps {
        Ok(left + right + delta / 15.0)
    } else {
        let left = recursive_adap(f.clone(), a, fa, m, fm, eps / 2.0, left, lm, flm)?;
        let right = recursive_adap(f.clone(), m, fm, b, fb, eps / 2.0, right, rm, frm)?;
        Ok(left + right)
    }
}

pub fn adaptive_simpson(
    f: &impl Fn(f64) -> Result<f64, EvalError>,
    a: f64,
    b: f64,
    eps: f64,
) -> Result<f64, IntegrationError> {
    if let (Ok(fa), Ok(fb)) = (f(a), f(b)) {
        let (m, fm, whole) = simpson_0(&f, a, fa, b, fb)?;
        let f_rc = Rc::new(f);
        recursive_adap(f_rc, a, fa, b, fb, eps, whole, m, fm)
    } else {
        Err(IntegrationError::DivergentIntegral)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    #[test]
    fn test_simpson_1() {
        // Test simpson_1 with a simple linear function: integral of f(x) = x from 0 to 1 is 0.5
        let result = simpson_1(&|x| Ok(x), 0.0, 1.0).unwrap();
        assert!((result - 0.5).abs() < 1e-5);

        // Test simpson_1 with a quadratic function: integral of f(x) = x^2 from 0 to 1 is 1/3
        let result = simpson_1(&|x| Ok(x * x), 0.0, 1.0).unwrap();
        assert!((result - (1.0 / 3.0)).abs() < 1e-5);

        // Test simpson_1 with a sine function: integral of sin(x) from 0 to pi is 2
        let result = simpson_1(&|x| Ok(x.sin()), 0.0, std::f64::consts::PI).unwrap();
        assert!((result - 2.0).abs() < 1e-5);
    }

    #[test]
    fn test_simpson_2() {
        // Test simpson_2 with a cubic function: integral of f(x) = x^3 from 0 to 1 is 0.25
        let result = simpson_2(&|x| Ok(x.powi(3)), 0.0, 1.0).unwrap();
        assert!((result - 0.25).abs() < 1e-5);

        // Test simpson_2 with an exponential function: integral of exp(x) from 0 to 1 is e - 1
        let result = simpson_2(&|x| Ok(x.exp()), 0.0, 1.0).unwrap();
        assert!((result - (std::f64::consts::E - 1.0)).abs() < 1e-5);

        // Test simpson_2 with a constant function: integral of f(x) = 5 from 0 to 2 is 10
        let result = simpson_2(&|_| Ok(5.0), 0.0, 2.0).unwrap();
        assert!((result - 10.0).abs() < 1e-5);
    }

    #[test]
    fn test_adaptive_simpson() {
        // Test adaptive_simpson with a linear function: integral of f(x) = x from 0 to 1 is 0.5
        let result = adaptive_simpson(&|x| Ok(x), 0.0, 1.0, 1e-5).unwrap();
        assert!((result - 0.5).abs() < 1e-5);

        // Test adaptive_simpson with a quadratic function: integral of f(x) = x^2 from 0 to 1 is 1/3
        let result = adaptive_simpson(&|x| Ok(x * x), 0.0, 1.0, 1e-5).unwrap();
        assert!((result - (1.0 / 3.0)).abs() < 1e-5);

        // Test adaptive_simpson with a sine function: integral of sin(x) from 0 to pi is 2
        let result = adaptive_simpson(&|x| Ok(x.sin()), 0.0, std::f64::consts::PI, 1e-5).unwrap();
        assert!((result - 2.0).abs() < 1e-5);

        // Test adaptive_simpson with an exponential function: integral of exp(x) from 0 to 1 is e - 1
        let result = adaptive_simpson(&|x| Ok(x.exp()), 0.0, 1.0, 1e-5).unwrap();
        assert!((result - (std::f64::consts::E - 1.0)).abs() < 1e-5);

        // Test adaptive_simpson with a constant function: integral of f(x) = 3 from 0 to 4 is 12
        let result = adaptive_simpson(&|_| Ok(3.0), 0.0, 4.0, 1e-5).unwrap();
        assert!((result - 12.0).abs() < 1e-5);
    }

    #[test]
    fn test_simpson_1_negative_interval() {
        // Test simpson_1 with reversed interval: integral of f(x) = x from 1 to 0 should be -0.5
        let result = simpson_1(&|x| Ok(x), 1.0, 0.0).unwrap();
        assert!((result + 0.5).abs() < 1e-5);
    }

    #[test]
    fn test_simpson_2_negative_interval() {
        // Test simpson_2 with reversed interval: integral of f(x) = x^3 from 1 to 0 should be -0.25
        let result = simpson_2(&|x| Ok(x.powi(3)), 1.0, 0.0).unwrap();
        assert!((result + 0.25).abs() < 1e-5);
    }

    #[test]
    fn test_adaptive_simpson_negative_interval() {
        // Test adaptive_simpson with reversed interval: integral of f(x) = x from 1 to 0 should be -0.5
        let result = adaptive_simpson(&|x| Ok(x), 1.0, 0.0, 1e-5).unwrap();
        assert!((result + 0.5).abs() < 1e-5);
    }

    #[test]
    fn test_simpson_1_small_interval() {
        // Test simpson_1 with a very small interval
        let result = simpson_1(&|x| Ok(x.sin()), 0.0, 1e-5).unwrap();
        let expected = 0.5 * 1e-5 * 1e-5; // sin(x) ~ x for small x, so integral ~ x^2/2
        assert!((result - expected).abs() < 1e-5);
    }

    #[test]
    fn test_simpson_2_small_interval() {
        // Test simpson_2 with a very small interval
        let result = simpson_2(&|x| Ok(x.cos()), 0.0, 1e-5).unwrap();
        let expected = 1e-5; // cos(x) ~ 1 for small x, so integral ~ x
        assert!((result - expected).abs() < 1e-5);
    }

    #[test]
    fn test_adaptive_simpson_small_interval() {
        // Test adaptive_simpson with a very small interval
        let result = adaptive_simpson(&|x| Ok(x.exp()), 0.0, 1e-5, 1e-12).unwrap();
        let expected = 1e-5 + 0.5 * 1e-12; // exp(x) ~ 1 + x for small x, so integral ~ x + x^2/2
        assert!((result - expected).abs() < 1e-5);
    }
}
