computorv1 1.0.11

use crate::aux::Polynomial;
use crate::solver::degree_2::degree_2;
use crate::solver::Roots;
use ft_lib::{
    ft_abs::ft_abs, ft_cbrt::ft_cbrt, ft_cos::ft_cos, ft_powi::ft_powi, ft_sqrt::ft_sqrt,
};
use std::f64::consts::PI;

// Degree 3 Polynomial (Cubic):
// a3 * x^3 + a2 * x^2 + a1 * x + a0 = 0, where a3 != 0
// Reduce to depressed cubic by Tschirnhaus transformation:
// Let x = t - a2 / (3 * a3)
// Resulting in: t^3 + p * t + q = 0
// where:
//   p = (3*a3*a1 - a2^2) / (3*a3^2)
//   q = (2*a2^3 - 9*a3*a2*a1 + 27*a3^2*a0) / (27*a3^3)
// Cardano’s Formula:
// t = cbrt(-q/2 + sqrt((q/2)^2 + (p/3)^3)) +
//     cbrt(-q/2 - sqrt((q/2)^2 + (p/3)^3))
// Then x = t - a2 / (3 * a3)
// Depending on the discriminant:
// Δ = (q/2)^2 + (p/3)^3
// - Δ > 0: one real root, two complex
// - Δ = 0: multiple real roots
// - Δ < 0: three real roots

pub fn degree_3(poly: Polynomial) -> Result<Roots, String> {
    let a: f64 = poly.terms.get(&3).copied().unwrap_or(0.0);
    let b: f64 = poly.terms.get(&2).copied().unwrap_or(0.0);
    let c: f64 = poly.terms.get(&1).copied().unwrap_or(0.0);
    let d: f64 = poly.terms.get(&0).copied().unwrap_or(0.0);

    if a == 0.0 {
        return degree_2(poly);
    }

    // Compute depressed cubic coefficients (t³ + p t + q = 0)
    let p: f64 = (3.0 * a * c - b * b) / (3.0 * a * a);
    let q: f64 = (2.0 * b * b * b - 9.0 * a * b * c + 27.0 * a * a * d) / (27.0 * a * a * a);

    // Discriminant Δ = (q/2)² + (p/3)³
    let delta: f64 = ft_powi(q / 2.0, 2) + ft_powi(p / 3.0, 3);
    // Shift to recover x from t (x = t - b/(3a))
    let shift: f64 = -b / (3.0 * a);

    let mut roots: Roots = Roots {
        real: Vec::new(),
        complex: Vec::new(),
    };

    if delta > 0.0 {
        // Case 1: Δ > 0 → One real root, two complex roots
        let sqrt_delta: f64 = ft_sqrt(delta);
        let u: f64 = ft_cbrt((-q / 2.0) + sqrt_delta);
        let v: f64 = ft_cbrt((-q / 2.0) - sqrt_delta);
        let real_root: f64 = u + v + shift;
        roots
            .real
            .push(format!("{:.5}", real_root).parse().unwrap());

        // Complex roots: -(u + v)/2 ± (u - v)√3 i / 2 + shift
        let real_part: f64 = -(u + v) / 2.0 + shift;
        let imag_part: f64 = (u - v) * (ft_sqrt(3.0) / 2.0);

        roots.complex.push((
            format!("{:.5}", real_part).parse().unwrap(),
            format!("{:.5}", imag_part).parse().unwrap(),
        ));
        roots.complex.push((
            format!("{:.5}", real_part).parse().unwrap(),
            format!("{:.5}", -imag_part).parse().unwrap(),
        ));
    } else if ft_abs(delta) < 1e-12 {
        // Case 2: Δ = 0 → Multiple real roots
        let u: f64 = ft_cbrt(-q / 2.0);
        let root1: f64 = 2.0 * u + shift;
        let root2: f64 = -u + shift;
        if ft_abs(root1 - root2) < 1e-12 {
            roots.real.push(format!("{:.5}", root1).parse().unwrap());
        } else {
            roots.real.push(format!("{:.5}", root1).parse().unwrap());
            roots.real.push(format!("{:.5}", root2).parse().unwrap());
        }
    } else {
        // all roots are real and distinct
        let r: f64 = (-p / 3.0).sqrt();
        let phi: f64 = ((-q / 2.0) / r.powi(3)).acos();
        let root1: f64 = 2.0 * r * ft_cos(phi / 3.0) + shift;
        let root2: f64 = 2.0 * r * ft_cos((phi + 2.0 * PI) / 3.0) + shift;
        let root3: f64 = 2.0 * r * ft_cos((phi + 4.0 * PI) / 3.0) + shift;
        roots.real.push(format!("{:.5}", root1).parse().unwrap());
        roots.real.push(format!("{:.5}", root2).parse().unwrap());
        roots.real.push(format!("{:.5}", root3).parse().unwrap());
    }
    Ok(roots)
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::parser::parse_equation;

    mod degree_3 {
        use super::*;

        #[test]
        fn test_one_real_root_with_parse() {
            let equation: &str = "X^3 = 1";
            let poly: Polynomial = parse_equation(equation).expect("Failed to parse equation");
            let result: Roots = degree_3(poly).unwrap();
            assert_eq!(result.real.len(), 1);
            assert!((result.real[0] - 1.0).abs() < 1e-6);
        }

        #[test]
        fn test_triple_root_with_parse() {
            let equation: &str = "X^3 - 3*X^2 + 3*X - 1 = 0";
            let poly: Polynomial = parse_equation(equation).expect("Failed to parse equation");
            let result: Roots = degree_3(poly).unwrap();
            assert_eq!(result.real.len(), 1);
            assert_eq!(result.real[0], 1.0);

            let equation: &str = "3X^3 - 16X^2 + 23X - 6 = 0";
            let poly: Polynomial = parse_equation(equation).expect("Failed to parse equation");
            let result: Roots = degree_3(poly).unwrap();
            assert_eq!(result.real.len(), 3);
            assert_eq!(result.real[0], 3.00000);
            assert_eq!(result.real[1], 0.33333);
            assert_eq!(result.real[2], 2.00000);
        }

        #[test]
        fn test_double_and_single_root_with_parse() {
            let equation: &str = "X^3 - 5*X^2 + 8*X - 4 = 0";
            let poly: Polynomial = parse_equation(equation).expect("Failed to parse equation");
            let result: Roots = degree_3(poly).unwrap();

            assert_eq!(result.real.len(), 2);
            assert_eq!(result.real[0], 1.00000);
            assert_eq!(result.real[1], 2.00000);
        }

        #[test]
        fn test_three_distinct_real_roots_with_parse() {
            let equation: &'static str = "X^3 - 7*X + 6 = 0";
            let poly: Polynomial = parse_equation(equation).expect("Failed to parse equation");
            let result: Roots = degree_3(poly).unwrap();

            assert_eq!(result.real.len(), 3);
            for root in [-3.0, 1.0, 2.0] {
                assert!(result.real.iter().any(|&x| (x - root).abs() < 1e-6));
            }
        }

        #[test]
        fn test_not_a_cubic_with_parse() {
            let equation: &str = "X^2 - 3*X + 2 = 0";
            let poly: Polynomial = parse_equation(equation).expect("Failed to parse equation");
            let result: Roots = degree_3(poly).unwrap();
            let expected: Vec<f64> = vec![2.0, 1.0];
            assert_eq!(result.real, expected);
        }

        #[test]
        fn test_not_a_cubic_all_0() {
            let equation: &str = "0X^3 + 0X^2 - 0*X + 0 = 0";
            let poly: Polynomial = parse_equation(equation).expect("Failed to parse equation");
            let result: Roots = degree_3(poly).unwrap();
            assert_eq!(result.real.len(), 0);
        }
    }

    mod degree_3_complex {
        use super::*;
        #[test]
        fn test_complex_roots() {
            let equation: &str = "X^3-3X^2+ 6X -4 = 0";
            let poly: Polynomial = parse_equation(equation).expect("Failed to parse equation");
            let result: Roots = degree_3(poly).unwrap();
            assert_eq!(result.real.len(), 1);
            assert_eq!(result.complex.len(), 2);
            assert_eq!(result.real[0], 1.0);
            assert_eq!(result.complex[0].0, 1.0);
            assert_eq!(result.complex[0].1, 1.73205);
            assert_eq!(result.complex[1].0, 1.0);
            assert_eq!(result.complex[1].1, -1.73205);
        }
        #[test]
        fn test_complex_roots_case1() {
            let equation: &str = "X^3 - 2X^2 + 4X -8 = 0";
            let poly: Polynomial = parse_equation(equation).expect("Failed to parse equation");
            let result: Roots = degree_3(poly).unwrap();
            assert_eq!(result.real.len(), 1);
            assert_eq!(result.complex.len(), 2);
            assert_eq!(result.real[0], 2.0);
            assert_eq!(result.complex[0].0, 0.0);
            assert_eq!(result.complex[0].1, 2.0);
            assert_eq!(result.complex[1].0, 0.0);
            assert_eq!(result.complex[1].1, -2.0);
        }
        #[test]
        fn test_complex_roots_case2() {
            let equation: &str = "X^3 + 6X^2 + 9X + 10 = 0";
            let poly: Polynomial = parse_equation(equation).expect("Failed to parse equation");
            let result: Roots = degree_3(poly).unwrap();
            assert_eq!(result.real.len(), 1);
            assert_eq!(result.complex.len(), 2);
            assert_eq!(result.real[0], -4.49203);

            assert_eq!(result.complex[0].0, -0.75398);
            assert_eq!(result.complex[0].1, 1.28751);
            assert_eq!(result.complex[1].0, -0.75398);
            assert_eq!(result.complex[1].1, -1.28751);
        }

        #[test]
        fn test_complex_roots_case3() {
            let equation: &str = "X^3 - 3X + 2 = 0";
            let poly: Polynomial = parse_equation(equation).expect("Failed to parse equation");
            let result: Roots = degree_3(poly).unwrap();
            assert_eq!(result.real.len(), 2);
            assert_eq!(result.complex.len(), 0);

            assert_eq!(result.real[0], -2.0);
            assert_eq!(result.real[1], 1.0);
        }

        #[test]
        fn test_simple_complex_roots() {
            let equation: &str = "X^3 + X = 0";
            let poly: Polynomial = parse_equation(equation).expect("Failed to parse equation");
            let result: Roots = degree_3(poly).unwrap();
            assert_eq!(result.real.len(), 1);
            assert_eq!(result.complex.len(), 2);
            assert_eq!(result.real[0], 0.0);

            assert_eq!(result.complex[0].0, 0.0);
            assert_eq!(result.complex[0].1, 1.0);
            assert_eq!(result.complex[1].0, 0.0);
            assert_eq!(result.complex[1].1, -1.0);
        }

        #[test]
        fn test_one_real_two_complex_roots() {
            let equation: &str = "X^3 - 2X^2 + 2X = 0";
            let poly: Polynomial = parse_equation(equation).expect("Failed to parse equation");
            let result: Roots = degree_3(poly).unwrap();
            assert_eq!(result.real.len(), 1);
            assert_eq!(result.complex.len(), 2);
            assert_eq!(result.real[0], 0.0);
            assert_eq!(result.complex[0], (1.0, 1.0));
            assert_eq!(result.complex[1], (1.0, -1.0));
        }

        #[test]
        fn test_complex_roots_with_constant() {
            let equation: &str = "X^3 + 8 = 0";
            let poly: Polynomial = parse_equation(equation).expect("Failed to parse equation");
            let result: Roots = degree_3(poly).unwrap();
            assert_eq!(result.real.len(), 1);
            assert_eq!(result.complex.len(), 2);
            assert_eq!(result.real[0], -2.0);
            assert_eq!(result.complex[0], (1.0, 1.73205));
            assert_eq!(result.complex[1], (1.0, -1.73205));
        }

        #[test]
        fn test_complex_roots_with_quadratic_term() {
            let equation: &str = "X^3 + 3X^2 + 4X + 2 = 0";
            let poly: Polynomial = parse_equation(equation).expect("Failed to parse equation");
            let result: Roots = degree_3(poly).unwrap();
            assert_eq!(result.real.len(), 1);
            assert_eq!(result.complex.len(), 2);
            assert_eq!(result.real[0], -1.0);
            assert_eq!(result.complex[0], (-1.0, 1.0));
            assert_eq!(result.complex[1], (-1.0, -1.0));
        }

        #[test]
        fn test_complex_roots_with_all_terms() {
            let equation: &str = "2X^3 - 4X^2 + 4X - 2 = 0";
            let poly: Polynomial = parse_equation(equation).expect("Failed to parse equation");
            let result: Roots = degree_3(poly).unwrap();
            assert_eq!(result.real.len(), 1);
            assert_eq!(result.complex.len(), 2);
            assert_eq!(result.real[0], 1.0);
            assert_eq!(result.complex[0], (0.5, 0.86603));
            assert_eq!(result.complex[1], (0.5, -0.86603));
        }
    }
}