BZOJ4818: [Sdoi2017]序列计数

BZOJ4818

容斥一下，变成任取不超过$m$的$n$个数且和为$p$的倍数的方案数-任取不超过$m$的$n$个合数且和为$p$的倍数的方案数。
令$f_{i,j}$表示$i$个数，模$p=j$的方案数。
容易写出方程：
$for(i=1……n)$
$for(j=0……p-1)$
$for(k=1……m)$
$f_{i,(j+k)Mod\space p}+=f_{i-1,j}$
复杂度上显然不允许
发现$m$其实并没有太大用，有用的是每个数$Mod\space m$的值。
那么就用$num_i$记录，所有数中$Mod \space p$为$i$的个数。
然后就可以优化成$O(n*p)$的$dp$
$for(i=1……n)$
$for(j=0……p-1)$
$f_{i,(j+p)Mod \space p}+=num_p*f_{i-1,j}$
复杂度还是太高了。
然后发现这个其实是一个线性齐次递推式。可以矩乘优化嘛。
线性筛素数的标记数组$Not_Prime$要开成$bool$的。。开成$int$CE了两次。。
然后$O(m+p^3*log_n)$$28s$给卡过去了。。可能是我打得太丑了$dsy$垫底$QAQ$
好像还可以再快一点。管他的。。 A <script type="math/tex" id="MathJax-Element-62">A</script>了就行了。。

【代码】

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <bitset>
#define N 20000005
#define Mod 20170408
#define INF 0x7fffffff
using namespace std;
typedef long long ll;

ll read()
{
    ll x=0,f=1;char ch=getchar();
    while(!isdigit(ch)){if(ch=='-') f=-1;ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
    return x*f;
}

int n,m,mod,ans;
int p[1270610],num[105];bool Not_Prime[N]={1,1};

class Matrix{
    public:
        int x,y;
        int v[105][105];
}a,b;

Matrix New(int x)
{
    Matrix rtn;rtn.x=mod,rtn.y=mod;
    for(int i=1;i<=mod;i++)
    for(int j=1;j<=mod;j++) rtn.v[i][j]=(i==j)?x:0;
    return rtn;
}

void Get_Prime()
{
    num[1]=1;
    for(register int i=2;i<=m;i++) {
        if(!Not_Prime[i]) p[++p[0]]=i;
        for(register int j=1;j<=p[0]&&i*p[j]<=m;j++) {
            Not_Prime[i*p[j]]=1;num[(i*p[j])%mod]++;
            if(i%p[j]==0) break;
        }
    }
}

Matrix operator *(Matrix A,Matrix B){
    Matrix rtn=New(0);
    for(int i=1;i<=A.x;i++)
    for(int j=1;j<=B.y;j++)
    for(int k=1;k<=A.y;k++)
        rtn.v[i][j]=(rtn.v[i][j]+1LL*A.v[i][k]*B.v[k][j]%Mod)%Mod;
    return rtn;
}

Matrix Qpow(Matrix X,int y)
{
    Matrix rtn=New(1);
    while(y) {
        if(y&1) rtn=rtn*X;
        X=X*X;y>>=1;
    }
    return rtn;
}

int Get_Ans()
{
    a.x=1;a.y=b.x=b.y=mod;
    for(int i=1;i<=mod;i++)
    for(int j=1;j<=mod;j++) 
        b.v[i][j]=num[(i-j+mod)%mod];
    a.v[1][1]=1;for(int i=2;i<=mod;i++) a.v[1][i]=0;
    a=a*Qpow(b,n);
    return a.v[1][1];
}

int main()
{
    n=read(),m=read(),mod=read();
    for(register int i=1;i<=m;i++) num[i%mod]++;
    ans+=Get_Ans();
    for(int i=0;i<mod;i++) num[i]=0;
    Get_Prime();
    ans-=Get_Ans();
    printf("%d\n",(ans+Mod)%Mod);
    return 0;
}