2019 ICPC Asia Nanchang Regional J. Summon(polya定理+矩阵快速幂优化dp)-CSDN博客

本文链接：https://blog.csdn.net/Code92007/article/details/106690998

本文探讨了一道复杂的算法题目，涉及项链染色方案计数问题。在考虑到旋转等价性和特定颜色序列限制的情况下，利用矩阵快速幂和数论技巧求解合法方案的数量。文章深入解析了解题思路，包括置换群理论的应用、矩阵乘法优化和倍增/块速递推技巧，并提供了两种实现代码示例。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题目

n(4<=n<=1e5)个珠子的项链，项链由四种颜色0123组成，

旋转相同时视为同种方案，

m(0<=m<=256)个限制，第i次给出四个整数a b c d(0<=a,b,c,d<=3)

表示顺时针看这串项链时，abcd段不能出现，

求合法的方案数%998244353

思路来源

https://www.cnblogs.com/zxytxdy/p/12582065.html

题解

poj2888原题既视感，高度相似的一个题，只是变二维为四维了

旋转步长为i时，每个置换里有n/gcd(i,n)，共gcd(i,n)个置换

即在环上每gcd(i,n)个就会遇到相同的颜色一次，

那在一个长度为n的环上无限制的走一圈，相当于同样地在一个长为gcd(i,n)的一段上走了n/gcd(i,n)次

而如果把走向下一段看成是走向这段的开头（因为颜色相同），可以认为是一个gcd(i,n)的环走了n/gcd(i,n)次

则对旋转步长为i时的答案求和时，i的贡献是f(gcd(i,n))的值，f()表示不考虑旋转时，禁止掉这些串的方案数

$\sum_{i=0}^{n-1}f(gcd(i,n))=\sum_{i=1}^{n}f(gcd(i,n))=\sum_{d=1}^{n}\sum_{i=1}^{n}f(d)*[gcd(i,n)==d]=\sum_{d=1}^{n}\sum_{\frac{i}{d}=1}^{\frac{n}{d}}f(d)*[gcd(\frac{i}{d},\frac{n}{d})==1]=\sum_{d=1}^{n}f(d)*\varphi (\frac{n}{d})$

考虑怎么求f(i)，设长度为i+2的且末三位数字为j,k,l的方案数为dp[i][j][k][l]

显然dp[i][j][k][l]可以向dp[i+1][k][l][m]转移，考虑j、k、l都在[0,3]之间，

把三个四进制的压维成一个64进制的，矩阵快速幂一波，

如果是一条链的话，转移i-2次即可

如果是一个环的话，在初始为s1s2s3的链上走i次，仍然要回到s1s2s3

因此，矩阵快速幂i次，最后统计系数矩阵中c[s1s2s3][s1s2s3]的即可（因为初始向量的值均为1）

复杂度 $O(64^3*logn*\sigma(n))$ ，其中因子个数 $\sigma(n)$ （枚举d）大致在 $\sqrt n$ 级别，

极限数据是83160(对应约数128)，4s难以通过

于是要卡卡常，

①矩阵乘法，c[i][k]=0的就没必要算下去了

②倍增：预处理系数矩阵的二进制次方数，这样在快速幂里x=x*x就不用次次现算了

③块速递推： $sz=\sqrt n$ ， $x=i*\sqrt n +j(i<\sqrt n, j<\sqrt n)$ ，

预处理 $i*\sqrt n$ 和 $j$ 的矩阵快速幂，复杂度 $O(64^3*\sqrt n)$

②/③选一即可，虽然理论复杂度③更优，但数据表现上②更优一些

题解中提到f(1)和f(2)要特判，但实际不必，因为1次幂的时候，只有形如111的能转移到111，2次幂同理

代码1（预处理倍增）

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
#define rep(i,a,b) for(int i=(a);i<=(b);++i)
#define sci(a) scanf("%d",&(a))
typedef long long ll;
ll modpow(ll x,ll n,ll mod){ll res=1;for(;n;n>>=1,x=x*x%mod)if(n&1)res=res*x%mod;return res;}
const int mod=998244353;
const int M=400,maxn=1e5+10;
bool ok[maxn];
int prime[maxn],phi[maxn],cnt;
void sieve(){
    phi[1]=1;
	for(ll i=2;i<maxn;++i){
		if(!ok[i]){
			prime[cnt++]=i;
			phi[i]=i-1;
		}
		for(int j=0;j<cnt;++j){
			if(i*prime[j]>=maxn)break;
			ok[i*prime[j]]=1;
			if(i%prime[j]==0){
				phi[i*prime[j]]=phi[i]*prime[j];//prime[j]是i的因子 prime[j]的素因子项包含在i的素因子项里
				break;
			}
			else phi[i*prime[j]]=phi[i]*(prime[j]-1);//prime[j]与i互质 phi[i*prime[j]=phi[i]*phi[prime[j]]
		}
	}
}
struct mat{
    static const int N=64;
    ll c[N][N];
    int m,n;
    mat(){
    	memset(c,0,sizeof(c));
    	m=n=N;
    }
    mat(int a,int b):m(a),n(b){
        memset(c,0,sizeof(c));
    }
    void clear(){
		memset(c,0,sizeof(c));
    }
    void E(){
        int mn=min(m,n);
        for(int i=0;i<mn;i++){
            c[i][i]=1;
        }
    }
    mat operator *(const mat& x){
        mat ans(m,x.n);
        for(int i=0;i<m;i++)
            for (int k=0;k<n;k++)
                {
                    if(!c[i][k])continue;
                    for (int j=0;j<x.n;j++)
                    (ans.c[i][j]+=c[i][k]*x.c[k][j]%mod)%=mod;//能不取模 尽量不取模
                }
                    //这里maxn=2 故不会超过ll 视具体情况 改变内部取模情况
        return ans;
    }
}bs,tmp,pbs[22];
mat modpow(mat x,ll n){//幂次一般为ll
        mat ans(x.m, x.m);
    	ans.E();
        int c=0;
        for(;n;n>>=1,c++){//x自乘部分可以预处理倍增，也可以用分块加速递推
            if(n&1)ans=ans*pbs[c];
        }
        return ans;
	}
int n,m,sz,v;
ll ans,now;
void init(){
    sieve();
    sci(n),sci(m);
    rep(i,0,63){
        rep(j,0,3){
            bs.c[i][(i%16)*4+j]=1;//201 -> 013
        }
    }
    rep(i,1,m){
        int x=0;
        rep(j,0,3){
            sci(v);
            x=x*4+v;
        }
        bs.c[x/4][x%64]=0;
    }
    sz=(int)(sqrt(n)+0.5);
    pbs[0]=bs;
    rep(i,1,19){
        pbs[i]=pbs[i-1]*pbs[i-1];
    }
}
int main(){
    init();
    rep(i,1,n){
        if(n%i)continue;
        tmp=modpow(bs,i);
        now=0;
        rep(j,0,63){
            now=(now+tmp.c[j][j])%mod;
        }
        now=now*phi[n/i]%mod;
        ans=(ans+now)%mod;
    }
    printf("%lld\n",ans*modpow(n,mod-2,mod)%mod);
    return 0;
}

代码2（分块加速）

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
#define rep(i,a,b) for(int i=(a);i<=(b);++i)
#define sci(a) scanf("%d",&(a))
typedef long long ll;
ll modpow(ll x,ll n,ll mod){ll res=1;for(;n;n>>=1,x=x*x%mod)if(n&1)res=res*x%mod;return res;}
const int mod=998244353;
const int M=400,maxn=1e5+10;
bool ok[maxn];
int prime[maxn],phi[maxn],cnt;
void sieve(){
    phi[1]=1;
	for(ll i=2;i<maxn;++i){
		if(!ok[i]){
			prime[cnt++]=i;
			phi[i]=i-1;
		}
		for(int j=0;j<cnt;++j){
			if(i*prime[j]>=maxn)break;
			ok[i*prime[j]]=1;
			if(i%prime[j]==0){
				phi[i*prime[j]]=phi[i]*prime[j];//prime[j]是i的因子 prime[j]的素因子项包含在i的素因子项里
				break;
			}
			else phi[i*prime[j]]=phi[i]*(prime[j]-1);//prime[j]与i互质 phi[i*prime[j]=phi[i]*phi[prime[j]]
		}
	}
}
struct mat{
    static const int N=64;
    ll c[N][N];
    int m,n;
    mat(){
    	memset(c,0,sizeof(c));
    	m=n=N;
    }
    mat(int a,int b):m(a),n(b){
        memset(c,0,sizeof(c));
    }
    void clear(){
		memset(c,0,sizeof(c));
    }
    void E(){
        int mn=min(m,n);
        for(int i=0;i<mn;i++){
            c[i][i]=1;
        }
    }
    mat operator *(const mat& x){
        mat ans(m,x.n);
        for(int i=0;i<m;i++)
            for (int k=0;k<n;k++)
                {
                    if(!c[i][k])continue;
                    for (int j=0;j<x.n;j++)
                    (ans.c[i][j]+=c[i][k]*x.c[k][j]%mod)%=mod;//能不取模 尽量不取模
                }
                    //这里maxn=2 故不会超过ll 视具体情况 改变内部取模情况
        return ans;
    }
    friend mat operator^(mat x,ll n){//幂次一般为ll
        mat ans(x.m, x.m);
    	ans.E();
        for(;n;n>>=1,x=x*x){//x自乘部分可以预处理倍增，也可以用分块加速递推
            if(n&1)ans=ans*x;
        }
        return ans;
	}
}bs,l[M],h[M],tmp;
int n,m,sz,v;
ll ans,now;
void init(){
    sieve();
    sci(n),sci(m);
    rep(i,0,63){
        rep(j,0,3){
            bs.c[i][(i%16)*4+j]=1;//201 -> 013
        }
    }
    rep(i,1,m){
        int x=0;
        rep(j,0,3){
            sci(v);
            x=x*4+v;
        }
        bs.c[x/4][x%64]=0;
    }
    sz=(int)(sqrt(n)+0.5);
    l[0].E();//low high 分块加速
    rep(i,1,sz){
        l[i]=l[i-1]*bs;
    }
    h[0].E();
    rep(i,1,sz){
        h[i]=h[i-1]*l[sz];
    }
}
int main(){
    init();
    rep(i,1,n){
        if(n%i)continue;
        tmp=h[i/sz]*l[i%sz];
        now=0;
        rep(j,0,63){
            now=(now+tmp.c[j][j])%mod;
        }
        now=now*phi[n/i]%mod;
        ans=(ans+now)%mod;
    }
    printf("%lld\n",ans*modpow(n,mod-2,mod)%mod);
    return 0;
}