[SDOI2017]序列计数

题目描述

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

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

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

输入输出格式

输入格式：

 

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

 

输出格式：

 

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

 

输入输出样例

输入样例#1： 复制 
3 5 3
输出样例#1： 复制 
33

说明

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

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

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

对100\%100%的数据，1\leq n \leq 10^9,1\leq m \leq 2\times 10^7,1\leq p\leq 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的数量

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

于是用矩阵快速幂就行

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

  1 #include<iostream>
  2 #include<cstdio>
  3 #include<cstring>
  4 #include<algorithm>
  5 using namespace std;
  6 typedef long long lol;
  7 struct Matrix
  8 {
  9   lol a[201][201];
 10 }pre,ans1,ans2,Mat1,Mat2;
 11 lol n,m,p,Mod=20170408;
 12 long long cnt1[201],cnt2[201];
 13 lol tot,pri[25000001];
 14 bool vis[25000001];
 15 Matrix operator*(const Matrix a,const Matrix b)
 16 {
 17   lol i,j,k;
 18   Matrix res;
 19   memset(res.a,0,sizeof(res.a));
 20   for (k=0;k<p;k++)
 21   for (i=0;i<p;i++)
 22     if (a.a[i][k])
 23     {
 24       for (j=0;j<p;j++)
 25     {
 26       res.a[i][j]+=a.a[i][k]*b.a[k][j];
 27       res.a[i][j]%=Mod;
 28     }
 29     }
 30   return res;
 31 }
 32 Matrix qpow1(lol y)
 33 {lol i;
 34   Matrix res;
 35   memset(res.a,0,sizeof(res.a));
 36   for (i=0;i<p;i++)
 37     res.a[i][i]=1;
 38   while (y)
 39     {
 40       if (y&1) res=res*Mat1;
 41       Mat1=Mat1*Mat1;
 42       y=y/2;
 43     }
 44   return res;
 45 }
 46 Matrix qpow2(lol y)
 47 {lol i;
 48   Matrix res;
 49   memset(res.a,0,sizeof(res.a));
 50   for (i=0;i<p;i++)
 51     res.a[i][i]=1;
 52   while (y)
 53     {
 54       if (y&1) res=res*Mat2;
 55       Mat2=Mat2*Mat2;
 56       y=y/2;
 57     }
 58   return res;
 59 }
 60 int main()
 61 {lol i,j;
 62   cin>>n>>m>>p;
 63   cnt1[1]++;
 64   for (i=2;i<=m;i++)
 65     {
 66       cnt1[i%p]++,cnt1[i%p]%=Mod;
 67       if (vis[i]==0)
 68     {
 69       ++tot;
 70       pri[tot]=i;
 71     }
 72       for (j=1;j<=tot;j++)
 73     {
 74       if (i*pri[j]>m) break;
 75       vis[i*pri[j]]=1;
 76       if (i%pri[j]==0) break;
 77     }
 78     }
 79   cnt2[1]++;
 80   for(i=2;i<=m;i++)
 81     if (vis[i]) cnt2[i%p]++,cnt2[i%p]%=Mod;
 82   for (i=0;i<p;i++)
 83     {
 84       for (j=0;j<p;j++)
 85     {
 86       Mat1.a[i][(i+j)%p]+=cnt1[j]%Mod;
 87       Mat1.a[i][(i+j)%p]%=Mod;
 88     }
 89     }
 90   for (i=0;i<p;i++)
 91   pre.a[0][i]=cnt1[i];
 92   ans1=qpow1(n-1);
 93     ans1=pre*ans1;
 94     for (i=0;i<p;i++)
 95     {
 96       for (j=0;j<p;j++)
 97     {
 98       Mat2.a[i][(i+j)%p]+=cnt2[j]%Mod;
 99       Mat2.a[i][(i+j)%p]%=Mod;
100     }
101     }
102     for (i=0;i<p;i++)
103     pre.a[0][i]=cnt2[i];
104   ans2=qpow2(n-1);
105     ans2=pre*ans2;
106     cout<<(ans1.a[0][0]-ans2.a[0][0]+Mod)%Mod;
107 }

 

posted @ 2017-10-28 20:27  Z-Y-Y-S  阅读(735)  评论(0)  编辑  收藏  举报