題意

Alice想要得到一個長度為n的序列，序列中的數都是不超過m的正整數，而且這n個數的和是p的倍數。
Alice還希望，這n個數中，至少有一個數是質數。
Alice想知道，有多少個序列滿足她的要求。
對100%的資料，1≤n≤109,1≤m≤2×107,1≤p≤100

分析

一眼容斥，用所有方案減去不含質數的方案。
設f[i,j]表示序列前i個數模p的餘數為j時的方案數。f[i,j]=∑f[i−1,(j−k)modp]
構建矩陣然後快速冪即可。
求不含質數的方案同理。
在一開始建矩陣的時候我的複雜度是mp，發現會超時。但是注意到矩陣的每一列都是迴圈同構的，於是就可以先O(m)處理好一列，然後再構建其他列即可。

程式碼

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;

typedef long long LL;

const int M=20000005;
const int P=105;
const int MOD=20170408;

int n,m,p,prime[M/10],tot,f[P];
struct arr{int a[P][P];}ans,a;
bool not_prime[M];

void
 
 get_prime(int n)
{
    not_prime[1]=1;
    for (int i=2;i<=n;i++)
    {
        if (!not_prime[i]) prime[++tot]=i;
        for (int j=1;j<=tot&&i*prime[j]<=n;j++)
        {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0) break;
        }
    }
}

void mul(arr &c,arr a,arr b)
{
    memset
 
(c.a,0,sizeof(c.a));
    for (int i=0;i<p;i++)
        for (int j=0;j<p;j++)
            for (int k=0;k<p;k++)
                c.a[i][j]=(c.a[i][j]+(LL)a.a[i][k]*b.a[k][j]%MOD)%MOD;
}

arr ksm(arr x,int y)
{
    memset(ans.a,0,sizeof(ans.a));
    for (int i=0;i<p;i++) ans.a[i][i]=1;
    while (y)
    {
        if (y&1) mul(ans,ans,x);
        mul(x,x,x);y>>=1;
    }
    return ans;
}

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

int solve2()
{
    memset(f,0,sizeof(f));
    for (int i=1;i<=m;i++) if (not_prime[i]) f[i%p]++;
    memset(a.a,0,sizeof(a.a));
    for (int j=1;j<=m;j++)
        if (not_prime[j]) a.a[(-j%p+p)%p][0]++;
    for (int i=1;i<p;i++)
        for (int j=0;j<p;j++)
            a.a[j][i]=a.a[(j-1+p)%p][i-1];
    a=ksm(a,n-1);
    int ans=0;
    for (int i=0;i<p;i++) ans=(ans+(LL)f[i]*a.a[i][0]%MOD)%MOD;
    return ans;
}

int main()
{
    scanf("%d%d%d",&n,&m,&p);
    get_prime(m);
    int x=solve1(),y=solve2();
    printf("%d",(x-y+MOD)%MOD);
    return 0;
}