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

// IntCombinatorics.hpp
// Combinatorial functions (factorial, binomial coefficient, etc.) for Int

#ifndef SANGI_INT_COMBINATORICS_HPP
#define SANGI_INT_COMBINATORICS_HPP

#include "IntBase.hpp"

namespace sangi {

/**
 * @brief Combinatorial functions
 *
 * Provides functions used in combinatorics, such as factorial and binomial
 * coefficient.
 */
class IntCombinatorics {
public:
    /**
     * @brief Factorial n!
     *
     * Computes n! = n * (n-1) * ... * 2 * 1.
     * Based on benchmarks (2026-02-14) the optimal algorithm is selected:
     *   - N <= 100: table lookup (O(1))
     *   - 100 < N <= 228: simple loop (sequential multiplication starting from the table)
     *   - N > 228: sieve method (prime factorization + tournament multiplication)
     * Note: the tournament method is worse than the sieve method across the
     * whole range, so it is not selected automatically.
     * Details: TODO/NOTES_FACTORIAL_BENCHMARK.md
     *
     * @param n Non-negative integer
     * @return The value of n!
     *
     * Special cases:
     *   factorial(0) = 1 (by definition)
     *   factorial(1) = 1
     *   factorial(n < 0) = NaN (factorial of a negative number is undefined)
     *   factorial(NaN) = NaN
     *   factorial(infinity) = infinity
     */
    static Int factorial(const Int& n);

    /**
     * @brief Double factorial n!!
     *
     * Computes n!! = n * (n-2) * (n-4) * ...
     * - n is even: n * (n-2) * ... * 4 * 2
     * - n is odd:  n * (n-2) * ... * 3 * 1
     *
     * @param n Non-negative integer
     * @return The value of n!!
     *
     * Special cases:
     *   doubleFactorial(0) = 1 (by definition)
     *   doubleFactorial(1) = 1
     *   doubleFactorial(-1) = 1 (by definition)
     *   doubleFactorial(n < -1) = NaN
     *   doubleFactorial(NaN) = NaN
     *   doubleFactorial(infinity) = infinity
     *
     * Examples:
     *   6!! = 6 * 4 * 2 = 48
     *   7!! = 7 * 5 * 3 * 1 = 105
     */
    static Int doubleFactorial(const Int& n);

    /**
     * @brief Binomial coefficient C(n, k) = n! / (k! * (n-k)!)
     *
     * Computes the number of combinations "choose k from n".
     *
     * @param n Non-negative integer
     * @param k Non-negative integer (0 <= k <= n)
     * @return The value of C(n, k)
     *
     * Special cases:
     *   binomial(n, 0) = 1 (one way to pick nothing)
     *   binomial(n, n) = 1 (one way to pick everything)
     *   binomial(n, k) = 0 (when k > n)
     *   binomial(n, k < 0) = NaN
     *   binomial(n < 0, k) = NaN
     *   binomial(NaN, k) = NaN
     *
     * Optimizations:
     *   - Use C(n, k) = C(n, n-k) to reduce the amount of computation
     *   - Alternate multiplications and divisions to avoid overflow
     *
     * Examples:
     *   binomial(5, 2) = 5!/(2!*3!) = 10
     *   binomial(10, 3) = 10!/(3!*7!) = 120
     */
    static Int binomial(const Int& n, const Int& k);

    /**
     * @brief Fibonacci number F(n)
     *
     * Computed via fast doubling, using O(log n) multiplications:
     *   F(2k)   = F(k) * (2*F(k+1) - F(k))
     *   F(2k+1) = F(k)^2 + F(k+1)^2
     *
     * @param n Non-negative integer
     * @return F(n)
     *
     * Special cases:
     *   fibonacci(0) = 0, fibonacci(1) = 1
     *   fibonacci(n < 0) = NaN
     */
    static Int fibonacci(const Int& n);

    /**
     * @brief Lucas number L(n)
     *
     * Uses the relation L(n) = F(n-1) + F(n+1) = 2*F(n+1) - F(n).
     *
     * @param n Non-negative integer
     * @return L(n)
     *
     * Special cases:
     *   lucas(0) = 2, lucas(1) = 1
     *   lucas(n < 0) = NaN
     */
    static Int lucas(const Int& n);

    /**
     * @brief Primorial primorial(n) = product of all primes <= n
     *
     * Equivalent to GMP's mpz_primorial_ui.
     * Enumerates primes up to n using the Sieve of Eratosthenes, then
     * multiplies them using a binary-splitting product.
     *
     * @param n Non-negative integer (int)
     * @return Product of all primes <= n. n < 2 -> 1.
     *
     * Example: primorial(10) = 2 * 3 * 5 * 7 = 210
     */
    static Int primorial(int n);

    /**
     * @brief m-step factorial mfac(n, m) = n * (n-m) * (n-2m) * ...
     *
     * Equivalent to GMP's mpz_mfac_uiui.
     * Reduces to n! when m=1, and to n!! (doubleFactorial) when m=2.
     *
     * @param n Non-negative integer
     * @param m Positive integer (step width)
     * @return n * (n-m) * (n-2m) * ... (while each factor is at least 1)
     *
     * Example: mfac(10, 3) = 10 * 7 * 4 * 1 = 280
     */
    static Int mfac(int n, int m);

    // Individual algorithms (exposed for tests and benchmarks)

    /** Table lookup: N <= 100 */
    static const Int& factorialTable(int n);

    /** Simple loop: sequential multiplication starting from the table (100!) */
    static Int factorialSimple(int n);

    /** Tournament method: split out powers of two + balanced edge-x-edge multiplication (for benchmarks; not used by the automatic selection) */
    static Int factorialTournament(int n);

    /** Sieve method: Eratosthenes-style prime factorization + tournament multiplication */
    static Int factorialSieve(int n);

private:
    static constexpr int TABLE_MAX = 100;
    static constexpr int SIMPLE_MAX = 228;
};

} // namespace sangi

#endif // SANGI_INT_COMBINATORICS_HPP