P10594 BZOJ2445 最大团

· · 题解

或许更好的阅读体验。

题意：

一个 n 个点的无向图被叫做是一个好图，当且仅当该图的任意一个极大连通子图拥有相同的点数，并且任意一个极大连通子图都是完全图。

现有 m 种颜色和所有含有 n 个点且节点有标号的好图。我们需要将每个好图都染上一种颜色，两个不同的好图可以染相同颜色，求方案数对 10^9 - 401 取模的结果。

n, m \le 2 \times 10^9

思路：

考虑每个子图的大小 d，那么相当于将 n 个点划分为 \frac{n}{d} 个大小为 d 的组的方案数，一般来说是：

\binom{n}{d, d, \cdots, d} = \frac{n!}{(d!)^{\frac{n}{d}}}

但是这样我们是钦定了分组的顺序的，而染色显然是不区分的，于是多算了 (\frac{n}{d})! 的贡献。

所以我们的答案是：

m^{\sum\limits_{d \mid n} \frac{n!}{(d!)^{\frac{n}{d}} (\frac{n}{d})!}}

考虑费马小定理，我们只需要算 \sum\limits_{d \mid n} \frac{n!}{(d!)^{\frac{n}{d}} (\frac{n}{d})!} \bmod {10^9 - 402} 的值，发现 10^9 - 402 = 2 \times 13 \times 5281 \times 7283，于是可以对每个质因子算出模它后的值，最后 CRT 合并即可。

考虑如何算 \frac{n!}{(d!)^{\frac{n}{d}} (\frac{n}{d})!} \bmod p（p 为质数），首先朴素算不太可能，考虑 exLucas 的思路，计算 \frac{(n!)_p}{(d!)_p^{\frac{n}{d}} (\frac{n}{d}!)_p} \cdot p^{v_p(n!) - \frac{n}{d}v_p(d) - v_p(\frac{n}{d})} 即可。

时间复杂度为 O(\sqrt{n} \log n)。

完整代码：

#include<bits/stdc++.h>
#define ls(k) k << 1
#define rs(k) k << 1 | 1
#define fi first
#define se second
#define popcnt(x) __builtin_popcount(x)
#define open(s1, s2) freopen(s1, "r", stdin), freopen(s2, "w", stdout);
using namespace std;
typedef __int128 __;
typedef long double lb;
typedef double db;
typedef unsigned long long ull;
typedef long long ll;
bool Begin;
inline ll read(){
    ll x = 0, f = 1;
    char c = getchar();
    while(c < '0' || c > '9'){
        if(c == '-')
            f = -1;
        c = getchar();
    }
    while(c >= '0' && c <= '9'){
        x = (x << 1) + (x << 3) + (c ^ 48);
        c = getchar();
    }
    return x * f;
}
inline void write(ll x){
    if(x < 0){
        putchar('-');
        x = -x;
    }
    if(x > 9)
        write(x / 10);
    putchar(x % 10 + '0');
}
namespace CRTll{
    inline ll X(ll x, ll y, ll p) {
        return (__) x * y % p; 
    }
    inline ll exgcd(ll a, ll b, ll & x, ll & y) {
        if (b == 0) {
            x = 1;
            y = 0;
            return a;
        }
        ll x1, y1, d;
        d = exgcd(b, a % b, x1, y1);
        x = y1, y = x1 - a / b * y1;
        return d;
    }
    inline ll solve(ll n, vector<ll> R, vector<ll> M) {
        ll d = 0, r1 = 0, m1 = 0, r2 = 0, m2 = 0, p = 0, q = 0;
        for (int i = 0; i < n; i++) {
            if (!i)
                m1 = M[i], r1 = R[i];
            else {
                m2 = M[i], r2 = R[i];
                d = exgcd(m1, m2, p, q);
                p = X(p, ((r2 - r1 % m2 + m2) % m2) / d, m2 / d);
                p = (p % (m2 / d) + m2 / d) % (m2 / d);
                r1 = m1 * p + r1;
                m1 = m1 / d * m2;
            }
        }
        return (r1 % m1 + m1) % m1;
    }
};
namespace CRT{
    inline int X(int x, int y, int p) {
        return 1ll * x * y % p; 
    }
    inline int exgcd(int a, int b, int & x, int & y) {
        if (b == 0) {
            x = 1;
            y = 0;
            return a;
        }
        int x1, y1, d;
        d = exgcd(b, a % b, x1, y1);
        x = y1, y = x1 - a / b * y1;
        return d;
    }
    inline int solve(int n, vector<int> R, vector<int> M) {
        int d = 0, r1 = 0, m1 = 0, r2 = 0, m2 = 0, p = 0, q = 0;
        for (int i = 0; i < n; i++) {
            if (!i)
                m1 = M[i], r1 = R[i];
            else {
                m2 = M[i], r2 = R[i];
                d = exgcd(m1, m2, p, q);
                p = X(p, ((r2 - r1 % m2 + m2) % m2) / d, m2 / d);
                p = (p % (m2 / d) + m2 / d) % (m2 / d);
                r1 = m1 * p + r1;
                m1 = m1 / d * m2;
            }
        }
        return (r1 % m1 + m1) % m1;
    }
};
namespace Lucas{
    const int N = 1e5 + 10;
    int mod;
    int fac[N], inv[N];
    inline int qpow(int a, int b, int mod){
        int ans = 1;
        while(b){
            if(b & 1)
                ans = 1ll * ans * a % mod;
            a = 1ll * a * a % mod;
            b >>= 1;
        }
        return ans;
    }
    inline void init(int _mod){
        mod = _mod;
        fac[0] = 1;
        for(int i = 1; i < mod; ++i)
            fac[i] = 1ll * fac[i - 1] * i % mod;
        inv[mod - 1] = qpow(fac[mod - 1], mod - 2, mod);
        for(int i = mod - 2; i >= 0; --i)
            inv[i] = 1ll * inv[i + 1] * (i + 1) % mod;
    }
    inline int C(int n, int m){
        if(n < m)
            return 0;
        return 1ll * fac[n] * inv[m] % mod * inv[n - m] % mod;
    }
    inline int lucas(ll n, ll m){
        if(!m)
            return 1;
        return 1ll * lucas(n / mod, m / mod) * C(n % mod, m % mod) % mod;
    }
};
namespace exLucas{
    const int N = 1e4 + 10;
    int now[4], p, pk;
    int pm[4][N], id[N];
    vector<pair<int, int>> prime;
    vector<int> R, M;
    inline int exgcd(int a, int b, int & x, int & y) {
        if (b == 0) {
            x = 1;
            y = 0;
            return a;
        }
        int x1, y1, d;
        d = exgcd(b, a % b, x1, y1);
        x = y1, y = x1 - a / b * y1;
        return d;
    }
    inline int inv(int a, int m){
        int x, y;
        exgcd(a, m, x, y);
        return (x % m + m) % m;
    }
    inline int qpow(int a, ll b, int mod){
        int ans = 1;
        while(b){
            if(b & 1)
                ans = 1ll * ans * a % mod;
            a = 1ll * a * a % mod;
            b >>= 1;
        }
        return ans;
    }
    inline void init(int _p, int _pk, int t){
        p = _p, pk = _pk;
        id[p] = t;
        pm[t][0] = 1;
        for(int i = 1; i <= pk; ++i){
            if(i % p == 0){
                pm[t][i] = pm[t][i - 1];
                continue;
            }
            pm[t][i] = 1ll * pm[t][i - 1] * i % pk;
        }
        now[t] = pm[t][pk];
    }
    inline void init(int p){
        prime.clear();
        for(int i = 2; i * i <= p; ++i){
            if(p % i == 0){
                int num = 1;
                while(p % i == 0){
                    p /= i;
                    num *= i;
                }
                prime.push_back({i, num});
            }
        }
        if(p > 1)
            prime.push_back({p, p});
        int cnt = 0;
        for(auto t : prime)
          init(t.fi, t.se, cnt++);
    }
    inline ll vp(ll n, ll p){
        ll sum = 0;
        while(n){
            n /= p;
            sum += n;
        }
        return sum;
    }
    inline ll np(ll n, ll p, ll pk){
        if(!n)
            return 1;
        int sum = 1ll * pm[id[p]][n % pk] * np(n / p, p, pk) % pk;
        if(now[id[p]] == 1)
            return sum;
        if((n / pk) & 1)
            return sum ? (pk - sum) : sum;
        return sum;
    }
    inline int binom(ll n, ll m){
        R.clear(), M.clear();
        for(auto t : prime){
            p = t.fi, pk = t.se;
            ll sum = vp(n, p) - vp(m, p) - vp(n - m, p);
            if(sum >= (pk / p)){
                R.push_back(0);
                M.push_back(pk);
                continue;
            }
            int ans = 1ll * np(n, p, pk) * inv(np(m, p, pk), pk) % pk * inv(np(n - m, p, pk), pk) % pk * qpow(p, sum, pk) % pk;
            R.push_back(ans);
            M.push_back(pk);
        }
        return CRT::solve(prime.size(), R, M);
    }
    inline int get(ll n, ll d){
        R.clear(), M.clear();
        for(auto t : prime){
            p = t.fi, pk = t.se;
            ll sum = vp(n, p) - (n / d) * vp(d, p) - vp(n / d, p);
            if(sum >= (pk / p)){
                R.push_back(0);
                M.push_back(pk);
                continue;               
            }
            int ans = 1ll * np(n, p, pk) * inv(np(n / d, p, pk), pk) % pk * qpow(inv(np(d, p, pk), pk), n / d, pk) % pk;
            R.push_back(ans);
            M.push_back(pk);
        }
        return CRT::solve(prime.size(), R, M);
    }
};
const int mod = 1e9 - 401;
int T, n, m, ans;
inline void solve(){
    ans = 0;
    n = read(), m = read();
    for(ll d = 1; d * d <= n; ++d){
        if(n % d == 0){
            ans = (ans + exLucas::get(n, d)) % (mod - 1);
            if(d * d != n)
              ans = (ans + exLucas::get(n, n / d)) % (mod - 1);
        }
    }
    write(exLucas::qpow(m, ans, mod));
    putchar('\n');
}
int main(){
    exLucas::init(mod - 1);
    T = read();
    while(T--)
      solve();
    return 0;
}