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

// PuiseuxSeries.hpp — Puiseux 級数 (CAS-7)
//
// Puiseux 級数は分数べきを許す形式的べき級数。
// f(x) = Σ c_k * x^(k/d)  (d = 分岐次数, k は整数)
//
// 代数曲線の特異点解析、Newton-Puiseux アルゴリズムで使用。

#ifndef CALX_MATH_CORE_PUISEUXSERIES_HPP
#define CALX_MATH_CORE_PUISEUXSERIES_HPP

#include <cmath>
#include <map>

namespace calx {

    // ================================================================
    // PuiseuxSeries — 分数べき級数
    // ================================================================

    template<typename T>
    class PuiseuxSeries {
        std::map<int, T> coeffs_;  // 指数の分子 → 係数 (x^(key/denom_) の係数)
        int denom_;                // 共通分母 (分岐次数)
        int order_;                // 打ち切り次数 (分子単位)

    public:
        /// コンストラクタ
        /// @param denom  分岐次数 (共通分母, 既定: 1 = 通常のべき級数)
        /// @param order  打ち切り次数 (order/denom まで保持, 既定: 10)
        explicit PuiseuxSeries(int denom = 1, int order = 10)
            : denom_(denom > 0 ? denom : 1)
            , order_(order * denom_) {}

        // ────────────────────────────────────────
        // 係数アクセス
        // ────────────────────────────────────────

        /// x^(expNum/denom) の係数を参照 (読み書き可)
        T& operator[](int expNum) {
            return coeffs_[expNum];
        }

        /// x^(expNum/denom) の係数を取得 (存在しなければ 0)
        T coefficient(int expNum) const {
            auto it = coeffs_.find(expNum);
            return (it != coeffs_.end()) ? it->second : T(0);
        }

        /// 分岐次数を返す
        int ramification() const { return denom_; }

        /// 打ち切り次数 (分子単位) を返す
        int truncationOrder() const { return order_; }

        /// 非零項数を返す
        size_t numTerms() const { return coeffs_.size(); }

        // ────────────────────────────────────────
        // 加算
        // ────────────────────────────────────────

        PuiseuxSeries operator+(const PuiseuxSeries& other) const {
            // 共通分母を LCM に統一
            int d = lcm(denom_, other.denom_);
            int scaleA = d / denom_;
            int scaleB = d / other.denom_;
            int newOrder = std::min(order_ * scaleA, other.order_ * scaleB);

            PuiseuxSeries result(1, 1);
            result.denom_ = d;
            result.order_ = newOrder;

            for (const auto& [k, c] : coeffs_) {
                int newK = k * scaleA;
                if (newK <= newOrder)
                    result.coeffs_[newK] += c;
            }
            for (const auto& [k, c] : other.coeffs_) {
                int newK = k * scaleB;
                if (newK <= newOrder)
                    result.coeffs_[newK] += c;
            }

            result.removeZeros();
            return result;
        }

        PuiseuxSeries& operator+=(const PuiseuxSeries& other) {
            *this = *this + other;
            return *this;
        }

        // ────────────────────────────────────────
        // 乗算 (畳み込み)
        // ────────────────────────────────────────

        PuiseuxSeries operator*(const PuiseuxSeries& other) const {
            int d = lcm(denom_, other.denom_);
            int scaleA = d / denom_;
            int scaleB = d / other.denom_;
            int newOrder = std::min(order_ * scaleA, other.order_ * scaleB);

            PuiseuxSeries result(1, 1);
            result.denom_ = d;
            result.order_ = newOrder;

            for (const auto& [ka, ca] : coeffs_) {
                int ka2 = ka * scaleA;
                for (const auto& [kb, cb] : other.coeffs_) {
                    int kb2 = kb * scaleB;
                    int ksum = ka2 + kb2;
                    if (ksum <= newOrder) {
                        result.coeffs_[ksum] += ca * cb;
                    }
                }
            }

            result.removeZeros();
            return result;
        }

        PuiseuxSeries& operator*=(const PuiseuxSeries& other) {
            *this = *this * other;
            return *this;
        }

        // ────────────────────────────────────────
        // 数値評価
        // ────────────────────────────────────────

        /// x における級数の値を計算
        T evaluate(T x) const {
            T sum = T(0);
            for (const auto& [k, c] : coeffs_) {
                double exponent = static_cast<double>(k) / static_cast<double>(denom_);
                T xpow = static_cast<T>(std::pow(static_cast<double>(x), exponent));
                sum += c * xpow;
            }
            return sum;
        }

    private:
        // ほぼ 0 の係数を除去
        void removeZeros() {
            for (auto it = coeffs_.begin(); it != coeffs_.end(); ) {
                if (std::abs(static_cast<double>(it->second)) < 1e-15)
                    it = coeffs_.erase(it);
                else
                    ++it;
            }
        }

        // 最小公倍数
        static int lcm(int a, int b) {
            return a / gcd(a, b) * b;
        }

        // 最大公約数
        static int gcd(int a, int b) {
            a = std::abs(a);
            b = std::abs(b);
            while (b) { int t = b; b = a % b; a = t; }
            return a;
        }
    };

}  // namespace calx

#endif  // CALX_MATH_CORE_PUISEUXSERIES_HPP
