【大意】
给定一个玩具，其周围有 (n) 个小珠子，中间有 (1) 个大珠子。要求用 (k) 个颜色对所有珠子着色，要求任意相邻小珠子颜色互不相同，大珠子与任意小珠子颜色互不相同。问旋转同构的方案数。
【分析】
优先枚举中心大珠子的着色方案，为 (k) 种。等价于用 ((k-1)) 个颜色对剩余 (n) 个小珠子着色的方案数。
显然，绕圆心均匀旋转 (0)~((n-1)) 个小珠子的操作 ({ rot_0, rot_1, rot_2, cdots, rot_{n-1} }) 构成一个置换群 (G) 。
考虑使用 Burniside 引理求解： (displaystyle kcdot{1over |G|}sum_{gin G}|X^g|={kover n}sum_{i=0}^{n-1}f(i, k-1))
其中，(f(i, k-1)) 表示所有相邻颜色不同的着色方案中，绕圆心旋转 (i) 个小珠子的操作 (rot_i) 不动点数
不难想到，所有着色方案中，若是 (rot_i) 的不动点，则对于任意序号 (a) 的小珠子，它的颜色必须和序号为 ((a+i)) 的小珠子相同
因此，(rot_i) 操作将小珠子的序号划分为 (gcd(i, n)) 个等价类，而同一等价类的染色方案必须相同
当 (gcd(i, n)=1) 时，由于相邻小珠子染色不同，方案数为 (0)
当 (gcd(i, n)=d) 时，用 (h(d, k-1)) 表示其染色方案数
代入原式得到：(displaystyle {kover n}sum_{i=0}^{n-1}h(gcd(i, n), k-1)={kover n}sum_{i=1}^n h(gcd(i, n), k-1))
采用莫比乌斯反演：(displaystyle {kover n}sum_{i=1}^n h(gcd(i, n), k-1)={kover n}sum_{dmid n\d>1}h(d, k-1)sum_{i=1}^n [gcd(i, n)=d]={kover n}sum_{dmid n\d>1}h(d, k-1)cdot oldsymbol varphi({nover d}))
现考虑如何计算 (h(d, k),d>1) ，表示 (n) 被划分为 (d) 个等价类，用 (k) 个颜色进行着色，使得相邻等价类颜色不同的旋转同构方案数
考虑将 (d) 个等价类看成 (d) 元环上的不同点，则 (h(d, k)) 等价于求解用 (k) 个颜色对 (d) 元环进行着色，使得相邻点颜色不同的旋转不同构方案数
故 (h(2, k)=k(k-1), h(3, k)=k(k-1)(k-2))
当 (d>3) 时，考虑破环成链，考虑链上第一个点和倒数第二个点颜色是否相同：
-
若相同，则最后一个点的着色方案数为 ((k-1)) 。前 ((d-1)) 个点由于第一个和倒数第二个颜色相同（等价于一个点），着色方案为 (h(d-2, k)) 。由乘法原理，方案数为 ((k-1)cdot h(d-2, k))
-
若不同，则最后一个点的着色方案数为 ((k-2)) 。前 ((d-1)) 个点由于第一个和倒数第二个颜色不同，可以直接相连，故着色方案数为 (h(d-1, k)) 。由乘法原理，方案数为 ((k-2)cdot h(d-1, k))
故根据加法原理，(h(d, k)=(k-1)cdot h(d-2, k)+(k-2)cdot h(d-1, k))
特征方程为 (x^2+(1-k)x+(2-k)=[x-(k-1)]cdot [x-(-1)]=0) ，得 (h(d, k)=C_1cdot (k-1)^d+C_2cdot (-1)^d)
代入 (h(2, k)) 与 (h(3, k)) 得 (h(d, k)=(k-1)^d+(k-1)cdot (-1)^d)
由于 (oldsymbol varphi({nover d})) 均为 (n) 的因数，因此本人通过 divit 预处理 (n) 的所有质因数，优化求解复杂度
【代码】
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MOD=1e9+7;
inline ll fpow(ll a,ll x) { ll ans=1; for(;x;x>>=1,a=a*a%MOD) if(x&1) ans=ans*a%MOD; return ans; }
inline int add(int a, int b) { return (a+=b)>=MOD?a-MOD:a; }
inline int dis(int a, int b) { return (a-=b)<0?a+MOD:a; }
const int Lim=4e4, MAXN=Lim+10;
int fc[MAXN], prime[MAXN], cntprime;
inline void sieve(){
for(int i=2;i<=Lim;++i){
if(!fc[i]) fc[i]=prime[++cntprime]=i;
for(int j=1;j<=cntprime;++j)
if(prime[j]>fc[i]||prime[j]*i>Lim) break;
else fc[ prime[j]*i ]=prime[j];
}
}
int n, k;
vector<int> divi;
inline void divit(int n){
divi.clear();
for(int i=2, I=sqrt(n)+1e-6;i<=I;++i) if(n%i==0) {
while(n%i==0) n/=i;
divi.push_back(i);
}
if(n>1)
divi.push_back(n);
}
inline int phi(int n){
int res=n;
for(auto p : divi)
if(n%p==0)
res=res/p*(p-1);
return res;
}
inline int f(int n, int k) {
return n&1?dis( fpow(k-1, n) , k-1):add( fpow(k-1, n) , k-1);
}
inline int ans(int n, int k){
int res=f(n, k-1);
divit(n);
for(int i=2, j, I=sqrt(n)+1e-6;i<=I;++i) if(n%i==0) {
j=n/i;
res=add(res, (ll)f(i, k-1)*phi(j)%MOD );
if(i!=j) res=add(res, (ll)f(j, k-1)*phi(i)%MOD );
}
return fpow(n, MOD-2)*k%MOD*res%MOD;
}
int main(){
ios::sync_with_stdio(0);
cin.tie(0); cout.tie(0);
sieve();
while(cin>>n>>k)
cout<<ans(n, k)<<"
";
cout.flush();
return 0;
}