「CQOI 2018」交错序列

Problem

Description

我们称一个 01 序列为"交错序列",当且仅当序列中不存在两个相邻的 1
000,001,101 都是交错序列,而 110 则不是。
对于一个交错序列,记其中 0 出现的次数为 x 1 出现的次数为 y
给定参数 a,b ,定义一个"交错序列"的特征值为 x^{a}y^{b}

这里规定任何整数的 0 次幂都等于 1 ,包括 0^{0}=1

求所有长度为 n 的交错序列的特征值的和。

Constraints

1\le n\le 10^{7},0\le a,b\le 45

Solution

Analysis

由于 x^{a}y^{b} 中有两个变量较难处理,考虑将 x 消去。
容易发现有 x+y=n x=n-y ,则

x^{a}y^{b}=y^{b}\sum_{i=0}^{a}(-1)^{a-i}\binom{a}{i}n^{i}y^{a-i}=\sum_{i=0}^{a}(-1)^{a-i}\binom{a}{i}n^{i}y^{a+b-i}

答案即为:

\sum_{p}\sum_{i=0}^{a}(-1)^{a-i}\binom{a}{i}n^{i}y^{a+b-i}=\sum_{i=0}^{a}\left((-1)^{a-i}\binom{a}{i}n^{i}\sum_{p}y^{a+b-i}\right)

其中 (-1)^{a-i}\begin{pmatrix}a\\i\end{pmatrix}n^{i} 可以 O\left(a^{2}\right) 处理后 O(1) 求出。

考虑如何求出 \sum_{p}y^{a+b-i}
f\left(len,i,v\right) 为所有长为 len 且以 v 结尾的交错序列的 y^{i} 之和。
f\left(len,i,0\right) 可以从 f\left(len-1,i,0\right) f\left(len-1,i,1\right) 转移而来且 1 的出现次数不变,即:

f\left(len,i,0\right)=f\left(len-1,i,0\right)+f\left(len-1,i,1\right)

f\left(len,i,1\right) 只能从 f\left(len-1,i,0\right) 转移而来且 1 的出现次数增加 1 ,即:

\begin{aligned} f\left(len,i,1\right)&=\sum\left(y'+1\right)^{i}\\ &=\sum\sum_{j=0}^{i}\binom{i}{j}\left(y'\right)^{j}\\ &=\sum_{j=0}^{i}\binom{i}{j}\sum\left(y'\right)^{j}\\ &=\sum_{j=0}^{i}\binom{i}{j}f\left(len-1,j,0\right) \end{aligned}

初始状态为 f\left(0,0,0\right)=1 ,可以用矩阵进行加速。
将转移矩阵填上后矩阵快速幂即可。

时间复杂度 O\left(\left(a+b\right)^{3}\log{n}\right)

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#include<cstdio>
#include<cstring>

using i64=long long;
int p;
template<typename _1,typename _2>inline void inc(_1&x,const _2&y){(x>=p-y)?(x+=y-p):(x+=y);}
template<typename _1,typename _2>inline void dec(_1&x,const _2&y){(x-=y);(x<0)&&(x+=p);}
template<typename T>inline T Plus(const T&x,const T&y){return(x>=p-y)?(x+y-p):(x+y);}

#define Maxn 165
const i64 Limit=(1ULL<<63)-1ULL;

template<typename T>struct Matrix{
int R,C;
T v[Maxn][Maxn];
Matrix():R(0),C(0){memset(v,0,sizeof v);}
Matrix(const int&_R,const int&_C):R(_R),C(_C){memset(v,0,sizeof v);}

inline T*operator[](const int&idx)const{return(T*)v[idx];}

inline Matrix<T>operator*(const Matrix<T>&rhs){
Matrix<T>Ans(R,rhs.C);
for(int i=0;i<R;++i)
for(int j=0;j<rhs.C;++j){
i64 cur=0;
for(int k=0;k<C;++k){
cur+=(i64)v[i][k]*rhs[k][j];
if(Limit<=cur)cur%=p;
}Ans[i][j]=cur%p;
}
return Ans;
}

};

int C[Maxn][Maxn];
i64 pown[Maxn];

int main(){
int n,a,b;
scanf("%d%d%d%d",&n,&a,&b,&p);
int c=a+b+1;
int d=c<<1;

for(int i=C[0][0]=1;i<c;++i)
for(int j=C[i][0]=C[i][i]=1;j<i;++j)
C[i][j]=Plus(C[i-1][j-1],C[i-1][j]);

for(int i=pown[0]=1;i<=a;++i)
pown[i]=(i64)pown[i-1]*n%p;

Matrix<int>trans(d,d);
for(int i=0;i<c;++i){
trans[i][i]=trans[i+c][i]=1;
for(int j=i;j<c;++j)
trans[i][j+c]=C[j][i];
}

Matrix<int>coef(1,d);
for(coef[0][0]=1;n;n>>=1,trans=trans*trans)
if(n&1)coef=coef*trans;

int Ans=0;
void(*opt[2])(int&,const int&)={inc,dec};
for(int i=0;i<=a;++i){
i64 coe=C[a][i]*pown[i]%p;
int Fac=Plus(coef[0][a+b-i],coef[0][a+b-i+c]);
opt[(a-i)&1](Ans,coe*Fac%p);
}

return(printf("%d",Ans),0);
}