bzoj 4818: [Sdoi2017]序列计数（DP+矩阵快速幂）

4818: [Sdoi2017]序列计数

Time Limit: 30 Sec  Memory Limit: 128 MB
Submit: 769  Solved: 463
[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

以下假设p=4，n和m未知：

dp[n][k]表示前n个数之和对p取模为k的总情况数，x[k]表示对p取模为k的数的个数（当然<=m）

那么可以列出矩阵：

中间那个矩阵的大小是p²的，所以用矩阵快速幂的话复杂度为O(p^3logn)

可是Alice还希望这n个数中有质数，这样的话还要再算一个不出现质数的dp然后相减

可以啊

#pragma comment(linker, "/STACK:102400000,102400000")#include<stdio.h>#include<string.h>#define mod 20170408#define LL long longLL p, x[115];int flag[20000010];typedef struct{LL i;LL a[105][105];void init(){memset(a, 0, sizeof(a));}void unit(){memset(a, 0, sizeof(a));for(i=1;i<=p;i++)a[i][i] = 1;}}Matrix;Matrix Xpow(Matrix p1, Matrix p2){Matrix pe;LL i, j, k;pe.init();for(i=1;i<=p;i++){for(j=1;j<=p;j++){for(k=1;k<=p;k++)pe.a[i][j] = (pe.a[i][j]+p1.a[i][k]*p2.a[k][j])%mod;}}return pe;}Matrix Powto(Matrix p, LL k){Matrix E;E.unit();while(k){if(k%2)E = Xpow(E, p);p = Xpow(p, p);k /= 2;}return E;}int main(void){Matrix Jz;LL n, m, i, j, ans;flag[0] = flag[1] = 1;for(i=2;i<=20000000;i++){if(flag[i])continue;for(j=i*i;j<=20000000;j+=i)flag[j] = 1;}scanf("%lld%lld%lld", &n, &m, &p);Jz.init();for(i=0;i<=p-1;i++){for(j=i;j<=m;j+=p)x[i]++;}x[0]--;for(i=1;i<=p;i++){for(j=1;j<=p;j++)Jz.a[i][j] = x[((i-j)+p)%p];}Jz = Powto(Jz, n);ans = Jz.a[1][1];Jz.init();memset(x, 0, sizeof(x));for(i=0;i<=p-1;i++){for(j=i;j<=m;j+=p)x[i] += flag[j];}x[0]--;for(i=1;i<=p;i++){for(j=1;j<=p;j++)Jz.a[i][j] = x[((i-j)+p)%p];}Jz = Powto(Jz, n);ans = (ans-Jz.a[1][1]+mod)%mod;printf("%lld\n", ans);return 0;}