[SDOI2017]序列计数（矩阵乘法）_矩阵列计数-优快云博客

本文链接：https://blog.youkuaiyun.com/zhangjianjunab/article/details/109059235

本文详细介绍了如何利用数学和编程技巧，特别是矩阵快速幂算法，来解决关于质数和无质数计数的问题。通过容斥原理简化思考，并给出了一种优化的解决方案，降低时间复杂度到O(p^3 log n)。同时，文章探讨了在无法立即找到最优解时，如何从复杂解法入手逐步优化。

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

题目

做法

~~这都没想到~~

首先，我们需要考虑有质数的，和没质数的情况。

但是我们发现，有质数的情况不可能用来推导没质数的情况，且总情况减去没有质数的情况就是有质数的情况，且有质数的情况需要用没质数的情况推导而成。

所以直接用处理无质数的情况和总情况（事实上，总情况总比无质数情况好处理，而无质数情况比有质数情况好处理），然后容斥求得即可，以此减少思考难度以及思考总量（事实上，用容斥一般可以减少常数甚至是复杂度，本题减少的是常数）。

现在说总情况的统计方法：

考虑DP， $f [i] [j]$ ，表示第 $i$ 个位置和模 $p$ 为 $j$ 的情况个数，然后通过转移方程便可以得到，设 $cnt_{i}$ 为不超过 $m$ 的且模 $p$ 等于 $i$ 的正整数个数，即： $f[i+1][(k+j)\%p]+=f[i][j]*cnt_k$ ，不难发现，这个转移可以写成矩阵的形式：
$\begin{bmatrix} cnt_0 &cnt_{p-1} & cnt_{p-2} & ...\\ cnt_1 & cnt_0 & cnt_{p-1} & ...\\ cnt_2 & cnt_1 & cnt_0 & ...\\ ... & ... & ... & ... \end{bmatrix}$

把 $f [i]$ 看成 $1 * p$ 的矩阵，那么 $f [i]$ 转移到 $f [i + 1]$ 就是乘这个矩阵。

直接矩阵快速幂加速即可。

至于没有质数的情况，一个很简单的方法就是把 $cnt_{i}$ 的定义加上不是质数。

如果你是在不想用容斥，也可以直接把DP设成 $f [i] [j] [0 / 1]$ ，分别表示有质数和没有质数，然后化成二维的 $D P$ ： $f [i] [j / j + p]$ ，但是事实上非常的没有必要，因为矩阵长度乘二，相当于常数乘 $8$ ，大概率过不了。

代码

时间复杂度： $O(p^3\log{n})$

#include<cstdio>
#include<cstring>
#define  N  21000000
#define  SN  3100000
using  namespace  std;
typedef  long  long  LL;
int  li[SN],top;
bool  v[N];
int  n,m,p;
LL  cn1[110],cn2[110];
LL  mod=20170408;
void  XXS()
{
	cn1[1]=cn2[1]=1;
	for(int  i=2;i<=m;i++)
	{
		cn1[i%p]++;
		if(!v[i])li[++top]=i;
		else  cn2[i%p]++;
		for(int  j=1;j<=top  &&  i*li[j]<=m;j++)
		{
			v[li[j]*i]=1;
			if(i%li[j]==0)break;
		}
	}
}
struct  node
{
	LL  a[110][110];
	node(){memset(a,0,sizeof(a));}
};
inline  node  operator*(node  x,node  y)
{
	node  z;
	for(int  i=0;i<p;i++)
	{
		for(int  j=0;j<p;j++)
		{
			for(int  k=0;k<p;k++)z.a[i][j]+=x.a[i][k]*y.a[k][j]%mod,z.a[i][j]%=mod;
		}
	}
	return  z;
}
inline  node  ksm(node  x,int  y)
{
	node  ans;
	for(int  i=0;i<p;i++)ans.a[i][i]=1;
	while(y)
	{
		if(y&1)ans=ans*x;
		x=x*x;y>>=1;
	}
	return  ans;
}
int  main()
{
	scanf("%d%d%d",&n,&m,&p);
	XXS();
	LL  ans=0;node  x;
	for(int  i=0;i<p;i++)
	{
		for(int  j=0;j<p;j++)x.a[i][(j+i)%p]=cn1[j];
	}
	x=ksm(x,n);
	ans=x.a[0][0];
	for(int  i=0;i<p;i++)
	{
		for(int  j=0;j<p;j++)x.a[i][(j+i)%p]=cn2[j];
	}
	x=ksm(x,n);
	ans-=x.a[0][0];
	ans=(ans+mod)%mod;
	printf("%lld\n",ans);
	return  0;
}