【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;
}