拉格朗日插值与牛顿插值

n+e posted @ 2015年6月30日 22:32 in Interesting with tags 多项式插值 , 1798 阅读

BC#45 1004  &&  HDU 5275

首先声明 $[bool]$ 表示当 $bool$ 为真时,值为 $1$,否则为 $0$

拉格朗日插值:

假设有三个点 $(x_1,y_1),(x_2,y_2),(x_3,y_3)$

我们的目的是让最终的多项式变成:当 $x=x_i$ 时 $y=y_i$

构造一个多项式使得当 $x=x_1$ 时的值为 $1$,当$x=x_2,x=x_3$ 时的值为 $0$

这个多项式长这样$$\frac{(x-x_2)\cdot (x-x_3)}{(x_1-x_2)\cdot (x_1-x_3)}$$

同理,构造 $[x=x_2]$ $$\frac{(x-x_1)\cdot (x-x_3)}{(x_2-x_1)\cdot (x_2-x_3)}$$

$[x=x_3]$ $$\frac{(x-x_1)\cdot (x-x_2)}{(x_3-x_1)\cdot (x_3-x_2)}$$

显然,对于最终的多项式,当 $x=x_1$ 时 $y=y_1$,就是 $y_1\cdot [x=x_1]$

还有两个条件,写出来就是 $y_2\cdot [x=x_2]$,$y_3\cdot [x=x_3]$

所以最终的多项式为$$\begin{array} \\f(x)&=&y_1\cdot [x=x_1]+y_2\cdot [x=x_2]+y_3\cdot [x=x_3]\\ &=&y_1\cdot \frac{(x-x_2)\cdot (x-x_3)}{(x_1-x_2)\cdot (x_1-x_3)}+y_2\cdot \frac{(x-x_1)\cdot (x-x_3)}{(x_2-x_1)\cdot (x_2-x_3)}+y_3\cdot \frac{(x-x_1)\cdot (x-x_2)}{(x_3-x_1)\cdot (x_3-x_2)}\\ \end{array}$$

推广到 $N$ 个点就是$$f(x)=\sum_{i=1}^N \prod_{j\ne i} \frac{x-x_j}{x_i-x_j}$$

把 $x$ 代进去就是答案。

 

至于这道题,首先先把1~250000的逆元线性求出来,然后开一个 $N*N$ 的数组记录一下前缀乘积 $\prod_{j\ne i}(x_i-x_j)$ 以及它的逆元,然后对于每一个询问,花费 $O(N)$ 的代价计算 $(x-x_j)$ 的前缀积和后缀积,然后就可以单次 $O(N)$ 得到答案啦~

#include<cstdio>
#define p 1000000007
#define N 3010
int n,x[N],y[N],q,l,r,x0,sum[N][N],isum[N][N],inv[250010],L[N],R[N],ans;
#define getinv(x) (x>=0?inv[x]:p-inv[-(x)])
int main(){
	scanf("%d",&n);inv[1]=1;
	for(int i=1;i<=n;i++)scanf("%d%d",x+i,y+i);
	for(int i=2;i<=250000;i++)inv[i]=1LL*(p-p/i)*inv[p%i]%p;
	for(int i=1;i<=n;i++){
		sum[i][0]=isum[i][0]=1;
		for(int j=1;j<i;j++)sum[i][j]=1LL*sum[i][j-1]*(p+x[i]-x[j])%p,isum[i][j]=1LL*isum[i][j-1]*getinv(x[i]-x[j])%p;
		sum[i][i]=sum[i][i-1],isum[i][i]=isum[i][i-1];
		for(int j=i+1;j<=n;j++)sum[i][j]=1LL*sum[i][j-1]*(p+x[i]-x[j])%p,isum[i][j]=1LL*isum[i][j-1]*getinv(x[i]-x[j])%p;
	}
	for(scanf("%d",&q);q--;printf("%d\n",ans)){
		scanf("%d%d%d",&l,&r,&x0);L[l-1]=R[r+1]=1;ans=0;
		for(int i=l;i<=r;i++)L[i]=1LL*L[i-1]*(p+x0-x[i])%p;
		for(int i=r;i>=l;i--)R[i]=1LL*R[i+1]*(p+x0-x[i])%p;
		for(int i=l;i<=r;i++)ans=(ans+1LL*sum[i][l-1]*isum[i][r]%p*L[i-1]%p*R[i+1]%p*y[i])%p;
	}
}

牛顿插值:

任何一个不高于 $n$ 次多项式,都可以表示成函数

$1,x-x_0,(x-x_0)(x-x_1),\cdots ,(x-x_0)(x-x_1)\cdots (x-x_{n-1})$

的线性组合。

$f(x)=c_0+c_1\cdot (x-x_0)+c_2\cdot (x-x_0)(x-x_1)+\cdots +c_n\cdot (x-x_0)(x-x_1)\cdots (x-x_{n-1})$

假设有 $3$ 个点 $(x_0,y_0),(x_1,y_1),(x_2,y_2)$

最后的多项式是 $f(x)=c_0+c_1\cdot (x-x_0)+c_2\cdot (x-x_0)(x-x_1)$

那么当 $x=x_0$ 的时候,$f(x)=y_0$,所以 $c_0=y_0$

当 $x=x_1$ 的时候,$f(x)=y_1$,所以

$$c_0+c_1\cdot (x_1-x_0)=y_1$$

所以 $$c_1=\frac{y_0-y_1}{x_0-x_1}$$

当 $x=x_2$ 的时候,$f(x)=y_2$,所以$$c_0+c_1\cdot (x_1-x_0)+c_2\cdot (x_2-x_0)(x_2-x_1)=y_2$$

所以 $$c_2=\frac{\frac{y_0-y_1}{x_0-x_1}-\frac{y_1-y_2}{x_1-x_2}}{x_0-x_2}$$

于是这个多项式就是$$f(x)=y_0+\frac{y_0-y_1}{x_0-x_1}\cdot (x-x_0)+\frac{\frac{y_0-y_1}{x_0-x_1}-\frac{y_1-y_2}{x_1-x_2}}{x_0-x_2}\cdot (x-x_0)(x-x_1)$$

推广一下就是下面这样:

定义 $f[x_i,x_j]=\frac{f(i)-f(j)}{x_i-x_j} (i\ne j,x_i\ne x_j)$ 为 $f(x)$ 在 $x_i,x_j$ 处的一阶差商。

$f[x_i,x_j,x_k]=\frac{f[x_i,x_j]-f[x_j,x_k]}{x_i-x_k}$ 为 $f(x)$ 在 $x_i,x_j,x_k$ 处的二阶差商。

$f[x_0,x_1,\cdots x_n]=\frac{f[x_0,x_1,\cdots ,x_{n-1}]-f[x_1,x_2,\cdots ,x_n]}{x_0-x_n}$ 为 $f(x)$ 在 $x_0,x_1,\cdots ,x_n$ 处的 $n$ 阶差商。

$x$ $f(x)$ 一阶差商 二阶差商 三阶差商
$x_0$ $f(x_0)$      
$x_1$ $f(x_1)$ $f[x_0,x_1]$    
$x_2$ $f(x_2)$ $f[x_1,x_2]$ $f[x_0,x_1,x_2]$  
$x_3$ $f(x_3)$ $f[x_2,x_3]$ $f[x_1,x_2,x_3]$ $f[x_0,x_1,x_2,x_3]$

然后最终的式子就是$$f(x)=f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+\cdots +f[x_0,x_1,\cdots ,x_{n}](x-x_0)(x-x_1)\cdots (x-x_{n-1})$$

就是上面这张表格对角线上的数即为牛顿插值多项式的系数。

于是处理出来这张表格就好了。

#include<cstdio>
#define N 3010
#define p 1000000007
int n,q,x[N],f[N][N],inv[250010];
#define getinv(x) (x>=0?inv[x]:p-inv[-(x)])
int main(){
	scanf("%d",&n);inv[1]=1;
	for(int i=2;i<=250000;i++)inv[i]=1LL*(p-p/i)*inv[p%i]%p;
	for(int i=1;i<=n;i++)scanf("%d%d",x+i,f[i]+1);
	for(int i=2;i<=n;i++)
	for(int j=2,k=i-1;j<=i;j++,k--)f[i][j]=1LL*(p+f[i-1][j-1]-f[i][j-1])*getinv(x[k]-x[i])%p;
	scanf("%d",&q);
	for(int l,r,x0,ans;q--;printf("%d\n",ans)){
		scanf("%d%d%d",&l,&r,&x0);ans=0;
		for(int i=l,j=1,sum=1;i<=r;i++,j++)ans=(ans+1LL*f[i][j]*sum)%p,sum=1LL*sum*(p+x0-x[i])%p;
	}
}
---By n+e
boardmodelpaper.com 说:
2024年1月06日 13:42

Board Model Papers 2024 Download with Suggestions for 10th Class Textbooks 2024 Pdf Download and SSLC New Syllabus Sample Question Paper 2024 and different types of model papers boardmodelpaper.com and question papers for following the website and Arts, Science, Commerce Stream Subject Wise Solved Question Bank for Hindi & English Medium Students with Exam Pattern & Blueprint and subject Wise with 11th & 12th Question Bank 2024 for General & Vocational Course Languages & Subjects Important Question for the above link.


登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter