/ 数列

数列求解的一些思路

组合数学里面数列是一个非常重要的概念。事实上,如果掌握得熟练的话,一类数列递推求解对高中生来说应该是十分easy的。

比如说,

an=3an1+2a_n=3a_{n-1}+2

常见做法有两种:an3n=an13n1+23n\frac{a_n}{3^n}=\frac{a_{n-1}}{3^{n-1}}+\frac2{3^n}和特征根法an1=3(an11)a_n-1=3(a_{n-1}-1)。这两种做法都十分地妙,但是老师从来没有讲过为什么会这么做,以至于我现在碰到这样的题还只是知其然不知其所以然。

本文将讨论一些更普适的做法以及程序的实现,并尽可能地提供这些做法的motivation。

生成函数法

为什么上来就开讲生成函数,因为生成函数实际上是另一种描述递推式的方式。对数列的操作可以十分简洁地用生成函数这种记号来表示,而多项式又有很多近乎甚至某些地方超越实数的美好性质,在计算的时候通常可以将我们肉眼无法发现的一些联系用无脑推导得出。

举个刚才那个例子,定义{an}\{a_n\}的生成函数为A(x)=iaixiA(x)=\sum_{i}a_ix^i,且a0=2a_0=2,则A(x)=3xA(x)+21xA(x)=3xA(x)+\frac2{1-x},解得

A(x)=2(1x)(13x)=2(ixi)(i(3x)i)\begin{aligned} &A(x)\\ =&\frac2{(1-x)(1-3x)}\\ =&2\left(\sum_ix^i\right)\left(\sum_i(3x)^i\right)\\ \end{aligned}

容易看出数列的第nn项就为i=0n3i\sum_{i=0}^n3^i

又如卡特兰数,它的递推式为Cn=i=0n1CiCn1iC_n=\sum_{i=0}^{n-1}C_iC_{n-1-i}。如何求解它的通项?有一个绝妙的组合解释能说明Cn=1n+1(2nn)C_n=\frac1{n+1}\binom{2n}{n},但是作为NOIP级别的选手,我表示只会暴力。

定义C(x)C(x)为卡特兰数的生成函数。观察卷积形式的递推式,可以得到方程C(x)=xC2(x)+1C(x)=xC^2(x)+1。解一下,就得到了C(x)=1±14x2xC(x)=\frac{1\pm\sqrt{1-4x}}{2x}。这两个解只有一个是正确的,代入x=0x=0,发现应该取负号。

接下来我要不加解释地给出一系列推导。

C(x)=114x2x=1i=0+(12i)(4x)i2x=i=1+(4x)i(j=0i1(12j))/i!2x=i=1+(2x)i(j=0i1(2j1))/i!2x=i=1+(2x)i(j=1i1(2j1))/i!2x=i=1+(2x)i(2(i1))!/(j=1i12j)/i!2x=i=1+2xi(2(i1))!/((i1)!i!)2x=i=1+2xii(2(i1)i1)2x=i=0+xii+1(2ii)\begin{matrix} C(x)&=&\frac{1-\sqrt{1-4x}}{2x}\\ &=&\frac{1-\sum_{i=0}^{+\infty}\binom{\frac12}{i}(-4x)^i}{2x}\\ &=&\frac{-\sum_{i=1}^{+\infty}(-4x)^i\left(\prod_{j=0}^{i-1}(\frac12-j)\right)/{i!}}{2x}\\ &=&\frac{-\sum_{i=1}^{+\infty}(2x)^i\left(\prod_{j=0}^{i-1}(2j-1)\right)/{i!}}{2x}\\ &=&\frac{\sum_{i=1}^{+\infty}(2x)^i\left(\prod_{j=1}^{i-1}(2j-1)\right)/{i!}}{2x}\\ &=&\frac{\sum_{i=1}^{+\infty}(2x)^i(2(i-1))!/(\prod_{j=1}^{i-1}2j)/{i!}}{2x}\\ &=&\frac{\sum_{i=1}^{+\infty}2x^i(2(i-1))!/((i-1)!i!)}{2x}\\ &=&\frac{\sum_{i=1}^{+\infty}\frac{2x^i}i\binom{2(i-1)}{i-1}}{2x}\\ &=&\sum_{i=0}^{+\infty}\frac{x^i}{i+1}\binom{2i}{i} \end{matrix}

就得到了相同的结论。虽然用牛顿二项式定理展开略显麻烦,但有一道天津OI的结论题就是这么做的,大部分人打表找规律就过了……

由于生成函数的系数的求解在算法界有大名鼎鼎的FFT作为靠山,很多时候直接从生成函数得到数列各项的值是可取的,尤其当我们只需要106\leq10^6项的时候。比如说,伯努利数的预处理就常常用xex1\frac{x}{e^x-1}这个生成函数解决。

在15年金策的论文里面有很多方法帮助我们求套着初等函数的生成函数的各项系数,甚至有ln(1F(x))\ln(1-F(x))eF(x)e^{F(x)}的求法;如果一个生成函数能表示成一个方程的解,可以用牛顿迭代法几乎不变复杂度地求解;求解复合函数,可以用多项式快速点值,然后IDFT插值;对于复合函数求复合逆,更是有拉格朗日反演。

这样一来,许多能用简单的生成函数形式表示的数列都可以得到比较小的项了。

较大的项怎么办?可能会达到101810^{18}级别。

线性性

因为项比较大,我们能求解的问题的范围也必然会相对应地缩小,因为项较小的问题可以规约成项较大的问题。所以我们要对问题做出一定的限制,比如说数列存在递推式,比如说递推式是线性的,比如说递推式是常系数的,比如说递推式是齐次的。

随手找一个满足以上所有限制的例子:an=an1+an2a_n=a_{n-1}+a_{n-2}。没错,这就是斐波那契数列。

我们求斐波那契数列第n(n107)n(n\leq10^7)项是很快的,但是对于nn更大的时候这条递归式的方法就不太行得通了。

在这里再不加解释地插入另一个看似无关也实际无关的问题。这个问题的引入要感谢《组合数学》,它的这一节具有非常启发式的思维,不愧是给人看的书。

求解一个十分基本的一阶微分方程:y=yy=y'。即找出一个求导等于自身的函数。

我还记得一个很冷很冷的笑话:天堂里,一个精神病人学了微积分后,老是对别人说:“我微分你,我微分你。”他认为,只要微分足够多次,就可以把别人削没。直到有一天牛顿来了,牛顿默默地对他说了一句:“我是exe^x。”

的确,exe^x发明出来仿佛就是为了能让它求导等于自己。从泰勒展开的过程来看,exe^x系函数似乎也是唯一解了。

另外,y=0y=0也是一个解。这个显然。

那答案是两个吗?非也。还有y=kexy=ke^x。由于求导的线性性,只要微分方程齐次,F(x)F(x)满足条件,kF(x)kF(x)也必然满足条件。但是微分方程的求解通常是伴随着初始值的,所以要根据初始值来确定最终kk的值。

这似乎给我们了一个思路。我们先找出几个满足这个最初的方程的解,我们尝试把几个解乘上系数相加,得到一个新的解,因为方程是线性的,而新的解又是线性组合出来的,所以这个新的解必然满足方程。但是我们可以通过调整这几个系数,来使解满足初值。

举另外一个例子。解二阶微分方程x=kxmx''=\frac{-kx}m。这其实是我自己推导简谐运动方程,把小球挂在一个劲度系数为kk的弹簧上,小球质量为mm

这回在解上带了一个系数,我们已经不能简单地用exe^x了。但是我们可以做若干合理的猜测,猜测它在eqte^{qt}这个空间内有解。

代入x=eqtx=e^{qt},得q2=kmq^2=\frac{-k}mk,m>0k,m>0,所以无解?怕什么虚数!q=±kmiq=\pm\sqrt{\frac km}i。代回去,就可以得到解的形式为:aekmit+bekmitae^{\sqrt{\frac km}it}+be^{-\sqrt{\frac km}it}。虽然结果不可信,但不妨代入一个初始值试试。假设初始位置在x=0x=0处,初始速度x=1x'=1,加速度由于方程已经限制过了,这里再限制没有用。得方程组:

{aekmi0+bekmi0=0kmaekmi0bkmekmi0=1\left\{\begin{matrix} ae^{\sqrt{\frac km}i\cdot0}+be^{-\sqrt{\frac km}i\cdot0}=0\\ \sqrt{\frac km}ae^{\sqrt{\frac km}i\cdot0}-b\sqrt{\frac km}e^{-\sqrt{\frac km}i\cdot0}=1 \end{matrix}\right.

解得a=12mk,b=12mka=\frac12\sqrt{\frac mk},b=-\frac12\sqrt{\frac mk}。故最终x=12mk(ekmitekmit)=mksin(kmt)x=\frac12\sqrt{\frac mk}\left(e^{\sqrt{\frac km}it}-e^{-\sqrt{\frac km}it}\right)=\sqrt{\frac mk}\sin\left(\sqrt{\frac km}t\right)。这正说明了简谐运动的图像是正弦型函数!而且解唯一,我们又得到了一个解,所以这就是那个解!

继续讲下去就扯远了。上面的一大段内容,只是想说明我们可以从一些特殊的但比较容易得到的解出发,再根据线性性得到更多的解乃至通解。

回归到斐波那契数列上。斐波那契数列的增长非常快,它事实上代表了辗转相除法的最坏的情况。我们可以合理地推测,kqnkq^n形式是斐波那契数列递推式的一种解。

kqnkq^n代入递推式,kqn=kqn1+kqn2kq^n=kq^{n-1}+kq^{n-2} \Rightarrow q2=q+1q^2=q+1 \Rightarrow q=1±52q=\frac{1\pm\sqrt5}2,真的是有解的!但是这个解不是整数,可信吗?

在刚才解二阶微分的时候,我们无意引入了虚数,但是最后虚数又奇迹般的消失了。现在,如果我们做同样的事情,也能得到一样的结果。虽然底不是整数,但是只要初始值是整数,二项式展开之后发现一定还是整数。

多项式

对于一个常系数齐次线性递推an=i=1kciania_n=\sum_{i=1}^kc_ia_{n-i},对应的方程就是qni=1kciqni=0q^n-\sum_{i=1}^kc_iq^{n-i}=0。如果我们想求一个常系数齐次线性递推数列的第nn项,有两步:

  1. 求出根q1,q2,...,qkq_1,q_2,...,q_k
  2. 根据初始项求出每个根应该乘的系数bib_i,则答案为i=1kbiqin\sum_{i=1}^kb_iq_i^n

这看起来非常美好。快速幂一上问题就解决了。然而,kk次方程的kk个根我们并不一定能解出来。事实上,有人证明了大于等于五次的方程是没有通解的,这限制了我们的能力。真的就没有办法了吗?

解析几何里有着很妙的一招,叫“设而不求”。假装得到了这个根,然后用它进行运算。在碰到这个根满足的方程的时候,我们就能进行等价替换。

在这里我们也先搁着那些根q1,q2,...,qkq_1,q_2,...,q_k。我们想求qinq_i^n。可是我们连求qiq_i都做不到,可能求出qinq_i^n吗?不要忘了,我们是有一些初始条件的。a1=i=1kbiqia_1=\sum_{i=1}^kb_iq_ia2=i=1kbiqi2a_2=\sum_{i=1}^kb_iq_i^2……我们也不是直接要求qinq_i^n,而是他们线性组合的结果,这些组合的方式和系数和前kk项一模一样。这是否能为我们提供方便?

考虑一下另外一个问题。一个平面上的点,可以用直角坐标描述,也可以用极坐标描述,但无论怎么样,我们都需要两个数;空间中的点,可以用空间直角坐标描述,也可以用球坐标描述,还可以用柱坐标描述,但无论怎么样,我们都需要三个数。一个必须用三个数确定的东西,我们至少得用三个数才能还原出那个信息。类似的例子还有点值和多项式的关系,一个nn次多项式,至少要取nn个点值就才能描述。

但是这个关系是双向的。我们至少需要nn个数来描述一个nn元信息,也就意味着我们至多需要nn个数来描述它。如果在二次函数上面取4个点,有一个点是无用的;试着对这四个点做拉格朗日插值,得到的,还是一个二次函数。

在这个递推数列中,每一项只依赖于前kk项。确定下前kk项,我们就确定下了所有的项。这是一个kk元信息的体系,当我们只需要kk个信息来描述这一切的时候,为什么会需要求qinq_i^n这么大一个东西?

于是,接下来的思考,似乎就是如何简化信息到之和前kk项有关了。我们知道一个方程qki=1kciqki=0q^k-\sum_{i=1}^kc_iq^{k-i}=0(这里及以后用qq代表对于每一个qiq_i),我们要利用这个信息。qnq^n可以减掉若干个qki=1kciqkiq^k-\sum_{i=1}^kc_iq^{k-i},它的值并不会发生变化。如果我们把qnq^n减掉一些之后能够让qq的次数降下来,那就可以继续利用已有的前kk项信息进行计算。

稍微了解多了些多项式的人都知道,这实际上是在取模。

现在我们得到了qnmod(qki=1kciqki)q^n\mathrm{mod}(q^k-\sum_{i=1}^kc_iq^{k-i}),然后怎么办?

i=0k1biqin\sum_{i=0}^{k-1}b_iq_i^n

=i=0k1bij=0k1djqij=\sum_{i=0}^{k-1}b_i\sum_{j=0}^{k-1}d_jq_i^j

=j=0k1dji=0k1biqij=\sum_{j=0}^{k-1}d_j\sum_{i=0}^{k-1}b_iq_i^j

=j=0k1djaj=\sum_{j=0}^{k-1}d_ja_j

So easy!

矩阵

因为这东西是线性的,线性的,线性的,所以一个常见的做法就是矩阵优化转移。

譬如说,对于斐波那契数列,

[1110][an1an2]=[anan1]\left[\begin{matrix}1&1\\1&0\end{matrix}\right]\left[\begin{matrix}a_{n-1}\\a_{n-2}\end{matrix}\right]=\left[\begin{matrix}a_n\\a_{n-1}\end{matrix}\right]

利用矩阵乘法的结合性,我们可以对转移矩阵做一个快速幂,最后再乘上初始向量。

利用矩阵乘法的结合性,能让我们本来要做很多次转移的一些DP,还有常系数线性递推,省了一些转移。但是有时先对转移矩阵做乘法的代价高到不可忽略。

这个时候的有一些可用的策略:

  • 放弃套路,另找策略,从复杂度较低的乘法处入手。譬如说,放弃矩阵带来的优势,改用暴力转移;或者先用矩阵预处理一部分,最后再暴力转移。
  • 可以想办法它的一些行为来描述该矩阵。Codeforces有一道叫Matrix God的题,让你验证两个矩阵相乘是否等于另一个。正解是,随机几个向量,用两边去乘一下看是否相等。
  • 使矩阵更加可乘。

我们来讨论一下第三种。

对角化

这里讲的不是Cantor的对角线法。

什么矩阵乘起来比较快?

对角矩阵。

对角矩阵乘普通矩阵O(n2)\mathrm O(n^2),对角矩阵乘对角矩阵O(n)\mathrm O(n)。如果要算对角矩阵的快速幂,直接每个位置原地快速幂就好了。

这么方便的性质,当然不是每个矩阵都能拥有。对于一个更加普通的矩阵AA,我们想得到它的幂,而且要快,怎么办?

引入对角化:找一个普通矩阵PP,一个对角矩阵Λ\Lambda,让P1ΛP=AP^{-1}\Lambda P=A。虽然在求幂的时候还是有一次的普通矩阵乘普通矩阵,但是绝大多数运算已经被简化了!

特征值和特征向量

现在来找PPΛ\Lambda吧。变一下形,AP=PΛAP=P\Lambda。由于Λ\Lambda是一个对角矩阵,PΛP\Lambda就是把PP的第一列乘上Λ11\Lambda_{11},第二行乘上Λ22\Lambda_{22}……由此看出,AAPP的变换也只能这么做。

一列一列考虑,对于第jj列,APj=ΛjjPjAP_{* j}=\Lambda_{jj}P_{* j}。因此只要找到一堆列向量VV和一个常数λ\lambda满足AV=λVAV=\lambda V就可以了。我们还不能只找一个,我们总共要nn个线性无关的解,否则PP不可逆。

为了找这个东西,搞出了一门学问。这个λ\lambda就叫AA的一个特征值,VV就是AAλ\lambda下的一个特征向量。

AV=λVAV=\lambda V

(λIA)V=0(\lambda I-A)V=0

我们想要的VV肯定不是零向量。所以λI−A=0。(感谢Owen指正,这是显然不对也不必要的。)

通常情况下,满秩的矩阵性质更有得讨论。但是在这里,我们却希望λIA=0\left|\lambda I-A\right|=0

为什么呢?如果我们把VV看成[x1x2xn]\left[\begin{matrix}x_1\\x_2\\\vdots\\x_n\end{matrix}\right],这就变成了一个常数项全部是00nn元一次方程组。而λIA0\left|\lambda I-A\right|\neq0即满秩的话,所有的变量都不得不取00,这是我们不希望看到的。于是,令其不再满秩,方程组也有了无数个非零解。

但因此我们收获了一个方程,这使得λ\lambda能解了。λIA=0\left|\lambda I-A\right|=0就叫矩阵AA的特征方程,λIA\left|\lambda I-A\right|AA的特征多项式。这些“特征”在不同的领域有不同的解释、含义和应用,却能用统一的语言来描述,这大概就是数学之美吧。

不妨做个验证。我们的确可以举出各种各样的例子,但我想举一个特殊的例子,来说明特征方程和对角化是多少有些关系的。

当我们得到AA对角化之后的一个矩阵Λ\Lambda,如果我们想把它再对角化,可以得到从AA出发得不到的结果吗?肯定是不能。但是从特征多项式出发,我们能得到相同的结论吗?如果Λ\Lambda不能得到更多的解,是否意味着Λ\Lambda的特征方程和AA是一样的?

而的确也是一样的。定义矩阵AA与矩阵BB相似当且仅当存在矩阵PP使得A=P1BPA=P^{-1}BP。矩阵相似关系是一种等价关系,即满足自反、对称性和传递性。为什么把两个看起来完全不一样的矩阵叫做“相似矩阵”?然而它们只是形式上看起来不一样,它们的很多性质是一样的。如果知道AB=AB\left|A\right|\cdot\left|B\right|=\left|AB\right|,就很容易证明出它们有相同的行列式。更进一步,它们拥有相同的特征多项式和特征方程。这就说明了,特征方程的确是矩阵的特征,矩阵对角化的所有解,都蕴含在这个特征中。

这样我们的很多矩阵都可对角化了。但是这和上面递推方程一样绕不开一个问题,那就是解高次方程。

Cayley-Hamiton定理

我们在上面面对一个高次方程的时候,并没有直接求出高次方程的解,而是根据解的性质解决问题。这能帮助我们绕开解方程这个难题。

但是特征方程的λ\lambda是什么意思?这东西除了能作为矩阵对角化后的对角线,又对解决线性数列递推有什么帮助?

Cayley-Hamilton定理:λIAλ=A=0\left|\lambda I-A\right|_ {\lambda=A}=0

如果我们想求一个矩阵AAnn次幂,我们就可以先求出λnmodλIA\lambda^n\mathrm{mod}\left|\lambda I-A\right|,然后再把AA代入。

如果试着把常系数齐次线性递推的转移矩阵手动展开出特征多项式,会发现得到的多项式和上面的方法得到的是完全一样的。那为什么要这样一个方法呢?因为这比上面的方法更强大。它可以求解非齐次常系数线性递推,以及优化理论上所有矩阵转移,优化程度因矩阵而异。

如果我们真的无法手算特征多项式,那还可以用求行列式的单位根点值然后IDFT出来。反正只要知道一定的初始值,就能推出来了。不过绝望到用O(n4)\mathrm O(n^4)的复杂度的时候……要不就是可以用暴力转移,要不就要写高精度了。

速证

介绍一些线性代数中的基本知识:

  • aija_{ij}的余子式MijM_{ij}:把第ii行和第jj列删去,剩下的四个小矩阵再拼起来形成的一个小一号的矩阵;
  • 拉普拉斯展开:求解矩阵行列式的一种方法,对某一列或某一行进行,如对第ii行:A=j(1)i+jaijMij\left|A\right|=\sum_ j(-1)^{i+j}a_ {ij}\left|M_ {ij}\right|
  • 伴随矩阵$A^* 满足a^* _ {ji}=(-1)^{i+j}\left|M_ {ij}\right|$:注意有转置!
  • AA=AA=AIA^* A=AA^* =\left|A\right|I(待证iji\neq j的情况);

AAnn阶方阵,设λIA\lambda I-A的伴随矩阵为BB

因为λIA\lambda I-A是关于λ\lambda的不超过nn次的多项式,所以BB的每一项都是一个不超过n1n-1次的多项式;故设B=i=0n1λiBiB=\sum_ {i=0}^{n-1}\lambda^iB_ i

由伴随矩阵性质得(λIA)B=λIAI=(i=0naiλi)I(\lambda I-A)B=\left|\lambda I-A\right|I=\left(\sum_ {i=0}^na_ i\lambda^i\right)I

又有

(λIA)B=i=0n1(λi+1BiλiABi)(\lambda I-A)B=\sum_ {i=0}^{n-1}(\lambda^{i+1}B_ i-\lambda^iAB_ i)

=λnBn1+i=1n1λi(Bi1ABi)AB0=\lambda^nB_{n-1}+\sum_{i=1}^{n-1}\lambda^i\left(B_ {i-1}-AB_i\right)-AB_0

对比系数会得到一个方程组aiI=Bi1ABia_ iI=B_ {i-1}-AB_ i,发现后面那两个东西形式看起来非常可裂项!

aiAi=AiBi1Ai+1Bia_ iA^i=A^iB_ {i-1}-A^{i+1}B_ i

把所有的方程加起来,原定理得证:

i=0naiAi=0\sum_ {i=0}^na_ iA^i=0

小结

本文从生成函数入手,讨论了若干种数列的算法,并引入了若干线性代数的概念,并利用线性代数里面我不会证明的一个定理,尝试解决一类数列递推问题。

因此与其称作“数列求解思路”,不如把这篇文章看为我线性代数的入门文章。