// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later
//
// example_rationalfunction.cpp
// Rational functions P(x)/Q(x): construction with auto-reduction, arithmetic,
// differentiation, composition, and partial-fraction decomposition.

#include <math/core/RationalFunction.hpp>
#include <math/core/RationalFunction_partialFractions.hpp>
#include <iomanip>
#include <iostream>

using namespace sangi;

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

    std::cout << "=== construction & auto-reduction ===\n";
    // (x^2 - 1) / (x - 1) cancels (x-1) -> x + 1
    RationalFunction<double> f({-1, 0, 1}, {-1, 1});
    std::cout << "  (x^2-1)/(x-1) = " << f.toString() << "\n";
    std::cout << "  f(3)          = " << f(3.0) << "\n\n";

    std::cout << "=== arithmetic (field) ===\n";
    RationalFunction<double> a({1}, {1, 1});      // 1/(1+x)
    RationalFunction<double> b({0, 1}, {1, 1});   // x/(1+x)
    std::cout << "  a = 1/(1+x), b = x/(1+x)\n";
    std::cout << "  a + b = " << (a + b).toString() << "\n";   // 1
    std::cout << "  a * b = " << (a * b).toString() << "\n";
    std::cout << "  a / b = " << (a / b).toString() << "\n\n";

    std::cout << "=== calculus ===\n";
    std::cout << "  d/dx a = " << a.derivative().toString() << "\n"; // -1/(1+x)^2
    std::cout << "  1 / a  = " << a.reciprocal().toString() << "\n"; // 1+x
    std::cout << "  a^2    = " << a.pow(2).toString() << "\n\n";

    std::cout << "=== composition ===\n";
    auto h = compose(a, b);  // a(b(x)) = (1+x)/(1+2x)
    std::cout << "  a(b(x)) = " << h.toString() << "\n";
    std::cout << "  a(b(1)) = " << h(1.0) << "\n\n";

    std::cout << "=== partial fractions ===\n";
    // 1 / (x^2 - 1) = 0.5/(x-1) - 0.5/(x+1)
    RationalFunction<double> r({1}, {-1, 0, 1});
    auto pfe = partialFractions(r);
    if (pfe.success) {
        std::cout << "  1/(x^2-1) = " << toString(pfe) << "\n";
        std::cout << "  #terms    = " << pfe.terms.size() << "\n";
    }
    auto ps = poles(r);
    std::cout << "  poles     = {";
    for (std::size_t i = 0; i < ps.size(); ++i)
        std::cout << (i ? ", " : "") << ps[i].re;
    std::cout << "}\n";

    return 0;
}
