2019 ICPC Nanchang Summon
题目链接
https://vjudge.net/contest/418216#problem/J
这道题也是\(Burnside\)计数问题，我之前也写过一篇讲\(Burnside\)的博客，链接如下
https://www.cnblogs.com/ssdfzhyf/p/14533929.html
这道题说你需要拿\(0,1,2,3\)四种宝石排列成一个长度\(n\)的环形阵，环形阵无首尾，若一个阵旋转一定角度后和另一个阵重合，则视为同一个阵。
有\(m\)种非法子串，每个子串长度都是\(4\)。假如你的阵里面出现了一种子串，那么就是非法的。
求有多少种合法的阵，其中\(0\le m\le 256,4\le n\le 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},k\in[0,n)\cap N\}\)
那么根据\(Burnside\)定理
答案就是\(\frac{1}{n}\sum_{k=0}^{n-1}\sum_{a\in S(n)}[f_k(a)=a]\)
表达式\([f_k(a)=a]\)，等价于
\([\forall i\in [0,n)\cap N,a_{(i+k)mod\,n}=a_i]\)
那么我们可以把\({0,1,2,...,n-1}\)划分成若干等价类
\(\{\{x|x=(i+kj)mod\,n,j\in N\}|i\in[0,n)\cap N\}\)
对于\(i\)所在等价类里面的元素
\(\{x|x=(i+kj)mod\,n,j\in N\}\)
\(x\equiv i+kj(mod\,n)\)
根据贝祖定理，该方程有解的等价条件为\(x\equiv 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}\)
我们现在需要对这个数列进行计数
如果\(t\ge 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}\)
之后枚举\(0\le i,j< 64\)，查询\(ji\)(这是一个长度为\(6\)的字符串)中是否存在非法子串。
如果存在，则说明将该序列首尾相接后会出现非法子串，忽略即可。
如果不存在，则加上\(P[i][j]\)
如果\(t=3\)，要求有多少个满足首尾详解后不存在非法子串的数列\(a=a_0,a_1,a_2,...,a_0,a_1,a_2\)
由于题目规定\(n\ge4\)，且\(t=gcd(k,n)\)，那么\(n\ge 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.\(t\ge4\)
每个情况写一个函数进行计算
对于第\(4\)种，肯定是最复杂的
第\(4\)种涉及一个基础的矩阵\(ATM\)，它可以预处理。
然后还要预处理出\(P=ATM^{t-3}\)中，哪些转移是合法的，它也可以预处理
每次计算\(ATM^{t-3}\)，需要使用矩阵快速幂计算，每次矩阵乘法的计算次数为\(64^3\approx 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;
}
