bzoj 4818 [Sdoi2017]序列计数（简单容斥+快速幂加速dp）

来源：互联网发布：龙腾世纪2 知乎编辑：程序博客网时间：2024/06/09 17:23

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

思路：至少有一个质数转换成全部方案-无质数方案。然后dp[i][j]：序列长度为i且这i个数字的加和%p=j的方案数。写出转移矩阵，快速幂即可。

#include <bits/stdc++.h>using namespace std;typedef long long LL;const int mod = 20170408;int n, m, p;const int maxn = 3e7 + 7;bool check[maxn];const int primeLimit = 1860000 + 5;int prime[primeLimit];int sieve(LL n){    memset(check, false, sizeof(check));//i为素数时， check[i]是false    check[0] = check[1] = true;    int tot = 0;    for(LL i = 2; i <= n; i++)    {        if(!check[i])   prime[tot++] = i;        for(LL j = 0; j < tot && i * prime[j] <= n; j++)        {            check[i * prime[j]] = true;            if(i % prime[j] == 0)   break;        }    }    return tot;}const int mSize = 105;struct Matrix{    int v[mSize][mSize];    friend Matrix operator* (const Matrix& a, const Matrix& b)    {        Matrix c;        memset(c.v, 0, sizeof(c.v));        for (int k = 0; k < mSize; k++)        {            for (int i = 0; i < mSize; i++)            {                if(a.v[i][k] == 0)  continue;                for (int j = 0; j < mSize; j++)                {                    if(b.v[k][j] == 0)  continue;                    c.v[i][j] = c.v[i][j] + 1LL * a.v[i][k] * b.v[k][j] % mod;                    c.v[i][j] %= (mod);                }            }        }        return c;    }    friend Matrix operator^ (Matrix x, LL n)    {        Matrix ans;        for (int i = 0; i < mSize; i++)            for (int j = 0; j < mSize; j++)                ans.v[i][j] = (i == j);        while (n)        {            if (n & 1)   ans = ans * x;            x = x * x;            n >>= 1;        }        return ans;    }};int b[200 + 5];int solve1(){    for(int i = 1; i <= m; i++) b[i % p] ++;    Matrix a;    memset(a.v, 0, sizeof(a.v));    for(int i = 1; i <= m; i++)    {        int num = (p - i % p) % p;        a.v[0][num]++;    }    for(int i = 1; i < p; i++)    {        for(int j = 0; j < p; j++)        {            a.v[i][j] = a.v[i - 1][(j - 1 + p) % p];        }    }    a = a ^ (n-1);    int ans = 0;    for(int i = 0; i < p; i++)    {        ans = ans + 1LL * b[i] * a.v[0][i] % mod;        ans %= mod;    }    return ans;}int solve2(){    memset(b, 0, sizeof(b));    for(int i = 1; i <= m; i++)    {        if(check[i]) b[i % p] ++;    }    Matrix a;    memset(a.v,0,sizeof(a.v));    for(int i = 1; i <= m; i++)    {        if(check[i])        {            int num = (p - i % p) % p;            a.v[0][num]++;        }    }    for(int i = 1; i < p; i++)    {        for(int j = 0; j < p; j++)        {            a.v[i][j] = a.v[i - 1][(j - 1 + p) % p];        }    }    a = a ^ (n-1);    int ans = 0;    for(int i = 0; i < p; i++)    {        ans = ans + 1LL * b[i] * a.v[0][i] % mod;        ans %= mod;    }    return ans;}int main(){    scanf("%d%d%d", &n, &m, &p);    sieve(m + 5);    printf("%d\n", (solve1() - solve2() + mod) % mod);    return 0;}

阅读全文

0 0