BZOJ5298 CQOI2018 交错序列 【DP+矩阵快速幂优化】
Description
我们称一个仅由0、1构成的序列为”交错序列”,当且仅当序列中没有相邻的1(可以有相邻的0)。例如,000,001,101,都是交错序列,而110则不是。对于一个长度为n的交错序列,统计其中0和1出现的次数,分别记为x和y。给定参数a、b,定义一个交错序列的特征值为 xayb x a y b 。注意这里规定任何整数的0次幂都等于1(包括 00=1 0 0 = 1 )。显然长度为n的交错序列可能有多个。我们想要知道,所有长度为n的交错序列的特征值的和,除以m的余数。(m是一个给定的质数)例如,全部长度为3的交错串为:000、001、010、100、101。当a=1,b=2时,可计算 31∗02+21∗12+21∗12+21∗12+11∗22=10 3 1 ∗ 0 2 + 2 1 ∗ 1 2 + 2 1 ∗ 1 2 + 2 1 ∗ 1 2 + 1 1 ∗ 2 2 = 10
Input
输入文件共一行,包含三个空格分开的整数n,a,b和m。
1≤n≤10000000,0≤a,b≤45,m<100000000
Output
输出文件共一行,为计算结果。
Sample Input
3 1 2 1009
Sample Output
10
首先我们看到这题,果断想到DP,但是
xayb
x
a
y
b
形式的答案十分不好维护,所以我们将它进行转化
∑xayb
∑
x
a
y
b
=∑xa(n−x)b
=
∑
x
a
(
n
−
x
)
b
=∑xa∑bi=0Cibni(−x)b−i
=
∑
x
a
∑
i
=
0
b
C
b
i
n
i
(
−
x
)
b
−
i
=∑bi=0Cibni(−1)b−i∑xa+b−i
=
∑
i
=
0
b
C
b
i
n
i
(
−
1
)
b
−
i
∑
x
a
+
b
−
i
现在我们只需要在优秀时间内算出
∑xa+b−i
∑
x
a
+
b
−
i
就好了,那么我们定义f[i][j][0/1]表示前i位的所有合法方案的0的个数的j次方之和,且第i位是0/1。
那么我们可以得到:
f[i][j][1]=f[i−1][j][0]
f
[
i
]
[
j
]
[
1
]
=
f
[
i
−
1
]
[
j
]
[
0
]
f[i][j][0]=∑jk−0Ckj∗(f[i−1][k][0]+f[i−1][k][1])
f
[
i
]
[
j
]
[
0
]
=
∑
k
−
0
j
C
j
k
∗
(
f
[
i
−
1
]
[
k
]
[
0
]
+
f
[
i
−
1
]
[
k
]
[
1
]
)
式子中的
Ckj
C
j
k
直接由二项式定理计算,然后我们很显然的发现DP方程是可以矩阵优化的。注意矩阵乘法的时候稍稍卡一下常数就好了。
时间
O((a+b)3∗log(n))
O
(
(
a
+
b
)
3
∗
l
o
g
(
n
)
)
空间
O((a+b)2)
O
(
(
a
+
b
)
2
)
#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define N 200
int read(){
int ans=0,w=1;
char c=getchar();
while(c!='-'&&!isdigit(c))c=getchar();
if(c=='-')w=-1,c=getchar();
while(isdigit(c))ans=ans*10+c-'0',c=getchar();
return ans*w;
}
int n,a,b,Mod;
LL J[N],inv[N];
struct Matrix{
int n;LL a[N][N];
Matrix(int n):n(n){
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
a[i][j]=0;
}
Matrix operator * (const Matrix &c)const{
Matrix ans=Matrix(n);
for(int i=0;i<n;i++)
for(int k=0;k<n;k++)if(a[i][k])
for(int j=0;j<n;j++)if(c.a[k][j])
ans.a[i][j]=(ans.a[i][j]+a[i][k]*c.a[k][j])%Mod;
return ans;
}
Matrix operator ^ (const int x)const{
Matrix a=*this,ans=Matrix(n);
for(int i=0;i<n;i++)ans.a[i][i]=1;
int b=x;
while(b){
if(b&1)ans=ans*a;
b>>=1;
a=a*a;
}
return ans;
}
};
LL C(int n,int m){
if(n<m)return 0;
return J[n]*inv[m]%Mod*inv[n-m]%Mod;
}
LL fast_pow(LL a,LL b){
LL ans=1;
while(b){
if(b&1)ans=ans*a%Mod;
b>>=1;
a=a*a%Mod;
}
return ans;
}
int main(){
n=read();a=read();b=read();Mod=read();
int MAX=a+b+1;
J[0]=1;
for(int i=1;i<=MAX;i++)J[i]=J[i-1]*i%Mod;
inv[MAX]=fast_pow(J[MAX],Mod-2);
for(int i=MAX-1;i>=0;i--)inv[i]=inv[i+1]*(i+1)%Mod;
Matrix move=Matrix(MAX*2);
for(int i=0;i<MAX;i++)move.a[i][i+MAX]=1;
for(int i=0;i<MAX;i++)
for(int j=0;j<=i;j++)
move.a[j][i]=move.a[j+MAX][i]=C(i,j);
move=move^n;
LL x=1,ans=0;
for(int i=0;i<=b;i++,x=x*n%Mod)
ans=(ans+C(b,i)*x%Mod*(move.a[0][a+b-i]+move.a[0][a+b-i+MAX])%Mod*(((b-i)&1)?(-1):(1)))%Mod;
printf("%lld",(ans+Mod)%Mod);
return 0;
}
多谢lyw大神指点 → → lyw大神