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

// RationalConstants.cpp
// Implementation of Bernoulli numbers and Stirling coefficients
//
// Improvements:
//   - Made thread-safe with std::mutex (replacing sangi's CxLock)
//   - camelCase naming convention
//   - C++23 style
//   - Bernoulli numbers: switched to the standard recurrence (the Sugi formula is inaccurate for K≥70)

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

namespace sangi {

//==========================================================================
// Bernoulli numbers (standard recurrence)
//==========================================================================
//
// B(n) = -1/(n+1) * Σ_{k=0}^{n-1} C(n+1, k) * B(k)
//
// Compute only even n = 2m (odd > 1 is 0).
// Using B(odd>1) = 0, accumulate only the even terms:
//   B(2m) = -1/(2m+1) * [B(0) + (2m+1)*B(1) + Σ_{j=1}^{m-1} C(2m+1, 2j)*B(2j)]
//
// Thanks to the cache, each B(2m) is computed only once.

Rational bernoulli(int n) {
    // Base cases (no cache needed)
    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);  // odd > 1 is 0

    // Thread-safe cache (filled sequentially)
    // 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];
    }

    // Sequentially compute the uncomputed entries: 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;      // target to compute: 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) — initial value 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)
            // Incremental computation (the division at each step is exact in integers)
            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 coefficients
//==========================================================================
//
// [Sources]
//   "Concerning Two Series for the Gamma Function" p618
//   "A New Derivation of Stirling's Approximation to n!" p828
//
// Recurrence:
//   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);

    // Thread-safe cache
    // c.size() indicates the computed range. If k < c.size(), it is already computed.
    // c[2ki-1] and c[2ki] are computed in pairs.
    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 already computed, return from the cache
    if (k < static_cast<int>(c.size())) {
        return c[k];
    }

    // The next pair index to compute
    // c[0] is already initialized. Pair ki=1 gives c[1],c[2], pair ki=2 gives c[3],c[4]...
    size_t kDone = (c.size() + 1) / 2;

    // Extend up to k+1 or k+2 (use an even size since computation is done in pairs)
    int newSize = (k & 1) ? k + 2 : k + 1;
    c.resize(newSize, Rational(0));

    for (size_t ki = kDone; ; ki++) {
        // Compute c[2*ki-1]
        Rational sum1(0);
        for (int i = 2; i <= static_cast<int>(2 * ki); i += 2) {
            Rational term = bernoulli(i);  // receive the return value directly (move)
            term *= c[2 * ki - i];
            term /= i;
            sum1 += term;
        }
        sum1 /= static_cast<int>(2 * ki - 1);
        c[2 * ki - 1] = std::move(sum1);

        // Compute 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 sangi