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

// RationalFunction_partialFractions.hpp
// 有理関数の部分分数展開 (Heaviside の方法)
// lib++20 有理式.h の部分分数Heaviside展開 からの移植 (2026-02-21)
//
// 使用法:
//   #include <math/core/RationalFunction.hpp>
//   #include <math/core/RationalFunction_partialFractions.hpp>

#pragma once

#include <math/core/RationalFunction.hpp>
#include <math/core/Complex.hpp>
#include <math/roots/polynomial_roots.hpp>
#include <vector>

namespace calx {

// ================================================================
// 部分分数展開の結果を保持する構造体
// ================================================================

/// 部分分数展開の各項: coefficient / (x - pole)
template<typename T>
struct PartialFraction {
    Complex<T> coefficient;  // 留数 (residue)
    Complex<T> pole;         // 極
};

/// 部分分数展開の結果
template<typename T>
struct PartialFractionExpansion {
    Polynomial<T> quotient;                     // 整式部分 (deg(num) >= deg(den) の場合)
    std::vector<PartialFraction<T>> terms;      // 部分分数の各項
    bool success = false;                       // 展開に成功したか
};

// ================================================================
// 部分分数展開 (Heaviside の方法)
// ================================================================

/// 有理関数 f = P(x)/Q(x) を部分分数展開する。
/// 前提条件: Q(x) の全根が単根 (重根は未サポート)
///
/// 結果: f(x) = quotient(x) + Σ A_k / (x - r_k)
/// ここで A_k = P(r_k) / Q'(r_k)  (留数公式)
///
/// T は浮動小数点型 (double, float) を想定
template<typename T>
[[nodiscard]] PartialFractionExpansion<T> partialFractions(
    const RationalFunction<T>& f,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000)
{
    static_assert(std::is_floating_point_v<T>,
                  "partialFractions requires floating-point coefficient type");

    PartialFractionExpansion<T> result;

    const auto& P = f.numerator();
    const auto& Q = f.denominator();

    if (Q.isZero()) {
        result.success = false;
        return result;
    }

    // deg(P) >= deg(Q) の場合、多項式除算で整式部分を取り出す
    Polynomial<T> remainder = P;
    if (P.degree() >= Q.degree()) {
        auto dr = P.divmod(Q);
        result.quotient = std::move(dr.quotient);
        remainder = std::move(dr.remainder);
    }

    // 分母が定数なら部分分数は不要
    if (Q.degree() <= 0) {
        result.success = true;
        return result;
    }

    // Q(x) の全根を求める
    auto roots = solvePolynomial(Q, eps, maxIter);

    // Q'(x) を計算
    Polynomial<T> Qprime = Q.derivative();

    // 各根 r_k について留数 A_k = remainder(r_k) / Q'(r_k) を計算
    result.terms.reserve(roots.size());
    for (const auto& r : roots) {
        Complex<T> numVal = remainder(r);
        Complex<T> denVal = Qprime(r);

        if (abs(denVal) < eps) {
            // Q'(r_k) ≈ 0 → 重根の可能性 → 未サポート
            result.success = false;
            return result;
        }

        PartialFraction<T> term;
        term.pole = r;
        term.coefficient = numVal / denVal;
        result.terms.push_back(term);
    }

    result.success = true;
    return result;
}

/// 部分分数展開の結果を文字列に変換
template<typename T>
[[nodiscard]] std::string toString(const PartialFractionExpansion<T>& pfe,
                                    const std::string& var = "x") {
    if (!pfe.success) return "(expansion failed)";

    std::ostringstream oss;
    bool first = true;

    // 整式部分
    if (!pfe.quotient.isZero()) {
        oss << pfe.quotient.toString(var);
        first = false;
    }

    // 各部分分数項
    for (const auto& t : pfe.terms) {
        if (!first) oss << " + ";

        oss << "(" << t.coefficient << ")/(" << var << " - (" << t.pole << "))";
        first = false;
    }

    if (first) oss << "0";
    return oss.str();
}

// ================================================================
// 極の列挙 (分母の根)
// ================================================================

/// 有理関数 f = P(x)/Q(x) の極 (Q(x) = 0 の根) を返す
template<typename T>
[[nodiscard]] std::vector<Complex<T>> poles(
    const RationalFunction<T>& f,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000)
{
    static_assert(std::is_floating_point_v<T>,
                  "poles requires floating-point coefficient type");
    const auto& Q = f.denominator();
    if (Q.degree() <= 0) return {};
    return solvePolynomial(Q, eps, maxIter);
}

// ================================================================
// 零点の列挙 (分子の根)
// ================================================================

/// 有理関数 f = P(x)/Q(x) の零点 (P(x) = 0 の根) を返す
template<typename T>
[[nodiscard]] std::vector<Complex<T>> zeros(
    const RationalFunction<T>& f,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000)
{
    static_assert(std::is_floating_point_v<T>,
                  "zeros requires floating-point coefficient type");
    const auto& P = f.numerator();
    if (P.degree() <= 0) return {};
    return solvePolynomial(P, eps, maxIter);
}

// ================================================================
// 留数 (指定した極での留数を計算)
// ================================================================

/// 有理関数 f = P(x)/Q(x) の単純極 pole における留数
/// Res(f, pole) = P(pole) / Q'(pole)
template<typename T>
[[nodiscard]] Complex<T> residue(
    const RationalFunction<T>& f,
    const Complex<T>& pole)
{
    static_assert(std::is_floating_point_v<T>,
                  "residue requires floating-point coefficient type");
    Complex<T> numVal = f.numerator()(pole);
    Complex<T> denDerVal = f.denominator().derivative()(pole);
    return numVal / denDerVal;
}

} // namespace calx
