利用简单线性表示优化线性递推

欢迎大家去看叉姐的论文啊~

欢迎大家去看这篇文章的姊妹篇：特征多项式优化线性递推

叉姐利用了特征多项式将线性递推从矩阵快速幂的 $O(k^3 \log n)$ 优化到了 $O(k^2 \log n)$ 。

但这方法对没有学过线性代数的同学极不友好（大雾）。我们尝试利用简单线性表示干同样的事情：将线性递推的复杂度从 $O(k^3 \log n)$ 优化到了 $O(k^2 \log n)$ 。

齐次线性递推

我们先来看一下齐次线性递推的概念。我们定义 $k$ 阶的齐次线性递推应该长这个样子：

对任意 $n \in \mathbb{N}$，有 $a_{n+k} = c_1 a_{n+k-1} + c_2 a_{n+k-2} + \cdots + c_k a_n$ ，其中 $c_1 , c_2 , \cdots , c_k \in \mathbb{N}$ 。

其实我们还希望 $c_k \neq 0$ ,不然这就会退化成一个 $k-1$ 阶的齐次线性递推。但在大多数情况下可以不考虑这种情况，它并不影响最终的时间复杂度。

对了，在本文中，若未特别提及线性递推均为 $k$ 阶齐次线性递推。

线性表示

如果你学过线性相关的话，就知道线性相关和线性表示其实是等价的。

当然，假装我们什么线性代数都没学过的话，线性表示讲述的是一个数和一些数的关系:如果这个数能被这些数加加减减得到，那么就称这个数能被这些数线性表示。

形式化地，对于一个数 $a_n$ 和一些数 $a_0,a_1,\cdots,a_{k-1}$ ，如果存在 $p_0,p_1,\cdots,p_{k-1}$ ，使得 $a_n = p_0a_0 + p_1a_1 + \cdots +p_{k-1}a_{k-1}$ ，则称 $a_n$ 能被 $a_0,a_1,…,a_{k-1}$ 线性表示。

Q: 为什么我们要花这么大力气来求线性表示
A: 因为这样我们就可以根据 $p_0,p_1,\cdots,p_{k-1}$ 线性时间内求出 $a_n$。

每一个 $a_n$ 都可以被线性表示

我们发现，对于每一个 $a_n$ ，他都可以被 $a_0,a_1,\cdots,a_{k-1}$ 线性表示出来。

先举个例子， 对于一个 $3$ 阶的齐次线性递推 $a_{n+3} = a_{n+2} + 2a_{n+1} - a_n $ 。我们有：

$a_0= 1a_0 + 0a_1 + 0a_2$ $a_1= 0a_0 + 1a_1 + 0a_2$ $a_2= 0a_0 + 0a_1 + 1a_2$ 我们尝试继续向下推： $a_3= -1a_0 + 2a_1 + 1a_2$ $a_4= -1a_1 + 2a_2 + 1a_3$ $= -1a_1 + 2a_2 + (-1a_0 + 2a_1 + 1a_2)$ $= -1a_0 + 1a_1 + 3a_2$ $a_5= -1a_2 + 2a_3 + 1a_4$ $= -1a_2 + 2(-1a_0 + 2a_1 + 1a_2) + (-1a_0 + 1a_1 + 3a_2)$ $= -3a_0 + 5a_1 + 4a_2$

显然我们可以一直推下去。换句话说，每一个 $a_n$ 都可以用 $a_0,a_1,a_2$ 的线性组合表示。

可加性

可加性讲的是：如果有

$a_n=p_0a_0+p_1a_1+\cdots+p_{k-1}a_{k-1}$

那么

$a_{n+x}=p_0a_x+p_1a_{x+1}+\cdots+p_{k-1}a_{x+k-1}\ ,\ x\in\mathbb{N}$

例如我们上面得到 $a_5 = -3a_0 + 5a_1 + 4a_2$ ，那么就有

$a_6= -3a_1 + 5a_2 + 4a_3$ $=-3a_1 + 5a_2 + 4(-1a_0 + 2a_1 + 1a_2)$ $=-4a_0+5a_1+9a_2$

可加性挺显然的，证明就留作读者自己思考吧。

这样递推显然方便很多。

优化线性递推

有了这两点显然我们就可以开始优化线性递推了。

在优化线性递推之前，我们先解决一个问题：

如果我们得到 $a_m,a_n$ 的线性表示，我们是否能得到 $a_{m+n}$ 的线性表示？如果能，时间复杂度是多少？

假设我们得到 $a_m=p_0a_0+p_1a_1+p_2a_2+\cdots+p_{k-1}a_{k-1}$ 和 $a_n=q_0a_0+q_1a_1+q_2a_2+\cdots+q_{k-1}a_{k-1}$。

我们先用一次可加性试试：

$a_{m+n}=p_0a_n+p_1a_{n+1}+p_2a_{n+2}+\cdots+p_{k-1}a_{n+k-1}$

我们再用几次可加性试试：

$a_{m+n}=p_0a_n+p_1a_{n+1}+p_2a_{n+2}+\cdots+p_{k-1}a_{n+k-1}$ $=p_0(q_0a_0+q_1a_1+q_2a_2+\cdots+q_{k-1}a_{k-1})$ $+p_1(q_0a_1+q_1a_2+q_2a_3+\cdots+q_{k-1}a_{k})$ $+p_2(q_0a_2+q_1a_3+q_2a_4+\cdots+q_{k-1}a_{k+1})$ $+\cdots+$ $p_{k-1}(q_0a_{k-1}+q_1a_{k}+q_2a_{k+1}+\cdots+q_{k-1}a_{2k-2})$

整理一下:

$a_{m+n}=\sum_{i=0}^{2k-2} \sum_{j=0}^{i} p_jq_{i-j} a_i$
这个式子算下来时间复杂度是 $O(k^2)$ 。

现在我们已经可以用 $a_0,a_1,a_2,\cdots,a_{2k-2}$ 线性表示 $a_{n+m}$ 了。

下一步就简单了，我们只需要知道 $a_{k},a_{k+1},\cdots,a_{2k-2}$ 的线性表示，带入进去就可以了，这部分是可以 $O(k^2)$ 预处理的。

这样算下来，已知 $a_m,a_n$ 求 $a_{m+n}$ 的总复杂度是 $O(k^2)$ 。

于是，我们利用快速幂一样算 $a_n$ ，复杂度为 $O(k^2 \log n)$ 。

补: 再次优化线性递推

相比大家都发现了我们整理下来的式子 $\sum_{i=0}^{2k-2} \sum_{j=0}^{i} p_jq_{i-j} a_i$ 其实是一个标准的卷积的形式。

那么我们根据 $a_m,a_n$ 的线性表示，利用FFT就可以轻松得到上面整理的式子。再利用多项式取模，就可以将 $a_{2k-2},a_{2k-3},\cdots,a_k$ 去掉。
总时间复杂度为 $O(k \log k \log n)$ 。

例：Project Euler 258 A lagged Fibonacci sequence

题意

定义递推数列 $\{g_n\}$ ，满足

• 初值: 对于 $n < 2000$,有 $g_n = 1$，
• 递推式：对于 $n \geq 2000$ ,有 $g_n=g_{n-2000}+g_{n-1999}$

求 $g_{10^{18}} \mod 20092010$。

分析

这有什么好分析的。就是两千阶的线性递推吗。直接上代码得了。

代码

#include<bits/stdc++.h> 

using namespace std;

const long long mod=20092010;
const int k=2000;

void mul(long long *p,long long *q,long long *c) //p,q为a_n,a_m的一种线性表示，c为系数数组，函数算出a_{m+n}的线性表示并存在数组p中
{
    long long tmp[2*k];
    memset(tmp,0,sizeof(tmp));
    for(int i=0;i<k;i++)
        for(int j=0;j<k;j++)
            tmp[i+j]=(tmp[i+j]+p[i]*q[j])%mod; //利用可加性将a_{m+n}利用(a_0,a_1,...,a_{2k-2})线性表示，并把系数存在tmp中
    for(int i=2*k-2;i>=k;i--) //利用a_n=c_0a_{n-1}+c_1a_{n-2}+...+c_ka_{n-k} 将 a_i 转化为 a_{i-1},a_{i-2},...,a_{i-k}
        if(tmp[i])
            for(int j=0;j<k;j++) //这样依次转化 a_{2k-2},a_{2k-3},...,a_k 就可以得到我们想要的线性表示
                if(c[j])
                    tmp[i-j-1]=(tmp[i-j-1]+tmp[i]*c[j])%mod;
    for(int i=0;i<k;i++)
        p[i]=tmp[i];
}

long long solve(long long *a,long long *c,long long n) //a为初始值数组 c为系数数组 和上文线性递推定义对应
{
//  if(k==1) return pw(c[0],n-1)*a[0]%mod;
    if(n<k) return a[n]; //如果已知就直接返回
    long long tmp[k],res[k]; //初始化快速幂，tmp存 a_{2^x} ，res存答案
    memset(tmp,0,sizeof(tmp));memset(res,0,sizeof(res));
    tmp[1]=res[0]=1;
    while(n>0)
    {
        if(n&1) mul(res,tmp,c);
        mul(tmp,tmp,c);
        n>>=1;
    }
    long long ans=0;
    for(int i=0;i<k;i++)
        ans=(ans+a[i]*res[i])%mod;
    return ans;
}

int main()
{
    long long a[k],c[k];
    fill(a,a+k,1);
    memset(c,0,sizeof(c));
    c[1998]=c[1999]=1;
    cout<<solve(a,c,1e18)<<endl;
    return 0;
}