题目链接
https://vjudge.net/contest/418216#problem/J
这道题也是(Burnside)计数问题，我之前也写过一篇讲(Burnside)的博客，链接如下
https://www.cnblogs.com/ssdfzhyf/p/14533929.html
这道题说你需要拿(0,1,2,3)四种宝石排列成一个长度(n)的环形阵，环形阵无首尾，若一个阵旋转一定角度后和另一个阵重合，则视为同一个阵。
有(m)种非法子串，每个子串长度都是(4)。假如你的阵里面出现了一种子串，那么就是非法的。
求有多少种合法的阵，其中(0le mle 256,4le nle 100000)
(Burnside)问题需要定义一个集合和对应的置换群
设一个长度为(n)的数列的集合(S(n)={a|a=a_0,a_1,...,a_{n-1},将a首尾连成环形后不出现非法的子串})
注意(a=0,1,2,3)与(b=1,2,3,0)是两个不同的元素，因为数列是有序的
设(S(n))上的置换(G)群
(G={f_k|b=f_k(a):b_i=a_{(i+k)mod\,n},kin[0,n)cap N})
那么根据(Burnside)定理
答案就是(frac{1}{n}sum_{k=0}^{n-1}sum_{ain S(n)}[f_k(a)=a])
表达式([f_k(a)=a])，等价于
([forall iin [0,n)cap N,a_{(i+k)mod\,n}=a_i])
那么我们可以把({0,1,2,...,n-1})划分成若干等价类
({{x|x=(i+kj)mod\,n,jin N}|iin[0,n)cap N})
对于(i)所在等价类里面的元素
({x|x=(i+kj)mod\,n,jin N})
(xequiv i+kj(mod\,n))
根据贝祖定理，该方程有解的等价条件为(xequiv i (mod\,gcd(k,n)))
那么共有(gcd(n,k))个等价类，每个等价类都有(frac{n}{gcd(k,n)})个元素
故对应着循环节为(gcd(k,n))的一个序列
不妨记(t=gcd(k,n))
那么满足条件的序列形如(a=a_0,a_1,...,a_{t-1},a_0,a_1,...,a_{t-1},...,a_0,a_1,...,a_{t-1})
我们现在需要对这个数列进行计数
如果(tge 4)，那么序列(a)首尾相连后不出现非法子串等价于(a_0,a_1,...,a_{t-1})首尾相连后不出现非法子串
那么我们可以统计集合(S(t))中有几个元素即可
(S(t))内的元素也是数列，是有顺序的，我们可以利用字符串自动机的思想进行统计
我们先计算(T(t)={a|a=a_0,a_1,...,a_{t-1},a中不出现非法的子串})的元素数量
我们知道非法字符串都是(4)位的，故我们可以把状态设置为长度为(3)的字符串，共(64)个状态
由于每个状态对应的长度都是(3)，转移的初始状态有(3)个字符，每一步转移给这个字符串末尾添加一个新字符，一共转移(t-3)次。
比如非法字符串只有(0123)一种，那么对于状态(012)，下一个状态可以是(120,121,122)
通过枚举状态和转移时新加入的字符，可以得到一个(64*64)的矩阵(ATM)
(ATM[i][j]=1)就表示从状态(i)添加一个字符可以到达状态(j)，且不新增违规字符串
(ATM[i][j]=0)则相反
计算出(P=ATM^{t-3})
(P[i][j])就表示如果这个字符串的长度为(3)的前缀对应状态(i)，且这个字符串的长度为(3)的后缀对应状态(j)，在考虑中间没有违规字符串的前提下，有多少种合法的序列
这样，将(P[i][j])直接求和即可算出不考虑数列首尾相接的情况，有多少个合法的串
如果考虑首尾相接，那么仍先计算(P=ATM^{t-3})
之后枚举(0le i,j< 64)，查询(ji)(这是一个长度为(6)的字符串)中是否存在非法子串。
如果存在，则说明将该序列首尾相接后会出现非法子串，忽略即可。
如果不存在，则加上(P[i][j])
如果(t=3)，要求有多少个满足首尾详解后不存在非法子串的数列(a=a_0,a_1,a_2,...,a_0,a_1,a_2)
由于题目规定(nge4)，且(t=gcd(k,n))，那么(nge 6)
那么它等价于(a_0,a_1,a_2,a_0)与(a_1,a_2,a_0,a_1)与(a_2,a_0,a_1,a_2)都不是非法子串
我们可以构造串(a_0,a_1,a_2,a_0,a_1,a_2)，在不考虑首尾相接的情况下有多少种合法的情况，枚举即可
如果(t=2)，要求有多少个满足首尾详解后不存在非法子串的数列(a=a_0,a_1,...,a_0,a_1)
那么它等价于(a_0,a_1,a_0,a_1)与(a_1,a_0,a_1,a_0)都不是非法子串
我们可以构造(a_0,a_1,a_0,a_1,a_0)，在不考虑首尾相接的情况下有多少种合法的情况，枚举即可
在确定(t)后，以上的方案数简记为(f(t))
那么答案为(frac{1}{n}sum_{k=0}^{n-1}f(gcd(k,n))=frac{1}{n}sum_{t|n}f(t)sum_{k=0}^{n-1}[t=gcd(k,n)]=frac{1}{n}sum_{t|n}f(t)varphi(frac{n}{t}))
编程实现
我们要计算(frac{1}{n}sum_{t|n}f(t)varphi(frac{n}{t}))，首先需要枚举(n)的所有因数，求出它们对应的欧拉函数，以及(f(t))
在求(f(t))时，要分(4)种情况讨论
1.(t=1)
2.(t=2)
3.(t=3)
4.(tge4)
每个情况写一个函数进行计算
对于第(4)种，肯定是最复杂的
第(4)种涉及一个基础的矩阵(ATM)，它可以预处理。
然后还要预处理出(P=ATM^{t-3})中，哪些转移是合法的，它也可以预处理
每次计算(ATM^{t-3})，需要使用矩阵快速幂计算，每次矩阵乘法的计算次数为(64^3approx 2.6e5)。
经过实际检测，当(n=98280)时，需要进行(1743)次的矩阵乘法运算，这是上限
那么最终的计算次数约(4.5e8)，在(4s)内应该可以通过
放代码
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#include<stack>
#include<vector>
using namespace std;
typedef long long ll;
const ll mod=998244353;
int prime[100010],first[100010],phi[100010];
bool isp[100010];
int euler(int n)//计算欧拉函数
{
int cnt=0;
isp[1]=0;
phi[1]=1;
for(int i=2;i<=n;i++)isp[i]=1;
for(int i=2;i<=n;i++)
{
if(isp[i])
{
prime[++cnt]=i;
first[i]=i;
phi[i]=i-1;
}
for(int j=1;j<=cnt&&prime[j]*i<=n;j++)
{
isp[i*prime[j]]=0;
first[i*prime[j]]=prime[j];
if(i%prime[j])
{
phi[i*prime[j]]=phi[i]*(prime[j]-1);
}
else
{
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
}
}
}
struct mat
{
vector<ll>a[65];
int n,m;
void init(int n_,int m_)//矩阵的初始化
{
n=n_;m=m_;
for(int i=0;i<n;i++)
{
a[i].resize(m);
for(int j=0;j<m;j++)a[i][j]=0;
}
}
void I(int n_)//生成单位矩阵
{
init(n_,n_);
for(int i=0;i<n;i++)a[i][i]=1;
}
mat operator *(const mat &b)const//模意义下矩阵乘法
{
mat res;
res.init(n,b.m);
for(int i=0;i<res.n;i++)
{
for(int j=0;j<res.m;j++)
{
for(int k=0;k<m;k++)
{
res.a[i][j]=(res.a[i][j]+a[i][k]*b.a[k][j])%mod;
}
}
}
return res;
}
}ATM;
mat qpow(mat A,int n)//矩阵快速幂
{
mat ans;
ans.I(A.n);
while(n)
{
if(n&1)ans=ans*A;
A=A*A;
n>>=1;
}
return ans;
}
int str[260][6];//非法子串
int list[200];//因数表
//pattern[i]=1就代表子串i是一个非法子串，i是一个8bit的整数，比如0123即1*16+2*4+3=27
bool pattern[260];
bool check(int a[],int len)//检查a[0]~a[len-1]中间有没有非法子串
{
int now=0;
for(int i=0;i<len;i++)
{
now=(now*4+a[i])%256;
if(i>=3&&pattern[now])return 1;
}
return 0;
}
int f_1()
{
int a[4];
int ans=0;
for(a[0]=0;a[0]<4;a[0]++)
{
a[3]=a[2]=a[1]=a[0];
if(!check(a,4))ans++;
}
return ans;
}
int f_2()
{
int a[6];
int ans=0;
for(a[0]=0;a[0]<4;a[0]++)
{
for(a[1]=0;a[1]<4;a[1]++)
{
a[4]=a[2]=a[0];
a[5]=a[3]=a[1];
if(!check(a,6))ans++;
}
}
return ans;
}
int f_3()
{
int a[6];
int ans=0;
for(a[0]=0;a[0]<4;a[0]++)
{
for(a[1]=0;a[1]<4;a[1]++)
{
for(a[2]=0;a[2]<4;a[2]++)
{
a[3]=a[0];
a[4]=a[1];
a[5]=a[2];
if(!check(a,6))ans++;
}
}
}
return ans;
}
ll trans[70][70];//trans[i][j]表示P[i][j]的转移是合法的
ll f(int d)
{
mat P;
P=qpow(ATM,d-3);//计算P
ll ans=0;
for(int i=0;i<64;i++)
for(int j=0;j<64;j++)
ans=(ans+P.a[i][j]*trans[i][j])%mod;
return ans;
}
ll qpow(ll a,ll n)
{
ll ans=1;
while(n)
{
if(n&1)ans=ans*a%mod;
a=a*a%mod;
n>>=1;
}
return ans;
}
ll inv(ll x)
{
return qpow(x,mod-2);
}
void solve(int n)
{
int sqr=sqrt(n);
int tot=0;
for(int i=1;i<=sqr;i++)//把因数放到列表里面
{
if(n%i==0)list[++tot]=i;
}
for(int i=sqr;i>0;i--)
{
if(n%i==0&&n/i!=sqr)list[++tot]=n/i;
}
ll ans=0;
for(int k=1;k<=tot;k++)
{
int i=list[k];
ll func=0;
if(i==1)func=f_1();
else if(i==2)func=f_2();
else if(i==3)func=f_3();
else func=f(i);
ans=(ans+func*phi[n/i])%mod;//求和
}
ans=ans*inv(n)%mod;
printf("%lld",ans);
}
int main()
{
euler(100000);
int n,m;
scanf("%d%d",&n,&m);
for(int i=1;i<=m;i++)
{
for(int j=1;j<=4;j++)scanf("%d",&str[i][j]);
int tmp=0;
for(int j=1;j<=4;j++)tmp=tmp*4+str[i][j];
pattern[tmp]=1;
}
ATM.init(64,64);
for(int i=0;i<64;i++)//建立一步转移矩阵
{
for(int ch=0;ch<4;ch++)
{
int now=i*4+ch;
if(pattern[now])continue;
else ATM.a[i][now%64]++;
}
}
for(int i=0;i<64;i++)//检测P[i][j]是否合法
{
for(int j=0;j<64;j++)
{
int a[6];
int num=j*64+i;
for(int k=5;k>=0;k--)
{
a[k]=num%4;
num>>=2;
}
if(!check(a,6))trans[i][j]++;
}
}
solve(n);
return 0;
}