【动态规划20】bzoj4818[sdoi2017]序列计数（dp+矩阵快速幂）

题目描述

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

输入输出格式

一行三个数，n,m,p。
1<=n<=10^9,1<=m<=2×10^7,1<=p<=100
一行一个数，满足Alice的要求的序列数量，答案对20170408取模。

题目显然容斥原理，用所有方案-不含质数的方案为所求。
首先，有一个非常菜的dp方程。
f[i][j]表示前i个数%p==j的序列个数
$f[i][j]+=f[i-1][(j-k)$%$p]$$(1<=k<=m)$
那么就是一个显然的矩阵乘法。
首先f[i][]只与f[i-1][]相关，那么我们可以将原来的i维消掉。
那么之后的转移就是将所有$f[(j-k)$%$p]$$(1<=k<=m)$转移到f[j]
那矩阵就很好构建了，从1到m枚举k，矩阵的第j行的(j-k)%p列的值就加一，很好理解。（对于质数在枚举过程中特殊判断即可）
但是若是每一行都从1到m枚举k，时间复杂度上是不允许的。
但是很显然，我们可以意识到矩阵的第j+1行的(j+1-k)%p列，实际上就是j行的前一列（直接写出来有点像废话..），所以直接暴力从矩阵的上一行转移，这样子最后就是一个p*p的矩阵。
直接快速幂搞就完事了。

#include<bits/stdc++.h>
#define fer(i,j,n) for(int i=j;i<=n;i++)
#define far(i,j,n) for(int i=j;i>=n;i--)
#define ll long long
const int maxn=20000010;
const int INF=1e9+7;
const int mod=20170408;
using namespace std;
/*----------------------------------------------------------------------------*/
inline ll read()
{
    char ls;ll x=0,sng=1;
    for(;ls<'0'||ls>'9';ls=getchar())if(ls=='-')sng=-1;
    for(;ls>='0'&&ls<='9';ls=getchar())x=x*10+ls-'0';
    return x*sng;
}
/*----------------------------------------------------------------------------*/
int n,m,p,cnt;
int prime[maxn],f[110];
bool flag[maxn];
struct kaga
{
    ll v[110][110];
    kaga friend operator *(kaga a,kaga b)
    {
        kaga c;
        fer(i,0,p-1)
            fer(j,0,p-1)
            {
                c.v[i][j]=0;
                fer(k,0,p-1)
                c.v[i][j]+=a.v[i][k]*b.v[k][j]%mod;
            }
        return c;
    }
    kaga friend operator ^(kaga a,ll k)
    {
        kaga c;
        fer(i,0,p-1)
            fer(j,0,p-1)
                if(i==j)c.v[i][j]=1;
                else c.v[i][j]=0;
        for(;k;k>>=1,a=a*a)
            if(k&1)c=c*a;
        return c;
    }
    void friend print(kaga a)
    {
        fer(i,0,p-1)
        {
            fer(j,0,p-1)
                cout<<a.v[i][j]<<" ";
            cout<<endl; 
        }   
    }
}a;
void Prime(int n)
{
    memset(flag,0,sizeof(flag));
    flag[1]=1;
    cnt=0;
    fer(i,2,n)
    {
        if(!flag[i])prime[++cnt]=i;
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++)
        {
            flag[i*prime[j]]=1;
            if(!(i%prime[j]))break; 
        }
    }
}
int solve1()
{
    fer(i,1,m)f[i%p]++;
    fer(i,1,m)a.v[(-i%p+p)%p][0]++;
    fer(i,1,p-1)
        fer(j,0,p-1)
            a.v[j][i]=a.v[(j-1+p)%p][i-1];
    a=a^(n-1);
    int ans=0;
    fer(i,0,p-1)ans=(ans+(ll)f[i]*a.v[i][0]%mod)%mod;
    return ans;
}
int solve2()
{
    memset(f,0,sizeof(f));
    fer(i,1,m)if(flag[i])f[i%p]++;
    memset(a.v,0,sizeof(a.v));
    fer(i,1,m)if(flag[i])a.v[(-i%p+p)%p][0]++;
    fer(i,1,p-1)
        fer(j,0,p-1)
            a.v[j][i]=a.v[(j-1+p)%p][i-1];
    a=a^(n-1);
    int ans=0;
    fer(i,0,p-1)ans=(ans+(ll)f[i]*a.v[i][0]%mod)%mod;
    return ans;
}
int main()
{
    n=read();m=read();p=read();
    Prime(m);
    cout<<(solve1()-solve2()+mod)%mod;
}