[SDOI 2017] 序列计数

[题目链接]

        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;
}