4818: [Sdoi2017]序列计数

4818: [Sdoi2017]序列计数

Time Limit: 30 Sec Memory Limit: 128 MB
Submit: 396 Solved: 267
[Submit][Status][Discuss]
Description

Alice想要得到一个长度为n的序列，序列中的数都是不超过m的正整数，而且这n个数的和是p的倍数。Alice还希望
，这n个数中，至少有一个数是质数。Alice想知道，有多少个序列满足她的要求。
Input

一行三个数，n,m,p。
1<=n<=10^9,1<=m<=2×10^7,1<=p<=100
Output

一行一个数，满足Alice的要求的序列数量，答案对20170408取模。
Sample Input

3 5 3
Sample Output

33
HINT

Source

鸣谢infinityedge上传

[Submit][Status][Discuss]

注意到题目中的p很小，那说明这么大的m是没用的
用线性筛预处理m以内的素数
枚举每个数字，mod p之后记录到一个数组里面
记f[i][j]为前i个数字的和为j的方案数
转移显然可以矩阵乘法
再统计一个g[i][j]含义相同，但是强制每个数字都是非素数
f[n][0]−g[n][0]就是答案了

#include<iostream>#include<cstdio>#include<cstring>#include<vector>#include<queue>#include<algorithm>#include<cmath>#include<stack>using namespace std;const int N = 100;const int maxm = 2E7 + 233;const int maxn = 1E6 + 3E5 + 233;typedef long long LL;typedef unsigned int u32;const LL mo = 20170408;int n,m,p,tot,c[N],d[N],pri[maxn];u32 mi[32],not_pri[maxm / 32];inline void Mark(int k) {not_pri[k / 32] |= mi[k % 32];}inline bool Search(int k) {return not_pri[k / 32] & mi[k % 32];}inline int Mul(const LL &x,const LL &y) {return x * y % mo;}inline int Add(const int &x,const int &y) {return x + y < mo ? x + y : x + y - mo;}inline int Dec(const int &x,const int &y) {return x - y >= 0 ? x - y : x - y + mo;}struct data{    LL a[N][N]; data(){memset(a,0,sizeof(a));}    data operator * (const data &b)    {        data c;        for (int k = 0; k < p; k++)            for (int i = 0; i < p; i++)                for (int j = 0; j < p; j++)                    c.a[i][j] += a[i][k] * b.a[k][j];        for (int i = 0; i < p; i++)            for (int j = 0; j < p; j++)                c.a[i][j] %= mo;        return c;    }}f,g,F,G;void ksm(){    for (int i = 0; i < p; i++)        F.a[i][i] = G.a[i][i] = 1;    for (; n; n >>= 1)    {        if (n & 1) F = F * f,G = G * g;        f = f * f; g = g * g;    }}int main(){    #ifdef DMC        freopen("DMC.txt","r",stdin);    #endif    mi[0] = 1; cin >> n >> m >> p;    for (int i = 1; i < 32; i++) mi[i] = mi[i - 1] << (u32)(1);    for (int i = 2; i <= m; i++)    {        int now = Search(i); ++c[i % p];        if (!now) pri[++tot] = i; else ++d[i % p];        for (int j = 1; j <= tot; j++)        {            int Nex = i * pri[j]; if (Nex > m) break;            Mark(Nex); if (i % pri[j] == 0) break;        }    }    ++c[1 % p]; ++d[1 % p];    for (int i = 0; i < p; i++)        for (int j = 0; j < p; j++)        {            int k = (i + j) % p;            f.a[i][k] = Add(f.a[i][k],c[j]);            g.a[i][k] = Add(g.a[i][k],d[j]);        }    ksm(); cout << Dec(F.a[0][0],G.a[0][0]) << endl;    return 0;}