Description：
题解：
这种东西肯定是burnside引理：
(sum置换后不动点数 over |置换数|)
一般来说，是枚举置换(i)，则(对所有x，满足a[x+i]=a[i])，然后a还要满足题目条件，但是仔细想一想，设(d=gcd(i,n))，只要a[0..d-1]满足就好了，所以：
(Ans=sum_{d|n}f(d)*phi(n/d)，f(d)表示)不考虑循环同构时的答案。
然后考虑dp：
枚举0这一列的块是什么，然后设(dp[i][1..4][1..4][1..3][1..3])，表示确定了后(i)列，第一行的颜色，和长度，第二行的颜色和长度，1、2表示第0列的颜色， 3、4表示不同于1、2列的颜色。
然后暴写dp，我死了……
n<=3的时候要特判，那么n更大的时候要矩阵乘法。
事实上有矩阵乘法的做法可以改成BM，最短递推式大概不超过35项
Code：
#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, B = y; i <= B; i ++)
#define ff(i, x, y) for(int i = x, B = y; i < B; i ++)
#define fd(i, x, y) for(int i = x, B = y; i >= B; i --)
#define ll long long
#define pp printf
#define hh pp("
")
using namespace std;
const int mo = 1e9 + 7;
ll ksm(ll x, ll y) {
ll s = 1;
for(; y; y /= 2, x = x * x % mo)
if(y & 1) s = s * x % mo;
return s;
}
namespace cz {
ll a[7] = {0, 0, 0, 24, 312, 1320, 3720};
ll cz(ll x) {
if(x <= 6) return a[x];
ll s = 0;
fo(i, 0, 6) {
ll xs = 1;
fo(j, 0, 6) if(i != j) xs = xs * (x - j) % mo * ksm(i - j, mo - 2) % mo;
s += xs * a[i];
}
return s % mo;
}
}
const int N = 200;
ll num, n, k, a[N];
ll f[N][5][5][4][4];
ll dp(int n, int cnt, int a0, int a1, int b0, int b1) {
ll m = k - cnt;
memset(f, 0, sizeof f);
f[n][1][cnt][a0][b0] = 1;
fd(i, n, 1) {
int c0 = 0, c1 = 0, c3 = 0, c4 = 0;
if(i - 1 < a0 || i - 1 >= n - a1) c0 = 1;
if(i - 1 < b0 || i - 1 >= n - b1) c1 = cnt;
if(i - 1 == a0 || i - 1 == n - a1 - 1) c3 = 1;
if(i - 1 == b0 || i - 1 == n - b1 - 1) c4 = cnt;
fo(x, 1, 4) fo(y, 1, 4) fo(lx, 1, 3) fo(ly, 1, 3) if(f[i][x][y][lx][ly]) {
ll F = f[i][x][y][lx][ly] % mo;
fo(u, 1, 3) {
ll xs = 1;
if(c0 && u != c0) xs = 0;
if(c1 && u != c1) xs = 0;
if(c3 && u == c3) xs = 0;
if(c4 && u == c4) xs = 0;
if(u == 2 && cnt == 1) xs = 0;
if(u <= 2 && (x == u || y == u)) xs = 0;
if(!xs) continue;
if(u == 3) {
xs *= m;
if(x >= 3) xs --;
if(y >= 3) xs --;
if(x == 3 && y == 3) xs ++;
}
xs = max(xs, 0ll);
f[i - 1][u][u][1][1] += xs * F % mo;
}
fo(u, 1, 3) fo(v, 1, 4) if(u != v && v != 3) {
ll xs = 1;
if(c0 && u != c0) xs = 0;
if(c1 && v != c1) xs = 0;
if(c3 && u == c3) xs = 0;
if(c4 && v == c4) xs = 0;
if((u == 2 || v == 2) && cnt == 1) xs = 0;
if(!xs) continue;
fo(r1, 0, 1) fo(r2, 0, 1) {
if(r1 && (u != x || lx == 3)) continue;
if(r2 && (v != y || ly == 3)) continue;
if(!r1 && u <= 2 && u == x) continue;
if(!r2 && v <= 2 && v == y) continue;
if(!r1 && !r2 && x != y) continue;
if(x == y && (r1 || r2)) continue;
int nlx = r1 ? lx + 1 : 1, nly = r2 ? ly + 1 : 1;
xs = 1;
if(u <= 2 || v <= 2 || r1 || r2) {
if(r1 && r2) {
xs = x != y;
} else
if(!r1 && !r2) {
xs *= (u <= 2 || r1) ? 1 : (m - (x >= 3));
xs *= (v <= 2 || r2) ? 1 : (m - (y >= 3));
} else {
if(r2) {
xs *= (u <= 2) ? 1 : (m - (x >= 3) - (y >= 3 && v >= 3));
} else {
xs *= (v <= 2) ? 1 : (m - (x >= 3 && u >= 3) - (y >= 3));
}
}
} else {
if(x <= 2 && y <= 2) xs *= m * (m - 1); else
if(x > 2 && y > 2) {
if(x == y) xs *= (m > 2) ? (m - 1) * (m - 2) : 0; else
xs *= (m > 3) ? (m - 2) * (m - 3) : 0;
} else {
xs *= (m > 1) ? (m - 1) * (m - 1) : 0;
}
}
xs = max(xs, 0ll);
xs %= mo;
f[i - 1][u][v][nlx][nly] += xs * F % mo;
}
}
}
}
ll ans = 0;
fo(x, 1, 4) fo(y, 1, 4) fo(lx, 1, 3) fo(ly, 1, 3) if(f[0][x][y][lx][ly]) {
ans = (ans + f[0][x][y][lx][ly]) % mo;
}
return ans;
}
void dp() {
for(int n = 4; n <= 70; n ++) {
a[n] = dp(n, 1, 1, 0, 1, 0) * k % mo;
fo(i, 1, 3) fo(j, 1, i) fo(u, 1, 3) fo(v, 1, u)
a[n] += dp(n, 2, j, i - j, v, u - v) * k % mo * (k - 1) % mo;
a[n] %= mo;
}
// pp("%lld
", a[5]);
a[1] = 0;
a[2] = (k * (k - 1) + k * (k - 1) % mo * (k - 2) * 2) % mo;
a[3] = cz :: cz(k);
}
int m; ll b[N];
namespace bm {
ll f[3][N], len[3], fail[3], delta[3];
void cp(int x, int y) {
fo(i, 1, len[x]) f[y][i] = f[x][i];
len[y] = len[x]; fail[y] = fail[x]; delta[y] = delta[x];
}
void solve(ll *a, int n) {
fo(i, 1, n) f[0][i] = 0; len[0] = 0;
fo(i, 1, n) {
ll tmp = a[i];
fo(j, 1, len[0]) tmp -= f[0][j] * a[i - j] % mo;
tmp = (tmp % mo + mo) % mo;
delta[0] = tmp;
if(!tmp) continue;
fail[0] = i;
if(!len[0]) {
cp(0, 1); f[0][len[0] = 1] = 0;
continue;
}
cp(0, 2);
ll mul = delta[0] * ksm(delta[1], mo - 2) % mo;
int st = i - fail[1]; len[0] = max(len[0], st + len[1]);
f[0][st] = (f[0][st] + mul) % mo;
fo(j, 1, len[1]) f[0][st + j] = (f[0][st + j] + (mo - mul) * f[1][j]) % mo;
if(len[2] - fail[2] < len[1] - fail[1]) cp(2, 1);
}
m = len[0] + 1; b[m] = 1;
fo(i, 1, len[0]) b[i] = -f[0][m - i];
}
}
ll c[N * 2], d[N * 2], e[N * 2];
void qmo(ll *a) {
fo(i, 0, 2 * m) a[i] %= mo;
fd(i, 2 * m, m) if(a[i]) {
ll v = a[i]; int c = i - m;
fo(j, 1, m) a[c + j] = (a[c + j] - v * b[j]) % mo;
}
}
ll solve(int n) {
if(n <= 35) return a[n];
fo(i, 0, 2 * m) c[i] = d[i] = 0;
c[0] = 1; d[1] = 1;
int y = n;
for(; y; y /= 2) {
if(y & 1) {
fo(i, 0, 2 * m) e[i] = c[i], c[i] = 0;
fo(i, 0, m) fo(j, 0, m) c[i + j] += e[i] * d[j] % mo;
qmo(c);
}
fo(i, 0, 2 * m) e[i] = d[i], d[i] = 0;
fo(i, 0, m) fo(j, 0, m) d[i + j] += e[i] * e[j] % mo;
qmo(d);
}
ll ans = 0;
fo(i, 1, m - 1) ans += c[i] * a[i] % mo;
return ans % mo;
}
ll phi(int n) {
ll s = n;
for(int i = 2; i * i <= n; i ++) if(n % i == 0) {
s = s / i * (i - 1);
for(; n % i == 0; n /= i);
}
if(n > 1) s = s / n * (n - 1);
return s;
}
int main() {
freopen("experience.in", "r", stdin);
freopen("experience.out", "w", stdout);
scanf("%lld %lld %lld", &num, &n, &k);
dp();
bm :: solve(a, 70);
ll ans = 0;
for(int d = 1; d * d <= n; d ++) if(n % d == 0) {
ans += solve(d) * phi(n / d) % mo;
if(d != n / d) ans += solve(n / d) * phi(d) % mo;
}
ans = (ans % mo + mo) % mo * ksm(n, mo - 2) % mo;
pp("%lld
", ans);
}