BZOJ 4818: [Sdoi2017]序列计数 （动态规划+矩阵乘法）

来源：http://blog.csdn.net/qq_33229466/article/details/70055284
题目：http://www.lydsy.com/JudgeOnline/problem.php?id=4818
分析

一眼容斥，用所有方案减去不含质数的方案。
设f[i,j]表示序列前i个数模p的余数为j时的方案数。f[i,j]=∑f[i−1,(j−k)modp]
构建矩阵然后快速幂即可。
求不含质数的方案同理。
在一开始建矩阵的时候我的复杂度是mp，发现会超时。但是注意到矩阵的每一列都是循环同构的，于是就可以先O(m)处理好一列，然后再构建其他列即可。
代码：

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;

typedef long long LL;

const int M=20000005;
const int P=105;
const int MOD=20170408;

int n,m,p,prime[M/10],tot,f[P];
struct arr{int a[P][P];}ans,a;
bool not_prime[M];

void get_prime(int n)
{
    not_prime[1]=1;
    for (int i=2;i<=n;i++)
    {
        if (!not_prime[i]) prime[++tot]=i;
        for (int j=1;j<=tot&&i*prime[j]<=n;j++)
        {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0) break;
        }
    }
}

void mul(arr &c,arr a,arr b)
{
    memset(c.a,0,sizeof(c.a));
    for (int i=0;i<p;i++)
        for (int j=0;j<p;j++)
            for (int k=0;k<p;k++)
                c.a[i][j]=(c.a[i][j]+(LL)a.a[i][k]*b.a[k][j]%MOD)%MOD;
}

arr ksm(arr x,int y)
{
    memset(ans.a,0,sizeof(ans.a));
    for (int i=0;i<p;i++) ans.a[i][i]=1;
    while (y)
    {
        if (y&1) mul(ans,ans,x);
        mul(x,x,x);y>>=1;
    }
    return ans;
}

int solve1()
{
    for (int i=1;i<=m;i++) f[i%p]++;
    for (int j=1;j<=m;j++) a.a[(-j%p+p)%p][0]++;
    for (int i=1;i<p;i++)
        for (int j=0;j<p;j++)
            a.a[j][i]=a.a[(j-1+p)%p][i-1];
    a=ksm(a,n-1);
    int ans=0;
    for (int i=0;i<p;i++) ans=(ans+(LL)f[i]*a.a[i][0]%MOD)%MOD;
    return ans;
}

int solve2()
{
    memset(f,0,sizeof(f));
    for (int i=1;i<=m;i++) if (not_prime[i]) f[i%p]++;
    memset(a.a,0,sizeof(a.a));
    for (int j=1;j<=m;j++)
        if (not_prime[j]) a.a[(-j%p+p)%p][0]++;
    for (int i=1;i<p;i++)
        for (int j=0;j<p;j++)
            a.a[j][i]=a.a[(j-1+p)%p][i-1];
    a=ksm(a,n-1);
    int ans=0;
    for (int i=0;i<p;i++) ans=(ans+(LL)f[i]*a.a[i][0]%MOD)%MOD;
    return ans;
}

int main()
{
    scanf("%d%d%d",&n,&m,&p);
    get_prime(m);
    int x=solve1(),y=solve2();
    printf("%d",(x-y+MOD)%MOD);
    return 0;
}