cp-library-cpp

This documentation is automatically generated by online-judge-tools/verification-helper

View the Project on GitHub suisen-cp/cp-library-cpp

:heavy_check_mark: test/src/sequence/stirling_number2_small_prime_mod/stirling_number_of_the_second_kind_small_p_large_n.test.cpp

Depends on

Code

#define PROBLEM "https://judge.yosupo.jp/problem/stirling_number_of_the_second_kind_small_p_large_n"

#include <iostream>
#include <atcoder/modint>

using mint = atcoder::modint;

#include "library/sequence/stirling_number2_small_prime_mod.hpp"

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int t, p;
    std::cin >> t >> p;
    mint::set_mod(p);

    suisen::StirlingNumber2SmallPrimeMod<mint> s2;

    while (t --> 0) {
        long long n, k;
        std::cin >> n >> k;
        std::cout << s2(n, k).val() << '\n';
    }

    return 0;
}
#line 1 "test/src/sequence/stirling_number2_small_prime_mod/stirling_number_of_the_second_kind_small_p_large_n.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/stirling_number_of_the_second_kind_small_p_large_n"

#include <iostream>
#include <atcoder/modint>

using mint = atcoder::modint;

#line 1 "library/sequence/stirling_number2_small_prime_mod.hpp"



#line 1 "library/sequence/stirling_number2.hpp"



#line 1 "library/math/factorial.hpp"



#include <cassert>
#include <vector>

namespace suisen {
    template <typename T, typename U = T>
    struct factorial {
        factorial() = default;
        factorial(int n) { ensure(n); }

        static void ensure(const int n) {
            int sz = _fac.size();
            if (n + 1 <= sz) return;
            int new_size = std::max(n + 1, sz * 2);
            _fac.resize(new_size), _fac_inv.resize(new_size);
            for (int i = sz; i < new_size; ++i) _fac[i] = _fac[i - 1] * i;
            _fac_inv[new_size - 1] = U(1) / _fac[new_size - 1];
            for (int i = new_size - 1; i > sz; --i) _fac_inv[i - 1] = _fac_inv[i] * i;
        }

        T fac(const int i) {
            ensure(i);
            return _fac[i];
        }
        T operator()(int i) {
            return fac(i);
        }
        U fac_inv(const int i) {
            ensure(i);
            return _fac_inv[i];
        }
        U binom(const int n, const int r) {
            if (n < 0 or r < 0 or n < r) return 0;
            ensure(n);
            return _fac[n] * _fac_inv[r] * _fac_inv[n - r];
        }
        template <typename ...Ds, std::enable_if_t<std::conjunction_v<std::is_integral<Ds>...>, std::nullptr_t> = nullptr>
        U polynom(const int n, const Ds& ...ds) {
            if (n < 0) return 0;
            ensure(n);
            int sumd = 0;
            U res = _fac[n];
            for (int d : { ds... }) {
                if (d < 0 or d > n) return 0;
                sumd += d;
                res *= _fac_inv[d];
            }
            if (sumd > n) return 0;
            res *= _fac_inv[n - sumd];
            return res;
        }
        U perm(const int n, const int r) {
            if (n < 0 or r < 0 or n < r) return 0;
            ensure(n);
            return _fac[n] * _fac_inv[n - r];
        }
    private:
        static std::vector<T> _fac;
        static std::vector<U> _fac_inv;
    };
    template <typename T, typename U>
    std::vector<T> factorial<T, U>::_fac{ 1 };
    template <typename T, typename U>
    std::vector<U> factorial<T, U>::_fac_inv{ 1 };
} // namespace suisen


#line 1 "library/sequence/powers.hpp"



#include <cstdint>
#line 1 "library/number/linear_sieve.hpp"



#line 5 "library/number/linear_sieve.hpp"
#include <numeric>
#line 7 "library/number/linear_sieve.hpp"

namespace suisen {
// referece: https://37zigen.com/linear-sieve/
class LinearSieve {
    public:
        LinearSieve(const int n) : _n(n), min_prime_factor(std::vector<int>(n + 1)) {
            std::iota(min_prime_factor.begin(), min_prime_factor.end(), 0);
            prime_list.reserve(_n / 20);
            for (int d = 2; d <= _n; ++d) {
                if (min_prime_factor[d] == d) prime_list.push_back(d);
                const int prime_max = std::min(min_prime_factor[d], _n / d);
                for (int prime : prime_list) {
                    if (prime > prime_max) break;
                    min_prime_factor[prime * d] = prime;
                }
            }
        }
        int prime_num() const noexcept { return prime_list.size(); }
        /**
         * Returns a vector of primes in [0, n].
         * It is guaranteed that the returned vector is sorted in ascending order.
         */
        const std::vector<int>& get_prime_list() const noexcept  {
            return prime_list;
        }
        const std::vector<int>& get_min_prime_factor() const noexcept { return min_prime_factor; }
        /**
         * Returns a vector of `{ prime, index }`.
         * It is guaranteed that the returned vector is sorted in ascending order.
         */
        std::vector<std::pair<int, int>> factorize(int n) const noexcept {
            assert(0 < n and n <= _n);
            std::vector<std::pair<int, int>> prime_powers;
            while (n > 1) {
                int p = min_prime_factor[n], c = 0;
                do { n /= p, ++c; } while (n % p == 0);
                prime_powers.emplace_back(p, c);
            }
            return prime_powers;
        }
    private:
        const int _n;
        std::vector<int> min_prime_factor;
        std::vector<int> prime_list;
};
} // namespace suisen


#line 6 "library/sequence/powers.hpp"

namespace suisen {
    // returns { 0^k, 1^k, ..., n^k }
    template <typename mint>
    std::vector<mint> powers(uint32_t n, uint64_t k) {
        const auto mpf = LinearSieve(n).get_min_prime_factor();
        std::vector<mint> res(n + 1);
        res[0] = k == 0;
        for (uint32_t i = 1; i <= n; ++i) res[i] = i == 1 ? 1 : uint32_t(mpf[i]) == i ? mint(i).pow(k) : res[mpf[i]] * res[i / mpf[i]];
        return res;
    }
} // namespace suisen


#line 6 "library/sequence/stirling_number2.hpp"

namespace suisen {
    /**
     * return:
     *   vector<mint> v s.t. v[i] = S2[n,i] for i=0,...,k
     * constraints:
     *   0 <= n <= 10^6
     */
    template <typename FPSType>
    std::vector<typename FPSType::value_type> stirling_number2(int n) {
        using mint = typename FPSType::value_type;
        std::vector<mint> pows = powers<mint>(n, n);
        factorial<mint> fac(n);
        FPSType a(n + 1), b(n + 1);
        for (int i = 0; i <= n; ++i) {
            a[i] = pows[i] * fac.fac_inv(i);
            b[i] = i & 1 ? -fac.fac_inv(i) : fac.fac_inv(i);
        }
        a *= b, a.cut(n + 1);
        return a;
    }
    template <typename mint>
    std::vector<std::vector<mint>> stirling_number2_table(int n) {
        std::vector dp(n + 1, std::vector<mint>{});
        for (int i = 0; i <= n; ++i) {
            dp[i].resize(i + 1);
            dp[i][0] = 0, dp[i][i] = 1;
            for (int j = 1; j < i; ++j) dp[i][j] = dp[i - 1][j - 1] + j * dp[i - 1][j];
        }
        return dp;
    }
} // namespace suisen


#line 1 "library/sequence/binomial_coefficient_small_prime_mod.hpp"



#line 5 "library/sequence/binomial_coefficient_small_prime_mod.hpp"

namespace suisen {
    template <typename mint>
    struct BinomialCoefficientSmallPrimeMod {
        mint operator()(long long n, long long r) const {
            return binom(n, r);
        }
        static mint binom(long long n, long long r) {
            factorial<mint> fac(mint::mod() - 1);

            if (r < 0 or n < r) return 0;
            r = std::min(r, n - r);
            // Lucas's theorem
            mint res = 1;
            while (r) {
                int ni = n % mint::mod(), ri = r % mint::mod();
                if (ni < ri) return 0;
                res *= fac.binom(ni, ri);
                n = n / mint::mod(), r = r / mint::mod();
            }
            return res;
        }
    };
} // namespace suisen



#line 6 "library/sequence/stirling_number2_small_prime_mod.hpp"

namespace suisen {
    template <typename mint>
    struct StirlingNumber2SmallPrimeMod {
        mint operator()(long long n, long long k) const {
            return s2(n, k);
        }
        static mint s2(long long n, long long k) {
            static const std::vector<std::vector<mint>> table = stirling_number2_table<mint>(mint::mod() - 1);
            static const BinomialCoefficientSmallPrimeMod<mint> binom{};
            static const int p = mint::mod();
            if (k < 0 or n < k) return 0;
            if (n == 0) return 1;
            long long c = k / p, d = k % p;
            long long a = (n - c - 1) / (p - 1), b = (n - c - 1) % (p - 1) + 1;
            if (b != p - 1) return d <= b ? binom(a, c) * table[b][d] : 0;
            else return d ? binom(a, c) * table[p - 1][d] : binom(a, c - 1);
        }
    };
} // namespace suisen



#line 9 "test/src/sequence/stirling_number2_small_prime_mod/stirling_number_of_the_second_kind_small_p_large_n.test.cpp"

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int t, p;
    std::cin >> t >> p;
    mint::set_mod(p);

    suisen::StirlingNumber2SmallPrimeMod<mint> s2;

    while (t --> 0) {
        long long n, k;
        std::cin >> n >> k;
        std::cout << s2(n, k).val() << '\n';
    }

    return 0;
}
Back to top page