// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// example_roots.cpp
// Find all roots of polynomials using various algorithms
//
// Low degree (1-4): robust closed-form formulas
// High degree (5+): Jenkins-Traub algorithm
//
// Build:
//   MSBuild or compile with /std:c++23 /I<sangi>/include

#include <math/core/Polynomial.hpp>
#include <math/core/Complex.hpp>
#include <math/roots/polynomial_roots.hpp>
#include <iostream>
#include <cmath>

using namespace sangi;

void print_roots(const std::vector<Complex<double>>& roots) {
    for (size_t i = 0; i < roots.size(); ++i) {
        std::cout << "  x" << i + 1 << " = ";
        if (std::abs(roots[i].im) < 1e-10) {
            std::cout << roots[i].re << std::endl;
        } else {
            std::cout << roots[i].re
                      << (roots[i].im >= 0 ? " + " : " - ")
                      << std::abs(roots[i].im) << "i" << std::endl;
        }
    }
}

int main() {
    std::cout << "=== sangi: Polynomial Root Finding ===" << std::endl;

    // --- Quadratic: x^2 - 5x + 6 = (x-2)(x-3) ---
    std::cout << "\n--- x^2 - 5x + 6 = 0 ---" << std::endl;
    Polynomial<double> p2({6, -5, 1});  // coefficients in ascending order
    auto roots2 = solvePolynomial(p2);
    print_roots(roots2);

    // --- Quadratic with complex roots: x^2 + 1 ---
    std::cout << "\n--- x^2 + 1 = 0 ---" << std::endl;
    Polynomial<double> p2c({1, 0, 1});
    auto roots2c = solvePolynomial(p2c);
    print_roots(roots2c);

    // --- Cubic: x^3 - 6x^2 + 11x - 6 = (x-1)(x-2)(x-3) ---
    std::cout << "\n--- x^3 - 6x^2 + 11x - 6 = 0 ---" << std::endl;
    Polynomial<double> p3({-6, 11, -6, 1});
    auto roots3 = solvePolynomial(p3);
    print_roots(roots3);

    // --- Quartic: x^4 - 10x^2 + 9 = (x-1)(x+1)(x-3)(x+3) ---
    std::cout << "\n--- x^4 - 10x^2 + 9 = 0 ---" << std::endl;
    Polynomial<double> p4({9, 0, -10, 0, 1});
    auto roots4 = solvePolynomial(p4);
    print_roots(roots4);

    // --- Degree 6: x^6 - 1 (6th roots of unity) ---
    std::cout << "\n--- x^6 - 1 = 0 (roots of unity, Jenkins-Traub) ---" << std::endl;
    Polynomial<double> p6({-1, 0, 0, 0, 0, 0, 1});
    auto roots6 = solvePolynomial(p6);
    print_roots(roots6);

    // Verify: each root should satisfy |x^6 - 1| < eps
    std::cout << "  Verification (|p(x)| for each root):" << std::endl;
    for (size_t i = 0; i < roots6.size(); ++i) {
        auto z = roots6[i];
        // Compute z^6
        Complex<double> z6(1.0, 0.0);
        for (int k = 0; k < 6; ++k) {
            double re = z6.re * z.re - z6.im * z.im;
            double im = z6.re * z.im + z6.im * z.re;
            z6 = Complex<double>(re, im);
        }
        double err = std::sqrt((z6.re - 1.0) * (z6.re - 1.0) + z6.im * z6.im);
        std::cout << "    |p(x" << i + 1 << ")| = " << err << std::endl;
    }

    return 0;
}