【BZOJ 4818】【SDOI 2017】序列计数

首先用一下容斥，变为所有方案减去没有质数的方案。
设$f[i][j]$表示前i个数和模p为j的方案数，显然答案就是$f[n][0]$。
递推式：$f[i][j]=\sum_{k=1}^{m}f[i-1][t]$，其中$t=((j-k)\%p+p)\%p$，也就是第$i-1$位取的数。
显然这个$n$太大了，考虑到这个方程是一个一阶线性方程，考虑构造矩阵。
下面是所有方案的矩阵：
$\begin{bmatrix}f[i][0]\\ f[i][1]\\ f[i][2]\\ :\\ f[i][p-1]\end{bmatrix}=\begin{bmatrix}1 & 1 & 1 & … & 1\\ 1 & 1 & 1 & … & 1\\ 1 & 1 & 1 & … & 1\\ : &: &: & … & :\\ 1 & 1 & 1 & … & 1\end{bmatrix}\begin{bmatrix}f[i-1][0]\\ f[i-1][1]\\ f[i-1][2]\\ :\\ f[i-1][p-1]\end{bmatrix}$
如果是没有质数的方案，就是改一下那个p*p的矩阵，使得质数的那些地方变成0就好，具体可以看程序。

#include<cmath>
#include<cstdio>
#include<vector>
#include<queue>
#include<cstring>
#include<iomanip>
#include<stdlib.h>
#include<iostream>
#include<algorithm>
#define ll long long
#define inf 1000000000
#define mo 20170408
#define M 20000005
#define P 105
#define fo(i,a,b) for(i=a;i<=b;i++)
#define fd(i,a,b) for(i=a;i>=b;i--)
using namespace std;
struct arr{int a[P][P];} a,res;
int not_prime[M],prime[M/10];
int f[P];
int n,m,p,ans,r1,r2;
void prm(int n)
{
    int i,j;
    not_prime[1] = 1;
    fo(i,2,n)
        {
            if (!not_prime[i]) prime[++prime[0]] = i;
            for (j = 1;j <= prime[0] && i*prime[j] <= n; j++)
                {
                    not_prime[i*prime[j]] = 1;
                    if (i % prime[j] == 0) break;
                }
        }
}
void multi(arr &x,arr y,arr z)
{
    int i,j,k;
    memset(x.a,0,sizeof(x.a));
    fo(i,0,p-1)
        fo(j,0,p-1)
            fo(k,0,p-1)
                x.a[i][j] = (x.a[i][j] + (ll)y.a[i][k] * z.a[k][j] % mo) % mo;
}

void quick_power(arr &x,int t)
{
    int i,j;
    memset(res.a,0,sizeof(res.a));
    fo(i,0,p-1) res.a[i][i] = 1;
    while (t)
        {
            if (t&1) multi(res,res,x);
            multi(x,x,x); t>>=1;
        }
    fo(i,0,p-1) fo(j,0,p-1) x.a[i][j] = res.a[i][j]; 
}
int solve1()
{
    int i,j; 
    memset(f,0,sizeof(f)); memset(a.a,0,sizeof(a.a));
    fo(i,1,m) f[i%p]++;
    fo(i,1,m) a.a[0][(-i%p+p)%p]++;
    fo(i,1,p-1)
        fo(j,0,p-1)
            a.a[i][j] = a.a[i-1][(j-1+p)%p];
    quick_power(a,n-1);
    int ans = 0;
    fo(i,0,p-1) ans = (ans + (ll)a.a[0][i]*f[i]%mo)%mo;
    return ans;
}
int solve2()
{
    int i,j; 
    memset(f,0,sizeof(f)); memset(a.a,0,sizeof(a.a));
    fo(i,1,m) if (not_prime[i]) f[i%p]++;
    fo(i,1,m) if (not_prime[i]) a.a[0][(-i%p+p)%p]++;
    fo(i,1,p-1)
        fo(j,0,p-1)
            a.a[i][j] = a.a[i-1][(j-1+p)%p];
    quick_power(a,n-1);
    int ans = 0;
    fo(i,0,p-1) ans = (ans + (ll)a.a[0][i]*f[i]%mo)%mo;
    return ans;
}
int main()
{
    scanf("%d%d%d",&n,&m,&p);
    prm(m);
    r1 = solve1();
    r2 = solve2();
    //cout<<r1<<" "<<r2<<endl;
    printf("%d\n",(r1-r2+mo)%mo);
}