1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
#pragma once

#include <cassert><--- Include file:  not found. Please note: Cppcheck does not need standard library headers to get proper results.
#include <valarray><--- Include file:  not found. Please note: Cppcheck does not need standard library headers to get proper results.

#include "snippet/aliases.hpp"

#include "internal/types.hpp"

#include "numeric/modular/modint.hpp"


namespace uni {


namespace internal {


template<class T>
struct factorial_base {
    using value_type = T;
    using size_type = internal::size_t;

    static constexpr value_type calc(const size_type& n) noexcept(NO_EXCEPT) {
        assert(n >= 0);<--- Unsigned positive
        value_type ans = 1;
        FOR(k, 1, n) ans *= k;
        return ans;
    }
};

} // namespace internal



template<class>
struct factorial {};


template<std::integral T>
struct factorial<T> : internal::factorial_base<T> {
    using value_type = T;
    using size_type = internal::size_t;

  private:
    size_type _n;
    std::valarray<value_type> _fact;

  public:
    factorial(const size_type n) noexcept(NO_EXCEPT) : _n(n), _fact(n + 1) {
        this->_fact[0] = this->_fact[1] = 1;
        FOR(i, 2, n) this->_fact[i] = this->_fact[i - 1] * i;
    }

    inline value_type fact(const size_type k) noexcept(NO_EXCEPT) {
        assert(0 <= k and k <= this->_n);<--- Unsigned positive
        return this->_fact[k];
    }
    inline value_type operator()(const size_type k) noexcept(NO_EXCEPT) {
        return this->fact(k);
    }

    auto _debug() const { return this->_fact; }
};


template<uni::internal::modint_family T>
struct factorial<T> : internal::factorial_base<T> {
    using value_type = T;
    using size_type = internal::size_t;

  private:
    size_type _n;
    std::valarray<value_type> _fact, _ifact, _inv;

  public:
    factorial(const size_type n) noexcept(NO_EXCEPT) : _n(n), _fact(n + 1), _ifact(n + 1), _inv(n + 1) {
        constexpr auto P = T::mod();

        this->_fact[0] = this->_fact[1] = 1;
        this->_ifact[0] = this->_ifact[1] = 1;
        this->_inv[1] = 1;

        FOR(i, 2, n) {
            this->_inv[i] = -this->_inv[P % i] * (P / i);

            this->_fact[i] = this->_fact[i - 1] * i;
            this->_ifact[i] = this->_ifact[i - 1] * this->_inv[i];
        }
    }

    inline value_type fact(const size_type k) noexcept(NO_EXCEPT) {
        assert(0 <= k and k <= this->_n);<--- Unsigned positive
        return this->_fact[k];
    }
    inline value_type operator()(const size_type k) noexcept(NO_EXCEPT) {
        return this->fact(k);
    }

    inline value_type ifact(const size_type k) noexcept(NO_EXCEPT) {
        assert(0 <= k and k <= this->_n);<--- Unsigned positive
        return this->_ifact[k];
    }

    inline value_type inv(const size_type k) noexcept(NO_EXCEPT) {
        assert(0 <= k and k <= this->_n);<--- Unsigned positive
        return this->_inv[k];
    }

    inline value_type bimom(const size_type n, const size_type r) noexcept(NO_EXCEPT) {
        assert(0 <= r and r <= n and n <= this->_n);<--- Unsigned positive
        return this->fact(n) * this->ifact(r) * this->ifact(n - r);
    }

    auto _debug() const { return this->_fact; }
};


} // namespace uni