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

// ContinuedFraction.hpp
//
// 連分数ツール (Continued Fraction Tools)
//
// 連分数の評価・展開・各種変換を提供する。
// 特殊関数の高精度計算の内部ツールとしても使用。
//
// 主な API:
//   evalContinuedFraction      — 一般連分数の評価 (修正 Lentz 法)
//   toContinuedFraction        — 実数 → 単純連分数展開 [a0; a1, a2, ...]
//   convergents                — 収束子列 p_n/q_n の計算
//   bestRationalApprox         — 最良有理近似 (連分数経由)
//   fromContinuedFraction      — 単純連分数 → 実数値
//   centredContinuedFraction   — 中心連分数展開 (最近整数)
//   lurothExpansion            — Lüroth 展開
//   engelExpansion             — Engel 展開

#ifndef CALX_CONTINUED_FRACTION_HPP
#define CALX_CONTINUED_FRACTION_HPP

#include <math/core/common.hpp>
#include <vector>
#include <cmath>
#include <stdexcept>
#include <limits>
#include <utility>
#include <cstdint>
#include <functional>

namespace calx {

// =========================================================================
// 一般連分数の評価 (修正 Lentz 法)
// =========================================================================

/**
 * @brief 一般連分数の評価 (Modified Lentz Algorithm)
 *
 * f = b_0 + a_1 / (b_1 + a_2 / (b_2 + ...))
 *
 * Lentz-Thompson-Barnett の修正アルゴリズム:
 *   C_0 = f_0 = b_0 (b_0=0 なら tiny)
 *   D_0 = 0
 *   C_j = b_j + a_j / C_{j-1}
 *   D_j = b_j + a_j · D_{j-1}  → D_j = 1/D_j
 *   Δ_j = C_j · D_j
 *   f_j = f_{j-1} · Δ_j
 *
 * @tparam T 浮動小数点型
 * @param a_func  a_j を返す関数 (j = 1, 2, ...)
 * @param b_func  b_j を返す関数 (j = 0, 1, 2, ...)
 * @param maxIter 最大反復数
 * @param tol     収束判定の相対誤差
 * @return 連分数の値
 */
template<typename T>
[[nodiscard]]
T evalContinuedFraction(
    std::function<T(std::size_t)> a_func,
    std::function<T(std::size_t)> b_func,
    std::size_t maxIter = 500,
    T tol = std::numeric_limits<T>::epsilon() * T(10))
{
    const T tiny = std::numeric_limits<T>::min() * T(1e4);

    T f = b_func(0);
    if (std::abs(f) < tiny) f = tiny;

    T C = f;
    T D = T(0);

    for (std::size_t j = 1; j <= maxIter; ++j) {
        T a = a_func(j);
        T b = b_func(j);

        D = b + a * D;
        if (std::abs(D) < tiny) D = tiny;
        D = T(1) / D;

        C = b + a / C;
        if (std::abs(C) < tiny) C = tiny;

        T delta = C * D;
        f *= delta;

        if (std::abs(delta - T(1)) < tol) break;
    }

    return f;
}

// =========================================================================
// 単純連分数展開 (Simple Continued Fraction Expansion)
// =========================================================================

/**
 * @brief 実数を単純連分数 [a0; a1, a2, ...] に展開する
 *
 * x = a0 + 1/(a1 + 1/(a2 + ...))
 * 各 a_i は整数 (a0 は任意、a1 以降は正整数)。
 *
 * @param x     展開する実数
 * @param maxTerms 最大項数
 * @param tol   残差の打ち切り閾値
 * @return 部分商の列 [a0, a1, a2, ...]
 */
template<typename T>
[[nodiscard]]
std::vector<int64_t> toContinuedFraction(T x, std::size_t maxTerms = 50,
                                          T tol = T(1e-12))
{
    std::vector<int64_t> cf;

    for (std::size_t i = 0; i < maxTerms; ++i) {
        int64_t a = static_cast<int64_t>(std::floor(x));
        cf.push_back(a);
        T frac = x - T(a);
        if (std::abs(frac) < tol) break;
        if (std::abs(frac) < std::numeric_limits<T>::epsilon() * std::abs(x) * T(100))
            break;
        x = T(1) / frac;
    }

    return cf;
}

// =========================================================================
// 収束子 (Convergents)
// =========================================================================

/**
 * @brief 収束子 p_n/q_n の列を計算する
 *
 * 連分数 [a0; a1, a2, ...] の第 n 収束子は:
 *   p_{-1}=1, p_0=a_0, p_n = a_n·p_{n-1} + p_{n-2}
 *   q_{-1}=0, q_0=1,   q_n = a_n·q_{n-1} + q_{n-2}
 *
 * @param cf 部分商の列
 * @return (p_n, q_n) のペア列
 */
[[nodiscard]]
inline std::vector<std::pair<int64_t, int64_t>>
convergents(const std::vector<int64_t>& cf)
{
    if (cf.empty()) return {};

    std::vector<std::pair<int64_t, int64_t>> conv;

    int64_t p_prev2 = 1, q_prev2 = 0;  // p_{-1}, q_{-1}
    int64_t p_prev1 = cf[0], q_prev1 = 1;  // p_0, q_0

    conv.emplace_back(p_prev1, q_prev1);

    for (std::size_t n = 1; n < cf.size(); ++n) {
        int64_t a = cf[n];
        int64_t p = a * p_prev1 + p_prev2;
        int64_t q = a * q_prev1 + q_prev2;
        conv.emplace_back(p, q);

        p_prev2 = p_prev1;
        q_prev2 = q_prev1;
        p_prev1 = p;
        q_prev1 = q;
    }

    return conv;
}

// =========================================================================
// 最良有理近似 (Best Rational Approximation)
// =========================================================================

/**
 * @brief 分母が maxDenom 以下の最良有理近似を返す
 *
 * 連分数の収束子は最良有理近似 (best rational approximation) を与える。
 * |x - p/q| < 1/(q·q_{n+1}) を満たす唯一の分数。
 *
 * @param x         近似する実数
 * @param maxDenom  分母の上限
 * @return (p, q) 最良近似分数
 */
template<typename T>
[[nodiscard]]
std::pair<int64_t, int64_t> bestRationalApprox(T x, int64_t maxDenom = 1000000)
{
    auto cf = toContinuedFraction(x, 80);
    auto conv = convergents(cf);

    // 最後の収束子で q ≤ maxDenom を満たすものを返す
    std::pair<int64_t, int64_t> best = {0, 1};
    for (auto& [p, q] : conv) {
        if (std::abs(q) <= maxDenom) {
            best = {p, q};
        } else {
            break;
        }
    }
    return best;
}

// =========================================================================
// 連分数 → 実数値 (fromContinuedFraction)
// =========================================================================

/**
 * @brief 単純連分数 [a0; a1, a2, ...] から実数値を計算する
 *
 * 後ろから畳み込み: x = a_{n-1}, x = a_{n-2} + 1/x, ...
 *
 * @param cf 部分商の列
 * @return 連分数の値
 */
template<typename T>
[[nodiscard]]
T fromContinuedFraction(const std::vector<int64_t>& cf)
{
    if (cf.empty()) return T(0);
    if (cf.size() == 1) return T(cf[0]);

    T val = T(cf.back());
    for (std::size_t i = cf.size() - 1; i > 0; --i) {
        val = T(cf[i - 1]) + T(1) / val;
    }
    return val;
}

// =========================================================================
// 中心連分数展開 (Centred / Nearest-Integer Continued Fraction)
// =========================================================================

/**
 * @brief 中心連分数展開 (最近整数型)
 *
 * 通常の連分数では a_n = floor(x_n) だが、
 * 中心連分数では a_n = round(x_n) (最近整数) を使う。
 * 結果として |x_n - a_n| ≤ 1/2 となり、残差が小さく収束が速い。
 * 部分商は正にも負にもなりうる (|a_n| ≥ 2 for n ≥ 1)。
 *
 * @param x     展開する実数
 * @param maxTerms 最大項数
 * @param tol   打ち切り閾値
 * @return 部分商の列
 */
template<typename T>
[[nodiscard]]
std::vector<int64_t> centredContinuedFraction(T x, std::size_t maxTerms = 50,
                                               T tol = T(1e-12))
{
    std::vector<int64_t> cf;

    for (std::size_t i = 0; i < maxTerms; ++i) {
        int64_t a = static_cast<int64_t>(std::round(x));
        cf.push_back(a);
        T frac = x - T(a);
        if (std::abs(frac) < tol) break;
        if (std::abs(frac) < std::numeric_limits<T>::epsilon() * std::abs(x) * T(100))
            break;
        x = T(1) / frac;
    }

    return cf;
}

// =========================================================================
// Lüroth 展開 (Lüroth Expansion)
// =========================================================================

/**
 * @brief Lüroth 展開
 *
 * x ∈ (0, 1] に対して:
 *   x = 1/d_1 + 1/(d_1(d_1-1)) · (1/d_2 + 1/(d_2(d_2-1)) · (...))
 * ただし d_n ≥ 2 (正整数)。
 *
 * アルゴリズム: d_n = ceil(1/x_n), x_{n+1} = d_n(d_n-1)·x_n - (d_n-1)
 *
 * @param x     展開する実数 (0 < x ≤ 1)
 * @param maxTerms 最大項数
 * @return Lüroth 数列 [d_1, d_2, ...]
 */
template<typename T>
[[nodiscard]]
std::vector<int64_t> lurothExpansion(T x, std::size_t maxTerms = 50)
{
    if (x <= T(0) || x > T(1))
        throw std::invalid_argument("lurothExpansion: x must be in (0, 1]");

    std::vector<int64_t> digits;
    T tol = std::numeric_limits<T>::epsilon() * T(100);

    for (std::size_t i = 0; i < maxTerms; ++i) {
        if (std::abs(x) < tol) break;

        int64_t d = static_cast<int64_t>(std::ceil(T(1) / x));
        if (d < 2) d = 2;
        digits.push_back(d);

        // x_{n+1} = d(d-1)x - (d-1)
        x = T(d) * T(d - 1) * x - T(d - 1);

        // 数値誤差で負になったらクリップ
        if (x < T(0)) x = T(0);
        if (x > T(1)) x = T(1);
    }

    return digits;
}

/**
 * @brief Lüroth 展開から実数値を復元する
 *
 * 後ろから畳み込み:
 *   x = 0, x = (x + d_n - 1) / (d_n * (d_n - 1))  を繰り返す
 *   最後に x = (x + d_1 - 1) / (d_1 * (d_1 - 1)) ... いや
 *   正しくは: 値 = 1/d_1 + 1/(d_1(d_1-1)) · (1/d_2 + ...)
 *   後ろから: val = 0, val = 1/d_n + val/(d_n(d_n-1)) ... いや
 *   val = 0
 *   for i = n-1 downto 0:
 *     val = (1 + val) / (d_i * (d_i - 1)) ... ← 要確認
 *
 * 正確には: x = Σ_{k=1}^{n} 1/(d_1(d_1-1)·d_2(d_2-1)·...·d_{k-1}(d_{k-1}-1)·d_k)
 * 後ろから: val = 1/d_n, val = 1/d_{n-1} + val/(d_{n-1}(d_{n-1}-1))
 */
template<typename T>
[[nodiscard]]
T fromLurothExpansion(const std::vector<int64_t>& digits)
{
    if (digits.empty()) return T(0);

    T val = T(0);
    for (std::size_t i = digits.size(); i > 0; --i) {
        int64_t d = digits[i - 1];
        val = T(1) / T(d) + val / (T(d) * T(d - 1));
    }
    return val;
}

// =========================================================================
// Engel 展開 (Engel Expansion)
// =========================================================================

/**
 * @brief Engel 展開
 *
 * x ∈ (0, ∞) に対して:
 *   x = 1/a_1 + 1/(a_1·a_2) + 1/(a_1·a_2·a_3) + ...
 *       = 1/a_1 · (1 + 1/a_2 · (1 + 1/a_3 · (...)))
 * ただし a_1 ≤ a_2 ≤ a_3 ≤ ... (非減少正整数列)。
 *
 * アルゴリズム: a_n = ceil(1/x_n), x_{n+1} = a_n·x_n - 1
 *
 * @param x     展開する正の実数
 * @param maxTerms 最大項数
 * @return Engel 数列 [a_1, a_2, ...]
 */
template<typename T>
[[nodiscard]]
std::vector<int64_t> engelExpansion(T x, std::size_t maxTerms = 50)
{
    if (x <= T(0))
        throw std::invalid_argument("engelExpansion: x must be positive");

    std::vector<int64_t> digits;
    T tol = std::numeric_limits<T>::epsilon() * T(100);

    for (std::size_t i = 0; i < maxTerms; ++i) {
        if (std::abs(x) < tol) break;

        int64_t a = static_cast<int64_t>(std::ceil(T(1) / x));
        if (a < 1) a = 1;
        digits.push_back(a);

        // x_{n+1} = a·x - 1
        x = T(a) * x - T(1);

        if (x < T(0)) x = T(0);
    }

    return digits;
}

/**
 * @brief Engel 展開から実数値を復元する
 *
 * x = 1/a_1 · (1 + 1/a_2 · (1 + 1/a_3 · (...)))
 * 後ろから: val = 0, val = (1 + val) / a_n, ..., val = (1 + val) / a_1
 */
template<typename T>
[[nodiscard]]
T fromEngelExpansion(const std::vector<int64_t>& digits)
{
    if (digits.empty()) return T(0);

    T val = T(0);
    for (std::size_t i = digits.size(); i > 0; --i) {
        val = (T(1) + val) / T(digits[i - 1]);
    }
    return val;
}

} // namespace calx

#endif // CALX_CONTINUED_FRACTION_HPP
