bzoj4818 SDOI2017 序列计数

题目链接

思路

先考虑暴力(dp)，(f[i][j])表示前(i)个数，数字之和模(P)余(j)的方案数。
我们先不考虑必须有质数这个条件，先统计出全部方案。然后再减去没有质数的方案就行了。
那么就有(f[i + 1][(j + k) \% p] += f[i][j](1le k le m))
然后发现这个其实并不需要(O(m))的转移，因为((j + k) \% p)对于每个余数来说，最少有(lfloorfrac{m}{p} floor)个。然后有(m\%p)个数有(lfloorfrac{m}{p} floor + 1)个.只要对这(m\%p)个数单独处理一下就行了。
然后计算没有质数的方案数，全部用合数转移即可。
然后发现这个是可以矩阵优化的。
转移矩阵的第(i)行第(j)列表示从(1 + i)到(m+i)中有多少个数字模(p)余(j)。同样可以计算一下，然后(O(p^2))的列出矩阵。
对于全部用合数转移的情况，转移矩阵可以(O(pfrac{m}{lnm}))转移,亲测吸氧可过

代码

/*
* @Author: wxyww
* @Date:   2019-02-13 16:22:09
* @Last Modified time: 2019-02-13 20:51:05
*/
#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<queue>
#include<vector>
#include<ctime>
using namespace std;
typedef long long ll;
// #define int ll
const int mod = 20170408;
ll read() {
    ll x=0,f=1;char c=getchar();
    while(c<'0'||c>'9') {
        if(c=='-') f=-1;
        c=getchar();
    }
    while(c>='0'&&c<='9') {
        x=x*10+c-'0';
        c=getchar();
    }
    return x*f;
}
int n,m,p;
namespace BF2 {
    const int N = 110;
    struct node {
        int a[N][N],n,m;
        node() {
            memset(a,0,sizeof(a));n = 0,m = 0;
        }
        node(int nn,int mm) {
            memset(a,0,sizeof(a));n = nn;m = mm;
        }
        node(int nn) {
            n = nn;
            memset(a,0,sizeof(a));
            for(int i = 0;i <= nn;++i) a[i][i] = 1;
        }
    };
    node operator * (const node &x,const node &y) {
        int n = x.n,m = y.m;
        node c(n,m);
        int k = x.m;
        for(int i = 0;i <= k;++i) {
            for(int j = 0;j <= n;++j) {
                for(int z = 0;z <= m;++z) {
                    c.a[j][z] += 1ll * x.a[j][i] * y.a[i][z] % mod;
                    c.a[j][z] %= mod;
                }
            }
        }
        return c;
    }
    node operator ^ (node x,int y) {
        node ans(x.n);
        for(;y;y >>= 1,x = x * x)
            if(y & 1) ans = ans * x;
        return ans;
    }
    int dis[2000003],vis[20000003];
    int tot,tmp[2000003];
    int js;
    void Eur() {
        vis[1] = 1;
        for(int i = 2;i <= m;++i) {
            if(!vis[i]) dis[++js] = i;
            for(int j = 1;1ll * dis[j] * i <= m;++j) {
                vis[dis[j] * i] = 1;
                if(dis[j] % i == 0) break;
            }
        }
    }
    void Main() {
        node C(p - 1,p - 1);
        int k = m % p;
        for(int i = 0;i < p;++i) {
            for(int j = i + 1;j <= i + k;++j) {
                C.a[i][j % p]++;
            }
        }
        for(int i = 0;i < p;++i) for(int j = 0;j < p;++j) C.a[i][j] += m / p;
        node ans(0,p - 1);
        ans.a[0][0] = 1;
        ans = ans * (C ^ n);
        int anss = ans.a[0][0];
        memset(ans.a,0,sizeof(ans.a));
        memset(C.a,0,sizeof(C.a));
        ans.a[0][0] = 1;
        Eur();
        for(int i = 0;i < p;++i) {
            for(int j = i + 1;j <= i + k;++j) {
                C.a[i][j % p]++;
            }
        }
        for(int i = 0;i < p;++i) for(int j = 0;j < p;++j) C.a[i][j] += m / p;
        for(int i = 0;i < p;++i)
            for(int j = 1;j <= js;++j)
                C.a[i][(i + dis[j]) % p]--;
        ans = ans * (C ^ n);
        cout<<(anss - ans.a[0][0] + mod) % mod;
    }

}
signed main() {
    n = read(),m = read(),p = read();
    BF2::Main();
    return 0;
}