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
[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 long
LL 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;
}