結果

問題 No.502 階乗を計算するだけ
ユーザー suisensuisen
提出日時 2024-01-30 21:10:33
言語 C++17
(gcc 12.3.0 + boost 1.83.0)
結果
AC  
実行時間 770 ms / 1,000 ms
コード長 12,959 bytes
コンパイル時間 2,875 ms
コンパイル使用メモリ 132,800 KB
実行使用メモリ 80,768 KB
最終ジャッジ日時 2024-09-28 10:17:30
合計ジャッジ時間 11,979 ms
ジャッジサーバーID
（参考情報） 		judge4 / judge1
このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 2 ms
6,812 KB
testcase_01 AC 2 ms
6,944 KB
testcase_02 AC 2 ms
6,944 KB
testcase_03 AC 2 ms
6,944 KB
testcase_04 AC 2 ms
6,940 KB
testcase_05 AC 2 ms
6,940 KB
testcase_06 AC 2 ms
6,944 KB
testcase_07 AC 2 ms
6,940 KB
testcase_08 AC 2 ms
6,944 KB
testcase_09 AC 2 ms
6,944 KB
testcase_10 AC 2 ms
6,940 KB
testcase_11 AC 2 ms
6,944 KB
testcase_12 AC 2 ms
6,940 KB
testcase_13 AC 2 ms
6,940 KB
testcase_14 AC 2 ms
6,944 KB
testcase_15 AC 2 ms
6,940 KB
testcase_16 AC 2 ms
6,940 KB
testcase_17 AC 2 ms
6,940 KB
testcase_18 AC 2 ms
6,940 KB
testcase_19 AC 2 ms
6,940 KB
testcase_20 AC 2 ms
6,940 KB
testcase_21 AC 2 ms
6,944 KB
testcase_22 AC 15 ms
10,560 KB
testcase_23 AC 7 ms
6,940 KB
testcase_24 AC 11 ms
8,092 KB
testcase_25 AC 3 ms
6,940 KB
testcase_26 AC 8 ms
6,940 KB
testcase_27 AC 6 ms
6,940 KB
testcase_28 AC 7 ms
6,940 KB
testcase_29 AC 4 ms
6,940 KB
testcase_30 AC 15 ms
10,084 KB
testcase_31 AC 9 ms
6,940 KB
testcase_32 AC 765 ms
80,640 KB
testcase_33 AC 762 ms
80,680 KB
testcase_34 AC 760 ms
80,760 KB
testcase_35 AC 770 ms
80,636 KB
testcase_36 AC 762 ms
80,728 KB
testcase_37 AC 759 ms
80,768 KB
testcase_38 AC 757 ms
80,768 KB
testcase_39 AC 763 ms
80,636 KB
testcase_40 AC 763 ms
80,636 KB
testcase_41 AC 765 ms
80,764 KB
testcase_42 AC 2 ms
5,376 KB
testcase_43 AC 2 ms
5,376 KB
testcase_44 AC 2 ms
5,376 KB
testcase_45 AC 2 ms
5,376 KB
testcase_46 AC 2 ms
5,376 KB
testcase_47 AC 2 ms
5,376 KB
testcase_48 AC 2 ms
5,376 KB
testcase_49 AC 2 ms
5,376 KB
testcase_50 AC 2 ms
5,376 KB
testcase_51 AC 2 ms
5,376 KB
権限があれば一括ダウンロードができます

ソースコード

diff # 

#define PROBLEM "https://yukicoder.me/problems/no/502"

#include <iostream>

#include <atcoder/modint>

using mint = atcoder::modint1000000007;

#include <utility>

#include <vector>
#include <atcoder/convolution>

#include <cassert>

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

namespace suisen {
    template <typename mint, typename Convolve,
        std::enable_if_t<std::is_invocable_r_v<std::vector<mint>, Convolve, std::vector<mint>, std::vector<mint>>, std::nullptr_t> = nullptr>
    std::vector<mint> shift_of_sampling_points(const std::vector<mint>& ys, mint t, int m, const Convolve &convolve) {
        const int n = ys.size();
        factorial<mint> fac(std::max(n, m));

        std::vector<mint> b = [&] {
            std::vector<mint> f(n), g(n);
            for (int i = 0; i < n; ++i) {
                f[i] = ys[i] * fac.fac_inv(i);
                g[i] = (i & 1 ? -1 : 1) * fac.fac_inv(i);
            }
            std::vector<mint> b = convolve(f, g);
            b.resize(n);
            return b;
        }();
        std::vector<mint> e = [&] {
            std::vector<mint> c(n);
            mint prd = 1;
            std::reverse(b.begin(), b.end());
            for (int i = 0; i < n; ++i) {
                b[i] *= fac.fac(n - i - 1);
                c[i] = prd * fac.fac_inv(i);
                prd *= t - i;
            }
            std::vector<mint> e = convolve(b, c);
            e.resize(n);
            return e;
        }();
        std::reverse(e.begin(), e.end());
        for (int i = 0; i < n; ++i) {
            e[i] *= fac.fac_inv(i);
        }

        std::vector<mint> f(m);
        for (int i = 0; i < m; ++i) f[i] = fac.fac_inv(i);
        std::vector<mint> res = convolve(e, f);
        res.resize(m);
        for (int i = 0; i < m; ++i) res[i] *= fac.fac(i);
        return res;
    }

    template <typename mint>
    std::vector<mint> shift_of_sampling_points(const std::vector<mint>& ys, mint t, int m) {
        auto convolve = [&](const std::vector<mint> &f, const std::vector<mint> &g) { return atcoder::convolution(f, g); };
        return shift_of_sampling_points(ys, t, m, convolve);
    }
} // namespace suisen

#include <atcoder/convolution>
#include <iostream>

namespace suisen::internal {
    template <typename T, typename R = T>
    std::vector<R> convolution_naive(const std::vector<T>& a, const std::vector<T>& b) {
        const int n = a.size(), m = b.size();
        std::vector<R> c(n + m - 1);
        if (n < m) {
            for (int j = 0; j < m; j++) for (int i = 0; i < n; i++) c[i + j] += R(a[i]) * b[j];
        } else {
            for (int i = 0; i < n; i++) for (int j = 0; j < m; j++) c[i + j] += R(a[i]) * b[j];
        }
        return c;
    }
} // namespace suisen

namespace suisen {
    template <typename mint, atcoder::internal::is_modint_t<mint>* = nullptr>
    std::vector<mint> arbitrary_mod_convolution(const std::vector<mint>& a, const std::vector<mint>& b) {
        int n = int(a.size()), m = int(b.size());

        if constexpr (atcoder::internal::is_static_modint<mint>::value) {
            int maxz = 1;
            while (not ((mint::mod() - 1) & maxz)) maxz <<= 1;
            int z = 1;
            while (z < n + m - 1) z <<= 1;
            if (z <= maxz) return atcoder::convolution<mint>(a, b);
        }

        if (n == 0 or m == 0) return {};
        if (std::min(n, m) <= 120) return internal::convolution_naive(a, b);

        static constexpr long long MOD1 = 754974721;  // 2^24
        static constexpr long long MOD2 = 167772161;  // 2^25
        static constexpr long long MOD3 = 469762049;  // 2^26
        static constexpr long long M1M2 = MOD1 * MOD2;
        static constexpr long long INV_M1_MOD2 = atcoder::internal::inv_gcd(MOD1, MOD2).second;
        static constexpr long long INV_M1M2_MOD3 = atcoder::internal::inv_gcd(M1M2, MOD3).second;

        std::vector<int> a2(n), b2(m);
        for (int i = 0; i < n; ++i) a2[i] = a[i].val();
        for (int i = 0; i < m; ++i) b2[i] = b[i].val();

        auto c1 = atcoder::convolution<MOD1>(a2, b2);
        auto c2 = atcoder::convolution<MOD2>(a2, b2);
        auto c3 = atcoder::convolution<MOD3>(a2, b2);

        const long long m1m2 = mint(M1M2).val();
        std::vector<mint> c(n + m - 1);
        for (int i = 0; i < n + m - 1; ++i) {
            // Garner's Algorithm
            // X = x1 + x2 * m1 + x3 * m1 * m2
            // x1 = c1[i], x2 = (c2[i] - x1) / m1 (mod m2), x3 = (c3[i] - x1 - x2 * m1) / m2 (mod m3)
            long long x1 = c1[i];
            long long x2 = (atcoder::static_modint<MOD2>(c2[i] - x1) * INV_M1_MOD2).val();
            long long x3 = (atcoder::static_modint<MOD3>(c3[i] - x1 - x2 * MOD1) * INV_M1M2_MOD3).val();
            c[i] = x1 + x2 * MOD1 + x3 * m1m2;
        }
        return c;
    }

    std::vector<__uint128_t> convolution_int(const std::vector<int> &a, const std::vector<int> &b) {
        int n = int(a.size()), m = int(b.size());

        auto check_nonnegative = [](int e) { return e >= 0; };
        assert(std::all_of(a.begin(), a.end(), check_nonnegative));
        assert(std::all_of(b.begin(), b.end(), check_nonnegative));

        if (n == 0 or m == 0) return {};
        if (std::min(n, m) <= 120) return internal::convolution_naive<int, __uint128_t>(a, b);

        static constexpr long long MOD1 = 754974721;  // 2^24
        static constexpr long long MOD2 = 167772161;  // 2^25
        static constexpr long long MOD3 = 469762049;  // 2^26
        static constexpr long long M1M2 = MOD1 * MOD2;
        static constexpr long long INV_M1_MOD2 = atcoder::internal::inv_gcd(MOD1, MOD2).second;
        static constexpr long long INV_M1M2_MOD3 = atcoder::internal::inv_gcd(M1M2, MOD3).second;

        auto c1 = atcoder::convolution<MOD1>(a, b);
        auto c2 = atcoder::convolution<MOD2>(a, b);
        auto c3 = atcoder::convolution<MOD3>(a, b);
        std::vector<__uint128_t> c(n + m - 1);
        for (int i = 0; i < n + m - 1; ++i) {
            // Garner's Algorithm
            // X = x1 + x2 * m1 + x3 * m1 * m2
            // x1 = c1[i], x2 = (c2[i] - x1) / m1 (mod m2), x3 = (c3[i] - x1 - x2 * m1) / m2 (mod m3)
            int x1 = c1[i];
            int x2 = (atcoder::static_modint<MOD2>(c2[i] - x1) * INV_M1_MOD2).val();
            int x3 = (atcoder::static_modint<MOD3>(c3[i] - x1 - x2 * MOD1) * INV_M1M2_MOD3).val();
            c[i] = x1 + x2 * MOD1 + __uint128_t(x3) * M1M2;
        }
        return c;
    }
} // namespace suisen

namespace suisen {
    // mod must be a prime number
    template <typename mint,
        std::enable_if_t<atcoder::internal::is_static_modint<mint>::value, std::nullptr_t> = nullptr>
    struct factorial_large {
        using value_type = mint;

        static constexpr int LOG_BLOCK_SIZE = 9;
        static constexpr int BLOCK_SIZE = 1 << LOG_BLOCK_SIZE;
        static constexpr int BLOCK_NUM = value_type::mod() >> LOG_BLOCK_SIZE;

        static inline int threshold = 2000000;

        static_assert(atcoder::internal::is_prime_constexpr(mint::mod()));

        static value_type fac(long long n) {
            if (n >= mint::mod()) return 0;
            return n <= threshold ? factorial<mint>{}.fac(n) : _large_fac(n);
        }
        static value_type fac_inv(long long n) {
            assert(n < (long long) mint::mod());
            return n <= threshold ? factorial<mint>{}.fac_inv(n) : _large_fac(n).inv();
        }
        static value_type binom(int n, int r) {
            if (r < 0 or r > n) return 0;
            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>
        static value_type polynom(const int n, const Ds& ...ds) {
            if (n < 0) return 0;
            long long sumd = 0;
            value_type 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;
        }
        static value_type perm(int n, int r) {
            if (r < 0 or r > n) return 0;
            return fac(n) * fac_inv(n - r);
        }
    private:
        static inline std::vector<value_type> _block_fac{};

        static void _build() {
            if (_block_fac.size()) {
                return;
            }
            std::vector<value_type> f{ 1 };
            f.reserve(BLOCK_SIZE);
            for (int i = 0; i < LOG_BLOCK_SIZE; ++i) {
                std::vector<value_type> g = shift_of_sampling_points<value_type>(f, 1 << i, 3 << i, arbitrary_mod_convolution<value_type>);
                const auto get = [&](int j) { return j < (1 << i) ? f[j] : g[j - (1 << i)]; };
                f.resize(2 << i);
                for (int j = 0; j < 2 << i; ++j) {
                    f[j] = get(2 * j) * get(2 * j + 1) * ((2 * j + 1) << i);
                }
            }
            // f_B(x) = (x+1) * ... * (x+B-1)
            if (BLOCK_NUM > BLOCK_SIZE) {
                std::vector<value_type> g = shift_of_sampling_points<value_type>(f, BLOCK_SIZE, BLOCK_NUM - BLOCK_SIZE, arbitrary_mod_convolution<value_type>);
                std::move(g.begin(), g.end(), std::back_inserter(f));
            } else {
                f.resize(BLOCK_NUM);
            }
            for (int i = 0; i < BLOCK_NUM; ++i) {
                f[i] *= value_type(i + 1) * BLOCK_SIZE;
            }
            // f[i] = (i*B + 1) * ... * (i*B + B)

            f.insert(f.begin(), 1);
            for (int i = 1; i <= BLOCK_NUM; ++i) {
                f[i] *= f[i - 1];
            }
            _block_fac = std::move(f);
        }

        static value_type _large_fac(int n) {
            _build();
            value_type res;
            int q = n / BLOCK_SIZE, r = n % BLOCK_SIZE;
            if (2 * r <= BLOCK_SIZE) {
                res = _block_fac[q];
                for (int i = 0; i < r; ++i) {
                    res *= value_type::raw(n - i);
                }
            } else if (q != factorial_large::BLOCK_NUM) {
                res = _block_fac[q + 1];
                value_type den = 1;
                for (int i = 1; i <= BLOCK_SIZE - r; ++i) {
                    den *= value_type::raw(n + i);
                }
                res /= den;
            } else {
                // Wilson's theorem
                res = value_type::mod() - 1;
                value_type den = 1;
                for (int i = value_type::mod() - 1; i > n; --i) {
                    den *= value_type::raw(i);
                }
                res /= den;
            }
            return res;
        }
    };
} // namespace suisen

int main() {
    long long n;
    std::cin >> n;
    std::cout << suisen::factorial_large<mint>{}.fac(n).val() << '\n';
}
0