Index
拉格朗日插值
说在前面
写FFT的时候发现要先置拉格朗日插值和卷积,写拉格朗日的时候发现要先写高斯消元?!
呼……
慢慢来吧。
时间真的不多了呢。
拉格朗日插值是拿来干什么的
插值,指的就是将一个多项式从点值表达变换到系数表达。
Q:多项式是什么?
A:???????
Q:这两个表达法是什么?
A
点值表达:一个n次多项式里有n+1个常数,如果我们知道这个多项式里n+1个点的坐标(代入x求出y),我们就可以建立一个n+1元1次方程组,解出n+1个常数。
系数表达:就是我们熟悉的那种
y
=
f
(
x
)
=
a
n
x
n
+
a
n
−
1
x
n
−
1
+
⋯
+
a
1
x
1
+
a
0
x
0
y=f(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x^1+a_0x^0
y=f(x)=anxn+an−1xn−1+⋯+a1x1+a0x0。
Q:什么叫变换?
A:我们要通过已知的n个点
(
x
i
,
y
i
)
(x_i,y_i)
(xi,yi)求出一个n-1次的多项式
f
(
x
)
f(x)
f(x)使之满足
∀
i
∈
[
1
,
n
]
&
i
∈
N
+
,
有
f
(
x
i
)
=
y
i
\forall ~ i \in [1,n] \And i ~\in N_+,有f(x_i)=y_i
∀ i∈[1,n]&i ∈N+,有f(xi)=yi。
当然,这n个x互不相同,如果有x相同而这两者y相同,那么就是浪费名额;如果有x相同而这两者y不同,那么这个东西就不满足函数的定义了,怎么求函数or多项式?
生活经验告诉我们,用两点可以确定一条线,三点确定一条抛物线或者一个圆,就是靠这些点的坐标建立若干个方程并解出多项式的各个系数。
Q: 为什么用n个点求出的是n-1次的多项式?
A: n-1次的多项式 f ( x ) = a 0 x 0 + a 1 x 1 + a 2 x 2 + ⋯ + a n − 1 x n − 1 f(x)=a_0x^0+a_1x^1+a_2x^2+\cdots+a_{n-1}x^{n-1} f(x)=a0x0+a1x1+a2x2+⋯+an−1xn−1,有 a 0 , a 1 , a 2 , ⋯ , a n − 1 {a_0,a_1,a_2,\cdots,a_{n-1}} a0,a1,a2,⋯,an−1共n个系数,根据常识(?),求解n个未知数需要n个方程,于是我们需要n个点的坐标。
那么怎么操作呢?
当然可以用高斯消元
Θ
(
n
3
)
\Theta(n^3)
Θ(n3)正面解决啦。
但是呢,贪心的数学家们找出了更快的方案,最后取名叫拉格朗日插值。
……什么?牛顿插值? 那是什么?不清楚!
格雷戈里:那我呢?!
拉格朗日插值的流程
我们当然可以写一个大大的分段函数:
f
(
x
)
=
{
y
1
i
f
x
=
x
1
y
2
i
f
x
=
x
2
y
3
i
f
x
=
x
3
⋮
y
n
i
f
x
=
x
n
f(x)= \begin{cases} y_1\quad if ~ x=x_1\\ y_2\quad if ~ x=x_2\\ y_3\quad if ~ x=x_3\\ \vdots \\ y_n\quad if ~ x=x_n \end{cases}
f(x)=⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧y1if x=x1y2if x=x2y3if x=x3⋮ynif x=xn
问题解决!
这哪里解决了!
这个函数①很耍赖并且②只适用于这n个点。
①:什么是耍赖呢?
我们希望得到的多项式肯定不是一个跟表格没差的巨型分段函数啊喂,至少好歹只用一个函数吧!
于是我们试着把这个巨型分段函数压扁压成一个函数。那么这个函数至少需要满足当x取xi时,f(xi)=yi,于是想到引入一个bool量,得到:
f
(
x
)
=
∑
1
≤
i
≤
n
[
x
=
x
i
]
y
i
f(x)=\sum_{1~\le ~ i ~ \le ~ n}[x=x_i]y_i
f(x)=1 ≤ i ≤ n∑[x=xi]yi
只有当x等于某个特定的xi时,bool量取1,让yi进入和中,其余时候取0,那么得到的和j就会等于yi。
好像大概解决了。
可是数学家们叫起来了,这不是正经的数学公式!数学公式哪来的bool!
好吧,我们现在来设法换掉这个bool值。
我们重新定义上文的
f
(
x
)
f(x)
f(x):
令
f
(
x
)
=
∑
1
≤
i
≤
n
l
i
(
x
)
⋅
y
i
,
其
中
,
l
i
(
x
)
=
∏
1
≤
j
≤
n
&
&
i
≠
j
x
−
x
j
x
i
−
x
j
=
x
−
x
1
x
i
−
x
1
×
x
−
x
2
x
i
−
x
2
×
⋯
x
−
x
i
−
1
x
i
−
x
i
−
1
×
x
−
x
i
+
1
x
i
−
x
i
+
1
×
⋯
×
x
−
x
n
x
i
−
x
n
\begin{aligned} &令f(x)=\sum_{1\le i\le n} l_i(x)\cdot y_i ,\\ &其中,l_i(x)=\prod_{1\le j \le n ~ \And\And ~ i ~\neq ~ j} \frac{x-x_j}{x_i-x_j} \\ &=\frac{x-x_1}{x_i-x_1}\times \frac{x-x_2}{x_i-x_2}\times \cdots \frac{x-x_{i-1}}{x_i-x_{i-1}}\times\frac{x-x_{i+1}}{x_i-x_{i+1}}\times\cdots\times\frac{x-x_n}{x_i-x_n} \end{aligned}
令f(x)=1≤i≤n∑li(x)⋅yi,其中,li(x)=1≤j≤n && i = j∏xi−xjx−xj=xi−x1x−x1×xi−x2x−x2×⋯xi−xi−1x−xi−1×xi−xi+1x−xi+1×⋯×xi−xnx−xn
人们巧妙地构造出这个叫做拉格朗日基本多项式,或者叫做插值基函数的 l i ( x ) l_i(x) li(x)来换掉了这个bool值。
我们来看看这个新的函数,假设我们现在求
f
(
x
1
)
f(x_1)
f(x1),结果就应该是
y
1
y_1
y1,为了方便,令
n
=
3
n=3
n=3,那么就有:
f
(
x
1
)
=
l
1
(
x
1
)
⋅
y
1
+
l
2
(
x
1
)
⋅
y
2
+
l
3
(
x
1
)
⋅
y
3
=
x
1
−
x
2
x
1
−
x
2
×
x
1
−
x
3
x
1
−
x
3
×
y
1
+
x
1
−
x
1
x
2
−
x
1
×
x
1
−
x
3
x
2
−
x
3
×
y
2
+
x
1
−
x
1
x
3
−
x
1
×
x
1
−
x
2
x
3
−
x
2
×
y
3
=
y
1
\begin{aligned} f(x_1)&=l_1(x_1)\cdot y_1+l_2(x_1)\cdot y_2+l_3(x_1)\cdot y_3 \\ &=\frac{x_1-x_2}{x_1-x_2}\times \frac{x_1-x_3}{x_1-x_3}\times y_1+\frac{x_1-x_1}{x_2-x_1}\times \frac{x_1-x_3}{x_2-x_3}\times y_2+\frac{x_1-x_1}{x_3-x_1}\times \frac{x_1-x_2}{x_3-x_2}\times y_3 \\ &=y_1 \end{aligned}
f(x1)=l1(x1)⋅y1+l2(x1)⋅y2+l3(x1)⋅y3=x1−x2x1−x2×x1−x3x1−x3×y1+x2−x1x1−x1×x2−x3x1−x3×y2+x3−x1x1−x1×x3−x2x1−x2×y3=y1
怎么样怎么样。
这个式子满足了我们最开始的要求(通过已知的
n
个
点
(
x
i
,
y
i
)
n个点(x_i,y_i)
n个点(xi,yi)求出一个
n
−
1
n-1
n−1次的多项式
f
(
x
)
f(x)
f(x),使之满足
∀
i
∈
[
1
,
n
]
&
i
∈
N
+
,
有
f
(
x
i
)
=
y
i
\forall ~ i \in [1,n] \And i ~\in N_+,有f(x_i)=y_i
∀ i∈[1,n]&i ∈N+,有f(xi)=yi)。
数学家们看到这一长串的求和求积,满意了。
现在我们解决子问题②:让它不只适用于这n个点。
…………
已经搞定了。
为什么说人们巧妙地构造了这个插值基函数呢?
因为刚刚构造出来的这个东西真的大致对任意
x
∈
C
x\in C
x∈C成立。
为什么呢?
因为,刚才这个东西就是一个n-1次的多项式,我们完美完成了插值的任务。
耶!
结论:
拉格朗日插值就是把一堆堆点代入这个式子:
f
(
x
)
=
∑
1
≤
i
≤
n
∏
1
≤
j
≤
n
&
&
i
≠
j
x
−
x
j
x
i
−
x
j
⋅
y
i
f(x)=\sum_{1\le i\le n}~ ~ \prod_{1\le j \le n ~ \And\And ~ i ~\neq ~ j} ~~ \frac{x-x_j}{x_i-x_j}\cdot y_i
f(x)=1≤i≤n∑ 1≤j≤n && i = j∏ xi−xjx−xj⋅yi
拉格朗日插值的问题——龙格现象
但是,实际在投入生产运用中时,很好用但是也会有BUG。
当次数过大的时候,就可能会出现较大误差,这被称为拉格朗日插值的龙格现象,也就是拉格朗日插值的龙格现象。
解决办法就是在次数很大时,分段用较低次数的函数组合起来。
(听上去如此轻易,但实际上怎么做呢,我还真不知道)
好像是因为,你拿很高的次数(比方说20个点)来逼近一个较低次数的函数(比方说2次方),显然是不好的。
所以求插值的时候一定要算好要整几个点。
拉格朗日插值在信息学竞赛中的使用
可以用来给点求点。
举例举例:
尺八
Time Limit: 1000ms, Memory Limit: 8 MB
题目描述
行走在峡谷中,隐隐听见飘来了 Yasuo 悠扬的箫声…
Yasuo 随身携带的古乐器‘ 尺八’ 上有 5 个孔(-.-||),于是他让你帮他求一
求
∑
1
≤
i
≤
n
i
5
\sum_{1\leq i \leq n}i^5
1≤i≤n∑i5
由于和的值可能太大,只要求输出答案%998244353 的值
输入格式
首先一个正整数 T 表示数据组数.
接下来 T 行每行一个正整数 n
输出格式
T 行,每行一个正整数表示该组数据的答案
样例输入
3
1
2
3
样例输出
1
33
276
数据规模与约定
对于 20% 的数据,
n
≤
1000
,
T
=
1000
n \leq 1000 , T = 1000
n≤1000,T=1000。
对于 100% 的数据,
n
n
n在
l
o
n
g
l
o
n
g
long long
longlong 范围内,
T
≤
1
0
3
T \leq 10^3
T≤103
这是Ark在noiP前夕的胡扯 胡策 互测模拟赛上出的题的第一题,本来是希望同学们:
(摘自题解)
我不会告诉你我写的20分算法都爆炸了
Freopen/Kyle/wk大佬说,什么五次方啊,n次方都没有问题!
#include<cstdio>
#include<cstring>
#include<algorithm>
#define LL long long
#define mod 998244353
using namespace std;
LL n;
LL a[20],b[20];
inline int Pow(int base,int k)
{
int ret = 1;
for(;k;k>>=1,base=1ll*base*base%mod) if(k&1) ret = 1ll * ret * base % mod;
return ret;
}
int main()
{
freopen("flute.in","r",stdin);
freopen("flute.out","w",stdout);
for(int i=1;i<=19;i++) a[i] = i * i * i * i * i + a[i-1];
for(int i=0;i<=7;i++)
{
b[i] = a[1] % mod;
for(int j=1;j<=18;j++) a[j] = a[j + 1] - a[j];
}
int T;
for(scanf("%d",&T);T--;)
{
scanf("%lld",&n);
LL C = 1 , ans = 0;
for(int i=0;i<=7;i++)
{
ans = (ans + C * b[i]) % mod;
C = C * (n % mod - i - 1) % mod * Pow(i+1,mod-2) % mod;
}
printf("%lld\n",(ans+mod)%mod);
}
}
就是先在这个函数里求出七个点,然后代入。
代码实现
就看wk的代码好不好
下面的代码中,我们在函数
y
=
f
(
x
)
=
x
2
+
2
x
+
3
y=f(x)=x^2+2x+3
y=f(x)=x2+2x+3上取3/5/10个点,然后再取10个点验证一下。
代码的流程是:①读入点(我的代码是在上述二次函数上取点)并生成函数②代入点并求值(取10个x分别代入两个函数,比对结果)
来看吧。
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
struct Point
{
int x,y;
void get(int xx)
{
x=xx,y=xx*xx+3*xx+2;
}
}given[1000];
int n=3;
void point_get()
{
for(int i=1;i<=n;i++) given[i].get(i);
}
double calculate_l(int x,int i)
{
double ans=1;
for(int j=1;j<=n;j++)
if(i!=j) ans*=1.0*(x-given[j].x)/(given[i].x-given[j].x);
return ans;
}
double Lagrange(int x)
{
double ans=0;
for(int i=1;i<=n;i++) ans+=calculate_l(x,i)*given[i].y;
return ans;
}
int main()
{
point_get();
printf("when n is %d:\n",n);
for(int x=500;x<=5000;x+=500)
printf("%d: %d %.4lf\n",x,x*x+3*x+2,Lagrange(x));
}
when n is 3:
500: 251502 251502.0000
1000: 1003002 1003002.0000
1500: 2254502 2254502.0000
2000: 4006002 4006002.0000
2500: 6257502 6257502.0000
3000: 9009002 9009002.0000
3500: 12260502 12260502.0000
4000: 16012002 16012002.0000
4500: 20263502 20263502.0000
5000: 25015002 25015002.0000
when n is 5:
500: 251502 251502.0000
1000: 1003002 1003002.0000
1500: 2254502 2254502.0000
2000: 4006002 4006002.0000
2500: 6257502 6257502.0000
3000: 9009002 9009002.0000
3500: 12260502 12260502.1250
4000: 16012002 16012002.0000
4500: 20263502 20263502.0000
5000: 25015002 25015002.0000
when n is 8:
500: 251502 251790.5000
1000: 1003002 995456.0000
1500: 2254502 2426880.0000
2000: 4006002 12820480.0000
2500: 6257502 -9076736.0000
3000: 9009002 4063232.0000
3500: 12260502 -35651584.0000
4000: 16012002 950009856.0000
4500: 20263502 -3588227072.0000
5000: 25015002 -6664749056.0000
when n is 10:
500: 251502 4751360.0000
1000: 1003002 -660602880.0000
1500: 2254502 -37849399296.0000
2000: 4006002 3443490029568.0000
2500: 6257502 14164802142208.0000
3000: 9009002 -81123342286848.0000
3500: 12260502 -86723979640832.0000
4000: 16012002 1112156011495424.0000
4500: 20263502 -8983009998929920.0000
5000: 25015002 -17222750137483264.0000
所以一定要算好求几个点啊。
//当然了,这个代码是可以更精简的,是为了顺应拉格朗日插值的公式的思路写成这个样子的。
优化:重心拉格朗日插值
请详见参看文献第二个。
参考文献
下面的每一篇都写的比我好:
https://www.zhihu.com/question/58333118/answer/262507694
https://www.cnblogs.com/ECJTUACM-873284962/p/6833391.html
https://blog.youkuaiyun.com/shener_m/article/details/81706358
机房电脑上不了维基百科啊啊啊啊,连萌娘百科都上不去了只能看mooncell度日了啊!
说在后面
逆转系列太好玩了。