P1025 · OFFICIAL SOLUTION

P1025 万象启星环官方题解

Gioush OJ · P1025 万象启星环

万象启星环

【题目简述】

求长度为 $n$ 的有序序列 $(x_1,x_2,\ldots,x_n)$ 的数量，满足 $1\le x_i\le m$ ， $\sum x_i$ 是 $p$ 的倍数，且至少有一个 $x_i$ 为质数。答案对 $998{,}244{,}353$ 取模。

【数据点 1∼41\sim41∼4】

【直接 DP】

直接统计“至少一个质数”不方便，因此改为“所有序列数”减去“不含质数的序列数”。令 $f_{i,r,0/1}$ 表示前 $i$ 个数已确定、和模 $p$ 为 $r$ 的方案数，最后一维区分可选全集与只选非质数。

【转移】

对每个 $j\in[1,m]$ 枚举加入的数，转移到 $(r+j)\bmod p$ 。复杂度为 $\mathcal{O}(nmp)$ 。

【数据点 5∼85\sim85∼8】

【重复转移】

滚动数组后，每一层都执行同一个线性变换。令列向量 $\bm{x}_i$ 记录各余数的方案数，若 $\bm{A}$ 表示一次转移，则 $\bm{x}_n=\bm{A}^n\bm{x}_0$ 。

【矩阵】

令

\operatorname{Count}_{0,r}=\left|\left\{x\mid 1\le x\le m,\ x\equiv r\pmod p\right\}\right|,

\operatorname{Count}_{1,r}=\left|\left\{x\mid 1\le x\le m,\ x\equiv r\pmod p,\ x\text{ 不是质数}\right\}\right|.

转移矩阵第 $j,k$ 项为 $\operatorname{Count}_{(j-k+p)\bmod p}$ 。分别对全集与非质数集合快速幂，答案是两者余数 $0$ 的方案数之差。

【正解】

【计算计数】

$\operatorname{Count}_{0,r}$ 可以用带余除法直接求出；再用线性筛找出 $[1,m]$ 中所有质数，从对应余数类扣除即可得到 $\operatorname{Count}_{1,r}$ 。

【解法】

构造两张 $p\times p$ 循环矩阵，分别快速幂 $n$ 次。初始向量只有余数 $0$ 的分量为 $1$ 。两次转移后取余数 $0$ 分量相减。

【复杂度分析】

线性筛为 $\mathcal{O}(m)$ ，矩阵快速幂为 $\mathcal{O}(p^3\log n)$ ，总复杂度为 $\mathcal{O}(m+p^3\log n)$ 。

【参考代码】

1#include <cstdio>2#include <cstring>3#include <algorithm>4#include <iostream>5using namespace std;6const int Mod=998244353;7inline int Add(int a,int b){return a+b>=Mod?a+b-Mod:a+b;}8inline int Mul(int a,int b){return 1ll*a*b%Mod;}9inline int Del(int a,int b){return a-b<0?a-b+Mod:a-b;}10int T,n,m,p;11int Prime[2000100],cnt=0,Count[2][110];12bool Tag[20000010];13void Tackle() {14    Count[0][1%p]++;Count[1][1%p]++;15    for (int i=2;i<=m;i++) {16        Count[0][i%p]++;17        if (!Tag[i]) {18            Prime[++cnt]=i;19        }20        else Count[1][i%p]++;21        for (int j=1;j<=cnt;j++) {22            if (1ll*i*Prime[j]>m) {break;}23            Tag[i*Prime[j]]=true;24            if (i%Prime[j]==0) {break;}25        }26    }27    return;28}29struct Matrix {30    int Num[110][110];31    void Init() {32        memset(Num,0,sizeof(Num));33    }34}I,Ans1,f,g,Ans2;35Matrix operator*(const Matrix &A,const Matrix &B) {36    Matrix C;37    for(int i=1;i<=p;i++) {38        for(int j=1;j<=p;j++) {39            C.Num[i][j]=0;40            for(int k=1;k<=p;k++) {41                C.Num[i][j]=Add(C.Num[i][j],Mul(A.Num[i][k],B.Num[k][j]));42            }43        }44    }45    return C;46}47void Print(Matrix A) {48    for(int i=1;i<=p;i++) {49        for(int j=1;j<=p;j++) {50            cout<<A.Num[i][j]<<" ";51        }52        puts("");53    }54    return;55}56Matrix operator^(Matrix A,int b) {57    Matrix C;58    C=I;59    while(b) {60        if (b&1) C=C*A;61        A=A*A;62        b>>=1;63    }64    return C;65}66Matrix operator-(const Matrix &A,const Matrix &B) {67    Matrix C;68    C.Init();69    for(int i=1;i<=p;i++) {70        for(int j=1;j<=p;j++) {71            C.Num[i][j]=Del(A.Num[i][j],B.Num[i][j]);72        }73    }74    return C;75}76void Solve() {77    scanf("%d%d%d",&n,&m,&p);78    cnt=0;79    memset(Count,0,sizeof(Count));80    memset(Tag,0,(m+1)*sizeof(Tag[0]));81    I.Init();82    Ans1.Init();Ans2.Init();83    for (int i=1;i<=p;i++) I.Num[i][i]=1;84    Tackle();85    f.Init();g.Init();86    Ans1.Num[p][1]=1;Ans2.Num[p][1]=1;87    for (int i=1;i<=p;i++) {88        for (int j=1;j<=p;j++) {89            f.Num[i][j]=Count[0][(j-i+p)%p];90        }91    }92    for (int i=1;i<=p;i++) {93        for (int j=1;j<=p;j++) {94            g.Num[i][j]=Count[1][(j-i+p)%p];95        }96    }97    Matrix Ans=(f^n)*Ans1-(g^n)*Ans2;98    printf("%d\n",Ans.Num[p][1]);99}100int main() {101    scanf("%d",&T);102    while(T--) Solve();103    return 0;104}