拉格朗日插值与牛顿插值

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

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

登录 *


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