PE 101 Optimum polynomial（拉格朗日插值）

柔情只为你懂 2022-07-14 13:38 197阅读 0赞

## Optimum polynomial ##

### Problem 101 ###

**If we are presented with the first k terms of a sequence it is impossible to say with certainty the value of the next term, as there are infinitely many polynomial functions that can model the sequence.**

**As an example, let us consider the sequence of cube numbers. This is defined by the generating function,  
un = n3: 1, 8, 27, 64, 125, 216, ...**

**Suppose we were only given the first two terms of this sequence. Working on the principle that "simple is best" we should assume a linear relationship and predict the next term to be 15 (common difference 7). Even if we were presented with the first three terms, by the same principle of simplicity, a quadratic relationship should be assumed.**

**We shall define OP(k, n) to be the nth term of the optimum polynomial generating function for the first k terms of a sequence. It should be clear that OP(k, n) will accurately generate the terms of the sequence for n ≤ k, and potentially the *first incorrect term* (FIT) will be OP(k,k\+1); in which case we shall call it a *bad OP* (BOP).**

**As a basis, if we were only given the first term of sequence, it would be most sensible to assume constancy; that is, for n ≥ 2, OP(1, n) =u1.**

**Hence we obtain the following OPs for the cubic sequence:**

**Clearly no BOPs exist for k ≥ 4.**

**By considering the sum of FITs generated by the BOPs (indicated in red above), we obtain 1 + 15 + 58 = 74.**

**Consider the following tenth degree polynomial generating function:**

**un = 1 − n \+ n2 − n3 \+ n4 − n5 \+ n6 − n7 \+ n8 − n9 \+ n10**

**Find the sum of FITs for the BOPs.**

**题意：**

**最优多项式**

**如果我们知道了一个数列的前k项，我们仍无法确定地给出下一项的值，因为有无穷个多项式生成函数都有可能是这个数列的模型。**

**例如，让我们考虑立方数的序列，它可以用如下函数生成，  
un = n3: 1, 8, 27, 64, 125, 216, …**

**如果我们只知道数列的前两项，秉承“简单至上”的原则，我们应当假定这个数列遵循线性关系，并且预测下一项为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) = u1。**

**由此，我们得到了立方数列的最优多项式如下：**

**显然，当k ≥ 4时不存在坏最优多项式。**

**所有坏最优多项式的第一个不正确项（用红色标示的数）之和为1 + 15 + 58 = 74。**

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

**un = 1 − n + n2 − n3 \+ n4 − n5 \+ n6 − n7 \+ n8 − n9 \+ n10**

**求其所有坏最优多项式的第一个不正确项之和。**

**题解：**

**因为我们知道了一个数列的前k项，但我们仍无法确定地给出下一项的值，因为有无穷个多项式生成函数都有可能是这个数列的模型。从这里，我们可以想到拉格朗日插值。什么是拉格朗日插值？我简单得解释一下。**

**我们都做过数学的规律题，就是给你一个数列，然后让你猜下一个数是什么。这题也是这样。**

**先写一下拉格朗日插值的定义：**

**对某个多项式[函数][Link 1]，已知有给定的k \+ 1个取值点：**

**![(x\_\{0\},y\_\{0\}),\\ldots ,(x\_\{k\},y\_\{k\})][x_0_y_0_ldots _x_k_y_k]**

**其中![x\_\{j\}][x_j]对应着[自变量][Link 2]的位置，而![y\_\{j\}][y_j]对应着函数在这个位置的取值。**

**假设任意两个不同的x*j*都互不相同，那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：**

**![L(x):=\\sum \_\{\{j=0\}\}^\{\{k\}\}y\_\{j\}\\ell \_\{j\}(x)][L_x_sum _j_0_k_y_j_ell _j_x]**

**其中每个![\\ell \_\{j\}(x)][ell _j_x]为拉格朗日基本多项式（或称插值基函数），其表达式为：**

**![\\ell \_\{j\}(x):=\\prod \_\{\{i=0,\\,i\\neq j\}\}^\{\{k\}\}\{\\frac  \{x-x\_\{i\}\}\{x\_\{j\}-x\_\{i\}\}\}=\{\\frac  \{(x-x\_\{0\})\}\{(x\_\{j\}-x\_\{0\})\}\}\\cdots \{\\frac  \{(x-x\_\{\{j-1\}\})\}\{(x\_\{j\}-x\_\{\{j-1\}\})\}\}\{\\frac  \{(x-x\_\{\{j+1\}\})\}\{(x\_\{j\}-x\_\{\{j+1\}\})\}\}\\cdots \{\\frac  \{(x-x\_\{\{k\}\})\}\{(x\_\{j\}-x\_\{\{k\}\})\}\}.][ell _j_x_prod _i_0_i_neq j_k_frac _x-x_i_x_j_-x_i_frac _x-x_0_x_j_-x_0_cdots _frac _x-x_j-1_x_j_-x_j-1_frac _x-x_j_1_x_j_-x_j_1_cdots _frac _x-x_k_x_j_-x_k_.]**

**拉格朗日基本多项式![\\ell \_\{j\}(x)][ell _j_x]的特点是在![x\_\{j\}][x_j]上取值为1，在其它的点![x\_\{i\},\\,i\\neq j][x_i_i_neq j]上取值为0。**

**所谓的插值，就是“插”“值”，就是指找出一个通过给出离散数据点的函数。例如，数列中给出数据可以表示为在坐标系上的点，x坐标就是第几项，y坐标就是该项的值。然后画个图，找出这些数在哪个函数数，求出这个函数，这样就可以求出下一项了。先试一个简单的数列：1、8、27…那下一个是什么呢？**

**我们知道:**

**![gif.latex_x_0_3D1][], ![gif.latex_f_28x_0_29_3D1][]**

**![gif.latex_x_1_3D2][], ![gif.latex_f_28x_1_29_3D8][]**

**![gif.latex_x_2_3D3][] , ![gif.latex_f_28x_2_29_3D27][]**

**那么：![20161211192232345][]![20161211192302694][]**

**那么下一项就是将x=4代入上面的二元一次方程，就是58啦。**

**如果你发现这个答案其实是错的！正确的是y=x^3? 其实答案是64。那就是拉格朗日插值是错的咯？**

**再仔细看一下坑爹的高数课本，才发现原来是我们一直搞错了。如果我们给的是n个点，那么拉格朗日给出的函数将会是(n-1)次的。**

**也就是拉格朗日插值:n次多项式可以用n+1个点唯一确定.  
  
这不坑爹吗…用公式之前还得想清楚这个函数是几次的，而且如果是更高次数的还没办法加上点去求。  
  
这就意味着，就算是1、2、3、4、5、6…这样的数列，拉格朗日插值法在耗尽你大量的考试时间去求出通项公式以后，还会给出一个超级坑爹的答案！  
  
那么这个方法还有什么用！  
别急，前面的计算都是为后面做铺垫的。现在才是重点。**

**无论是分布得多么奇怪的点，拉格朗日插值法总能给出一条经过这些点的函数图象。也就是说，就算是1、2、3、4、5、6、1568这样明显不靠谱的答案也是“有规律的”。因为你总可以设一个六次多项式，找出这个数列的通项公式。  
所以说：  
1、3、5、7、9、1598，是对的  
3、1、4、1、5、9、2、999，是对的  
1，1，2，3，5，7，8989，是对的  
2，4，6，8，5，是对的**

**你想怎么写都是对的，即使看起来不靠谱，但其仍然存在内在关系与规律。**

**解释完啦！**

**插值公式是:**

**![20160326175633403][]**

**拉格朗日插值法可以帮助我们解决以下的问题：  
例如：已知x取值0，1，-1，2时， f\{x\}取值2，2，0，6 。  
求x=3时f\{x\}的值。**

#include<cstdio>
 #include<iostream>
 using namespace std;
 int xs[4]={0,1,-1,2};
 int ys[4]={2,2,0,6};
 //f(3)?
 
 int lagrange_Polynomial(int x, int xs[],int ys[]) 
 {
 int c1,c2; 
 int n=4;
 c1=1; 
 for(int i=0;i<n;i++)
 	{ 
 c2=ys[i]; 
 for(int j=0;j<n;j++)
 {
 if(j!=i) c2=c2*((x-xs[j])*1.0)/((xs[i]-xs[j])*1.0); 
 }
 	c1+=c2; 
 } 
 return c1;
 }
 int main()
 {
 	int val=lagrange_Polynomial(3,xs,ys);
 	cout<<val<<endl; //输出 ：20 
 	return 0;
 }

**同理，PE101这题就可以按相同的方法解出来啦。**

**如果你做出了PE101，你可以做另外一题：Codeforces 622F**

[Link 1]: http://zh.wikipedia.org/wiki/%E5%87%BD%E6%95%B0
[x_0_y_0_ldots _x_k_y_k]: /images/20220714/5fc22238c9464c938de1b46257dc20b4.png
[x_j]: http://upload.wikimedia.org/math/7/7/d/77dc7a56686018973d465a321f70b9a5.png
[Link 2]: http://zh.wikipedia.org/w/index.php?title=%E8%87%AA%E5%8F%98%E9%87%8F&action=edit&redlink=1
[y_j]: /images/20220714/d9f1620fe2e9415d805ea5785d0dd715.png
[L_x_sum _j_0_k_y_j_ell _j_x]: /images/20220714/eb8d088b1175407eb39f3573293d18c2.png
[ell _j_x]: /images/20220714/48ae2fd319044ff49b0799c73d34e0a0.png
[ell _j_x_prod _i_0_i_neq j_k_frac _x-x_i_x_j_-x_i_frac _x-x_0_x_j_-x_0_cdots _frac _x-x_j-1_x_j_-x_j-1_frac _x-x_j_1_x_j_-x_j_1_cdots _frac _x-x_k_x_j_-x_k_.]: /images/20220714/344b8cac27fa4d4683cd927839149f42.png
[x_i_i_neq j]: /images/20220714/90a7c5cc82a24adfbcc5fd5c280763e8.png
[gif.latex_x_0_3D1]: http://latex.codecogs.com/gif.latex?x_0%3D1
[gif.latex_f_28x_0_29_3D1]: /images/20220714/b69703c2f99f475e917ea419d0c6a709.png
[gif.latex_x_1_3D2]: http://latex.codecogs.com/gif.latex?x_1%3D2
[gif.latex_f_28x_1_29_3D8]: /images/20220714/7dfe1141f2ee4024af26b04b182c73a9.png
[gif.latex_x_2_3D3]: http://latex.codecogs.com/gif.latex?x_2%3D3
[gif.latex_f_28x_2_29_3D27]: /images/20220714/11a6deca3d68474fbe9a14cda0b994c8.png
[20161211192232345]: https://img-blog.csdn.net/20161211192232345
[20161211192302694]: /images/20220714/2c3b5519a8fb40e58034181d2fe506fa.png
[20160326175633403]: /images/20220714/9eb621202d5344b0ab57108f8b10791f.png