[SDOI2017]序列计数
题意
Alice 想要得到一个长度为 \(n\) 的正整数序列 \(\{a_i\}\)，满足：
- \(a_i\le m\)
- \(p|\sum a_i\)
- \(\exists i\in[1,n]\)，使 \(a_i\) 是质数
Alice 想知道，有多少个序列满足她的要求。
朴素
由于题目有【\(n\) 个数】【和】【有否质数】这三个主要限制，因此：
设 \(F(i,j,1/0)\) 表示前 \(i\) 个数，和为 \(j\)，是否有了一个质数。
列出朴素 dp 方程：
\[F(i,j,0)=\sum_{1\le k\le m}F(i-1,j-k,0)\times [k \text{ isn't prime}]\\
F(i,j,1)=\sum_{1\le k\le m}F(i-1,j-k,1)+F(i-1,j-k,0)\times [k\text{ is prime}]
\]
所求即为 \(\sum_{p|i,i\le nm}F(n,i,1)\)。
发现具体的【和】不重要，因为只关心【和 \(\bmod p\)】，所以 \(j,k\) 均可在 \(\pmod p\) 意义下枚举。
设 \(A_k,B_k,C_k\) 分别表示 \(\mod p=k\) 的质数、非质数、数的个数。
\[F(i,j,0)=\sum_{1\le k\le m}F(i-1,j-k,0)\times B_k\\
F(i,j,1)=\sum_{1\le k\le m}F(i-1,j-k,1)\times C_k+F(i-1,j-k,0)\times A_k
\]
所求即为 \(F(n,0,1)\)。
矩阵加速
每次转移在做一样的事情，所以矩阵快速幂可以用来加速转移。
dp 矩阵：
\[\left[
\begin{matrix}
F_{0,0}\ F_{1,0}\ F_{2,0}\ \cdots\ F_{p-1,0}\ F_{0,1}\ \cdots \ F_{p-1,1}
\end{matrix}
\right]\\
\]
转移矩阵：
| 编号 | 1 | 2 | p | 1 | 2 | p |
|---|---|---|---|---|---|---|
| 1 | A0 | A1 | … | B0 | B1 | … |
| 2 | Ap-1 | A0 | … | Bp-1 | B0 | … |
| … | Ap-1 | … | … | Bp-1 | … | |
| … | … | … | … | … | … | |
| A3 | … | … | … | … | … | |
| A2 | A3 | … | … | … | … | |
| p | A1 | A2 | … | B1 | B2 | … |
| p+1 | C0 | C1 | … | 0 | … | … |
| Cp-1 | C0 | … | 0 | … | … | |
| … | Cp-1 | … | … | … | … | |
| … | … | … | … | … | … | |
| … | … | … | … | … | … | |
| … | … | … | 0 | … | … | |
| 2p | C1 | C2 | … | 0 | … | … |
这样就可以利用 \(O(p^3\log n)\) 时间得到，其中 \(\log n\) 为快速幂，\(p^3\) 为两个 \((2p)^2\) 的矩阵的乘法。可以通过此题。
#include <bits/stdc++.h>
using namespace std;
const int N=205,M=2e7+5,mod=20170408;
int n,m,p,tot,a[N][N],f[N],prime[M],A[N],B[N],C[N];
bool v[M];
void mul(int f[N],int a[N][N]){
int c[N];memset(c,0,sizeof(c));
for(int i=1;i<=2*p;i++)
for(int j=1;j<=2*p;j++)
c[i]+=1ll*f[j]*a[j][i]%mod,c[i]%=mod;
memcpy(f,c,sizeof(c));
}
void mul2(int a[N][N],int b[N][N]){
int c[N][N];memset(c,0,sizeof(c));
for(int i=1;i<=2*p;i++)
for(int j=1;j<=2*p;j++)
for(int k=1;k<=2*p;k++)
c[i][j]+=1ll*a[i][k]*b[k][j]%mod,c[i][j]%=mod;
memcpy(a,c,sizeof(c));
}
void qp(int a[N][N],int b){
int c[N][N];memset(c,0,sizeof(c));
for(int i=1;i<=2*p;i++)c[i][i]=1;
while(b){
if(b&1)mul2(c,a);
b>>=1,mul2(a,a);
}
memcpy(a,c,sizeof(c));
}
int main(){
cin>>n>>m>>p;
v[1]=1;
for(int i=2;i<=m;i++){
if(!v[i])prime[++tot]=i;
for(int j=1;j<=m&&prime[j]*i<=m;j++){
v[prime[j]*i]=1;
if(i%prime[j]==0)break;
}
}
for(int i=1;i<=m;i++)A[i%p]+=!v[i],B[i%p]+=v[i],C[i%p]++;
for(int i=1;i<=p;i++)for(int j=1;j<=p;j++)a[i][j]=B[(p-i+j)%p];
for(int i=1;i<=p;i++)for(int j=1;j<=p;j++)a[i][j+p]=A[(p-i+j)%p];
for(int i=1;i<=p;i++)for(int j=1;j<=p;j++)a[i+p][j+p]=C[(p-i+j)%p];
f[1]=1;
qp(a,n);
mul(f,a);
cout<<f[p+1];
}
进一步优化
\(p\) 再大一点就过不了了，但其实完全可以再优化。
观察转移矩阵，其实我们是知道的，第 \(a_{1,2}...a_{p,2}\) 其实就是 \(a_{1,1}...a_{p,1}\) 向下平移一格（环状平移）得到的。其它三个矩形区域也是类似。即便在做了快速幂之后，这种格局也是不会改变的。考虑将这些很类似的列中只取一个 \(p\times 1\) 的列作为代表进行计算这样只用 \(O(p^2)\)，并再用 \(O(p^2)\) 进行复制、平移。具体地，按原来矩乘方式计算结果矩阵中 \((1,1)...(p,1);(p+1,1)...(2p,1);(1,p+1)...(p,p+1);(p+1,p+1)...(2p,p+1)\) 这些元素的值，然后按照规则填充空白部分。复杂度 \(O(p^2\log n)\)。更加可以通过。
#include <bits/stdc++.h>
using namespace std;
const int N=205,M=2e7+5,mod=20170408;
int n,m,p,tot,a[N][N],f[N],prime[M],A[N],B[N],C[N];
bool v[M];
void mul(int f[N],int a[N][N]){
int c[N];memset(c,0,sizeof(c));
for(int i=1;i<=2*p;i++)
for(int j=1;j<=2*p;j++)
c[i]+=1ll*f[j]*a[j][i]%mod,c[i]%=mod;
memcpy(f,c,sizeof(c));
}
void mul2(int a[N][N],int b[N][N]){
int c[N][N];memset(c,0,sizeof(c));
for(int i=1;i<=2*p;i++)
for(int j=1;j<=2*p;j++)
c[i][1]+=1ll*a[i][j]*b[j][1]%mod,c[i][1]%=mod;
for(int i=1;i<=2*p;i++)
for(int j=1;j<=2*p;j++)
c[i][p+1]+=1ll*a[i][j]*b[j][p+1]%mod,c[i][p+1]%=mod;
for(int j=2;j<=p;j++){
for(int i=1;i<=p;i++)c[i][j]=c[(i-2+p)%p+1][j-1];
for(int i=1;i<=p;i++)c[i+p][j]=c[(i-2+p)%p+1+p][j-1];
}
for(int j=p+2;j<=2*p;j++){
for(int i=1;i<=p;i++)c[i][j]=c[(i-2+p)%p+1][j-1];
for(int i=1;i<=p;i++)c[i+p][j]=c[(i-2+p)%p+1+p][j-1];
}
memcpy(a,c,sizeof(c));
}
void qp(int a[N][N],int b){
int c[N][N];memset(c,0,sizeof(c));
for(int i=1;i<=2*p;i++)c[i][i]=1;
while(b){
if(b&1)mul2(c,a);
b>>=1,mul2(a,a);
}
memcpy(a,c,sizeof(c));
}
int main(){
cin>>n>>m>>p;
v[1]=1;
for(int i=2;i<=m;i++){
if(!v[i])prime[++tot]=i;
for(int j=1;j<=m&&prime[j]*i<=m;j++){
v[prime[j]*i]=1;
if(i%prime[j]==0)break;
}
}
for(int i=1;i<=m;i++)A[i%p]+=!v[i],B[i%p]+=v[i],C[i%p]++;
for(int i=1;i<=p;i++)for(int j=1;j<=p;j++)a[i][j]=B[(p-i+j)%p];
for(int i=1;i<=p;i++)for(int j=1;j<=p;j++)a[i][j+p]=A[(p-i+j)%p];
for(int i=1;i<=p;i++)for(int j=1;j<=p;j++)a[i+p][j+p]=C[(p-i+j)%p];
f[1]=1;
qp(a,n);
mul(f,a);
cout<<f[p+1];
}
