101. 最优多项式(Optimum polynomial)

如果我们知道了一个数列的前\(k\)项,我们仍无法确定地给出下一项的值,因为有无穷个多项式生成函数都有可能是这个数列的模型。例如,让我们考虑立方数的序列,它可以用如下函数生成\(u_n=n^3:1,8,27,64,125,216,\cdots\)

如果我们只知道数列的前两项,秉承“简单至上”的原则,我们应当假定这个数列遵循线性关系,并且预测下一项为\(15\)(公差为\(7\))。如果我们知道了数列的前三项,根据同样的原则,我们也应当首先假定数列遵循二次函数关系。

给定数列的前\(k\)项,定义\(OP(k, n)\)是由最优多项式给出的第\(n\)项的值。显然\(OP(k, n)\)可以精确地给出\(n ≤ k\)的那些项,而可能的第一个不正确项(First Incorrect Term,简记为FIT)将会是\(OP(k, k+1)\)。如果事实的确如此,我们称这个多项式为坏最优多项式(Bad OP,简记为BOP)。

在最基本的情况下,如果我们只得到了数列的第一项,我们应当假定数列为常数,也就是说,对于\(n ≥ 2,OP(1,n) = u_1\)。由此,我们得到了立方数列的最优多项式如下:

$$ \begin{matrix} op\left( 1,n \right) =1& &1,{\color{red} 1},1,1,\cdots\\ op\left( 2,n \right) =7n-6& &1,8,{\color{red} 15},\cdots\\ op\left( 3,n \right) =6n^2-11n+6& &1,8,27,{\color{red} 58},\cdots\\ op\left( 4,n \right) =n^3& &1,8,27,64,125,\cdots\\ \end{matrix} $$
显然,当\(k ≥ 4\)时不存在坏最优多项式。所有坏最优多项式的第一个不正确项(用红色标示的数)之和为\(1 + 15 + 58 = 74\)

考虑下面这个十阶多项式生成函数:

$$ u_n = 1 − n + n^2 − n^3 + n^4 − n^5 + n^6 − n^7 + n^8 − n^9 + n^{10} $$
求其所有坏最优多项式的第一个不正确项之和。

分析:这道题涉及的主要知识点是多项式插值,关于多项式插值,网上有很多很好的介绍文章,我这里就不详细说明了。简单来说,多项式插值就是根据给它的数据点,找到一条次数尽可能低的多项式曲线,使得曲线穿过每一个数据点。如在题目给出的立方数的例子中,给出前两个立方数,也就是两个数据点,尽可能低次数的多项式就是一次多项式,也就是一条直线;同理,给出三个数据点,就可以用二次多项式;一般地,给定\(k\)个数据点,可以\(k-1\)次多项式进行插值,只要数据点两两相异,这个多项式必然存在且唯一。

为了找到这个多项式,可以使用多种插值方法,常用的包括拉格朗日插值、牛顿插值等等。我在这道题目中使用了牛顿插值法,主要是考虑到题目中给出的特殊条件可以把一般的牛顿插值法大大简化,从而使得我们只需通过计算多阶差分的形式就可以求出问题的解。比如对于给它的离散数据点\((x_1,y_1),(x_2,y_2),\cdots(x_k,y_k),\cdots\),我们可以计算它们的前向差分和后向差分,一阶前向差分\(\Delta y_k=y_{k+1}-y_k\),一阶后向差分\(\nabla y_k=y_k-y_{k-1}\);我们还可以计算差分的差分即二阶差分,二阶前向差分\(\Delta ^2y_k=\Delta y_{k+1}-\Delta y_k\),二阶后向差分\(\nabla ^2y_k=\nabla y_k-\nabla y_{k-1}\)。对于题目中给出的立方数,取出它的前四个数,分别计算其一阶、二阶和三阶差分,如下:

img

上图中第一行为立方数,第二行为一阶差分,第三行为二阶差分,第四行为三阶差分,根据差分的定义,这些数都可以很容易的计算出来。通过上图我们可以观察到,如果我们把曲线内的数字相加,如\(8+7=15\)正好等于\(op(2,3)\),而\(27+19+12=58\)正好等于\(op(3,4)\),依次类推,显然通过上图我们可以计算出这些坏多项式的第一个不正确数字之和,只需要我们把四个立方数及其各阶差分加起来就可以了。下面我将简要说明一下这种算法的原理:

一、差分与差商

根据题目的定义,\(op(k,n)\)表示使用数列的前\(k\)项生成的最优多项式给出的第\(n\)项的值,那么第一个不正确项即为\(op(k,k+1)\),如\(op(2,3)=15,op(3,4)=58\)等等。为了求出第一个不正确项,我们需要先求出对应的多项式,再把\((k+1)\)的值代入进去。为了求这个多项式,我们需要用到牛顿插值法。为了介绍牛顿插值法,我们还需要介绍一个差商或者叫均差的概念。差商可以基于差分来计算,如同样对于离散数据点\((x_1,y_1),(x_2,y_2),\cdots(x_k,y_k),\cdots\),它的一阶差商可以记为\([y_k,y_{k-1}]\),则有:

$$ \left[ y_k,y_{k+1} \right] =\frac{y_{k+1}-y_k}{x_{k+1}-x_k}=\frac{\Delta y_k}{x_{k+1}-x_k}=\frac{\nabla y_{k+1}}{x_{k+1}-x_k} $$

对于题目中的问题,其自变量项实际上为自然数列\(\{1,2,3,4,\cdots\}\),显然\((x_{k+1}-x_k)\)恒为一,则上式可以进一步化简为\(\left[ y_k,y_{k+1} \right] =\Delta y_k=\nabla y_{k+1}\),也即在这种特殊情况下,一阶差分和一阶差商是相同的。与二阶差分类似,二阶差商定义为差商的差商,则有:

$$ \left[ y_k,y_{k+1},y_{k+2} \right] =\frac{\left[ y_{k+1},y_{k+2} \right] -\left[ y_k,y_{k+1} \right]}{x_{k+2}-x_k}=\frac{\Delta y_{k+1}-\Delta y_k}{2}=\frac{\Delta ^2y_k}{2} $$

依次类推,\(n\)阶差商为:

$$ \left[ y_k,y_{k+1},\cdots ,y_{k+n} \right] =\frac{\Delta ^ny_{k+1}-\Delta ^ny_k}{n!}=\frac{\Delta ^ny_k}{n!} $$

也即\(y_k\)处的\(n\)阶差商等于其对应的\(n\)阶差分除以\(n!\),如对于立方数的前四项,我们可以计算其前三阶差商如下:

img

对比一下这里的差商和之前的差分,可以发现一阶差分和一阶差商是一样的,二阶差商等于二阶差分除以\(2!\)也就是二,三阶差商等于三阶差分除以\(3!\)也就是六,和我们公式推导的一致。

二、牛顿插值法

在计算出数列的各阶差商后,我们就可以用牛顿插值法来计算其最优多项式,对于题目中自变量为为自然数列的情况,使用数列的前\(k\)项,其对应的最优多项式为:

$$ op\left( k,x \right) =\left[ y_1 \right] +\left[ y_1,y_2 \right] \left( x-1 \right) +\left[ y_1,y_2,y_3 \right] \left( x-1 \right) \left( x-2 \right) \cdots +\left[ y_1,\cdots ,y_k \right] \left( x-1 \right) \cdots \left( x-k+1 \right) $$

如我们取出前三个数据点\(\left( 1,1 \right) ,\left( 2,8 \right) ,\left( 3,27 \right)\),则把相应的差商代入,有:

$$ op\left( 3,x \right) =1+7\left( x-1 \right) +6\left( x-1 \right) \left( x-2 \right) =6x^2-11x+6 $$

与题目中给出的二次多项式相符。除开上面这种形式以外,我们还可以把数据点的顺序颠倒过来,可以形成一样的多项式:

$$ \begin{align*} op\left( k,x \right) =\left[ y_k \right] &+\left[ y_{k-1},y_k \right] \left( x-k \right) \\&+\left[ y_{k-2},y_{k-1},y_k \right] \left( x-k \right) \left( x-k+1 \right)+ \cdots \\&+\left[ y_1,y_2,\cdots ,y_k \right] \left( x-k \right) \cdots \left( x-2 \right) \end{align*} $$

如同样取出前三项,最优多项式也可以写成:

$$ op\left( 3,x \right) =27+19\left( x-3 \right) +6\left( x-3 \right) \left( x-2 \right) =6x^2-11x+6 $$

根据上面推导的差分与差商的关系,最优多项式可进一步表示为:

$$ op\left( k,x \right) =\left[ y_k \right] +\Delta y_{k-1}\left( x-k \right) +\frac{\Delta ^2y_{k-2}}{2!}\left( x-k \right) \left( x-k+1 \right) \cdots +\frac{\Delta ^{k-1}y_1}{\left( k-1 \right) !}\left( x-k \right) \cdots \left( x-2 \right) $$

上式可以进一步表示成为:

$$ op\left( k,x \right) =\sum_{i=1}^k{\frac{\Delta ^{i-1}y_{k-i+1}}{\left( i-1 \right) !}}\prod_{j=0}^{i-2}{\left( x-k+j \right)} $$

题目要求最优多项式的第一个不正确项,即求\(op(k,k+1)\),则有:

$$ op\left( k,k+1 \right) =\sum_{i=1}^k{\frac{\Delta ^{i-1}y_{k-i+1}}{\left( i-1 \right) !}}\prod_{j=0}^{i-2}{\left( j+1 \right)}=\sum_{i=1}^k{\frac{\Delta ^{i-1}y_{k-i+1}}{\left( i-1 \right) !}}\left( i-1 \right) !=\sum_{i=1}^k{\Delta ^{i-1}y_{k-i+1}} $$

可以发现,上式中的阶乘项可以相互约掉,从而只保留了差分项,如我要求\(op(3,4)\),有:

$$ op\left( 3,4 \right) =\sum_{i=1}^3{\Delta ^{i-1}y_{k-i+1}}=\Delta ^0y_3+\Delta ^1y_2+\Delta ^2y_1=27+19+12=58 $$

可以看到,这里显示的三项正是我们在计算差分的那个图形显示的三项,就此也证明了我们算法的正确性。

三、实现

根据以上原理,我们只需要计算出题目给出的十次多项式的前十项,然后计算这十项的逐阶差分,直到只有一个差分项为止,然后我们把这些差分项逐一加总,即为题目所求。需要注意的,题目中给出的十次多项式是一个首项为\(1\),公比为\(-n\)的对比数列,则我们可以用等比数列公式进行求和,有:

$$ 1 − n + n^2 − n^3 + n^4 − n^5 + n^6 − n^7 + n^8 − n^9 + n^{10}=\frac{n^{11}+1}{n+1} $$

如此,可以进一步简化计算。代码如下:

# time cost = 22.6 µs ± 87 ns

def main():
    res = 0
    tdp = lambda x: (x**11+1)//(x+1)
    arr = [tdp(x) for x in range(1,11)]
    diff = lambda arr:[x-y for x,y in zip(arr[1:],arr[:-1])]
    for _ in range(10):
        res += sum(arr)
        arr = diff(arr)
    return res