洛谷 P3702 [SDOI2017]序列计数 dp+快速幂

题目描述

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

Alice还希望，这 $n$ 个数中，至少有一个数是质数。

Alice想知道，有多少个序列满足她的要求。

输入输出格式

输入格式：
一行三个数， $n,m,p$。

输出格式：
一行一个数，满足Alice的要求的序列数量，答案对 $20170408$ 取模。

输入输出样例

输入样例#1：
3 5 3
输出样例#1：
33
说明

对 20\% 的数据， $1≤n,m≤100$

对 50\% 的数据， $1≤m≤100$

对 80\% 的数据，$1≤m≤10^6$

对 100\% 的数据， $1≤n≤10^9 ,1≤m≤2×10^7,1≤p≤100$

时间限制：3s

空间限制：128MB

分析：
显然可以用所有序列减去全不为质数序列。
设$f[i][j]$为$i$个数，模$p$等于$j$的方案数，有

$$f[x+y][i]=\sum_{i=j+k}f[x][j]*f[y][k]$$
显然是一个卷积形式，而且每次$x$可以翻倍，$p$很小，直接暴力卷积然后快速幂就好了。

代码：

#include <iostream>
#include <cstdio>
#include <cmath>
#define LL long long

const LL maxn=1007;
const LL mod=20170408; 

using namespace std;

LL n,m,cnt,p;
LL f[maxn],g[maxn],f1[maxn],g1[maxn],a[maxn],b[maxn];
LL prime[50007];
int not_prime[20000007];

void getprime(LL n)
{
    for (LL i=2;i<=n;i++)
    {
        if (!not_prime[i])
        {
            prime[++cnt]=i;
        }
        for (LL j=1;j<=cnt;j++)
        {
            if (i*prime[j]>n) break;
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0) break;
        }
    }
}

void calc(LL *x,LL *y,LL *z,LL n,LL m)
{
    for (LL i=0;i<n;i++) a[i]=x[i];
    for (LL i=0;i<m;i++) b[i]=y[i];
    for (LL i=0;i<n+m;i++) z[i]=0;
    for (LL i=0;i<n;i++)
    {
        for (LL j=0;j<m;j++)
        {
            z[(i+j)%p]=(z[(i+j)%p]+a[i]*b[j]%mod)%mod;
        }
    }
}

void solve(LL n,LL m)
{
    if (n==1) return;
    solve(n/2,m);
    calc(f,f,f,m,m);
    calc(g,g,g,m,m);
    if (n%2)
    {
        calc(f,f1,f,m,m);
        calc(g,g1,g,m,m);
    }
}

int main()
{
    scanf("%lld%lld%lld",&n,&m,&p);
    getprime(m);
    LL j=1;
    for (LL i=1;i<=m;i++)
    {
        while (prime[j]<i) j++;
        if (i!=prime[j])
        {
            g[i%p]++;
            g1[i%p]++;
        }
        f[i%p]++;
        f1[i%p]++;
    }   
    solve(n,p);
    printf("%lld",(f[0]+mod-g[0])%mod);
}