LibreOJ #2002. 「SDOI2017」序列计数

/*
    LibreOJ #2002. 「SDOI2017」序列计数
    线性筛 + 矩阵优化dp

    先构造出全部情况的矩阵
    用矩阵快速幂计算答案
    再构造出全不是质数的矩阵
    计算出答案
    前一个答案减后一个答案即可
*/
#include <cstdio>
#include <iostream>
#include <cstring>

const int BUF = 12312312;
char Buf[BUF], *buf = Buf;

#define Mod 20170408

inline void read (int &now)
{
    for (now = 0; !isdigit (*buf); ++ buf);
    for (; isdigit (*buf); now = now * 10 + *buf - ‘0‘, ++ buf);
}

#define Max 20000007
#define L 100
struct Matrix
{
    int c[L][L];
    Matrix () { memset (c, 0, sizeof c); }

    Matrix operator * (const Matrix &now) const
    {
        Matrix res;
        for (register int i = 0, j, k; i < L; ++ i)
            for (k = 0; k < L; ++ k)
                if (this->c[i][k])
                    for (j = 0; j < L; ++ j)
                        res.c[i][j] = (res.c[i][j] + 1LL * this->c[i][k] * now.c[k][j]) % Mod;
        return res;    
    }
};

int operator ^ (Matrix A, int p)
{
    Matrix res; res.c[0][0] = 1;
    for (; p; A = A * A, p >>= 1)
        if (p & 1)
            res = res * A;
    return res.c[0][0];
}

Matrix A;
int prime[Max];
bool is[Max];
int C;

void Euler (const int N)
{
    is[1] = true;
    for (register int i = 2, j; i <= N; ++ i)
    {
        if (!is[i]) prime[++ C] = i;
        for (j = 1; i * prime[j] <= N && j <= C; ++ j)
        {
            is[i * prime[j]] = true;
            if (i % prime[j] == 0) break;
        }
    }
}

int data[Max];

int Main ()
{
    fread (buf, 1, BUF, stdin);
    int N, M, K; read (N), read (M), read (K);
    Euler (Max - 7); register int i, j;

    for (i = 1; i <= M; ++ i) ++ data[i % K];
    for (i = 0; i < K; ++ i) 
        for (j = 0; j < K; ++ j)
            A.c[i][(i + j) % K] = data[j] % Mod;
    int Answer = A ^ N; memset (A.c, 0, sizeof A.c);
    memset (data, 0, sizeof data);
    for (i = 1; i <= M; ++ i)
        if (is[i]) ++ data[i % K];
    for (i = 0; i < K; ++ i)
        for (j = 0; j < K; ++ j)
            A.c[i][(i + j) % K] = data[j] % Mod;
    Answer -= A ^ N;

    printf ("%d", (Answer + Mod) % Mod);
    return 0;
}
int ZlycerQan = Main ();
int main (int argc, char *argv[]) {;}