Submission #74370337

---

Source Code Expand

// C++14

#ifndef use_int128
# define use_int128 1
#endif

#include <cstdint>
#include <type_traits>
#include <limits>
#include <iterator>
#include <numeric>
#include <functional>
#include <iostream>
#include <algorithm>

#define expect(cond) __builtin_expect(static_cast<bool>(cond), true)
#define unexpect(cond) __builtin_expect(static_cast<bool>(cond), false)

namespace modint_nm {

    using u32 = std::uint32_t;
    using u64 = std::uint64_t;
    using umx = std::uintmax_t;
    using i32 = std::int32_t;
    using imx = std::intmax_t;

    // for primitive root
    constexpr u32 phi(u32 n) noexcept {
        u32 ph = n;
        for (u32 i = 2; u64(i) * i <= n; ++i)if (n % i == 0u){ph = ph / i * (i - 1);while (n % i == 0u)n /= i;}
        return n > 1? ph / n * (n - 1) : ph;
    }
    constexpr bool is_prime(u32 v) noexcept {
        if (v <= 1)return 0;
        for (u32 i = 2; u64(i) * i <= v; ++i)if (v % i == 0)return 0;
        return 1;
    }

    template <u32 mod, typename type, typename = void> struct modval_t;
    template <u32 mod, typename type> struct modval_t <mod, type, std::enable_if_t<std::is_unsigned<type>::value> >
        { static constexpr u32 modval(type v) noexcept { return v % mod; } };
    template <u32 mod, typename type> struct modval_t <mod, type, std::enable_if_t<std::is_signed<type>::value> >
        { static constexpr u32 modval(type v) noexcept { v %= i32(mod);return v < 0? v + mod : v; } };
    template <u32 mod, typename type> struct modval_t <mod, type, std::enable_if_t<std::is_floating_point<type>::value> >
        { static constexpr u32 modval(type v) noexcept { return modval_t<mod, imx>::modval(imx(v)); } };
    #if use_int128
    template <u32 mod> struct modval_t <mod, __int128_t>
        { static constexpr u32 modval(__int128_t v) noexcept { v %= i32(mod);return v < 0? v + mod : v; } };
    template <u32 mod> struct modval_t <mod, __uint128_t>
        { static constexpr u32 modval(__uint128_t v) noexcept { return v % mod; } };
    #endif

    template <u32 mod> constexpr auto modsplus(u32 &a, u32 b) noexcept
        -> std::enable_if_t<(std::numeric_limits<u32>::max() / 2 + 1 < mod), void> { if (a >= mod - b)a -= mod - b;else a += b; }
    template <u32 mod> constexpr auto modsplus(u32 &a, u32 b) noexcept
        -> std::enable_if_t<std::numeric_limits<u32>::max() / 2 + 1 >= mod, void> { a += b;if (a >= mod)a -= mod; }
    template <u32 mod> constexpr auto modplus(u32 a, u32 b) noexcept
        -> std::enable_if_t<(std::numeric_limits<u32>::max() / 2 + 1 < mod), u32> { return a >= mod - b? a - (mod - b) : a + b; }
    template <u32 mod> constexpr auto modplus(u32 a, u32 b) noexcept
        -> std::enable_if_t<std::numeric_limits<u32>::max() / 2 + 1 >= mod, u32> { a += b;return a >= mod? a - mod : a; }

    template <u32 mod> constexpr void modssubt(u32 &a, u32 b) noexcept { if (a < b)a += mod - b; else a -= b; }
    template <u32 mod> constexpr u32 modsubt(u32 a, u32 b) noexcept { return a < b? a + (mod - b) : a - b; }

    template <u32 mod> constexpr auto modmult(u32 a, u32 b) noexcept
        -> std::enable_if_t<(u64(mod - 1) * (mod - 1) <= u32(-1)), u32> { return a * b % mod; }
    template <u32 mod> constexpr auto modmult(u32 a, u32 b) noexcept
        -> std::enable_if_t<(u64(mod - 1) * (mod - 1) > u32(-1)), u32> { return u64(a) * b % mod; }
    template <u32 mod> constexpr void modsmult(u32 &a, u32 b) noexcept { a = modmult<mod>(a, b); }

    template <u32 mod> constexpr u32 modfpow(u32 a, umx b) noexcept
        { u32 ret = 1;for (; b; b >>= 1, a = modmult<mod>(a, a))if (b & 1)ret = modmult<mod>(ret, a);return ret; }
    template <u32 mod> constexpr void modsfpow(u32 &a, umx b) noexcept { a = modfpow<mod>(a, b); }

    // 假定值在模范围内，减少一次取模操作
    struct in_mod_range_t { };
    constexpr in_mod_range_t in_mod_range {};

    struct modint_tag { };

    template <u32> struct modint;

    template <u32 mod_>
        struct modint {
            typedef modint_tag is_modint_t;
            static constexpr u32 mod = mod_;
            static constexpr bool mod_is_prime = is_prime(mod);
            u32 val;
            constexpr modint() noexcept : val() { }
            template <typename Val_t, typename = decltype(modval_t<mod, Val_t>::modval(std::declval<Val_t>())) >
                constexpr modint(Val_t Val) noexcept : val(modval_t<mod, Val_t>::modval(Val)) {}
            constexpr modint(const modint &) noexcept = default;
            constexpr modint(modint &&) noexcept = default;
            constexpr modint(in_mod_range_t, u32 v) noexcept : val{v} {}
            template <typename Val_t> constexpr modint &operator = (Val_t Val) noexcept
                { val = modval_t<mod, Val_t>::modval(Val);return *this; }
            constexpr modint &operator = (const modint &) noexcept = default;
            constexpr modint &operator = (modint &&) noexcept = default;
            constexpr modint &operator += (modint o) noexcept { modsplus<mod>(val, o.val);return *this; }
            constexpr modint &operator -= (modint o) noexcept { modssubt<mod>(val, o.val);return *this; }
            constexpr modint &operator *= (modint o) noexcept { modsmult<mod>(val, o.val);return *this; }
            constexpr modint &operator ++ () noexcept {if (unexpect(++val == mod))val = 0;return *this;}
            constexpr modint operator ++ (int) noexcept {modint R = *this;if (unexpect(++val == mod))val = 0;return R;}
            constexpr modint &operator -- () noexcept {if (unexpect(--val == u32(-1)))val = mod - 1;return *this;}
            constexpr modint operator -- (int) noexcept {modint R = *this;if (unexpect(--val == u32(-1)))val = mod - 1;return R;}
            constexpr modint operator - () const noexcept { return modint(in_mod_range, val? mod - val : 0u); }
            constexpr modint inv() const noexcept { return pow(*this, phi_mod - 1); }
            constexpr modint &operator /= (modint o) noexcept { modsmult<mod>(val, o.inv().val);return *this; }
            constexpr u32 value() const noexcept { return val; }
            explicit constexpr operator u32 () const noexcept { return val; }
            template <typename _Tp> constexpr operator _Tp () const noexcept { return (_Tp)val; }

            constexpr modint& plussubt(bool v, modint x) noexcept //+=(-1)^v*x
                { if (v)*this -= x; else *this += x; return *this; }
            template <typename _Tp, typename = std::enable_if_t<std::is_integral<_Tp>::value> >
                constexpr modint& plussubt(_Tp v, modint x) noexcept //+=(-1)^v*x
                    { if (v & 1)*this -= x; else *this += x; return *this; }

            static constexpr modint nan = modint(in_mod_range, mod);

            static constexpr u32 phi_mod = phi(mod);

        }; // __is_trivially_copyable(modint<...>) = true

    template <u32 mod>
        constexpr u32 modint<mod>::mod;
    template <u32 mod>
        constexpr modint<mod> modint<mod>::nan;
    template <u32 mod>
        constexpr bool modint<mod>::mod_is_prime;
    template <u32 mod>
        constexpr u32 modint<mod>::phi_mod;

    template <typename...> using __void_t = void;
    template <u32 mod, typename type> struct modval_t <mod, type,
            __void_t<std::enable_if_t<!std::is_arithmetic<type>::value>, decltype(type::operator modint<mod>()) > >
        { static constexpr u32 modval(type v) noexcept { return ((modint<mod>)v).value(); } };

    template <u32 mod> constexpr modint<mod> operator + (modint<mod> a, modint<mod> b) noexcept
        { return modint<mod>(in_mod_range, modplus<mod>(a.val, b.val)); }
    template <u32 mod, typename _T> constexpr modint<mod> operator + (_T a, modint<mod> b) noexcept { return modint<mod>(a) + b; }
    template <u32 mod, typename _T> constexpr modint<mod> operator + (modint<mod> a, _T b) noexcept { return a + modint<mod>(b); }

    template <u32 mod> constexpr modint<mod> operator - (modint<mod> a, modint<mod> b) noexcept
        { return modint<mod>(in_mod_range, modsubt<mod>(a.val, b.val)); }
    template <u32 mod, typename _T> constexpr modint<mod> operator - (_T a, modint<mod> b) noexcept { return modint<mod>(a) - b; }
    template <u32 mod, typename _T> constexpr modint<mod> operator - (modint<mod> a, _T b) noexcept { return a - modint<mod>(b); }

    template <u32 mod> constexpr modint<mod> operator * (modint<mod> a, modint<mod> b) noexcept
        { return modint<mod>(in_mod_range, modmult<mod>(a.val, b.val)); }
    template <u32 mod, typename _T> constexpr modint<mod> operator * (_T a, modint<mod> b) noexcept { return modint<mod>(a) * b; }
    template <u32 mod, typename _T> constexpr modint<mod> operator * (modint<mod> a, _T b) noexcept { return a * modint<mod>(b); }

    template <u32 mod> constexpr modint<mod> operator / (modint<mod> a, modint<mod> b) noexcept
        { return modint<mod>(in_mod_range, modmult<mod>(a.val, b.inv().val)); }
    template <u32 mod, typename _T> constexpr modint<mod> operator / (_T a, modint<mod> b) noexcept { return modint<mod>(a) / b; }
    template <u32 mod, typename _T> constexpr modint<mod> operator / (modint<mod> a, _T b) noexcept { return a / modint<mod>(b); }

    template <typename _Tp, typename = void> struct is_modint : std::false_type { };
    template <typename _Tp> struct is_modint <_Tp, typename _Tp::is_modint_t> : std::true_type { };

    namespace __modint_help {

        template <typename _T>
            constexpr auto __safe_l(_T a, u32 b) noexcept -> std::enable_if_t<std::is_unsigned<_T>::value, bool> {
                return a < b;
            }
        template <typename _T>
            constexpr auto __safe_l(_T a, u32 b) noexcept -> std::enable_if_t<std::is_signed<_T>::value, bool> {
                return a < _T(0)? true : std::make_unsigned_t<_T>(a) < b;
            }
        #if use_int128
        constexpr bool __safe_l(__uint128_t a, u32 b){
            return a < b;
        }
        constexpr bool __safe_l(__int128_t a, u32 b){
            return a < 0? true : __uint128_t(a) < b;
        }
        #endif

    }

    template <u32 mod, typename _T>
        constexpr std::enable_if_t<(mod > 2 && is_prime(mod)) &&
                std::is_convertible<_T, modint<mod - 1> >::value, modint<mod> >
            pow(modint<mod> a, _T b) noexcept
                { return modint<mod>(in_mod_range, modfpow<mod>(a.val, modint<mod - 1>(b).value())); }
    template <u32 mod, typename _T>
        constexpr std::enable_if_t<!(mod > 2 && is_prime(mod)) && std::is_unsigned<_T>::value, modint<mod> >
            pow(modint<mod> a, _T b) noexcept
                { constexpr u32 p = phi(mod); return modint<mod>(in_mod_range, modfpow<mod>(a.val, __modint_help::__safe_l(b, p)? b : b % p + p)); }
    template <u32 mod, typename _T>
        constexpr std::enable_if_t<!(mod > 2 && is_prime(mod)) && std::is_signed<_T>::value, modint<mod> >
            pow(modint<mod> a, _T b) noexcept {
                constexpr u32 p = phi(mod);
                if (b >= 0)return modint<mod>(in_mod_range, modfpow<mod>(a.val, __modint_help::__safe_l(b, p)? b : b % p + p));
                b = -b;
                umx T = umx(__modint_help::__safe_l(b, p)? b : b % p + p) * (p - 1);
                return modint<mod>(in_mod_range, modfpow<mod>(a.val, T < p? T : T % p + p));
            }

    template <u32 mod> constexpr modint<mod> fma(modint<mod> a, modint<mod> b, modint<mod> c) noexcept
        { return modint<mod>(in_mod_range, (u64(a.val) * b.val + c.val) % mod); }
    template <u32 mod> constexpr modint<mod> fma(modint<mod> a, modint<mod> b, u32 c) noexcept
        { return modint<mod>(in_mod_range, (u64(a.val) * b.val + c) % mod); }
    template <u32 mod> constexpr modint<mod> fma(in_mod_range_t, modint<mod> a, modint<mod> b, u64 c) noexcept //assume c + (mod-1)^2 < u64 max
        { return modint<mod>(in_mod_range, (u64(a.val) * b.val + c) % mod); }

    template <u32 mod> constexpr modint<mod> fmam(modint<mod> a, modint<mod> b, modint<mod> c, modint<mod> d) noexcept
        { return modint<mod>(in_mod_range, (u64(a.val) * b.val + u64(c.val) * d.val) % mod); }

    template <u32 mod> constexpr modint<mod> fms(modint<mod> a, modint<mod> b, modint<mod> c) noexcept
        { return modint<mod>(in_mod_range, (u64(a.val) * b.val + mod - c.val) % mod); }

    template <u32 mod> constexpr modint<mod> fmsm(modint<mod> a, modint<mod> b, modint<mod> c, modint<mod> d) noexcept
        { constexpr u64 mm = u64(mod) * mod;return modint<mod>(in_mod_range, (u64(a.val) * b.val + mm - u64(c.val) * d.val) % mod); }

    template <u32 mod> // a += b * c
        constexpr void splusmult(modint<mod> &a, modint<mod> b, modint<mod> c) noexcept { a = fma(b, c, a); }

    template <u32 mod> struct modint_partial_order
        { constexpr bool operator () (modint<mod> a, modint<mod> b) const noexcept { return a.val < b.val; } };// 偏序关系

    template <u32 mod> constexpr bool operator == (modint<mod> a, modint<mod> b) noexcept { return a.val == b.val; }
    template <u32 mod, typename _T> constexpr bool operator == (_T a, modint<mod> b) noexcept { return modint<mod>(a) == b; }
    template <u32 mod, typename _T> constexpr bool operator == (modint<mod> a, _T b) noexcept { return a == modint<mod>(b); }
    template <u32 mod> constexpr bool operator != (modint<mod> a, modint<mod> b) noexcept { return a.val != b.val; }
    template <u32 mod, typename _T> constexpr bool operator != (_T a, modint<mod> b) noexcept { return modint<mod>(a) != b; }
    template <u32 mod, typename _T> constexpr bool operator != (modint<mod> a, _T b) noexcept { return a != modint<mod>(b); }

    template <typename iter_t>
        using is_modint_iter = std::enable_if_t<std::is_base_of<std::random_access_iterator_tag,
            typename std::iterator_traits<iter_t>::iterator_category>::value, typename std::iterator_traits<iter_t>::value_type::is_modint_t>;

    template <typename T>
        struct __idnty {
            using type = T;
        };

    template <typename iter_t, typename = is_modint_iter<iter_t>, u32 mod = std::iterator_traits<iter_t>::value_type::mod>
        constexpr void calc_power(u32 n, iter_t pw, typename __idnty<modint<mod> >::type b){//let pw_i=b^i, 0<=i<n
            if (n == 0)return;
            pw[0] = 1;
            for (u32 i = 1; i < n; ++i)pw[i] = pw[i - 1] * b;
        }

    template <typename iter_t, typename = is_modint_iter<iter_t>, u32 mod = std::iterator_traits<iter_t>::value_type::mod>
        constexpr void sieve_power(u32 n, iter_t pw, u32 e){//let pw_i=i^e, 0<=i<n
            if (n == 0)return;
            if (n == 1){
                pw[0] = 0;
                return;
            }
            bool *npr = new bool[n]{};
            u32 *prm = new u32[n <= 100? n - (n / 2) + 1 : 2 * n / std::__lg(n)]{};
            u32 l = 0;
            pw[0] = 0;
            pw[1] = 1;
            for (u32 i = 2; i < n; ++i){
                if (!npr[i]){
                    pw[i] = pow(modint<mod>(i), e);
                    prm[l++] = i;
                }
                for (u32 j = 0; j < l; ++j){
                    if (u64(i) * prm[j] >= n)break;
                    npr[i * prm[j]] = 1;
                    pw[i * prm[j]] = pw[i] * pw[prm[j]];
                    if (i % prm[j] == 0)break;
                }
            }
            delete [] npr; delete [] prm;
        }

    template <typename iter_t, typename = is_modint_iter<iter_t> >
        constexpr void calc_inv(iter_t bg, iter_t ed){//let bg[i] = bg[i]^{-1}
            constexpr u32 mod = std::iterator_traits<iter_t>::value_type::mod;
            if (bg == ed)return;
            u32 n = ed - bg;
            modint<mod> *s = new modint<mod>[n + 1];
            modint<mod> *sv = new modint<mod>[n + 1];
            s[0] = 1;
            for (u32 i = 1; i <= n; ++i)s[i] = s[i - 1] * bg[i - 1];
            sv[n] = s[n].inv();
            for (u32 i = n; i; --i)sv[i - 1] = sv[i] * bg[i - 1];
            for (u32 i = 1; i <= n; ++i)bg[i - 1] = sv[i] * s[i - 1];
            delete []s; delete []sv;
        }

    template <u32 mod> constexpr bool isnan(modint<mod> a) noexcept { return a == modint<mod>::nan; }

    // 加速求和
    template <u32 mod>
        struct accumulate_temporary_result {
            umx val;
            static constexpr umx limit = std::numeric_limits<umx>::max() - mod;
            static constexpr umx limit2 = std::numeric_limits<umx>::max() - umx(mod) * mod;
            constexpr accumulate_temporary_result() noexcept : val() {}
            constexpr accumulate_temporary_result(modint<mod> v) noexcept : val(v.value()) {}
            constexpr accumulate_temporary_result(umx v) noexcept : val(v) {}
            constexpr accumulate_temporary_result(const accumulate_temporary_result &r) noexcept : val(r.val) {}
            constexpr accumulate_temporary_result(const accumulate_temporary_result &&r) noexcept : val(r.val) {}
            template <typename _Val> constexpr accumulate_temporary_result(_Val v) noexcept : val(modval_t<mod, _Val>::modval(v)) {}
            template <typename _Val> constexpr accumulate_temporary_result &operator += (umx v) noexcept
                { if (unexpect(val >= limit))val %= mod; val += modval_t<mod, _Val>::modval(v); return *this; }
            constexpr accumulate_temporary_result &operator += (modint<mod> v) noexcept
                { if (unexpect(val >= limit))val %= mod; val += v.value(); return *this; }
            constexpr void plus_modint_prod(modint<mod> a, modint<mod> b) noexcept
                { if (unexpect(val >= limit2))val %= mod; val += umx(a.value()) * b.value(); }// += a * b
            constexpr umx value() const noexcept { return val; }
            constexpr modint<mod> mod_val() const noexcept { return modint<mod>(val); }
            constexpr operator modint<mod> () const noexcept { return modint<mod>(val); }
        };

    // 多项相加优化
    template <typename iter_t, typename iter_value_t = typename std::iterator_traits<iter_t>::value_type,
            u32 mod = iter_value_t::mod>
        constexpr modint<mod> mod_plus(iter_t bg, iter_t ed){
            typedef modint<mod> type;
            u64 sum = std::accumulate(bg, ed, u64(), [](u64 a, iter_value_t v){
                return a + type(v).value();
            });
            return type(sum);
        }
    template <typename iter_t, u32 mod = std::iterator_traits<iter_t>::value_type::mod>
        constexpr modint<mod> mod_plus(iter_t bg, u32 n){
            typedef modint<mod> type;
            if (n < 4){
                switch (n){
                    case 0: return type();
                    case 1: return *bg;
                    case 2: return *bg + *(bg + 1);
                    case 3: return *bg + *(bg + 1) + *(bg + 2);
                    default : __builtin_unreachable();
                }
            }
            return mod_plus(bg, bg + n);
        }

    // 多项相乘
    template <typename iter_t, u32 mod = std::iterator_traits<iter_t>::value_type::mod>
        constexpr modint<mod> mod_mult(iter_t bg, iter_t ed){
            typedef modint<mod> type;
            return std::accumulate(bg, ed, type(1), std::multiplies<type>{});
        }

    template <typename iter1_t, typename iter2_t,
            u32 mod = std::iterator_traits<iter1_t>::value_type::mod,
            u32 mod2 = std::iterator_traits<iter2_t>::value_type::mod,
            typename = std::enable_if_t<mod == mod2> >
        constexpr modint<mod> mod_inner_product(iter1_t bg1, iter1_t ed1, iter2_t bg2, modint<mod> init = 0){
            accumulate_temporary_result<mod> atr(init);
            for (; bg1 != ed1; ++bg1, ++bg2)atr.plus_modint_prod(*bg1, *bg2);
            return atr.mod_val();
        }

    inline namespace predefine_modint {
        using modp0 = modint<998244353>;
        using modp1 = modint<1000000007>;
        using modp2 = modint<167772161>;
        using modp3 = modint<469762049>;
        using modp4 = modint<754974721>;
        using modp5 = modint<1004535809>;
    }

}

#include <functional>

namespace std {

    using u32 = std::uint32_t;
    using u64 = std::uint64_t;

    template <u32 mod> inline ostream &operator << (ostream &os, modint_nm::modint<mod> v){ return os << v.value(); }
    template <u32 mod> inline istream &operator >> (istream &is, modint_nm::modint<mod> &v){ u64 x;is >> x;v = x;return is; }

    template <u32 mod>
        struct numeric_limits <modint_nm::modint<mod> > : numeric_limits<u32> { };

    template <u32 mod>
        struct hash<modint_nm::modint<mod> > {
            inline size_t operator () (const modint_nm::modint<mod> &v) const {
                return hash<u32>{}(v.value());
            }
        };

    template <u32 mod> constexpr modint_nm::modint<mod> fma(modint_nm::modint<mod> a,
            modint_nm::modint<mod> b, modint_nm::modint<mod> c) noexcept
        { return modint_nm::fma(a, b, c); }
    template <u32 mod> constexpr modint_nm::modint<mod> fma(modint_nm::modint<mod> a,
            modint_nm::modint<mod> b, u32 c) noexcept
        { return modint_nm::fma(a, b, c); }

}

/**
 * @todo
 * 1. 加入左移和右移支持
 * 2. 加入对卡特兰数等常见数列的支持
 * 3. 部分函数加入参数为常量时的优化
 * 4. 加入模数非质数时的除法
 */

#include <cstdint>
#include <algorithm>

namespace fact {

    using u32 = std::uint32_t;
    using u64 = std::uint64_t;

    struct default_mod {
        u32 mod;
        constexpr u32 operator () (u32 a, u32 b) noexcept {
            return u64(a) * b % mod;
        }
    };

    /// @brief 对于所有 0<=i<n，令 fc[i] = i! \bmod mod
    /// 传入的 mul 需要实现快速模乘
    template <typename val_t, typename modmul>
    constexpr void calfc(u32 n, val_t *fc, u32 mod, modmul &&mul){
        if (n == 0)return;
        if (mod == 0)throw "mod = 0.";
        if (mod == 1){
            std::fill_n(fc, n, val_t(0));
            return;
        }
        fc[0] = 1;
        for (u32 i = 1; i < n; ++i)fc[i] = mul(fc[i - 1], i);
    }
    template <typename val_t>
    constexpr void calfc(u32 n, val_t *fc, u32 mod){
        calfc(n, fc, mod, default_mod{mod});
    }

    /// @brief 对于所有 0<=i<n，令 fc[i] = i^{-1} \bmod mod，跳过 i\bmod mod = 0 的位置
    /// @pre mod 为质数
    /// 传入的 mul 需要实现快速模乘
    template <typename val_t, typename modmul>
    constexpr void caliv(u32 n, val_t *iv, u32 mod, modmul &&mul) noexcept {
        if (n == 0)return;
        if (n > mod){
            caliv(mod, iv, mod, (modmul&&)mul);
            u32 i = mod;
            for (; i + mod < n; i += mod)
                std::copy_n(iv + 1, mod - 1, iv + i + 1);
                    // 以 mod 为周期填充，即 [+1,+mod) 拷贝到 [+1+mod,2mod), [+1+2mod,3mod), ...
            std::copy_n(iv + 1, n - i - 1, iv + i + 1);
            return;
        }
        // n <= mod
        iv[1] = 1;
        for (u32 i = 2; i < n; ++i)iv[i] = mul(iv[mod % i], mod - mod / i);
    }
    template <typename val_t>
    constexpr void caliv(u32 n, val_t *iv, u32 mod) noexcept {
        caliv(n, iv, mod, default_mod{mod});
    }

    /// @brief 对于所有 0<=i<n，令 ivfc[i] = i!^{-1} \bmod mod
    /// @pre n <= mod，mod 为质数
    /// 传入的 mul 需要实现快速模乘
    template <typename val_t, typename modmul>
    constexpr void calivfc(u32 n, val_t *ivfc, u32 mod, modmul &&mul) noexcept {
        if (n == 0)return;
        caliv(n, ivfc, mod, (modmul&&)mul);
        ivfc[0] = 1;
        for (u32 i = 1; i < n; ++i)ivfc[i] = mul(ivfc[i - 1], ivfc[i]);
    }
    template <typename val_t>
    constexpr void calivfc(u32 n, val_t *ivfc, u32 mod) noexcept {
        calivfc(n, ivfc, mod, default_mod{mod});
    }

    /// @brief 求 n! \bmod mod
    /// 保证 O(n)
    /// 传入的 mul 需要实现快速模乘
    template <typename modmul>
    constexpr u32 factorial(u32 n, u32 mod, modmul &&mul) noexcept {
        if (__builtin_expect(mod == 1, false))return 0;
        if (__builtin_expect(n >= mod, false))return 0;
        u32 ret = 1;
        while (n)ret = mul(ret, n--);
        return ret;
    }
    constexpr u32 factorial(u32 n, u32 mod) noexcept {
        return factorial(n, mod, default_mod{mod});
    }

    template <typename val_t, typename modmul = default_mod>
    struct binom_base {
        modmul mul;
        const val_t *fc; // 要求构造时指向一块已经计算好阶乘的空间
        const val_t *ivfc; // 要求构造时指向一块已经计算好阶乘逆的空间
        constexpr val_t fact(u32 n) const noexcept {
            return fc[n];
        }
        constexpr val_t invfact(u32 n) const noexcept {
            return ivfc[n];
        }
    };

    struct in_range_t { } in_range;

    using i32 = std::int32_t;

    template <typename val_t, typename modmul = default_mod>
    struct combination {
        binom_base<val_t, modmul> base;
        template <typename U, typename V> constexpr auto mul(U a, V b) const noexcept { return base.mul(a, b); }
        constexpr u32 operator () (i32 n, i32 m) const noexcept
            { return ((n - m) | m) >> 31? 0 : u32(mul(base.fc[n], mul(base.ivfc[m], base.ivfc[n - m]))); }
        constexpr u32 operator () (in_range_t, i32 n, i32 m) const noexcept
            { return mul(base.fc[n], mul(base.ivfc[m], base.ivfc[n - m])); }
        constexpr u32 iv(i32 n, i32 m) const noexcept
            { return ((n - m) | m) >> 31? u32(-1) : u32(mul(base.ivfc[n], mul(base.fc[m], base.fc[n - m]))); }
        constexpr u32 iv(in_range_t, i32 n, i32 m) const noexcept
            { return mul(base.ivfc[n], mul(base.fc[m], base.fc[n - m])); }
    };

    template <typename val_t, typename modmul = default_mod>
    struct permutation {
        binom_base<val_t, modmul> base;
        template <typename U, typename V> constexpr auto mul(U a, V b) const noexcept { return base.mul(a, b); }
        constexpr u32 operator () (i32 n, i32 m) const noexcept
            { return ((n - m) | m) >> 31? 0 : u32(mul(base.fc[n], base.ivfc[n - m])); }
        constexpr u32 operator () (in_range_t, i32 n, i32 m) const noexcept
            { return mul(base.fc[n], base.ivfc[n - m]); }
        constexpr u32 iv(i32 n, i32 m) const noexcept
            { return ((n - m) | m) >> 31? u32(-1) : u32(mul(base.ivfc[n], base.fc[n - m])); }
        constexpr u32 iv(in_range_t, i32 n, i32 m) const noexcept
            { return mul(base.ivfc[n], base.fc[n - m]); }
    };

    template <typename val_t, typename modmul = default_mod>
    struct grid { // grid(n,m) : 网格 (0,0) 走到 (n,m) 的最短路径方案数，等于 C(n+m,n)
        binom_base<val_t, modmul> base;
        template <typename U, typename V> constexpr auto mul(U a, V b) const noexcept { return base.mul(a, b); }
        constexpr u32 operator () (i32 n, i32 m) const noexcept
            { return (n | m) >> 31? 0 : u32(mul(base.fc[n + m], mul(base.ivfc[n], base.ivfc[m]))); }
        constexpr u32 operator () (in_range_t, i32 n, i32 m) const noexcept
            { return mul(base.fc[n + m], mul(base.ivfc[n], base.ivfc[m])); }
        constexpr u32 iv(i32 n, i32 m) const noexcept
            { return (n | m) >> 31? u32(-1) : u32(mul(base.ivfc[n + m], mul(base.fc[n], base.fc[m]))); }
        constexpr u32 iv(in_range_t, i32 n, i32 m) const noexcept
            { return mul(base.ivfc[n + m], mul(base.fc[n], base.fc[m])); }
    };

    template <typename val_t, typename modmul = default_mod>
    struct split { // split(n,m) : 将 n 划分为 m 个非负整数的方案数，split(0,0)=1，其它 split(n,m)=C(n+m-1,m-1)
        binom_base<val_t, modmul> base;
        template <typename U, typename V> constexpr auto mul(U a, V b) const noexcept { return base.mul(a, b); }
        constexpr u32 operator () (i32 n, i32 m) const noexcept
            { return !m? u32(n == 0) : (n | (m - 1)) >> 31? 0 : u32(mul(base.fc[n + m - 1], mul(base.ivfc[n], base.ivfc[m - 1]))); }
        constexpr u32 operator () (in_range_t, i32 n, i32 m) const noexcept
            { return mul(base.fc[n + m - 1], mul(base.ivfc[n], base.ivfc[m - 1])); }
        constexpr u32 iv(i32 n, i32 m) const noexcept
            { return !m? n? u32(-1) : 0 : (n | (m - 1)) >> 31? u32(-1) : u32(mul(base.ivfc[n + m - 1], mul(base.fc[n], base.fc[m - 1]))); }
        constexpr u32 iv(in_range_t, i32 n, i32 m) const noexcept
            { return mul(base.ivfc[n + m - 1], mul(base.fc[n], base.fc[m - 1])); }
    };

    /// @brief 组合数上指标求和，求 \sum_{i=m}^n \binom im
    /// binom 实现了计算组合数
    template <typename _Binom>
        constexpr auto upper_prefix_sum(u32 n, u32 m, _Binom &&binom) noexcept -> decltype(binom(u32{}, u32{})) {
            if (n < m)return 0;
            return binom(n + 1, m + 1);
        }

}

#include <bits/stdc++.h>
using namespace std;
using mod = modint_nm::modp1;
int n, k, a[310], b[310];
bool cr[610][610]; // exists cross edge
int ce[610];
inline int cem(int l, int r){
    return ce[r - 1] - ce[l - 1];
}
mod fc[610], ivfc[610];
mod piv2[610];
inline mod mt(int c){ // match 2c elem into c pairs
    return fc[c + c] * piv2[c] * ivfc[c];
}
mod f[610][610]; // f_{l,r} : [l,r) is a block
int main(){
    ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    cin >> n >> k;
    fill_n(ce + 1, n + n, 1);
    for (int i = 1; i <= k; ++i)cin >> a[i] >> b[i], a[i] > b[i] && (swap(a[i], b[i]), 0), ce[a[i]] = ce[b[i]] = 0;
    for (int i = 1; i <= n + n; ++i)ce[i] += ce[i - 1];
    for (int l = 1; l <= n + n; ++l)
        for (int r = l + 1; r <= n + n; ++r)
            cr[l][r] = cr[r][l] = [&]{
                // clog << l << "~" << r << endl;
                for (int i = 1; i <= k; ++i){
                    if (l <= a[i] && b[i] < r)continue;
                    if (b[i] < l || a[i] >= r)continue;
                    if (a[i] < l && b[i] >= r)continue;
                    return 1;
                }
                return 0;
            }();
    // for (int l = 1; l <= n + n; ++l)
    //     for (int r = l + 1; r <= n + n; ++r)
    //         if (cr[l][r])clog << l << "-" << r << endl;
    fact::calfc(n + n + 1, fc, mod::mod, multiplies<>{});
    fact::calivfc(n + n + 1, ivfc, mod::mod, multiplies<>{});
    modint_nm::calc_power(n + n + 1, piv2, mod(1) / 2);
    for (int l = 1; l <= n + n; ++l)
        for (int r = l + 1; r <= n + n; ++r){
            if (cr[l][r] || (cem(l, r) & 1))continue;
            f[l][r] = mt(cem(l, r) >> 1);
            for (int k = l; k < r; ++k)
                if (~cem(k, r) & 1)
                    f[l][r] -= f[l][k] * mt(cem(k, r) >> 1);
            // clog << l << "," << r << ":" << f[l][r] << endl;
        }
    mod rs = mt(n - k);
    for (int l = 1; l <= n + n; ++l)
        for (int r = l + 1; r <= n + n; ++r){
            // [l,r) and [r,n]-[1,l)
            if (cr[l][r])continue; // exists cross
            if (cem(l, r) & 1)continue;
            // clog << l << "~" << r << ":" << f[l][r] << "*" << mt((n - k) - (cem(l, r) >> 1)) << endl;
            rs += f[l][r] * mt((n - k) - (cem(l, r) >> 1));
        }
    cout << rs << endl;
    return 0;
}

Submission Info

Judge Result