[线性代数学习笔记] 线性递推数列及 Berlekamp-Massey 算法的详细推导过程

线性递推数列

线性递推

对于无限数列 { 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\le n$ ，有 ${\sum }_{i=0}^{m-1}{a}_{n-i}{r}_i=0$ ，则称数列 $r$ 为数列 $a$ 的线性递归式。

若 $r_0=1$，称数列 $r$ 为数列 $a$ 的线性递推式。

称存在线性递推式的无限数列为线性递推数列。

$r_0=1$ 说明可以写成 $a_n=-(a_{n-1}\ast r_1+a_{n-2}\ast r_2+...+a_{n-m+1}\ast r_{m-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\le p\le n$ ，有 ${\sum }_{i=0}^{m-1}{a}_{p-i}{r}_i=0$ ，则称数列 $r$ 为数列 $a$ 的线性递归式。

若 $r_0=1$，称数列 $r$ 为数列 $a$ 的线性递推式。

称这个线性递推式的阶数为它的长度 $-1$。（因为下标是从 $0$ 开始的）

称数列 $a$ 阶数最小的线性递推式为数列 $a$ 的最短线性递推式。

生成函数

对于有限数列 { a 0 , . . . , a n − 1 } \{a_0,...,a_{n-1}\} {a0​,...,an−1​}，定义它的生成函数为多项式 $A(x)={\sum }_{i=0}^{n-1}{a}_i{x}^i$。

对于无限数列 { a 0 , a 1 , . . . } \{a_0,a_1,...\} {a0​,a1​,...}，定义它的生成函数为形式幂级数 $A(x)={\sum }_0^{\infty }{a}_i{x}^i$。

对于无限数列 { 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$，数列 $r$ 的生成函数为 $R$。

数列 $r$ 为数列 $a$ 的线性递归式等价于存在次数不超过 $m-2$ 次的多项式 $C$，满足 $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$，数列 $r$ 的生成函数为 $R$。

数列 $r$ 为数列 $a$ 的线性递归式等价于存在次数不超过 $m-2$ 次的多项式 $C$，满足 $AR+C\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

Berlekamp-Massey

注意：上面的都是下标从 $0$ 开始，这里的算法流程下标是从 $1$ 开始的。

Berlekamp-Massey 算法是用来在 $O(n^2)$ 的时间里求一个数列的最短线性递推式。

方法思想是增量法。

重申明确一下问题模型：给定 $n$ 个元素的数列 $a_1,...a_n$，求一个最短的数列 $r_1,...r_m$。要求满足 ${\forall }_{m<i\le n}{a}_i={\sum }_{j=1}^m{a}_{i-j}{r}_j$。要求在 $O(n^2)$ 时间内解决。

注意1：这个递推的下标可以看作是一种卷积形式，相加和为一个“定值”，这有利于你快速理解下面的一些下标变换与对应。

注意2：这个条件只需要 $>m$ 的下标做到，换言之如果被当成了 $r$ 中的一个“基”，是不需要满足这个式子的，因为会有下标为负的情况。

注意3：$a_i=a[i]=a(i)$ 只是表达方式不同。

假设递推式已经经过了 $t$ 次更新，第 $i$ 次更新的递推式记为 $r(i)$（这是一个多项式），长度为 $m(i)$。初始时，定义 $r(0)=\varnothing ,m(0)=0$。

一个一个考虑加入数列 $a$ 中元素，现在即将加入 $a_i$。假设现在递推式长度为 $m=m(t)$。

设 $\Delta (i)=a_i-{\sum }_{j=1}^{m}a[i-j]\cdot r(i)[j]$。

• $\Delta (i)=0$，证明目前递推式仍然满足条件，合法不用修改，继续考虑加入下一个。

• $\Delta (i)\ne 0$。设 $\text{fail}(t)=i$，表示 $r(t)$ 递推式第一次失效的位置在 $i$。

  • $t=0$。

    意味着 $a_i$ 前面的数都为 $0$。那么不管递推数列怎么配求和都是 $0$。显然只能将 $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$ 满足上面求和的式子了，它已经作为递推式的一个“基”了。

  • $t\ne 0$。

    考虑构造一个递推式 $R$，满足 ${\forall }_{|R|<k<i}{\sum }_{j=1}^{|R|}{a}_{k-j}{R}_j=0$，${\sum }_{j=1}^{|R|}{a}_{i-j}{R}_j=\Delta (i)$。

    寻找之前某次失效的递推式，$0\le p<t$，显然这个递推式的失效位置为 $\text{fail}(p)$。

    同时，设 $\omega =\frac{\Delta (i)}{\Delta (\text{fail}(p))}$，则 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$ 即可。（加法遵循多项式加法原则，即对位系数相加）

    ---

    中场休息——请确保上面的下标之间关系理清楚了，再继续往下看正确性

    ---

    为什么这样是对的？理论解释不清楚，不妨眼见为实，我们直接带进去算。

    • 对于 $a_i$。

      $R[1,i-\text{fail}(p)-1]$ 都为 $0$，对应 $a[i-1,\text{fail}(p)+1]$ 相乘和为 $0$。

      真正乘起来有值的是从 $R[i-\text{fail}(p)]\leftrightarrow a[\text{fail}(p)]$ 项开始的。

      并且后面对应相乘，$r(p)[1]\leftrightarrow a[\text{fail}(p)-1],r(p)[2]\leftrightarrow a[\text{fail(p)-2}]......$

      等等！！！

      这相乘求和不就是第 $p$ 个递推式失效的位置对应的表达式吗？转化一下即 $=a_{\text{fail(p)}}-\Delta (fail(p))$。

      $\Rightarrow \omega \cdot a_{\text{fail(p)}}-\omega \cdot (a_{\text{fail(p)}}-\Delta (\text{fail}(p)\left)\right)=\omega \ast \Delta (\text{fail}(p))=\Delta (i)$。

      发现 $R$ 把 $r(t)$ 失效的差距部分恰好补上了。

      由此可知 $r(t+1)$ 对 $a_i$ 成立。

    • 对于 $m\le k<i$。其实推导与上面差不多。

      $R[1,i-\text{fail}(p)-1]$ 都为 $0$，对应 $a[k-1,k-i+\text{fail}(p)+1]$ 相乘和为 $0$。

      真正乘起来有值的是从 $R[i-\text{fail}(p)]\leftrightarrow a[k-i+\text{fail}(p)]$ 项开始的。

      并且后面对应相乘，$r(p)[1]\leftrightarrow a[k-i+\text{fail}(p)-1],r(p)[2]\leftrightarrow a[k-i+\text{fail(p)-2}]......$

      转化一下即 $=a_{k-i+\text{fail}(p)}-\Delta (k-i+\text{fail}(p))$。

      因为 $k<i$，所以 $k-i+\text{fail}(p)<\text{fail}(p)$。

      还记得 $\text{fail}(p)$ 的定义吗？——是 $r(p)$ 第一个失效的位置。换言之，在此位置之前 $r(p)$ 都是成立的。

      所以 $\Delta (k-i+\text{fail}(p))=0$。

      $\Rightarrow \omega \cdot a_{k-i+\text{fail(p)}}-\omega \cdot (a_{k-i+\text{fail(p)}}-\Delta (k-i+\text{fail}(p)\left)\right)=\omega \ast \Delta (k-i+\text{fail}(p))=0$。

    综上，我们构造的这个 $R$ 本质上只起到了对 $a_i$ 补充 $\Delta (i)$ 的效果，对于其余 $k$ 的贡献都是 $0$。

    这是利用了 $r(p)$ 在 $1\sim \text{fail}(p)-1$ 都满足关系，而在 $\text{fail}(p)$ 相差 $\Delta (p)$ 的性质。

    ---

    此外我们还希望递推式的长度越短越好，也就是说 $\max (m(t),i-\text{fail}(p)+m(p))$ 最短。

    贪心地只需要动态维护最短的 $i-\text{fail}(p)+m(p)$，每次算出 $r(t+1)$ 时都与之前的 $r(p)$ 比一下谁更短即可。

一共递推 $O(n)$ 次，每次修改 $O(n)$ 次，时间复杂度为 $O(n^2)$。

代码实现中，$r$ 用的是 $\text{vector}$，下标从 $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的最小化 。

打表题也可以试着找一下规律，因为只能找线性递推的规律，所以大部分情况要灵活处理，比如说先差分再找规律之类的。