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

// factorial_utils.hpp
// Factorial-related utilities (header-only templates)
//
// Functions provided:
//   doubleFactorial(n)         — double factorial n!!
//   risingFactorial(x, n)     — rising factorial (Pochhammer symbol) (x)_n
//   fallingFactorial(x, n)    — falling factorial x^{(n)}
//   binomialCoefficient(n, k) — binomial coefficient C(n,k)
//
// Supported types:
//   integer types               — doubleFactorial, binomialCoefficient
//   float, double, long double  — risingFactorial, fallingFactorial

#ifndef SANGI_SPECIAL_FACTORIAL_UTILS_HPP
#define SANGI_SPECIAL_FACTORIAL_UTILS_HPP

#include <math/special/special_common.hpp>
#include <concepts>
#include <cmath>
#include <limits>
#include <stdexcept>

namespace sangi {
namespace special {

// ================================================================
// Double factorial n!! = n·(n-2)·(n-4)·...
// ================================================================

/// Integer version: 0!! = 1, 1!! = 1, (-1)!! = 1
[[nodiscard]] inline long long doubleFactorial(int n) {
    if (n <= 1) return 1;
    long long result = 1;
    for (int k = n; k >= 2; k -= 2) {
        result *= k;
    }
    return result;
}

// ================================================================
// Rising factorial (Pochhammer symbol) (x)_n = x·(x+1)·...·(x+n-1)
// ================================================================

/// Native floating-point types
template<IsNativeFloat T>
[[nodiscard]] T risingFactorial(T x, int n) {
    if (n == 0) return T(1);
    if (n < 0) {
        // (x)_{-n} is rarely used, but by definition it is Γ(x)/Γ(x+n)
        return std::numeric_limits<T>::quiet_NaN();
    }
    T result = T(1);
    for (int k = 0; k < n; k++) {
        result *= (x + static_cast<T>(k));
    }
    return result;
}

// Float type: to be added later as needed

// ================================================================
// Falling factorial x^{(n)} = x·(x-1)·...·(x-n+1)
// ================================================================

/// Native floating-point types
template<IsNativeFloat T>
[[nodiscard]] T fallingFactorial(T x, int n) {
    if (n == 0) return T(1);
    if (n < 0) {
        return std::numeric_limits<T>::quiet_NaN();
    }
    T result = T(1);
    for (int k = 0; k < n; k++) {
        result *= (x - static_cast<T>(k));
    }
    return result;
}

// Float type: to be added later as needed

// ================================================================
// Binomial coefficient C(n, k) = n! / (k! · (n-k)!)
// ================================================================

/// Integer version (watch for overflow: n ≤ 62 or so is safe)
[[nodiscard]] inline long long binomialCoefficient(int n, int k) {
    if (k < 0 || k > n) return 0;
    if (k == 0 || k == n) return 1;
    if (k > n - k) k = n - k;  // exploit symmetry
    long long result = 1;
    for (int i = 0; i < k; i++) {
        result = result * (n - i) / (i + 1);
    }
    return result;
}

/// Floating-point version (handles large n, and real n too)
template<IsNativeFloat T>
[[nodiscard]] T binomialCoefficientReal(T n, int k) {
    if (k < 0) return T(0);
    if (k == 0) return T(1);
    T result = T(1);
    for (int i = 0; i < k; i++) {
        result *= (n - static_cast<T>(i)) / static_cast<T>(i + 1);
    }
    return result;
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_FACTORIAL_UTILS_HPP