这真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]);
}
}