题解：P10594 BZOJ2445 最大团

· · 题解

题面

更好的阅读体验

其实是参考了另一篇题解，但感觉那篇题解很多细节讲的不是很清楚，所以自己写了一篇题解加深印象。

我们假设有 k 种合法的 n 个点的图，现在要对每个图染一种颜色，求方案数。所以答案是 m^k。

现在我们要求 k 可以枚举每个子图的大小 d,那么将 n 个点每 d 个为一组进行分组的方案数为 {n\choose d,d,\dots}，但是这样我们是钦定了顺序的，而每个图内染一种颜色显然是不关心顺序的，所以每种方案都算了 (\frac nd)!，所以 k=\sum\limits_{d|n}\frac{n!}{(d!)^{\frac nd}(\frac nd)!}。

最终答案为 m^{\sum\limits_{d|n}\frac{n!}{(d!)^{\frac nd}(\frac nd)!}} \bmod 10^9-401。

跟据费马小定理，我们实际要求 m^{\sum\limits_{d|n}\frac{n!}{(d!)^{\frac nd}(\frac nd)!}\bmod 10^9-402}\bmod 10^9-401，关键是指数部分怎么求。首先对 10^9-401 分解质因数，10^9-402=2\times13\times5281\times7283，那么我们可以先对每一个质因数算出取模后的答案再用 CRT 合并。

于是我们转化为解决子问题：\sum\limits_{d|n}\frac{n!}{(d!)^{\frac nd}(\frac nd)!}\bmod p，这里问题在于 d!，(\frac nd)! 中可能含有因子 p 无法做逆元操作，于是我们用类似 ExLucas 的思路处理一下。

\frac{n!}{(d!)^{\frac nd}(\frac nd)!}\bmod p=\frac{\frac{n!}{p^i}}{(\frac{d!}{p^j})^{\frac nd}\frac{(\frac nd)!}{p^k}}\times p^{i-j\frac nd-k}\bmod p

那么我们现在需要求的就是形如 \frac{n!}{p^x}\bmod p，我们首先将 1~n 中 p 的倍数提一个 p 的因子,剩下的按模 p 的模数分组，最后有可能还会剩下几个无法分为完整的一组，即 n!=p^{\lfloor\frac np\rfloor}（\lfloor\frac np\rfloor）!\times(p-1)!^{\lfloor\frac np\rfloor}\times\prod\limits_{i=p\lfloor\frac np\rfloor+1}^{i<=n}i 到这一步我们就可以递归处理下去了。除去所以因子 p 后我们就可以求逆元了，那么这道题就做完了。

code

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const ll mod=1e9-401,N=1e4+5,p[5]={0,2,13,5281,7283},M=1e9-402;
ll qp(ll x,ll y,ll P)
{
    ll res=1;
    while(y)
    {
        if(y&1) (res*=x)%=P;
        (x*=x)%=P,y>>=1;
    }
    return res;
}
void exgcd(ll a,ll b,ll &x,ll &y){
    ll t;
    if(!b){x=1,y=0;return;}
    exgcd(b,a%b,x,y);
    t=x,x=y,y=t-a/b*y;
}
ll inv(ll a,ll P)
{
    ll x,y;
    exgcd(a,P,x,y);
    return (x+P)%P;
}
ll g(ll i,ll P){return (i<P)?0:(g(i/P,P)+i/P);}
ll n,m,ans[5],fac[N],res;
ll S(ll now, ll P)
{
    if(now<P) return fac[now]%P;
    ll res=1;
    for(ll i=(now/P)*P+1;i<=now;i++) (res*=i)%=P;
    return qp(fac[P-1]%P,now/P,P)*S(now/P,P)%P*res%P;
}
ll solve(ll now,ll d,ll P)
{
    if(g(now,P)-(now/d)*g(d,P)-g(now/d,P)) return 0;
    return S(now,P)*inv(qp(S(d,P),now/d,P),P)%P*inv(S(now/d,P),P)%P;
}
inline ll rd()
{
    char c;ll f=1;
    while(!isdigit(c=getchar())) if(c=='-')f=-1;
    ll x=c-'0';
    while(isdigit(c=getchar())) x=x*10+(c^48);
    return x*f;
}
int main()
{
    fac[0]=1;
    for(int i=1;i<p[4];i++) fac[i]=fac[i-1]*i%M;
    for(int T=rd();T--;)
    {
        n=rd(),m=rd(),res=0;
        for(int i=1;i<=4;i++)
        {
            ans[i]=0;
            for(ll d=1;d*d<=n;d++) if(n%d==0)
            {
                (ans[i]+=solve(n,d,p[i]))%=p[i];
                if(d*d!=n) (ans[i]+=solve(n,n/d,p[i]))%=p[i];   
            }
        }
        for(int i=1;i<=4;i++)
        {
            ll Mi=M/p[i];
            (res+=ans[i]*Mi%M*inv(Mi,p[i]))%=M;
        }
        cout<<qp(m,res,mod)<<endl;
    }
    return 0; 
}