「THUSCH 2017」如果奇迹有颜色（burnside引理+状压dp+打表+BM+常系数齐次线性递推）

https://loj.ac/problem/2981

​

这真TM是个防AK题。

​

考虑先套上burnside，枚举置换(i)，发现问题变成(gcd(i,n))长的环，染m种颜色，连续m个不会出现m个颜色的方案数，记(f(d))表示长为d的环的方案数，则(Ans=sum_{d|n}f(d)*phi(n/d))

​

现在就是怎么求一个(f(n))的问题了。

​

暴力：

状态压缩dp，记(dp[i][S])表示确定了前(i)个，第(i-m+2)到(i)个颜色状态是(S)，的方案数。

至于环的限制，只要(n)步之后(S)等于一开始的(S)就好了。

用矩阵乘法优化，时间复杂度：((m^{m-1})^3*log~n)

可以通过(m<=4)的数据获得55分。

​

BM：

当(m)确定时，(f[n])肯定是一个有递推式的序列，且由上一个dp知道递推式长度(le m^{m-1})。

这样就可以BM后常系数齐次线性递推。

那么问题就是怎么快速算(f[1..n])。

考虑不用矩阵乘法，直接dp：
先枚举一个(w)，表示前(w)个的颜色是(1,2,3,...,w)，第(w+1)个颜色(in[1,w])。

然后用刚才那个(dp[i][S])，在(n-w)步后，看看(S)和前w个是否冲突，不冲突，则(f[n]+=dp[n-w][S]*P_m^w)

这样算(f[1..n])的复杂度是(O(n*(m^{m-1})*m^2))，(m=7,n=1000)大概要算20s，本地算出，打表即可。

接着跑BM，实际上，当(m=7)时，递推式长度也才(410)。

这样就不用多项式取模，直接暴力取模的常系数齐次线性递推即可。

​

代码（打表后的BM和递推）：https://loj.ac/submission/782204

代码（打表程序）：

#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;
}

const int M = 117649;

int n, m, t, am[8];


int bz[M][8], to[M][8];

namespace sub1 {
	int cx[10];
	void build1() {
		ff(i, 0, t)	{
			ff(j, 0, m) {
				ff(k, 0, m) cx[k] = 0; cx[j] = 1;
				ff(k, 0, m - 1) cx[i / am[k] % m] = 1;
				int ky = 0;
				ff(k, 0, m) if(!cx[k]) ky = 1;
				bz[i][j] = ky;
				if(ky) to[i][j] = (i / m) + j * am[m - 2];
			}
		}
	}
}
	

int po[M][8];

namespace sub2 {
	int c[20], cx[10];
	int pd(int n) {
		fo(i, 1, n - m + 1) {
			fo(j, 0, m - 1) cx[j] = 0;
			fo(j, 0, m - 1) cx[c[i + j]] = 1;
			int ky = 0;
			fo(j ,0, m - 1) if(!cx[j]) ky = 1;
			if(!ky) return 0;
		}
		return 1;
	}
	void build2() {
		ff(s, 0, am[m - 1]) {
			fo(i, 1, m - 1) c[i] = s / am[i - 1] % m;
			fo(w, 1, m - 1) {
				fo(j, m, m + w - 1) c[j] = j - m;;
				po[s][w] = pd(m + w - 1);
			}
		}
	}
}

ll fac[10], nf[10];
ll P(int n, int m) { return fac[n] * nf[n - m] % mo;}

ll a[2005], f[2][M], o;

int main() {
	freopen("a.out", "w", stdout);
	for(m = 0; m <= 7; m ++) {
		fo(i, 0, 1000) a[i] = 0;
		am[0] = 1; fo(i, 1, m) am[i] = am[i - 1] * m;
		t = am[m - 1];
		sub1 :: build1();
		sub2 :: build2();
		fo(i, 0, m - 1) a[i] = ksm(m, i);
		fac[0] = 1; fo(i, 1, m) fac[i] = fac[i - 1] * i;
		fo(i, 0, m) nf[i] = ksm(fac[i], mo - 2);
		fo(w, 1, m - 1) {
			memset(f, 0, sizeof f);
			int st = 0;
			fo(i, 0, w - 1) st += am[(m - 1 - w) + i] * i;
			f[o][st] = 1;
			fo(p, w + 1, 1000) {
				memset(f[!o], 0, sizeof f[!o]);
				ff(i, 0, t) if(f[o][i]) {
					fo(j, 0, (p == w + 1 ? w - 1 : m - 1)) if(bz[i][j]) {
						f[!o][to[i][j]] = (f[!o][to[i][j]] + f[o][i]) % mo;
					}
				}
				o = !o;
				if(p >= m) {
					ff(i, 0, t) if(po[i][w]) {
						a[p] = (a[p] + P(m, w) * f[o][i]) % mo;
					}
				}
			}
		}
		pp("{");
		fo(i, 0, 1000) {
			pp("%lld", a[i]);
			if(i != 1000) pp(",");
		}
		pp("},");
		hh;
//		fo(i, 1, 100) pp("%lld ", a[i]);
	}
	
}