// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later
//
// example_doc_polynomial_samples.cpp
// Verification harness for samples in api/Polynomial.html (ja+en).

#include <math/core/Polynomial.hpp>
#include <math/core/Polynomial_factorization.hpp>
#include <math/core/mp/Rational.hpp>
#include <iomanip>
#include <iostream>

using namespace sangi;

template <typename T>
static void show_poly(const char* name, const Polynomial<T>& p) {
    std::cout << "  " << name << " = " << p << "  (deg=" << p.degree() << ")\n";
}

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

    // ---- ctor / accessors ---------------------------------------
    {
        std::cout << "[ctor]\n";
        Polynomial<double> p = {1.0, 2.0, 3.0};
        show_poly("p={1,2,3}", p);
        std::cout << "  leading=" << p.leadingCoefficient()
                  << " constant=" << p.constantTerm()
                  << " degree=" << p.degree() << "\n";
    }

    // ---- composition ---------------------------------------------
    {
        std::cout << "[compose p(q)]\n";
        Polynomial<double> p({1, 0, 1});   // 1 + x^2
        Polynomial<double> q({0, 1, 1});   // x + x^2
        auto r = p(q);                      // p(q(x)) = 1 + (x+x^2)^2
        show_poly("r", r);
    }

    // ---- main example: arithmetic + divmod + eval + d/i + gcd ---
    {
        std::cout << "[main example]\n";
        Polynomial<double> p = {1.0, -3.0, 2.0};   // 1 - 3x + 2x^2
        Polynomial<double> q = {1.0, 1.0};          // 1 + x
        show_poly("p", p);
        show_poly("q", q);
        show_poly("p*q", p * q);

        auto [quot, rem] = p.divmod(q);
        show_poly("p/q quot", quot);
        show_poly("p/q rem", rem);

        std::cout << "  p(2.0) = " << p(2.0) << "  (expected 3)\n";
        show_poly("p'", p.derivative());
        show_poly("integ p", p.integral());

        Polynomial<double> a = {-1.0, 0.0, 1.0};   // x^2 - 1
        Polynomial<double> b = {-1.0, 1.0};         // x - 1
        show_poly("gcd(x^2-1, x-1)", gcd(a, b));
    }

    // ---- discriminant of x^2 - 2 ---------------------------------
    {
        std::cout << "[discriminant]\n";
        Polynomial<int> f = {-2, 0, 1};
        std::cout << "  disc(x^2 - 2) = " << discriminant(f) << "  (expected 8)\n";
    }

    // ---- factorize x^4 - 1 ---------------------------------------
    {
        std::cout << "[factorize x^4-1]\n";
        Polynomial<int> f = {-1, 0, 0, 0, 1};
        auto factors = factorize(f);
        for (std::size_t i = 0; i < factors.size(); ++i) {
            std::cout << "  factor[" << i << "] = " << factors[i]
                      << "^" << factors.exponent(i) << "\n";
        }
    }

    // ---- orthogonal polynomials -----------------------------------
    {
        std::cout << "[orthogonal]\n";
        auto T5 = chebyshevT<double>(5);
        show_poly("T_5", T5);
        auto L3 = legendre<double>(3);
        show_poly("P_3", L3);
    }

    return 0;
}