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

[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
", (ans - A.a[0][0] + mod) % mod);
	return 0;
}

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