AC日记——「SDOI2017」序列计数 LibreOJ 2002

「SDOI2017」序列计数

思路：

　　矩阵快速幂；

代码：

#include <bits/stdc++.h>
using namespace std;
#define mod 20170408
#define ll long long
struct MatrixType {
    int n,m;
    ll ai[105][105];
    void mem(int n_,int m_)
    {
        n=n_,m=m_;
        for(int i=0;i<=n;i++) for(int v=0;v<=m;v++) ai[i][v]=0;
    }
    MatrixType operator*(const MatrixType pos)const
    {
        MatrixType res;res.mem(this->n,pos.m);
        for(int i=0;i<=res.n;i++)
            for(int v=0;v<=res.m;v++)
                for(int x=0;x<=this->m;x++)
                    res.ai[i][v]=(res.ai[i][v]+(this->ai[i][x]*pos.ai[x][v])%mod)%mod;
        return res;
    }
    void debug()
    {
        puts("");
        printf("%d %d
",n,m);
        for(int i=0;i<=2;i++)
        {
            for(int v=0;v<=2;v++) printf("%d ",this->ai[i][v]);
            putchar('
');
        }
        puts("");
    }
};
int n,m,p,pi[2000005],cnt,ai[1001],bi[1001];
bool if_p[20000005];
MatrixType poww(MatrixType pos,int mi)
{
    MatrixType res=pos;mi--;
    while(mi)
    {
        if(mi&1) res=res*pos;
        mi=mi>>1,pos=pos*pos;
    }
    return res;
}
void ouler(int lit)
{
    for(int i=2;i<=lit;i++)
    {
        if(!if_p[i]) pi[++cnt]=i;
        for(int j=1;pi[j]*i<=lit&&j<=cnt;j++)
        {
            if_p[i*pi[j]]=true;
            if(i%pi[j]==0) break;
        }
    }
}
int main()
{
    //freopen("ans3.txt","w",stdout);
    scanf("%d%d%d",&n,&m,&p),ouler(m);
    MatrixType sta1;sta1.mem(0,p-1);
    sta1.ai[0][0]=1;
    MatrixType pos1,pos2,ans1,ans2;
    pos1.mem(p-1,p-1),pos2.mem(p-1,p-1);
    for(int i=1;i<=m;i++) ai[i%p]++;
    for(int i=0;i<p;i++) bi[i]=ai[i];
    for(int i=1;i<=cnt;i++) bi[pi[i]%p]--;
    for(int i=0;i<p;i++)
    {
        for(int v=0;v<p;v++)
        {
            if(i>v) pos1.ai[i][v]=ai[i-v],pos2.ai[i][v]=bi[i-v];
            if(i==v) pos1.ai[i][v]=ai[0],pos2.ai[i][v]=bi[0];
            if(i<v) pos1.ai[i][v]=ai[p-v+i],pos2.ai[i][v]=bi[p-v+i];
        }
    }
    pos1=poww(pos1,n),pos2=poww(pos2,n);
    ans1=sta1*pos1,ans2=sta1*pos2;
    cout<<(ans1.ai[0][0]-ans2.ai[0][0]+mod)%mod;
    return 0;
}