【BZOJ 4818】 4818: [Sdoi2017]序列计数 （矩阵乘法、容斥计数）

4818: [Sdoi2017]序列计数

Time Limit: 30 Sec  Memory Limit: 128 MB
Submit: 560  Solved: 359

Description

Alice想要得到一个长度为n的序列，序列中的数都是不超过m的正整数，而且这n个数的和是p的倍数。Alice还希望

，这n个数中，至少有一个数是质数。Alice想知道，有多少个序列满足她的要求。

Input

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

1<=n<=10^9,1<=m<=2×10^7,1<=p<=100

Output

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

Sample Input

3 5 3

Sample Output

33

HINT

Source

鸣谢infinityedge上传

【分析】

　　f[i][j]表示填i个数，mod为j的方案数。

　　这个转移可以用矩阵加速啦，那么n很大也没关系了。

　　然后要至少有一个质数，就用总方案数-没有一个质数。

　　然后一开始跑20s，看到别人2s，然后觉得这个循环矩阵的矩阵乘法可以O(n^2)，改了就2s了。。

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
#define LL long long
#define Maxn 20000010
#define Mod 20170408

int pri[Maxn],pl;
bool vis[Maxn];

void init(int n)
{
    for(int i=1;i<=n;i++) vis[i]=0;
    pl=0;vis[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i]) pri[++pl]=i;
        for(int j=1;j<=pl;j++)
        {
            if(pri[j]*i>n) break;
            vis[i*pri[j]]=1;
            if(i%pri[j]==0) break;
        }
    }
}

int sm[110],p;

struct Matrix
{
    int w[110][110];
    Matrix() {memset(w,0,sizeof(w));}
    friend Matrix operator * (Matrix x,Matrix y)
    {
        Matrix ret;
        for(int k=0;k<p;k++)
        {
            // for(int i=0;i<p;i++)
             for(int j=0;j<p;j++)
             {
                 ret.w[0][j]=(ret.w[0][j]+1LL*x.w[0][k]*y.w[k][j]%Mod)%Mod;
             }
        }
        for(int i=1;i<p;i++)
        {
            for(int j=0;j<p;j++) ret.w[i][(j+i)%p]=ret.w[0][j];
        }
        return ret;
    }
    friend Matrix operator ^ (Matrix x,int b)
    {
        Matrix ret;
        for(int i=0;i<p;i++) ret.w[i][i]=1;
        while(b)
        {
            if(b&1) ret=ret*x;
            x=x*x;
            b>>=1;
        }
        return ret;
    }
};

int main()
{
    int n,m;
    scanf("%d%d%d",&n,&m,&p);
    init(m);
    for(int i=0;i<p;i++) sm[p]=0;
    Matrix nw;
    for(int i=1;i<=m;i++) sm[i%p]++;
    for(int i=0;i<p;i++) for(int j=0;j<p;j++) (nw.w[i][(i+j)%p]+=sm[j])%Mod;
    nw=nw^n;
    int ans=nw.w[0][0];
    
    for(int i=1;i<=m;i++) if(!vis[i]) sm[i%p]--;
    memset(nw.w,0,sizeof(nw.w));
    for(int i=0;i<p;i++) for(int j=0;j<p;j++) (nw.w[i][(i+j)%p]+=sm[j])%Mod;
    nw=nw^n;
    ans=((ans-nw.w[0][0])%Mod+Mod)%Mod;
    
    printf("%d
",ans);
    return 0;
}

2017-04-28 10:51:39