【NOI2019模拟2019.7.4】朝夕相处 （动态规划+BM）

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("\n")
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\n", 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\n", ans);
}