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

// PoissonSeries.hpp — Poisson 級数 (CAS-9)
//
// Poisson 級数は三角関数と多項式の積和で表現される級数。
// 天体力学の摂動展開などで用いられる。
//
// 各項: coeff * cos/sin(n1*θ1 + n2*θ2 + ...)
// 乗算には積和公式 (product-to-sum) を使用。

#ifndef CALX_MATH_CORE_POISSONSERIES_HPP
#define CALX_MATH_CORE_POISSONSERIES_HPP

#include <algorithm>
#include <cmath>
#include <vector>

namespace calx {

    // ================================================================
    // PoissonTerm — Poisson 級数の 1 項
    // ================================================================

    template<typename T>
    struct PoissonTerm {
        T coefficient;
        std::vector<int> frequencies;  // 整数線形結合 [n1, n2, ...]
        bool isSine;                   // true = sin, false = cos

        PoissonTerm() : coefficient{}, isSine(false) {}
        PoissonTerm(T coeff, std::vector<int> freq, bool sine)
            : coefficient(std::move(coeff))
            , frequencies(std::move(freq))
            , isSine(sine) {}
    };

    // ================================================================
    // PoissonSeries — Poisson 級数
    // ================================================================

    template<typename T>
    class PoissonSeries {
        std::vector<PoissonTerm<T>> terms_;

    public:
        PoissonSeries() = default;

        /// 項を追加する
        void addTerm(T coeff, std::vector<int> freq, bool isSine) {
            terms_.emplace_back(std::move(coeff), std::move(freq), isSine);
        }

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

        /// 項へのアクセス
        const PoissonTerm<T>& term(size_t i) const { return terms_[i]; }

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

        PoissonSeries operator+(const PoissonSeries& other) const {
            PoissonSeries result = *this;
            result.terms_.insert(result.terms_.end(),
                                 other.terms_.begin(), other.terms_.end());
            result.simplify();
            return result;
        }

        PoissonSeries& operator+=(const PoissonSeries& other) {
            terms_.insert(terms_.end(), other.terms_.begin(), other.terms_.end());
            simplify();
            return *this;
        }

        // ────────────────────────────────────────
        // スカラー乗算
        // ────────────────────────────────────────

        PoissonSeries operator*(T scalar) const {
            PoissonSeries result;
            result.terms_.reserve(terms_.size());
            for (const auto& t : terms_) {
                result.terms_.emplace_back(t.coefficient * scalar,
                                           t.frequencies, t.isSine);
            }
            return result;
        }

        // ────────────────────────────────────────
        // 級数同士の乗算 (積和公式)
        // ────────────────────────────────────────
        //
        // cos(a)*cos(b) = 0.5*(cos(a-b) + cos(a+b))
        // sin(a)*sin(b) = 0.5*(cos(a-b) - cos(a+b))
        // sin(a)*cos(b) = 0.5*(sin(a+b) + sin(a-b))
        // cos(a)*sin(b) = 0.5*(sin(a+b) - sin(a-b))

        PoissonSeries operator*(const PoissonSeries& other) const {
            PoissonSeries result;
            result.terms_.reserve(terms_.size() * other.terms_.size() * 2);

            for (const auto& a : terms_) {
                for (const auto& b : other.terms_) {
                    T halfCoeff = a.coefficient * b.coefficient * T(0.5);
                    auto sumFreq = addFreq(a.frequencies, b.frequencies);
                    auto difFreq = subFreq(a.frequencies, b.frequencies);

                    if (!a.isSine && !b.isSine) {
                        // cos * cos → 0.5*(cos(a-b) + cos(a+b))
                        result.addTerm(halfCoeff, difFreq, false);
                        result.addTerm(halfCoeff, sumFreq, false);
                    } else if (a.isSine && b.isSine) {
                        // sin * sin → 0.5*(cos(a-b) - cos(a+b))
                        result.addTerm(halfCoeff, difFreq, false);
                        result.addTerm(-halfCoeff, sumFreq, false);
                    } else if (a.isSine && !b.isSine) {
                        // sin * cos → 0.5*(sin(a+b) + sin(a-b))
                        result.addTerm(halfCoeff, sumFreq, true);
                        result.addTerm(halfCoeff, difFreq, true);
                    } else {
                        // cos * sin → 0.5*(sin(a+b) - sin(a-b))
                        result.addTerm(halfCoeff, sumFreq, true);
                        result.addTerm(-halfCoeff, difFreq, true);
                    }
                }
            }

            result.simplify();
            return result;
        }

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

        /// angles[i] = θ_i の値
        T evaluate(const std::vector<T>& angles) const {
            T sum = T(0);
            for (const auto& t : terms_) {
                T phase = T(0);
                size_t n = std::min(t.frequencies.size(), angles.size());
                for (size_t i = 0; i < n; ++i) {
                    phase += T(t.frequencies[i]) * angles[i];
                }
                if (t.isSine)
                    sum += t.coefficient * std::sin(phase);
                else
                    sum += t.coefficient * std::cos(phase);
            }
            return sum;
        }

        // ────────────────────────────────────────
        // 簡約化 (同類項の結合)
        // ────────────────────────────────────────

        void simplify() {
            if (terms_.empty()) return;

            // 同じ (frequencies, isSine) の項をまとめる
            std::vector<PoissonTerm<T>> merged;
            merged.reserve(terms_.size());

            for (auto& t : terms_) {
                bool found = false;
                for (auto& m : merged) {
                    if (m.isSine == t.isSine && m.frequencies == t.frequencies) {
                        m.coefficient += t.coefficient;
                        found = true;
                        break;
                    }
                }
                if (!found) {
                    merged.push_back(std::move(t));
                }
            }

            // 係数がほぼ 0 の項を除去
            terms_.clear();
            for (auto& m : merged) {
                if (std::abs(static_cast<double>(m.coefficient)) > 1e-15) {
                    terms_.push_back(std::move(m));
                }
            }
        }

    private:
        // 周波数ベクトルの加算
        static std::vector<int> addFreq(const std::vector<int>& a,
                                         const std::vector<int>& b) {
            size_t n = std::max(a.size(), b.size());
            std::vector<int> result(n, 0);
            for (size_t i = 0; i < a.size(); ++i) result[i] += a[i];
            for (size_t i = 0; i < b.size(); ++i) result[i] += b[i];
            return result;
        }

        // 周波数ベクトルの減算
        static std::vector<int> subFreq(const std::vector<int>& a,
                                         const std::vector<int>& b) {
            size_t n = std::max(a.size(), b.size());
            std::vector<int> result(n, 0);
            for (size_t i = 0; i < a.size(); ++i) result[i] += a[i];
            for (size_t i = 0; i < b.size(); ++i) result[i] -= b[i];
            return result;
        }
    };

}  // namespace calx

#endif  // CALX_MATH_CORE_POISSONSERIES_HPP
