[bzoj4818][Sdoi2017]序列计数_矩阵乘法_欧拉筛

Posted 2022-09-02 shurak

[Sdoi2017]序列计数

题目大意：https://www.lydsy.com/JudgeOnline/problem.php?id=4818.

---

题解：

首先列出来一个递推式子

$f[i][0]$表示$i$个任意数的答案。

$f[i][1]$表示$i$个合数的答案。

转移的时候发现可以用矩阵优化这个过程。

至于怎么把矩阵建出来，我们可以开个桶来解决这个问题。

代码：

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

const int mod = 20170408 ;

char *p1, *p2, buf[100000];

#define nc() (p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF : *p1 ++ )

int rd() 
	int x = 0;
	char c = nc();
	while (c < 48) 
		c = nc();
	
	while (c > 47) 
		x = (((x << 2) + x) << 1) + (c ^ 48), c = nc();
	
	return x;


int p, n, m;

struct Matr 
	int a[110][110];
	Matr() memset(a, 0, sizeof a);
	friend Matr operator * (const Matr &a, const Matr &b) 
		Matr re;
		for (int i = 0; i < p; i ++ ) 
			for (int j = 0; j < p; j ++ ) 
				for (int k = 0; k < p; k ++ ) 
					(re.a[i][j] += (ll)a.a[i][k] * b.a[k][j] % mod) %= mod;
				
			
		
		return re;
	
	friend Matr operator ^ (Matr x, int y) 
		Matr re;
		for (int i = 0; i < p; i ++ ) 
			re.a[i][i] = 1;
		
		while (y) 
			if (y & 1) 
				re = re * x;
			
			y >>= 1;
			x = x * x;
		
		return re;
	
M, A;

int bu[110];

bool vis[20000010];

int prime[20000010], cnt;

int main() 
	n = rd(), m = rd(), p = rd();	
	for (int i = 0; i < p; i ++ ) 
		bu[i] = m / p;
		if (i) 
			if (i <= m % p) 
				bu[i] ++ ;
			
		
	
	for (int i = 0; i < p; i ++ ) 
		for (int j = 0; j < p; j ++ ) 
			M.a[i][j] = bu[(j - i + p) % p];
		
	
	for (int i = 0; i < p; i ++ ) 
		A.a[0][i] = bu[i];
	
	A = A * (M ^ (n - 1));
	int ans = A.a[0][0];
	vis[1] = true;
	for (int i = 2; i <= m; i ++ ) 
		if (!vis[i]) 
			prime[ ++ cnt] = i;
		
		for (int j = 1; j <= cnt && (ll)i * prime[j] <= m; j ++ ) 
			vis[i * prime[j]] = true;
			if (i % prime[j] == 0) 
				break;
			
		
	
	for (int i = 1; i <= m; i ++ ) 
		if (!vis[i]) 
			bu[i % p] -- ;
		
	
	for (int i = 0; i < p; i ++ ) 
		A.a[0][i] = bu[i];
	
	for (int i = 0; i < p; i ++ ) 
		for (int j = 0; j < p; j ++ ) 
			M.a[i][j] = bu[(j - i + p) % p];
		
	
	A = A * (M ^ (n - 1));
	printf("%d\n", (ans - A.a[0][0] + mod) % mod);
	return 0;

小结：就是这种求存在的问题，可以转化成全部-不存在。