题面
其实是参考了另一篇题解，但感觉那篇题解很多细节讲的不是很清楚，所以自己写了一篇题解加深印象。
题意极其抽象，所以我们根据答案反推题面，这道题的意思应该是：列出所有合法的 \(n\) 个点的图后（假设有 \(k\) 种),接着对每个图染一种颜色，求方案数。所以答案是 \(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)!}} (\mod 10^9-401)\)
跟据费马小定理，我们实际要求 \(m^{\sum\limits_{d|n}\frac{n!}{(d!)^{\frac nd}(\frac nd)!}(\mod 10^9-402)}(\mod 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)!}(\mod p)\)，这里问题在于 \(d!\),\((\frac nd)!\) 中可能含有因子 \(p\) 无法做逆元操作，于是我们用类似exlucas的思路处理一下
那么我们现在需要求的就是形如 \(\frac{n!}{p^x}(\mod 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\) 后我们就可以求逆元了，那么这道题就做完了。
代码
#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;
}
