【bzoj4818】[Sdoi2017]序列计数

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

题解
f[i][j]表示前i个数和为j(mod p)的方案数。
显然可以用矩阵加速。
考虑如何至少选取一个质数，那就是减去一个质数都不选的方案数。

代码

#include<bits/stdc++.h>typedef long long ll;const int mod=20170408;const int N=250005;const double eps=0.00000001;using namespace std;inline int read(){    int x=0,f=1;char ch=getchar();    while (ch<'0'||ch>'9'){if (ch=='-')f=-1;ch=getchar();}    while (ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}    return x*f;}int n,m,p,a[105][105],f[105][105],g[105][105],tmp[105][105];int sum1[105],sum2[105],tot,pri[5000005];bool mark[20000005];void mul(int a[105][105],int b[105][105],int c[105][105]){    for (int i=1;i<=p;i++)        for (int j=1;j<=p;j++)        {            tmp[i][j]=0;            for (int k=1;k<=p;k++)                (tmp[i][j]+=(ll)a[i][k]*b[k][j]%mod)%=mod;        }    for (int i=1;i<=p;i++)        for (int j=1;j<=p;j++)            c[i][j]=tmp[i][j];}int main(){    n=read();m=read();p=read();    sum1[1%p]=sum2[1%p]=1;    for (int i=2;i<=m;i++)    {        sum1[i%p]++;sum2[i%p]++;        if (!mark[i]) pri[++tot]=i,sum2[i%p]--;        for (int j=1;j<=tot&&(ll)pri[j]*i<=m;j++)        {            mark[pri[j]*i]=1;            if (i%pri[j]==0) break;        }    }    for (int i=1;i<=p;i++)        for (int j=1;j<=p;j++)            a[i][j]=sum1[(i-j+p)%p];    for (int i=1;i<=p;i++) f[i][i]=1;    int N=n;    while (N)    {        if (N&1) mul(f,a,f);        mul(a,a,a);        N>>=1;    }    for (int i=1;i<=p;i++)        for (int j=1;j<=p;j++)            a[i][j]=sum2[(i-j+p)%p];    for (int i=1;i<=p;i++) g[i][i]=1;    N=n;    while (N)    {        if (N&1) mul(g,a,g);        mul(a,a,a);        N>>=1;    }    cout<<(f[1][1]-g[1][1]+mod)%mod;    return 0;}