[SDOI2017]序列计数

题目描述

Alice想要得到一个长度为nn的序列，序列中的数都是不超过mm的正整数，而且这nn个数的和是pp的倍数。

Alice还希望，这nn个数中，至少有一个数是质数。

Alice想知道，有多少个序列满足她的要求。

输入输出格式

输入格式：

一行三个数，n,m,pn,m,p。

输出格式：

一行一个数，满足Alice的要求的序列数量，答案对2017040820170408取模。

输入输出样例

输入样例#1： 复制

3 5 3

输出样例#1： 复制

33

说明

对20\%20%的数据，1leq n,mleq1001≤n,m≤100

对50\%50%的数据，1leq m leq 1001≤m≤100

对80\%80%的数据，1leq mleq 10^61≤m≤106

对100\%100%的数据，1leq n leq 10^9,1leq m leq 2 imes 10^7,1leq pleq 1001≤n≤109,1≤m≤2×107,1≤p≤100

时间限制：3s

空间限制：128MB

至少有一个素数的方案=所有方案-没有素数的方案

于是用容斥就变成了简单的dp，先讨论所有方案

令f[i][j]表示i个数，和%p为j的方案数

f[i][j]=∑f[i-1][(j-k+p)%p]*cnt[k]

cnt[k]是1~m中%p等于k的数量

发现显然上式可以写为矩阵

于是用矩阵快速幂就行

然后用欧拉筛把素数筛掉，再做一次

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long lol;
struct Matrix
{
  lol a[201][201];
}pre,ans1,ans2,Mat1,Mat2;
lol n,m,p,Mod=20170408;
long long cnt1[201],cnt2[201];
lol tot,pri[25000001];
bool vis[25000001];
Matrix operator*(const Matrix a,const Matrix b)
{
  lol i,j,k;
  Matrix res;
  memset(res.a,0,sizeof(res.a));
  for (k=0;k<p;k++)
  for (i=0;i<p;i++)
    if (a.a[i][k])
    {
      for (j=0;j<p;j++)
    {
      res.a[i][j]+=a.a[i][k]*b.a[k][j];
      res.a[i][j]%=Mod;
    }
    }
  return res;
}
Matrix qpow1(lol y)
{lol i;
  Matrix res;
  memset(res.a,0,sizeof(res.a));
  for (i=0;i<p;i++)
    res.a[i][i]=1;
  while (y)
    {
      if (y&1) res=res*Mat1;
      Mat1=Mat1*Mat1;
      y=y/2;
    }
  return res;
}
Matrix qpow2(lol y)
{lol i;
  Matrix res;
  memset(res.a,0,sizeof(res.a));
  for (i=0;i<p;i++)
    res.a[i][i]=1;
  while (y)
    {
      if (y&1) res=res*Mat2;
      Mat2=Mat2*Mat2;
      y=y/2;
    }
  return res;
}
int main()
{lol i,j;
  cin>>n>>m>>p;
  cnt1[1]++;
  for (i=2;i<=m;i++)
    {
      cnt1[i%p]++,cnt1[i%p]%=Mod;
      if (vis[i]==0)
    {
      ++tot;
      pri[tot]=i;
    }
      for (j=1;j<=tot;j++)
    {
      if (i*pri[j]>m) break;
      vis[i*pri[j]]=1;
      if (i%pri[j]==0) break;
    }
    }
  cnt2[1]++;
  for(i=2;i<=m;i++)
    if (vis[i]) cnt2[i%p]++,cnt2[i%p]%=Mod;
  for (i=0;i<p;i++)
    {
      for (j=0;j<p;j++)
    {
      Mat1.a[i][(i+j)%p]+=cnt1[j]%Mod;
      Mat1.a[i][(i+j)%p]%=Mod;
    }
    }
  for (i=0;i<p;i++)
  pre.a[0][i]=cnt1[i];
  ans1=qpow1(n-1);
    ans1=pre*ans1;
    for (i=0;i<p;i++)
    {
      for (j=0;j<p;j++)
    {
      Mat2.a[i][(i+j)%p]+=cnt2[j]%Mod;
      Mat2.a[i][(i+j)%p]%=Mod;
    }
    }
    for (i=0;i<p;i++)
    pre.a[0][i]=cnt2[i];
  ans2=qpow2(n-1);
    ans2=pre*ans2;
    cout<<(ans1.a[0][0]-ans2.a[0][0]+Mod)%Mod;
}