「THUSCH 2017」如果奇迹有颜色（burnside引理+状压dp+打表+BM+常系数齐次线性递推）
这真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("\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;
}
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]);
}
}
