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

// RationalConstants.cpp
// Bernoulli 数・Stirling 係数の実装
//
// 移植元: lib++20 Bernoulli数.cpp, Stirling係数.cpp
// 改良点:
//   - std::mutex でスレッドセーフ化 (lib++20 の CxLock 置換)
//   - camelCase 命名規則
//   - C++23 スタイル
//   - Bernoulli 数: 標準再帰式に変更 (Sugi 公式は K≥70 で不正確)

#include <math/core/mp/Rational/RationalConstants.hpp>
#include <vector>
#include <mutex>

namespace calx {

//==========================================================================
// Bernoulli 数 (標準再帰式)
//==========================================================================
//
// B(n) = -1/(n+1) * Σ_{k=0}^{n-1} C(n+1, k) * B(k)
//
// 偶数 n = 2m のみ計算 (奇数 > 1 は 0)。
// B(odd>1) = 0 を利用して偶数項のみ累積:
//   B(2m) = -1/(2m+1) * [B(0) + (2m+1)*B(1) + Σ_{j=1}^{m-1} C(2m+1, 2j)*B(2j)]
//
// キャッシュにより各 B(2m) は一度だけ計算される。

Rational bernoulli(int n) {
    // 基本ケース (キャッシュ不要)
    if (n < 0) return Rational(0);
    if (n == 0) return Rational(1);
    if (n == 1) return Rational(-1, 2);
    if (n & 1) return Rational(0);  // 奇数 > 1 は 0

    // スレッドセーフなキャッシュ (順次充填)
    // cache[i] = B(2*(i+1)) for i >= 0
    // cache[0] = B(2), cache[1] = B(4), ...
    static std::mutex mtx;
    static std::vector<Rational> cache;

    int m = n / 2;     // B(n) = B(2m)
    int idx = m - 1;   // cache index

    std::lock_guard<std::mutex> lock(mtx);

    if (idx < static_cast<int>(cache.size())) {
        return cache[idx];
    }

    // 未計算分を順次計算: B(2*(size+1)) .. B(2m)
    cache.reserve(m);

    for (int mi = static_cast<int>(cache.size()) + 1; mi <= m; mi++) {
        int K = 2 * mi;      // 計算対象: B(K)
        int K1 = K + 1;      // K + 1

        // B(K) = -1/K1 * [C(K1,0)*B(0) + C(K1,1)*B(1) + Σ C(K1,2j)*B(2j)]
        Rational sum(1);                 // C(K1, 0) * B(0) = 1
        sum += Rational(-K1, 2);         // C(K1, 1) * B(1) = K1 * (-1/2)

        Int binom(1);   // C(K1, 2j) — 初期値 C(K1, 0) = 1
        for (int j = 1; j < mi; j++) {
            // C(K1, 2j) = C(K1, 2j-2) * (K1-2j+2)*(K1-2j+1) / ((2j-1)*2j)
            // 逐次計算 (各ステップの除算は整数として正確)
            binom *= (K1 - 2*j + 2);
            binom /= (2*j - 1);           // → C(K1, 2j-1)
            binom *= (K1 - 2*j + 1);
            binom /= (2*j);               // → C(K1, 2j)

            sum += Rational(binom) * cache[j - 1];
        }

        sum /= -K1;
        cache.push_back(std::move(sum));
    }

    return cache[idx];
}

//==========================================================================
// Stirling 係数
//==========================================================================
//
// 【出典】
//   "Concerning Two Series for the Gamma Function" p618
//   "A New Derivation of Stirling's Approximation to n!" p828
//
// 漸化式:
//   c[0] = 1
//   c[2k-1] = (Σ_{i=2,4,...,2k} B(i) * c[2k-i] / i) / (2k-1)
//   c[2k]   = (Σ_{i=2,4,...,2k} B(i) * c[2k-(i-1)] / i) / (2k)

Rational stirlingCoefficient(int k) {
    if (k < 0) return Rational(0);

    // スレッドセーフなキャッシュ
    // c.size() が計算済みの範囲を示す。k < c.size() なら計算済み。
    // c[2ki-1] と c[2ki] はペアで計算される。
    static std::mutex mtx;
    static std::vector<Rational> c;
    static bool initialized = false;

    std::lock_guard<std::mutex> lock(mtx);

    if (!initialized) {
        c.resize(1);
        c[0] = Rational(1);
        initialized = true;
    }

    // 計算済みならキャッシュから返す
    if (k < static_cast<int>(c.size())) {
        return c[k];
    }

    // 次に計算すべきペアインデックス
    // c[0] は初期化済み。ペア ki=1 で c[1],c[2]、ペア ki=2 で c[3],c[4]...
    size_t kDone = (c.size() + 1) / 2;

    // k+1 か k+2 まで拡張 (ペアで計算するため偶数サイズにする)
    int newSize = (k & 1) ? k + 2 : k + 1;
    c.resize(newSize, Rational(0));

    for (size_t ki = kDone; ; ki++) {
        // c[2*ki-1] の計算
        Rational sum1(0);
        for (int i = 2; i <= static_cast<int>(2 * ki); i += 2) {
            Rational term = bernoulli(i);  // 戻り値を直接受け取り (move)
            term *= c[2 * ki - i];
            term /= i;
            sum1 += term;
        }
        sum1 /= static_cast<int>(2 * ki - 1);
        c[2 * ki - 1] = std::move(sum1);

        // c[2*ki] の計算
        Rational sum2(0);
        for (int i = 2; i <= static_cast<int>(2 * ki); i += 2) {
            Rational term = bernoulli(i);
            term *= c[2 * ki - (i - 1)];
            term /= i;
            sum2 += term;
        }
        sum2 /= static_cast<int>(2 * ki);
        c[2 * ki] = std::move(sum2);

        if (static_cast<int>(2 * ki) >= k) break;
    }

    return c[k];
}

} // namespace calx