2017.8.7 序列计数 思考记录

来源:互联网 发布:dvr监控软件下载 编辑:程序博客网 时间:2024/05/22 15:52

这个题真心是矩乘裸题,和上一个题基本一毛一样,稍微做一些矩乘题就可以轻松搞出、(然而我并不轻松)

看到倍数就应该想到余数转移,看到10^2就应该想到n^3的矩乘、所以直接全部-全合数两遍跑完即可


30s时限,卡得很死(bzoj慢的要死,不开o2



dev的调试坏了、重载乘号直接卡主不动崩溃了好几次、

然后这题卡10^7的int   所以欧拉筛和桶就共用一个了、

码:

#include<iostream>#include<cstdio>#include<cstring>using namespace std;#define M 20170408 #define N 20000009 int n,m,tot,tong[N],i,j,k,p,ans;bool he[N];void eular(){for(int i=2;i<=m;i++){if(!he[i]){tong[++tot]=i;}for(int j=1;i*tong[j]<=m&&j<=tot;j++){int k=tong[j]*i;he[k]=1;if(i%tong[j]==0)break;}}}struct jz{int p[105][105],c,k;}a,b,u,v;jz operator *(jz x,jz y)  {      int i,j,k;    for(i=0;i<=x.k;i++)    for(j=0;j<=y.c;j++)    u.p[i][j]=0;      u.k=x.k;      u.c=x.c;    for(i=0;i<=x.k;i++)    for(j=0;j<=x.c;j++)     for(k=0;k<=y.c;k++)      u.p[i][k]=(u.p[i][k]+1ll*x.p[i][j]*y.p[j][k]%M)%M;      return u;  }  jz operator ^(jz x,int ci)  {      int i,j,k;      v.k=x.k;      v.c=x.c;    for(i=0;i<=x.k;i++)    for(j=0;j<=x.c;j++)    if(i==j)v.p[i][i]=1;else v.p[i][j]=0;      while(ci)      {      if(ci&1)v=v*x;            ci>>=1;          x=x*x;            }      return v;  }  int main(){scanf("%d%d%d",&n,&m,&p);a.k=p-1;b.c=p-1;b.k=p-1;for(i=1;i<=m;i++){tong[i%p]++;}   for(i=0;i<p;i++)   for(j=0;j<p;j++)   {   b.p[(j+i)%p][i]+=tong[j];   }    n--;   for(i=1;i<=m;i++)   a.p[i%p][0]++;a=(b^n)*a;ans+=a.p[0][0];eular();memset(tong,0,sizeof(tong));for(i=0;i<p;i++)for(j=0;j<p;j++)b.p[i][j]=a.p[i][j]=0;for(i=1;i<=m;i++){if(he[i]||i==1)tong[i%p]++;}   for(i=0;i<p;i++)   for(j=0;j<p;j++)   {   b.p[(j+i)%p][i]+=tong[j];   }   for(i=1;i<=m;i++)   if(he[i]||i==1)a.p[i%p][0]++;a=(b^n)*a;ans-=a.p[0][0];ans=(ans%M+M)%M;printf("%d",ans);}