// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later
//
// example_doc_roots_samples.cpp
// Verification harness for samples written in api/Roots.html.

#include <math/roots/root_finding.hpp>
#include <math/roots/polynomial_roots.hpp>
#include <math/core/Polynomial.hpp>
#include <math/core/Complex.hpp>
#include <math/linalg/solvers.hpp>   // needed for Matrix<T>::solve member
#include <cmath>
#include <iomanip>
#include <iostream>

using namespace sangi;

int main() {
    std::cout << std::setprecision(15);

    // ---- sample 1: bisection sqrt(2) -----------------------------
    {
        std::cout << "[bisection sqrt(2)]\n";
        auto r = bisection<double>(
            [](double x) { return x * x - 2.0; }, 1.0, 2.0);
        std::cout << "  converged = " << std::boolalpha << r.converged
                  << "  iter = " << r.iterations
                  << "  root = " << (r.root ? *r.root : 0.0)
                  << "  err = " << r.error_estimate << "\n";
    }

    // ---- sample 2: brent cos(x) = x -------------------------------
    {
        std::cout << "[brent cos(x) - x = 0]\n";
        auto r = brent_method<double>(
            [](double x) { return std::cos(x) - x; }, 0.0, 1.0);
        std::cout << "  converged = " << std::boolalpha << r.converged
                  << "  iter = " << r.iterations
                  << "  root = " << (r.root ? *r.root : 0.0)
                  << "  err = " << r.error_estimate << "\n";
    }

    // ---- sample 4: newton_raphson sqrt(2) -------------------------
    {
        std::cout << "[newton_raphson sqrt(2)]\n";
        auto r = newton_raphson<double>(
            [](double x) { return x * x - 2.0; },
            [](double x) { return 2.0 * x; },
            1.0);
        std::cout << "  converged = " << std::boolalpha << r.converged
                  << "  iter = " << r.iterations
                  << "  root = " << (r.root ? *r.root : 0.0) << "\n";
    }

    // ---- sample 5: secant Dottie ----------------------------------
    {
        std::cout << "[secant cos(x)-x]\n";
        auto r = secant_method<double>(
            [](double x) { return x - std::cos(x); }, 0.0, 1.0);
        std::cout << "  converged = " << std::boolalpha << r.converged
                  << "  iter = " << r.iterations
                  << "  root = " << (r.root ? *r.root : 0.0) << "\n";
    }

    // ---- sample 6: halley cube root of 2 --------------------------
    {
        std::cout << "[halley cube root of 2]\n";
        auto r = halley_method<double>(
            [](double x) { return x * x * x - 2.0; },
            [](double x) { return 3.0 * x * x; },
            [](double x) { return 6.0 * x; },
            1.0);
        std::cout << "  converged = " << std::boolalpha << r.converged
                  << "  iter = " << r.iterations
                  << "  root = " << (r.root ? *r.root : 0.0) << "\n";
    }

    // ---- sample 7: solveQuadratic x^2 - 5x + 6 --------------------
    {
        std::cout << "[solveQuadratic x^2 - 5x + 6 = 0]\n";
        Polynomial<double> p({6.0, -5.0, 1.0});  // ascending
        auto roots = solveQuadratic(p);
        std::cout << "  roots =";
        for (auto& r : roots) std::cout << ' ' << r.re;
        std::cout << "\n";
    }

    // ---- sample 8: jenkinsTraub x^6 - 1 ---------------------------
    {
        std::cout << "[jenkinsTraub x^6 - 1 = 0]\n";
        Polynomial<double> p({-1.0, 0, 0, 0, 0, 0, 1.0});
        auto roots = jenkinsTraub(p);
        double max_residual = 0.0;
        for (auto& r : roots) {
            auto r2 = r * r;
            auto r6 = r2 * r2 * r2;
            auto resid_re = r6.re - 1.0;
            auto resid_im = r6.im;
            double m = std::sqrt(resid_re * resid_re + resid_im * resid_im);
            if (m > max_residual) max_residual = m;
        }
        std::cout << "  n_roots = " << roots.size()
                  << "  max |p(root)| = " << max_residual << "\n";
    }

    // ---- sample 3: newton_raphson_nd, unit circle ∩ y = x ----
    {
        std::cout << "[newton_raphson_nd unit circle ∩ y=x]\n";
        auto F = [](const Vector<double>& v) {
            Vector<double> r(2);
            r[0] = v[0] * v[0] + v[1] * v[1] - 1.0;
            r[1] = v[0] - v[1];
            return r;
        };
        auto J = [](const Vector<double>& v) {
            Matrix<double> m(2, 2);
            m(0, 0) = 2 * v[0]; m(0, 1) = 2 * v[1];
            m(1, 0) = 1.0;      m(1, 1) = -1.0;
            return m;
        };
        Vector<double> x0(2);
        x0[0] = 0.5; x0[1] = 0.5;
        auto r = newton_raphson_nd<double, Vector<double>, Matrix<double>>(F, J, x0);
        std::cout << "  converged = " << std::boolalpha << r.converged
                  << "  iter = " << r.iterations;
        if (r.root) {
            const auto& v = *r.root;
            std::cout << "  root = (" << v[0] << ", " << v[1] << ")";
        }
        std::cout << "\n";
    }
    return 0;
}