2017.8.7 序列计数 思考记录

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

看到倍数就应该想到余数转移，看到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
,i,j,k,p,ans;
bool he
;
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);
}