[题目链接]
https://www.lydsy.com/JudgeOnline/problem.php?id=4818
[算法]
考虑容斥 , 用有至少有一个质数的合法序列数 - 没有质数的合法序列数
这两个问题是等价的 , 为方便讨论 , 我们考虑前者该如何计算 :
用fi , j表示前i个数 , 模p余j的合法序列数
显然有fi , j = sigma{ fi - 1 , j - k }
矩阵优化即可
时间复杂度 : O(M + logN)
[代码]
#include<bits/stdc++.h> using namespace std; #ifndef LOCAL #define eprintf(...) fprintf(stderr , _VA_ARGS) #else #define eprintf(...) 42 #endif typedef long long ll; typedef long double ld; typedef vector< int > vi; typedef pair<int , int> pii; typedef pair<ll , int> pli; typedef pair<ll , ll> pll; typedef unsigned long long ull; #define mp make_pair #define fi first #define se second const int N = 2e7 + 10; const int P = 20170408; struct Tmatrix { int mat[110][110]; } a; int n , m , p , tot; int prime[N] , sa[110] , sb[110]; bool lab[N] , f[N]; template <typename T> inline void chkmax(T &x , T y) { x = max(x , y); } template <typename T> inline void chkmin(T &x , T y) { x = min(x , y); } template <typename T> inline void read(T &x) { T f = 1; x = 0; char c = getchar(); for (; !isdigit(c); c = getchar()) if (c == '-') f = -f; for (; isdigit(c); c = getchar()) x = (x << 3) + (x << 1) + c - '0'; x *= f; } inline void sieve( ) { for (int i = 2; i <= m; ++i) { if (!f[i]) { prime[++tot] = i; lab[i] = true; } for (int j = 1; j <= tot; ++j) { int tmp = i * prime[j]; if (tmp > m) break; f[tmp] = true; if (i % prime[j] == 0) break; } } } inline int sub(int x , int y) { x -= y; while (x < 0) x += P; return x; } inline void multipy(Tmatrix &a , Tmatrix b) { Tmatrix c; for (int i = 0; i < p; ++i) { for (int j = 0; j < p; ++j) { c.mat[i][j] = 0; } } for (int i = 0; i < p; ++i) { for (int j = 0; j < p; ++j) { for (int k = 0; k < p; ++k) { c.mat[i][j] = (c.mat[i][j] + 1ll * a.mat[i][k] * b.mat[k][j] % P) % P; } } } for (int i = 0; i < p; ++i) { for (int j = 0; j < p; ++j) { a.mat[i][j] = c.mat[i][j]; } } } inline void qpow(Tmatrix &a , int n) { Tmatrix b , res; for (int i = 0; i < p; ++i) { for (int j = 0; j < p; ++j) { res.mat[i][j] = (i == j); b.mat[i][j] = a.mat[i][j]; } } while (n > 0) { if (n & 1) multipy(res , b); multipy(b , b); n >>= 1; } for (int i = 0; i < p; ++i) { for (int j = 0; j < p; ++j) { a.mat[i][j] = res.mat[i][j]; } } } inline int calc(int n , int type) { for (int i = 0; i < p; ++i) { for (int j = 0; j < p; ++j) { int k = (i - j + p) % p; a.mat[i][j] = type == 0 ? sa[k] : sb[k]; } } qpow(a , n); return a.mat[0][0]; } int main() { read(n); read(m); read(p); sieve( ); for (int i = 1; i <= m; ++i) { sa[i % p] = (sa[i % p] + 1) % P; if (!lab[i]) sb[i % p] = (sb[i % p] + 1) % P; } printf("%d " , sub(calc(n , 0) , calc(n , 1))); return 0; }