线性递推数列
线性递推
对于无限数列 { a 0 , a 1 , . . . } \{a_0,a_1,...\} {a0,a1,...} 和有限非空数列 { r 0 , r 1 , . . . , r m − 1 } \{r_{0},r_1,...,r_{m-1}\} {r0,r1,...,rm−1} 。
若对于任意 m − 1 ≤ n m-1\le n m−1≤n ,有 ∑ i = 0 m − 1 a n − i r i = 0 \sum_{i=0}^{m-1}a_{n-i}r_i=0 ∑i=0m−1an−iri=0 ,则称数列 r r r 为数列 a a a 的线性递归式。
若 r 0 = 1 r_0=1 r0=1,称数列 r r r 为数列 a a a 的线性递推式。
称存在线性递推式的无限数列为线性递推数列。
r 0 = 1 r_0=1 r0=1 说明可以写成 a n = − ( a n − 1 ∗ r 1 + a n − 2 ∗ r 2 + . . . + a n − m + 1 ∗ r m − 1 ) a_n=-(a_{n-1}*r_1+a_{n-2}*r_2+...+a_{n-m+1}*r_{m-1}) an=−(an−1∗r1+an−2∗r2+...+an−m+1∗rm−1)。
对于有限数列 { a 0 , a 1 , . . . a n − 1 } \{a_0,a_1,...a_{n-1}\} {a0,a1,...an−1} 和有限非空数列 { r 0 , r 1 , . . . , r m − 1 } \{r_{0},r_1,...,r_{m-1}\} {r0,r1,...,rm−1} 。
若对于任意 m − 1 ≤ p ≤ n m-1\le p\le n m−1≤p≤n ,有 ∑ i = 0 m − 1 a p − i r i = 0 \sum_{i=0}^{m-1}a_{p-i}r_i=0 ∑i=0m−1ap−iri=0 ,则称数列 r r r 为数列 a a a 的线性递归式。
若 r 0 = 1 r_0=1 r0=1,称数列 r r r 为数列 a a a 的线性递推式。
称这个线性递推式的阶数为它的长度 − 1 -1 −1。(因为下标是从 0 0 0 开始的)
称数列 a a a 阶数最小的线性递推式为数列 a a a 的最短线性递推式。
生成函数
对于有限数列 { a 0 , . . . , a n − 1 } \{a_0,...,a_{n-1}\} {a0,...,an−1},定义它的生成函数为多项式 A ( x ) = ∑ i = 0 n − 1 a i x i A(x)=\sum_{i=0}^{n-1}a_ix^i A(x)=∑i=0n−1aixi。
对于无限数列 { a 0 , a 1 , . . . } \{a_0,a_1,...\} {a0,a1,...},定义它的生成函数为形式幂级数 A ( x ) = ∑ 0 ∞ a i x i A(x)=\sum_{0}^∞a_ix^i A(x)=∑0∞aixi。
对于无限数列 { a 0 , a 1 , . . . } \{a_0,a_1,...\} {a0,a1,...} 和有限非空数列 { r 0 , r 1 , . . . , r m − 1 } \{r_{0},r_1,...,r_{m-1}\} {r0,r1,...,rm−1} 。
设数列 a a a 的生成函数为 A A A,数列 r r r 的生成函数为 R R R。
数列 r r r 为数列 a a a 的线性递归式等价于存在次数不超过 m − 2 m-2 m−2 次的多项式 C C C,满足 A R + C = 0 AR+C=0 AR+C=0。
对于有限数列 { a 0 , a 1 , . . . , a n − 1 } \{a_0,a_1,...,a_{n-1}\} {a0,a1,...,an−1} 和有限非空数列 { r 0 , r 1 , . . . , r m − 1 } \{r_{0},r_1,...,r_{m-1}\} {r0,r1,...,rm−1} 。
设数列 a a a 的生成函数为 A A A,数列 r r r 的生成函数为 R R R。
数列 r r r 为数列 a a a 的线性递归式等价于存在次数不超过 m − 2 m-2 m−2 次的多项式 C C C,满足 A R + C ≡ 0 ( m o d x n ) AR+C\equiv 0\pmod {x^n} AR+C≡0(modxn)。
Berlekamp-Massey
注意:上面的都是下标从 0 0 0 开始,这里的算法流程下标是从 1 1 1 开始的。
Berlekamp-Massey
算法是用来在
O
(
n
2
)
O(n^2)
O(n2) 的时间里求一个数列的最短线性递推式。
方法思想是增量法。
重申明确一下问题模型:给定 n n n 个元素的数列 a 1 , . . . a n a_1,...a_n a1,...an,求一个最短的数列 r 1 , . . . r m r_1,...r_m r1,...rm。要求满足 ∀ m < i ≤ n a i = ∑ j = 1 m a i − j r j \forall_{m<i\le n}\ a_i=\sum_{j=1}^{m}a_{i-j}r_j ∀m<i≤n ai=∑j=1mai−jrj。要求在 O ( n 2 ) O(n^2) O(n2) 时间内解决。
注意1:这个递推的下标可以看作是一种卷积形式,相加和为一个“定值”,这有利于你快速理解下面的一些下标变换与对应。
注意2:这个条件只需要 > m >m >m 的下标做到,换言之如果被当成了 r r r 中的一个“基”,是不需要满足这个式子的,因为会有下标为负的情况。
注意3: a i = a [ i ] = a ( i ) a_i=a[i]=a(i) ai=a[i]=a(i) 只是表达方式不同。
假设递推式已经经过了 t t t 次更新,第 i i i 次更新的递推式记为 r ( i ) r(i) r(i)(这是一个多项式),长度为 m ( i ) m(i) m(i)。初始时,定义 r ( 0 ) = ∅ , m ( 0 ) = 0 r(0)=\empty,m(0)=0 r(0)=∅,m(0)=0。
一个一个考虑加入数列 a a a 中元素,现在即将加入 a i a_i ai。假设现在递推式长度为 m = m ( t ) m=m(t) m=m(t)。
设 Δ ( i ) = a i − ∑ j = 1 m a [ i − j ] ⋅ r ( i ) [ j ] \Delta(i)=a_i-\sum_{j=1}^ma[i-j]·r(i)[j] Δ(i)=ai−∑j=1ma[i−j]⋅r(i)[j]。
-
Δ ( i ) = 0 \Delta(i)=0 Δ(i)=0,证明目前递推式仍然满足条件,合法不用修改,继续考虑加入下一个。
-
Δ ( i ) ≠ 0 \Delta(i)\neq 0 Δ(i)=0。设 fail ( t ) = i \text{fail}(t)=i fail(t)=i,表示 r ( t ) r(t) r(t) 递推式第一次失效的位置在 i i i。
-
t = 0 t=0 t=0。
意味着 a i a_i ai 前面的数都为 0 0 0。那么不管递推数列怎么配求和都是 0 0 0。显然只能将 i i i 这个下标也被递推式包含。
所以新的递推式直接就是 r ( t + 1 ) = { 0 , 0 , . . . , 0 ⏟ i } r(t+1)=\{\underbrace{0,0,...,0}_i\} r(t+1)={i 0,0,...,0}。
这样就不需要 a i a_i ai 满足上面求和的式子了,它已经作为递推式的一个“基”了。
-
t ≠ 0 t\ne 0 t=0。
考虑构造一个递推式 R R R,满足 ∀ ∣ R ∣ < k < i ∑ j = 1 ∣ R ∣ a k − j R j = 0 \forall_{|R|<k<i}\ \sum_{j=1}^{|R|}a_{k-j}R_j=0 ∀∣R∣<k<i ∑j=1∣R∣ak−jRj=0, ∑ j = 1 ∣ R ∣ a i − j R j = Δ ( i ) \sum_{j=1}^{|R|}a_{i-j}R_j=\Delta(i) ∑j=1∣R∣ai−jRj=Δ(i)。
寻找之前某次失效的递推式, 0 ≤ p < t 0\le p<t 0≤p<t,显然这个递推式的失效位置为 fail ( p ) \text{fail}(p) fail(p)。
同时,设 ω = Δ ( i ) Δ ( fail ( p ) ) \omega=\frac{\Delta(i)}{\Delta\big(\text{fail}(p)\big)} ω=Δ(fail(p))Δ(i),则 R = { 0 , 0 , … , 0 ⏟ i − fail ( p ) − 1 , ω , − ω ⋅ r ( p ) [ 1 ] , − ω ⋅ r ( p ) [ 2 ] , . . . , − ω ⋅ r ( p ) [ ∣ r ( p ) ∣ ] } R=\{\underbrace{0,0,…,0}_{i-\text{fail}(p)-1},\omega,-\omega·r(p)[1],-\omega·r(p)[2],...,-\omega·r(p)\Big[\big|r(p)\big|\Big]\} R={i−fail(p)−1 0,0,…,0,ω,−ω⋅r(p)[1],−ω⋅r(p)[2],...,−ω⋅r(p)[∣∣r(p)∣∣]}。
令 r ( t + 1 ) = r ( t ) + R r(t+1)=r(t)+R r(t+1)=r(t)+R 即可。(加法遵循多项式加法原则,即对位系数相加)
中场休息——请确保上面的下标之间关系理清楚了,再继续往下看正确性
为什么这样是对的?理论解释不清楚,不妨眼见为实,我们直接带进去算。
-
对于 a i a_i ai。
R [ 1 , i − fail ( p ) − 1 ] R[1,i-\text{fail}(p)-1] R[1,i−fail(p)−1] 都为 0 0 0,对应 a [ i − 1 , fail ( p ) + 1 ] a[i-1,\text{fail}(p)+1] a[i−1,fail(p)+1] 相乘和为 0 0 0。
真正乘起来有值的是从 R [ i − fail ( p ) ] ↔ a [ fail ( p ) ] R\big[i-\text{fail}(p)\big]\leftrightarrow a\big[\text{fail}(p)\big] R[i−fail(p)]↔a[fail(p)] 项开始的。
并且后面对应相乘, r ( p ) [ 1 ] ↔ a [ fail ( p ) − 1 ] , r ( p ) [ 2 ] ↔ a [ fail(p)-2 ] . . . . . . r(p)[1]\leftrightarrow a[\text{fail}(p)-1],r(p)[2]\leftrightarrow a[\text{fail(p)-2}]...... r(p)[1]↔a[fail(p)−1],r(p)[2]↔a[fail(p)-2]......
等等!!!
这相乘求和不就是第 p p p 个递推式失效的位置对应的表达式吗?转化一下即 = a fail(p) − Δ ( f a i l ( p ) ) =a_{\text{fail(p)}}-\Delta\big(fail(p)\big) =afail(p)−Δ(fail(p))。
⇒ ω ⋅ a fail(p) − ω ⋅ ( a fail(p) − Δ ( fail ( p ) ) ) = ω ∗ Δ ( fail ( p ) ) = Δ ( i ) \Rightarrow \omega·a_{\text{fail(p)}}-\omega·\Big(a_{\text{fail(p)}}-\Delta\big(\text{fail}(p)\big)\Big)=\omega*\Delta\big(\text{fail}(p)\big)=\Delta(i) ⇒ω⋅afail(p)−ω⋅(afail(p)−Δ(fail(p)))=ω∗Δ(fail(p))=Δ(i)。
发现 R R R 把 r ( t ) r(t) r(t) 失效的差距部分恰好补上了。
由此可知 r ( t + 1 ) r(t+1) r(t+1) 对 a i a_i ai 成立。
-
对于 m ≤ k < i m\le k<i m≤k<i。其实推导与上面差不多。
R [ 1 , i − fail ( p ) − 1 ] R[1,i-\text{fail}(p)-1] R[1,i−fail(p)−1] 都为 0 0 0,对应 a [ k − 1 , k − i + fail ( p ) + 1 ] a[k-1,k-i+\text{fail}(p)+1] a[k−1,k−i+fail(p)+1] 相乘和为 0 0 0。
真正乘起来有值的是从 R [ i − fail ( p ) ] ↔ a [ k − i + fail ( p ) ] R\big[i-\text{fail}(p)\big]\leftrightarrow a\big[k-i+\text{fail}(p)\big] R[i−fail(p)]↔a[k−i+fail(p)] 项开始的。
并且后面对应相乘, r ( p ) [ 1 ] ↔ a [ k − i + fail ( p ) − 1 ] , r ( p ) [ 2 ] ↔ a [ k − i + fail(p)-2 ] . . . . . . r(p)[1]\leftrightarrow a[k-i+\text{fail}(p)-1],r(p)[2]\leftrightarrow a[k-i+\text{fail(p)-2}]...... r(p)[1]↔a[k−i+fail(p)−1],r(p)[2]↔a[k−i+fail(p)-2]......
转化一下即 = a k − i + fail ( p ) − Δ ( k − i + fail ( p ) ) =a_{k-i+\text{fail}(p)}-\Delta(k-i+\text{fail}(p)) =ak−i+fail(p)−Δ(k−i+fail(p))。
因为 k < i k<i k<i,所以 k − i + fail ( p ) < fail ( p ) k-i+\text{fail}(p)<\text{fail}(p) k−i+fail(p)<fail(p)。
还记得 fail ( p ) \text{fail}(p) fail(p) 的定义吗?——是 r ( p ) r(p) r(p) 第一个失效的位置。换言之,在此位置之前 r ( p ) r(p) r(p) 都是成立的。
所以 Δ ( k − i + fail ( p ) ) = 0 \Delta(k-i+\text{fail}(p))=0 Δ(k−i+fail(p))=0。
⇒ ω ⋅ a k − i + fail(p) − ω ⋅ ( a k − i + fail(p) − Δ ( k − i + fail ( p ) ) ) = ω ∗ Δ ( k − i + fail ( p ) ) = 0 \Rightarrow \omega·a_{k-i+\text{fail(p)}}-\omega·\Big(a_{k-i+\text{fail(p)}}-\Delta\big(k-i+\text{fail}(p)\big)\Big)=\omega*\Delta\big(k-i+\text{fail}(p)\big)=0 ⇒ω⋅ak−i+fail(p)−ω⋅(ak−i+fail(p)−Δ(k−i+fail(p)))=ω∗Δ(k−i+fail(p))=0。
综上,我们构造的这个 R R R 本质上只起到了对 a i a_i ai 补充 Δ ( i ) \Delta(i) Δ(i) 的效果,对于其余 k k k 的贡献都是 0 0 0。
这是利用了 r ( p ) r(p) r(p) 在 1 ∼ fail ( p ) − 1 1\sim \text{fail}(p)-1 1∼fail(p)−1 都满足关系,而在 fail ( p ) \text{fail}(p) fail(p) 相差 Δ ( p ) \Delta(p) Δ(p) 的性质。
此外我们还希望递推式的长度越短越好,也就是说 max ( m ( t ) , i − fail ( p ) + m ( p ) ) \max\big(m(t),i-\text{fail}(p)+m(p)\big) max(m(t),i−fail(p)+m(p)) 最短。
贪心地只需要动态维护最短的 i − fail ( p ) + m ( p ) i-\text{fail}(p)+m(p) i−fail(p)+m(p),每次算出 r ( t + 1 ) r(t+1) r(t+1) 时都与之前的 r ( p ) r(p) r(p) 比一下谁更短即可。
-
-
一共递推 O ( n ) O(n) O(n) 次,每次修改 O ( n ) O(n) O(n) 次,时间复杂度为 O ( n 2 ) O(n^2) O(n2)。
代码实现中, r r r 用的是 vector \text{vector} vector,下标从 0 0 0 开始,所以有些下标会与上面的推导略有差异。且不保证一定正确!!!
//std参考code
#include <bits/stdc++.h>
using namespace std;
#define mod 1000000007
#define int long long
#define maxn 1005
int n, cnt;
int a[maxn], fail[maxn], delta[maxn];
vector < int > r[maxn];
int qkpow( int x, int y ) {
int ans = 1;
while( y ) {
if( y & 1 ) ans = ans * x % mod;
x = x * x % mod;
y >>= 1;
}
return ans;
}
signed main() {
scanf( "%lld", &n );
for( int i = 1;i <= n;i ++ ) scanf( "%lld", &a[i] );
for( int i = 1;i <= n;i ++ )
if( ! cnt ) {
if( a[i] ) {
fail[cnt ++] = i;
delta[i] = a[i];
r[cnt].resize( i, 0 );
}
continue;
}
else {
fail[cnt] = i;
delta[i] = a[i];
for( int j = 0;j < r[cnt].size();j ++ )
( delta[i] -= a[i - j - 1] * r[cnt][j] ) %= mod;
if( ! delta[i] ) continue;
// int len = 0x7f7f7f7f, p;
int len = i - fail[cnt - 1] + r[cnt - 1].size(), p = cnt - 1;
for( int j = 0;j < cnt;j ++ )
if( i - fail[j] + r[j].size() < len )
len = i - fail[j] + r[j].size(), p = j;
int omega = delta[i] * qkpow( delta[fail[p]], mod - 2 ) % mod;
r[cnt + 1] = r[cnt];
cnt ++;
while( r[cnt].size() < len ) r[cnt].push_back( 0 );
( r[cnt][i - fail[p] - 1] += omega ) %= mod;
for( int j = 0;j < r[p].size();j ++ )
( r[cnt][i - fail[p] + j] -= omega * r[p][j] ) %= mod;
}
printf( "%lld\n", r[cnt].size() );
for( int i : r[cnt] ) printf( "%lld ", ( i + mod ) % mod );
return 0;
}
/*
Input 1
7
1 2 4 9 20 40 90
Output 1
4
0 0 10 0
Input 2
18
2 4 8 16 32 64 128 256 512 2 4 8 16 32 64 128 256 512
Output 2
0 0 0 0 0 0 0 0 1
*/
BM算法和线性递推的应用
当我们确信这个题是 矩阵加速 / 线性递推 的题且矩阵很难构造。
那么可以 暴力递推/打表 出超过两倍矩阵阶数的长度,然后 BM+线性递推 得到这个序列的任意项。
不知道矩阵阶数就 能找多长,找多长 了。
某些特殊情况下也可以代替 DFA的最小化 。
打表题也可以试着找一下规律,因为只能找线性递推的规律,所以大部分情况要灵活处理,比如说先差分再找规律之类的。