loj#2002. 「SDOI2017」序列计数(dp 矩阵乘法)
程序员文章站
2023-09-07 16:31:35
题意 "题目链接" Sol 质数的限制并没有什么卵用,直接容斥一下:答案 = 忽略质数总的方案 没有质数的方案 那么直接dp,设$f[i][j]$表示到第i个位置,当前和为j的方案数 $f[i + 1][(j + k) \% p] += f[i][j]$ 矩乘优化一下。 cpp include de ......
题意
sol
质数的限制并没有什么卵用,直接容斥一下:答案 = 忽略质数总的方案 - 没有质数的方案
那么直接dp,设\(f[i][j]\)表示到第i个位置,当前和为j的方案数
\(f[i + 1][(j + k) \% p] += f[i][j]\)
矩乘优化一下。
#include<bits/stdc++.h> #define ll long long using namespace std; const int maxn = 2e7 + 10, mod = 20170408, ss = 1e5 + 10; ll gg = 1e17; inline int read() { char c = getchar(); int x = 0, f = 1; while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();} while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar(); return x * f; } template<typename a, typename b> inline int add(a x, b y) { if(x + y < 0) return x + y + mod; else return x + y >= mod ? x + y - mod : x + y; } template<typename a, typename b> inline void add2(a &x, b y) { if(x + y < 0) x = x + y + mod; else x = (x + y >= mod ? x + y - mod : x + y); } template<typename a, typename b> inline int mul(a x, b y) { return 1ll * x * y % mod; } int n, m, p, lim;//1 - m, ºíêçpµä±¶êý int f[ss], vis[maxn], mu[maxn], prime[maxn], tot, cnt, num[ss], tim[ss], val[ss]; struct ma { int m[201][201]; ma() { memset(m, 0, sizeof(m)); } void init() { for(int i = 0; i <= lim; i++) m[i][i] = 1; } void clear() { memset(m, 0, sizeof(m)); } void print() { for(int i = 0; i <= lim; i++, puts("")) for(int j = 0; j <= lim; j++) printf("%d ", m[i][j]); } ma operator * (const ma &rhs) const { ma ans = {}; for(int i = 0; i <= lim; i++) for(int j = 0; j <= lim; j++) { __int128 tmp = 0; for(int k = 0; k <= lim; k++) { tmp += 1ll * m[i][k] * rhs.m[k][j]; } ans.m[i][j] = tmp % mod; } return ans; } }g; ma matrixpow(ma a, int p) { ma base; base.init(); while(p) { if(p & 1) base = base * a; a = a * a; p >>= 1; } return base; } void sieve(int n) { vis[1] = 1; mu[1] = 1; for(int i = 2; i <= n; i++) { if(!vis[i]) prime[++tot] = i, mu[i] = -1; for(int j = 1; j <= tot && i * prime[j] <= n; j++) { vis[i * prime[j]] = 1; if(i % prime[j]) mu[i * prime[j]] = -mu[i]; else {mu[i * prime[j]] = 0; break;} } } for(int i = 1; i <= n; i++) if(vis[i]) num[i % p]++; } int solve1() {//ºöêóöêêýµäïþöæ for(int i = 1; i <= m; i++) f[i % p]++; for(int j = 0; j < p; j++) { memset(tim, 0, sizeof(tim)); memset(val, 0, sizeof(val)); int step = m; for(int k = 1; k <= m; k++) { int nxt = (j + k) % p; if(tim[nxt]) {step = k - 1; break;} tim[nxt] = 1; val[nxt]++; } if(step) for(int k = 0; k <= lim; k++) g.m[k][j] = m / step * val[k]; for(int k = m / step * step + 1; k <= m; k++) g.m[(j + k) % p][j]++; } ma ans = matrixpow(g, n - 1); int out = 0; for(int i = 0; i <= lim; i++) add2(out, mul(ans.m[0][i], f[i])); return out; } int solve2() {//îþöêêý memset(f, 0, sizeof(f)); g.clear(); for(int i = 1; i <= m; i++) f[i % p] += (vis[i]); for(int j = 0; j < p; j++) for(int k = 0; k < p; k++) g.m[(j + k) % p][j] += num[k]; ma ans = matrixpow(g, n - 1); int out = 0; for(int i = 0; i <= lim; i++) add2(out, mul(ans.m[0][i], f[i])); return out; } int main() { n = read(); m = read(); lim = p = read(); sieve(m); cout << (solve1() - solve2() + mod) % mod; return 0; }