bzoj4818 [Sdoi2017]序列计数

Description

---

Alice想要得到一个长度为n的序列，序列中的数都是不超过m的正整数，而且这n个数的和是p的倍数。Alice还希望
，这n个数中，至少有一个数是质数。Alice想知道，有多少个序列满足她的要求。

一行三个数，n,m,p。
1<=n<=10^9,1<=m<=2×10^7,1<=p<=100

来自 https://www.lydsy.com/JudgeOnline/problem.php?id=4818

Solution

---

如果模数有原根的话是可以NTT+快速幂的，

先写出dp柿子，f[I,j]=Σf[i-1,j-k]。对于这种当前层只和上一层的状态有关、且n特别大的dp都可以套矩阵乘法。转移矩阵的求法可以考虑把f压成一维向量，转移的实质就是依靠矩阵中的系数，注意乘的方式即可

朴素地求转移矩阵是O(mp)的，氮素考虑到相邻两列只差了1，也就是他们是循环同构的，可以先求一列然后p^2推出整个转移矩阵

人丑自带大常数

Code

---

#include <stdio.h>
#include <string.h>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)
#define fill(x,t) memset(x,t,sizeof(x))

typedef long long LL;
const int N=20000005;
const int P=105;
const int MOD=20170408;

int prime[N/15],a[P],n,m,p;

bool not_prime[N+1];

struct mat {
    int a[P][P];
    mat operator *(mat b) const {
        mat c; fill(c.a,0);
        rep(i,0,p-1) rep(j,0,p-1) {
            rep(k,0,p-1) c.a[i][j]=(c.a[i][j]+(LL)a[i][k]*b.a[k][j]%MOD)%MOD;
        }
        return c;
    }
} tmp,B;

mat ksm(mat x,LL dep) {
    fill(tmp.a,0); rep(i,0,p-1) tmp.a[i][i]=1;
    while (dep) {
        if (dep&1) tmp=tmp*x;
        x=x*x; dep/=2;
    }
    return tmp;
}

void pre_work() {
    rep(i,2,N) {
        if (!not_prime[i]) {
            prime[++prime[0]]=i;
        }
        for (int j=1;i*prime[j]<=N&&j<=prime[0];j++) {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0) break;
        }
    }
    not_prime[1]=1;
}

LL solve1() {
    fill(a,0); fill(B.a,0);
    rep(i,1,m) a[i%p]++;
    rep(i,1,m) B.a[(p-i%p)%p][0]++;
    rep(i,1,p-1) {
        rep(j,0,p-1) {
            B.a[j][i]=B.a[(j+p-1)%p][i-1];
        }
    }
    B=ksm(B,n-1);
    LL ans=0;
    rep(i,0,p-1) ans=(ans+(LL)a[i]*B.a[i][0]%MOD)%MOD;
    return ans;
}

LL solve2() {
    fill(a,0); fill(B.a,0);
    rep(i,1,m) if (not_prime[i]) a[i%p]++;
    rep(i,1,m) if (not_prime[i]) B.a[(p-i%p)%p][0]++;
    rep(i,1,p-1) {
        rep(j,0,p-1) {
            B.a[j][i]=B.a[(j+p-1)%p][i-1];
        }
    }
    B=ksm(B,n-1);
    LL ans=0;
    rep(i,0,p-1) ans=(ans+(LL)a[i]*B.a[i][0]%MOD)%MOD;
    return ans;
}

int main(void) {
    pre_work();
    scanf("%d%d%d",&n,&m,&p);
    LL ans1=solve1();
    LL ans2=solve2();
    printf("%lld\n", (ans1-ans2+MOD)%MOD);
    return 0;
}