应该能想到polya定理，考虑一下polya定理的公式是怎么一个表达。
(f)是一个置换，(m(f))是置换的循环的数量。
因为有对限制，所以这就导致了这道题不能直接做，我们从原理入手。
考虑置换群(G={0,1,...,n-1})。我们知道对于一个置换(L)，有循环数量为(gcd(L,n))，而且第一个元素在第一个循环中，第二个元素在第二个循环中，...，第(gcd(L,n))个元素在第(gcd(L,n))个循环中，而且每个循环中的元素颜色相同。
所以问题可以抽象为一个链条长度为(gcd(L,n))，且在上面染色的计数问题。
接下来我们处理这样的一个问题。
如果考虑是两个限制的话，即((a,b))表示这两个颜色不能相邻，实际题目中是四个，先考虑简单的。
我们可以用一个矩阵(m(i,j)=1)表示(i,j)两色可以相邻，否则不行。
然后我们可以写一个(dp)方程，设(f(i,j))表示处理到前(i)位最后一位颜色是(j)的答案，转移方程：
当然要保证最后一位的颜色和第一位的颜色满足条件，所以我们可以取：
为答案，就相当于是枚举第一个元素的颜色(j)，你在尾部接一个(j)，到这里一定是满足条件的。
为什么说这个呢，在后面的推导中用到了这样的(dp)思想，但实际上用的是矩阵快速幂。
在离散数学中，如果用(0,1)矩阵表示无向图连通情况的话，(A^k(i,j))次方代表的就是一个点经过(k)条路（也就是(k)次转移，到了第(k+1)位）后(i)点到(j)点的方案数，又由于我们要起始元素颜色和终止元素颜色相同。
- 所以矩阵此时的对角线之和就是答案。
考虑一个四元组((a,b,c,d))。
如果第一位是(a)且后面跟着(b,c)，那么第二位就不能是(b)且后面跟着(c,d)。
这样就抽象到二维了，也就是上面的那个问题。
所以我们就设矩阵(m(i,j))表示(i)后面跟着(x_1,x_2)，(j)后面跟着(y_1,y_2)是否可行。
这样的(i,j)就分别表示一个三元组，跑矩阵快速幂就是答案。
设对于长度为(i)的链条，方案数为(S(i))。
那么我的答案就是：
枚举(gcd(L,n)=d)。
预处理矩阵的阶乘，在快速幂的时候不要频繁的调动(a=a*a\%mod)。
最后别忘了除以(n)，直接乘上(n)在(mod 998244353)意义下的逆元。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 998244353;
const int maxn = 1e5+10;
int primes[maxn], cnt;
int phi[maxn];
bool v[maxn];
void init(int n)
{
phi[1] = 1;
for(int i = 2; i <= n; i++)
{
if(!v[i])
{
primes[++cnt] = i;
phi[i] = i-1;
}
for(int j = 1; primes[j] <= n/i; j++)
{
v[primes[j]*i] = 1;
if(i%primes[j] == 0)
{
phi[primes[j]*i] = phi[i]*primes[j];
break;
}
else phi[primes[j]*i] = phi[i]*(primes[j]-1);
}
}
}
struct mat
{
ll a[67][67];
}base, pbase[25];
mat mul(mat a, mat b)
{
mat res;
memset(res.a, 0, sizeof res.a);
for(int i = 0; i < 64; i++)
for(int j = 0; j < 64; j++)
for(int k = 0; k < 64; k++)
res.a[i][j] = (res.a[i][j]%mod+a.a[i][k]%mod*b.a[k][j]%mod)%mod;
return res;
}
ll qmi(ll a, ll b)
{
ll res = 1;
while(b)
{
if(b&1) res = (res*a)%mod;
a = a*a%mod;
b >>= 1;
}
return res%mod;
}
mat qmi(int b)
{
mat res; memset(res.a, 0, sizeof res.a);
for(int i = 0; i < 64; i++) res.a[i][i] = 1;
int t = 0;
while(b)
{
if(b&1) res = mul(res, pbase[t]);
b >>= 1; t++;
}
return res;
}
int main()
{
init(maxn-5);
int n, m;
scanf("%d%d", &n, &m);
for(int i = 0; i < 64; i++)
for(int k = 0; k < 4; k++)
{
int nxt = (i*4+k)%64;
base.a[i][nxt] = 1;
}
while(m--)
{
int x = 0, t;
for(int i = 0; i < 4; i++)
{
scanf("%d", &t);
x = x*4+t;
}
base.a[x/4][x%64] = 0;
}
pbase[0] = base;
for(int i = 1; i < 20; i++) pbase[i] = mul(pbase[i-1], pbase[i-1]);
ll ans = 0;
for(int i = 1; i <= n; i++)
{
if(n%i) continue;
mat tt = qmi(i);
ll t = 0;
for(int j = 0; j < 64; j++)
t = (t+tt.a[j][j])%mod;
t = t*phi[n/i]%mod;
ans = (ans+t)%mod;
}
ans = ans%mod*qmi(n, mod-2)%mod;
cout << ans << endl;
return 0;
}