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

// LazyPowerSeries.hpp — 遅延評価べき級数 (CAS-8)
//
// 係数をオンデマンドで計算し、キャッシュする形式的べき級数。
// f(x) = Σ a_n * x^n  (a_n はジェネレータ関数で生成)
//
// 加算・乗算 (Cauchy 積)・微分・積分を遅延評価で実装。

#ifndef CALX_MATH_CORE_LAZYPOWERSERIES_HPP
#define CALX_MATH_CORE_LAZYPOWERSERIES_HPP

#include <cmath>
#include <functional>
#include <vector>

namespace calx {

    // ================================================================
    // LazyPowerSeries — 遅延評価べき級数
    // ================================================================

    template<typename T>
    class LazyPowerSeries {
        std::function<T(size_t)> generator_;
        mutable std::vector<T> cache_;

    public:
        /// ジェネレータ関数から構築。generator(n) = a_n
        explicit LazyPowerSeries(std::function<T(size_t)> gen)
            : generator_(std::move(gen)) {}

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

        /// n 次の係数 a_n を返す (遅延計算+キャッシュ)
        T coefficient(size_t n) const {
            if (n < cache_.size()) return cache_[n];
            // 途中の係数もすべてキャッシュする
            cache_.reserve(n + 1);
            for (size_t i = cache_.size(); i <= n; ++i) {
                cache_.push_back(generator_(i));
            }
            return cache_[n];
        }

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

        /// x における級数の値を最大 maxTerms 項で近似
        T evaluate(T x, size_t maxTerms = 50) const {
            T sum = T(0);
            T xpow = T(1);
            for (size_t n = 0; n < maxTerms; ++n) {
                T an = coefficient(n);
                sum += an * xpow;
                xpow *= x;
                // 収束判定: 項の寄与が十分小さければ打ち切り
                if (n > 0 && std::abs(static_cast<double>(an * xpow)) < 1e-16 * (1.0 + std::abs(static_cast<double>(sum)))) {
                    break;
                }
            }
            return sum;
        }

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

        LazyPowerSeries operator+(const LazyPowerSeries& other) const {
            auto genA = generator_;
            auto genB = other.generator_;
            return LazyPowerSeries([genA, genB](size_t n) -> T {
                return genA(n) + genB(n);
            });
        }

        // ────────────────────────────────────────
        // 乗算 (Cauchy 積)
        // ────────────────────────────────────────
        //
        // (f*g)_n = Σ_{k=0}^{n} f_k * g_{n-k}

        LazyPowerSeries operator*(const LazyPowerSeries& other) const {
            // shared_ptr でキャプチャ (コピー可能にするため)
            auto self = std::make_shared<LazyPowerSeries>(*this);
            auto oth = std::make_shared<LazyPowerSeries>(other);
            return LazyPowerSeries([self, oth](size_t n) -> T {
                T sum = T(0);
                for (size_t k = 0; k <= n; ++k) {
                    sum += self->coefficient(k) * oth->coefficient(n - k);
                }
                return sum;
            });
        }

        // ────────────────────────────────────────
        // 微分: (f')_n = (n+1) * f_{n+1}
        // ────────────────────────────────────────

        LazyPowerSeries derivative() const {
            auto gen = generator_;
            return LazyPowerSeries([gen](size_t n) -> T {
                return T(n + 1) * gen(n + 1);
            });
        }

        // ────────────────────────────────────────
        // 積分: (∫f)_n = f_{n-1} / n  (定数 c = 0)
        // ────────────────────────────────────────

        LazyPowerSeries integral() const {
            auto gen = generator_;
            return LazyPowerSeries([gen](size_t n) -> T {
                if (n == 0) return T(0);
                return gen(n - 1) / T(n);
            });
        }

        // ────────────────────────────────────────
        // 定番の級数
        // ────────────────────────────────────────

        /// e^x = Σ x^n / n!
        static LazyPowerSeries exp() {
            return LazyPowerSeries([](size_t n) -> T {
                T fact = T(1);
                for (size_t i = 2; i <= n; ++i)
                    fact *= T(i);
                return T(1) / fact;
            });
        }

        /// sin(x) = Σ (-1)^k x^(2k+1) / (2k+1)!
        static LazyPowerSeries sin() {
            return LazyPowerSeries([](size_t n) -> T {
                if (n % 2 == 0) return T(0);
                size_t k = n / 2;
                T fact = T(1);
                for (size_t i = 2; i <= n; ++i)
                    fact *= T(i);
                return ((k % 2 == 0) ? T(1) : T(-1)) / fact;
            });
        }

        /// cos(x) = Σ (-1)^k x^(2k) / (2k)!
        static LazyPowerSeries cos() {
            return LazyPowerSeries([](size_t n) -> T {
                if (n % 2 != 0) return T(0);
                size_t k = n / 2;
                T fact = T(1);
                for (size_t i = 2; i <= n; ++i)
                    fact *= T(i);
                return ((k % 2 == 0) ? T(1) : T(-1)) / fact;
            });
        }

        /// 1/(1-r*x) = Σ r^n x^n
        static LazyPowerSeries geometric(T r = T(1)) {
            return LazyPowerSeries([r](size_t n) -> T {
                T result = T(1);
                for (size_t i = 0; i < n; ++i)
                    result *= r;
                return result;
            });
        }
    };

}  // namespace calx

#endif  // CALX_MATH_CORE_LAZYPOWERSERIES_HPP
