bzoj4818 [SDOI2017] 序列计数

, ,

[SDOI2017] 序列计数

https://www.luogu.com.cn/problem/P3702

Solution

一种很简单的做法是全部的减去不包含质数的，用矩阵乘法维护，在这里不赘述，除去筛质数的复杂度是 $O(p^3\log n)$ 的。

观察我们的矩阵，发现是 $w_0,w_1,..,w_{p-1}$ 每次向右移动一位，移动 $p$ 次构成的矩阵。这个矩阵很特殊，再找性质。

回到最初始的 dp，$f_{i,j}\times w_k \to f_{i+1,(j+k)% p}$，可以看成 $f_{i,j}\times f_{1,k} \to f_{i+1,(j+k)% p}$，进而，可以看成 $f_{i,j}\times f_{k,l}\to f_{i+j,(k+l)% p}$。

这个东西显然有结合律，就像矩阵一样。所以也可以快速幂。直接做是 $O(p^2\log n)$ 的。

而且这个是卷积的形式，也可以做到 $O(p\log p\log n)$ 的复杂度。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
#include<bits/stdc++.h>
#define int long long
#define ull unsigned long long
#define maxn 205
#define put() putchar('\n')
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
using namespace std;
void read(int &x){
    int f=1;x=0;char c=getchar();
    while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}
    while (c>='0'&&c<='9') {x=x*10+c-'0';c=getchar();}
    x*=f;
}
namespace Debug{
	Tp void _debug(char* f,Ty t){cerr<<f<<'='<<t<<endl;}
	Ts void _debug(char* f,Ty x,Ar... y){while(*f!=',') cerr<<*f++;cerr<<'='<<x<<",";_debug(f+1,y...);}
	Tp ostream& operator<<(ostream& os,vector<Ty>& V){os<<"[";for(auto& vv:V) os<<vv<<",";os<<"]";return os;}
	#define gdb(...) _debug((char*)#__VA_ARGS__,__VA_ARGS__)
}using namespace Debug;
#define fi first
#define se second
#define mk make_pair
const int mod=20170408;
int power(int x,int y=mod-2) {
	int sum=1;
	while (y) {
		if (y&1) sum=sum*x%mod;
		x=x*x%mod;y>>=1;
	}
	return sum;
}
int n,m,p;
int f[maxn],g[maxn],F[maxn],G[maxn];
bitset<20000005>vis;
int c[maxn],prime[5000005],cnt;
void mul(int *a,int *b) {
	int i,j;
	for (i=0;i<p;i++) c[i]=0;
	for (i=0;i<p;i++)
		for (j=0;j<p;j++)
			c[(i+j)%p]=(c[(i+j)%p]+a[i]*b[j])%mod;
	for (i=0;i<p;i++) a[i]=c[i];
}
signed main(void){
	int i,j;
	read(n);read(m);read(p);
	vis[1]=1;f[1]++,g[1]++;
	for (i=2;i<=m;i++) {
		f[i%p]++;
		if (!vis[i]) prime[++cnt]=i;
		else g[i%p]++;
		for (j=1;j<=cnt&&prime[j]*i<=m;j++) {
			vis[i*prime[j]]=1;
			if (i%prime[j]==0) break;
		}
//		gdb(i,vis[i]);
	}
	F[0]=G[0]=1;
	while (n) {
		if (n&1) mul(F,f),mul(G,g);
		mul(f,f),mul(g,g),n>>=1;
	}
	printf("%lld",(F[0]-G[0]+mod)%mod);
	return 0;
}