知识点 - 多项式插值法

知识点 - 多项式插值法

解决问题类型：

• 已知f(0),f(1),f(2)…f(n)，求一个次数界为 $n$ 的多项式，满足这些取值。

• 求$n^k$ 的前缀和，或$n^k$ 的k阶前缀和公式（直接插值，不过会T。对于第一个问题，后面给出了组合数解法。）

公式blog

讲义

2.多项式插值算法

2.1 多项式插值的存在唯一性

多项式一直以来备受数学家们青睐，一方面它构造起来简单，另一方面它有非常美妙的性质，下面介绍多项式插值算法。

如果给定 $n$ 个横纵坐标分别互不相同的点 $(x_i,y_i),i=1,2,\dots ,n$，那么我们能否构造一个次数界为 $n$ 的多项式函数，使得它的函数图像恰好经过这 $n$ 个点？答案是肯定的，而且这个多项式函数是唯一的，证明如下：

设存在这样的一个多项式
$$f(x)=\sum _{i=0}^{n-1}{a}_i{x}^i$$
根据构造条件，有
$$\{\begin{array}{l}f(x_1)=y_1\\ f(x_2)=y_2\\ \dots \\ f(x_n)=y_n\end{array}$$
将上述线性方程组中的 $a_i,i=1,2,\dots ,n$ ，视为未知量，其系数矩阵的行列式 $A$ 恰好为范德蒙行列式，
$$V(x_1,x_2,\cdots ,x_n)=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ x_1& x_2& \cdots & x_n\\ x_1^2& x_2^2& \cdots & x_n^2\\ ⋮& ⋮& & ⋮\\ x_1^{n-1}& x_2^{n-1}& \cdots & x_n^{n-1}\end{array}\right]$$
故(数学归纳法证明)
$$det(A)=\prod _{1\le i<j\le n}(x_j-x_i)$$
又 $x_i互不相同,i=1,2,\dots ,n$，因此 $det(A)\ne 0$，方程组有唯一解

2.2 Lagrange多项式插值 O(N^2)

那么问题就来了，如何求这个多项式呢？利用高斯消元解上述线性方程组是一个办法，算法的复杂度为 $O(n^3)$ ，其实这个多项式可以直接被构造出来
$$f(x)=\sum _{i=1}^n\left(\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}\right)y_i$$
这就是Lagrange插值公式，不难验证其次数至多为 $n-1$ ，且满足上述线性方程组，因此这就是我们要求的多项式。

参考代码：O(n^2)

//输出在x处的取值
#include <bits/stdc++.h>
#define rep(i,a,b) for(ll i=a;i<=b;i++)
using namespace std;
typedef long long ll;
typedef long double type;
type lagrange(vector<type> x,vector<type> y,type X)
{
	int n=x.size()-1;
	type ans=0;
	rep(k,0,n)
	{
		type temp=y[k];
		rep(j,0,n)if(j!=k)temp*=(X-x[j])/(x[k]-x[j]);
		ans+=temp;
	}
	return ans;
}
int main()
{
	vector<type> x={0,1,2,3};
	vector<type> y={0,1,4,9};
	type X;
	while(cin>>X)cout<<lagrange(x,y,X)<<endl;
	return 0;
}

在ACM竞赛中，如果某个组合数学类题目刚好是输入一个整数 $n$ ，输出一个多项式函数的值 $f(n)$ ，那么上述算法只需要放入某几项即可，不需要推导复杂的公式，例如2018年icpc南京现场赛的G题，将上述算法的除法修改为乘法逆元即可，至于最终的公式是啥以及如何推导，本文暂不讨论。

参考代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mo=1e9+7;
ll fpow(ll a,ll b){
	ll ans=1;
	while(b>0){if(b&1)ans=ans*a%mo;b>>=1;a=a*a%mo;}
	return ans;
}
ll lagrange(vector<ll> x,vector<ll> y,ll X){
	auto p=y.begin();
	ll ans=0;
	for(auto k:x){
		ll a=*p++%mo,b=1;
		for(auto j:x)if(j!=k)a=(X-j)%mo*a%mo,b=(k-j)%mo*b%mo;
		ans=(ans+mo+a*fpow(b,mo-2)%mo)%mo;
	}
	return ans;
}
int main(){
	vector<ll> x={0,1,2,3,4};
	vector<ll> y={0,1,5,15,35};
	ll n;
	while(cin>>n)cout<<lagrange(x,y,n)<<endl;
	return 0;
}

其实在1.1的参考代码中已经给出了输出插值多项式的函数，可用如下方式调用

int main(){
	vector<type> x={frac(0,1),frac(1,1),frac(2,1),frac(3,1),frac(4,1)};
	vector<type> y={frac(0,1),frac(1,1),frac(5,1),frac(15,1),frac(35,1)};
	Poly f=Lagrange(x,y);
	cout<<f;
	ll X;
	while(cin>>X)cout<<f(frac(X,1))<<endl;
	return 0;
}

运行结果

(1/24)x^4+(1/4)x^3+(11/24)x^2+(1/4)x^1+(0)

2.2.1 x连续时的O(N)优化

当$x_i\in [0,n]$ 时，注意到下式的分母会变成阶乘
$$f(x)=\sum _{i=1}^n\left(\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}\right)y_i$$
前面是递减的连乘，直到$i==j$后面是带负号的,于是我们将分母写出来有
$$i!(-1{)}^{(n-i)}(n-i)!$$

对于分子可以预处理出连乘${\prod }_{j=1}^n(x-j)$，然后第 i 个的分子为 $S(i,x)=\frac{{\prod }_{j=1}^n(x-j)}{x-i}$ (实现用逆元或维护前后缀积)。

这样便$O(1)$地求出$\left({\prod }_{j\ne i}\frac{x-x_j}{x_i-x_j}\right)y_i$。总复杂度$O(N)$

南昌邀请赛的银牌题Polynomial： 题库链接

2.2.2 O(n)在线单点增加操作： 重心拉格朗日插值法

可以将拉格朗日基本多项式重新写为：（公式本地可见）

[外链图片转存失败(img-tO1Ptqxl-1565618397136)(http://upload.wikimedia.org/math/f/6/5/f65c55d7b4c330438e7a68482d1590e7.png)]

定义重心权[7][8]

[外链图片转存失败(img-8NiabXBb-1565618397137)(http://upload.wikimedia.org/math/f/5/b/f5b96ca21fad59d37cd3ea6a4c70246a.png)]

上面的表达式可以简化为：

[外链图片转存失败(img-rXSJuUtt-1565618397138)(http://upload.wikimedia.org/math/c/f/0/cf0ee67a5ae90253b230f700a86cb5c1.png)]

于是拉格朗日插值多项式变为：

[外链图片转存失败(img-V2KRbr1n-1565618397139)(http://upload.wikimedia.org/math/e/2/b/e2b06ed5610f4fcc046fef1b0a9f493e.png)]

即所谓的重心拉格朗日插值公式（第一型）或改进拉格朗日插值公式。它的优点是当插值点的个数增加一个时，将每个都除以[外链图片转存失败(img-iPPfzXMU-1565618397140)(http://upload.wikimedia.org/math/7/e/9/7e9644833883af0bcd5e70adb818536b.png)]，就可以得到新的重心权[外链图片转存失败(img-0ewyQlBX-1565618397141)(http://upload.wikimedia.org/math/e/9/2/e926afcf008b0e517216b346c6883710.png)]，计算复杂度为[外链图片转存失败(img-AW5mUFDb-1565618397142)(http://upload.wikimedia.org/math/d/4/3/d43308deb858cf186f862f2451af08b6.png)]，比重新计算每个基本多项式所需要的复杂度[外链图片转存失败(img-0FZGN8jw-1565618397142)(http://upload.wikimedia.org/math/5/6/d/56d1e7b3b345d81e79b9ff2a0023adb1.png)]降了一个量级。

2.3 Newton多项式插值

Lagrange插值算法对于ACM竞赛中的相关题目来说可能已经足够了，但Lagrange插值算法最初并不是为了这么用的，它的主要用途是构造一个多项式来逼近另外一个函数，例如我们用计算机可能没办法计算三角函数 $sin(x)$ 的精确值，但是如果已知其中某些点的值，就可以构造这样的一个多项式来逼近 $sin(x)$ ，就可以计算其近似值，误差即为 $sin(x)$ 的泰勒展开式中的Lagrange余项。对于复杂的函数，如果增加一个插值点，那么多项式就需要重新构造，求解单点处的值的复杂度为 $O(n^2)$ ，于是数学家们想出了另一个算法—Newton多项式插值

首先定义差商：

零阶差商
$$F(x_i)=y_i$$
$n$ 阶差商
$$F(x_1,x_2,\dots ,x_{n+1})=\frac{F(x_2,x_3,\dots ,x_{n+1})-F(x_1,x_2,\dots ,x_n)}{x_n-x_1}$$
那么
$$f(x)=\sum _{i=1}^n\left(F(x_1,\dots ,x_i)\prod _{j=1}^i(x-x_j)\right)$$
这样一来，这个算法就有了很好的继承性，每次添加一个插值点，复杂度为 $O(n)$ ，每次计算单点处的值，复杂度为 $O(n)$ （请参考秦九韶算法）。

参考代码（连续函数）

#include <bits/stdc++.h>
#define rep(i,a,b) for(ll i=a;i<=b;i++)
using namespace std;
typedef long long ll;
typedef long double type;
class NewtonPoly{
	public:
	type f[105],d[105],x[105];
	ll n=0;
	void add(type X,type Y){
		x[n]=X,f[n]=Y;
		rep(i,1,n)f[n-i]=(f[n-i+1]-f[n-i])/(x[n]-x[n-i]);
		d[n++]=f[0];
	}
	type cal(type X){
		type ans=0,t=1;
		rep(i,0,n-1)ans+=d[i]*t,t*=X-x[i];
		return ans;
	}
}P;
int main(){
	P.add(0,0);
	P.add(1,1);
	P.add(2,5);
	P.add(3,15);
	P.add(4,35);
	type x;
	while(cin>>x)cout<<P.cal(x)<<endl;
	return 0;
}

参考代码（离散函数，取余数，可用于ACM竞赛）

#include <bits/stdc++.h>
#define rep(i,a,b) for(ll i=a;i<=b;i++)
using namespace std;
typedef long long ll;
const ll mo=1e9+7;
ll fpow(ll a,ll b){
	ll ans=1;
	while(b>0){if(b&1)ans=ans*a%mo;b>>=1;a=a*a%mo;}
	return ans;
}
class NewtonPoly{
	public:
	ll f[105],d[105],x[105],n=0;
	void add(ll X,ll Y){
		x[n]=X,f[n]=Y%mo;
		rep(i,1,n)f[n-i]=(f[n-i+1]-f[n-i])%mo*fpow((x[n]-x[n-i])%mo,mo-2)%mo;
		d[n++]=f[0];
	}
	ll cal(ll X){
		ll ans=0,t=1;
		rep(i,0,n-1)ans=(ans+d[i]*t)%mo,t=(X-x[i])%mo*t%mo;
		return ans+mo*(ans<0);
	}
}P;
int main(){
	P.add(0,0);
	P.add(1,1);
	P.add(2,5);
	P.add(3,15);
	P.add(4,35);
	ll x;
	while(cin>>x)cout<<P.cal(x)<<endl;
	return 0;
}

2.5 高维插值整式

我们解决了一维的情况，二维的情况可由一维的Lagrange插值函数推广得来，更高维也类似
$$f(x,y)=\sum _{n=1}^N\sum _{m=1}^M\left(\prod _{i\ne n}\frac{x-x_i}{x_n-x_i}\right)\left(\prod _{j\ne m}\frac{y-y_j}{y_m-y_j}\right)z_{n,m}$$
参考代码（连续函数）

#include <bits/stdc++.h>
#define rep(i,a,b) for(ll i=a;i<=b;i++)
using namespace std;
typedef long long ll;
typedef long double type;
type lagrange2(vector<type> x,vector<type> y,vector<vector<type> > z,type X,type Y){
	int M=x.size()-1,N=y.size()-1;
	type ans=0;
	rep(m,0,M)rep(n,0,N){
		type t=z[m][n];
		rep(i,0,M)if(i!=m)t*=(X-x[i])/(x[m]-x[i]);
		rep(i,0,N)if(i!=n)t*=(Y-y[i])/(y[n]-y[i]);
		ans+=t;
	}
	return ans;
}
int main(){
	vector<type> x={1,2};
	vector<type> y={3,4};
	vector<vector<type> > z={{3,4},{6,8}};
	type X,Y;
	while(cin>>X>>Y)cout<<lagrange2(x,y,z,X,Y)<<endl;
	return 0;
}

参考代码（离散函数，取余数，可用于ACM竞赛）

#include <bits/stdc++.h>
#define rep(i,a,b) for(ll i=a;i<=b;i++)
using namespace std;
typedef long long ll;
const ll mo=1e9+7;
ll fpow(ll a,ll b){
	ll ans=1;
	while(b>0){if(b&1)ans=ans*a%mo;b>>=1;a=a*a%mo;}
	return ans;
}
ll lagrange2(vector<ll> x,vector<ll> y,vector<vector<ll> > z,ll X,ll Y){
	ll M=x.size()-1,N=y.size()-1,ans=0;
	rep(m,0,M)rep(n,0,N){
		ll a=z[m][n]%mo,b=1;
		rep(i,0,M)if(i!=m)a=(X-x[i])%mo*a%mo,b=(x[m]-x[i])%mo*b%mo;
		rep(i,0,N)if(i!=n)a=(Y-y[i])%mo*a%mo,b=(y[n]-y[i])%mo*b%mo;
		ans=(ans+a*fpow(b,mo-2)%mo)%mo;
	}
	return ans+mo*(ans<0);
}
int main(){
	vector<ll> x={0,1,2,3,4,5,6,7,8,9};
	vector<ll> y={0,1,2,3,4,5,6,7,8,9};
	vector<vector<ll> > z={
		{0,0,0,0,0,0,0,0,0,0},
		{0,1,2,3,4,5,6,7,8,9},
		{-2,2,6,10,14,18,22,26,30,34},
		{-10,0,10,20,30,40,50,60,70,80},
		{-30,-10,10,30,50,70,90,110,130,150},
		{-70,-35,0,35,70,105,140,175,210,245},
		{-140,-84,-28,28,84,140,196,252,308,364},
		{-252,-168,-84,0,84,168,252,336,420,504},
		{-420,-300,-180,-60,60,180,300,420,540,660},
		{-660,-495,-330,-165,0,165,330,495,660,825}
	};
	ll X,Y;
	while(cin>>X>>Y){
		cout<<lagrange2(x,y,z,X,Y)<<endl;
	}
	return 0;
}

1. 多项式的前缀和

1.1 $n^k$ 的前缀和

首先来看几个众所周知的公式
$$\sum _{i=1}^{n}1=n$$

$$\sum _{i=1}^{n}i=\frac{n(n+1)}{2}$$

$$\sum _{i=1}^n{i}^2=\frac{n(n+1)(2n+1)}{6}$$

用数学归纳法很容易验证上述公式的正确性，但是对于任何给定的非负整数 $k$ ，如何求出 ${\sum }_{i=1}^n{i}^k$ 呢？下面给出一种利用杨辉三角的计算方法。
$$\begin{array}{ccccc}1& 1& 1& 1& 1\\ 1& 2& 3& 4& 5\\ 1& 3& 6& 10& 15\\ 1& 4& 10& 20& 35\\ 1& 5& 15& 35& 70\end{array}$$

杨辉三角的自然数形式

$$\begin{array}{ccccc}{C}_0^0& C_1^1& C_2^2& C_3^3& C_4^4\\ C_1^0& C_2^1& C_3^2& C_4^3& C_5^4\\ C_2^0& C_3^1& C_4^2& C_5^3& C_6^4\\ C_3^0& C_4^1& C_5^2& C_6^3& C_7^4\\ C_4^0& C_5^1& C_6^2& C_7^3& C_8^4\end{array}$$

杨辉三角的组合数形式

不难发现，第$i$列的前$n$个数字之和刚好等于第$i+1$列第$n$个数字，即
$$C_k^k+C_{k+1}^k+\cdots +C_{k+n}^k=C_{k+n+1}^{k+1}$$
证 根据杨辉恒等式 $C_n^k=C_{n-1}^{k-1}+C_{n-1}^k$ ，$左边=C_k^k+C_{k+1}^k+\cdots +C_{k+n}^k=C_{k+1}^{k+1}+C_{k+1}^k+\cdots +C_{k+n}^k=C_{k+2}^{k+1}+\cdots +C_{k+n}^k=\cdots =C_{k+n+1}^{k+1}=右边$

换言之，除第一列外，每一列都是前一列的“前缀和”，而每一列的数字，都是 $C_n^k(k=列数-1)$ 的形式，其本质就是 $n$ 的多项式，例如：
$$第一列:C_n^0=1$$

$$第二列：C_n^1=\sum _{i=0}^{n-1}{C}_i^0=\sum _{i=0}^{n-1}1=n$$

$$第三列：C_n^2=\sum _{i=1}^{n-1}{C}_i^1=\sum _{i=1}^{n-1}i=\frac{n(n-1)}{2}$$

$$第四列：C_n^3=\sum _{i=2}^{n-1}{C}_i^2=\sum _{i=2}^{n-1}\frac{i(i-1)}{2}=\frac{n(n-1)(n-2)}{6}$$

根据第三列，得到
$$\sum _{i=1}^{n}i=\sum _{i=1}^{n-1}i+n=\frac{n(n+1)}{2}$$
根据第四列，得到
$$\sum _{i=2}^{n-1}\frac{i(i-1)}{2}=\frac{1}{2}\sum _{i=2}^{n-1}{i}^2-\frac{1}{2}\sum _{i=2}^{n-1}i=\frac{1}{2}\left(\sum _{i=1}^n{i}^2-1-n^2\right)-\frac{(n+1)(n-2)}{4}=\frac{n(n-1)(n-2)}{6}$$
进而得到
$$\sum _{i=1}^n{i}^2=\frac{n(n+1)(2n+1)}{6}$$
代码太长，令开一篇给出。

例题

南昌邀请赛的银牌题Polynomial：题库链接 题意：x取值连续的拉格朗日插值+前缀和的板题

代码

#include<bits/stdc++.h>
using namespace std;
const int mod=9999991;
const int N=1005;
long long a[N];
long long sum[N];
long long fac[N];
long long inv[N];
long long _inv[mod+5];
long long quickmod(long long a,long long b)
{
    long long ans=1;
    while(b)
    {
        if(b%2==1)
            ans=ans*a%mod;
        a=a*a%mod;
        b=b/2;
    }
    return ans;
}
int main()
{
    fac[0]=1;
    for(int i=1; i<N; i++)
        fac[i]=1ll*(fac[i-1])*i%mod;
    inv[N-1]=quickmod(fac[N-1],mod-2);
    for(int i=N-2; i>=0; i--)
        inv[i]=(1ll*inv[i+1]*(i+1))%mod;
    _inv[0]=_inv[1]=1;
    for(int i=2;i<mod+5;i++)
        _inv[i]=((mod-mod/i)*_inv[mod%i])%mod;
    int t;
    scanf("%d",&t);
    while(t--)
    {
        int n,m;
        scanf("%d%d",&n,&m);
        memset(a,0,sizeof(a));
        memset(sum,0,sizeof(sum));
        long long z=1;
        for(int i=0; i<=n; i++)
        {
            scanf("%d",&a[i]);
            a[i]=a[i]%mod;
            z=z*((n+1-i+mod)%mod)%mod;
            if(i!=0)
                sum[i]=(sum[i-1]+a[i])%mod;
            else
                sum[i]=a[i];
        }
        for(int i=0; i<=n; i++)
        {
            if((n-i)%2==1)
                a[n+1]=(a[n+1]-z%mod*_inv[n+1-i]%mod*inv[n-i]%mod*inv[i]%mod%mod*a[i]%mod+mod)%mod;
            else
                a[n+1]=(a[n+1]+z%mod*_inv[n+1-i]%mod*inv[n-i]%mod*inv[i]%mod%mod*a[i]%mod+mod)%mod;
        }
        sum[n+1]=(sum[n]+a[n+1])%mod;
        while(m--)
        {
            int l,r;
            scanf("%d%d",&l,&r);
            if(r<=n+1)
            {
                printf("%d\n",(sum[r]-sum[l-1]+mod)%mod);
                continue;
            }
            long long x=1,y=1;
            for(int i=0; i<=n+1; i++)
            {
                x=x*(r-i)%mod;
                y=y*(l-1-i)%mod;
            }
            long long ans=0;
            for(int i=0; i<=n+1; i++)
            {
                if((n+1-i)%2==1)
                    ans=(ans-x%mod*_inv[r-i]%mod*inv[n+1-i]%mod*inv[i]%mod%mod*sum[i]%mod+mod)%mod;
                else
                    ans=(ans+x%mod*_inv[r-i]%mod*inv[n+1-i]%mod*inv[i]%mod%mod*sum[i]%mod+mod)%mod;
            }
            if(l-1<=n+1)
                ans=(ans-sum[l-1]+mod)%mod;
            else
            {
                for(int i=0; i<=n+1; i++)
                {
                    if((n+1-i)%2==1)
                        ans=(ans+y%mod*_inv[l-1-i]%mod*inv[n+1-i]%mod*inv[i]%mod%mod*sum[i]%mod+mod)%mod;
                    else
                        ans=(ans-y%mod*_inv[l-1-i]%mod*inv[n+1-i]%mod*inv[i]%mod%mod*sum[i]%mod+mod)%mod;
                }
            }
            printf("%lld\n",ans);
        }
    }
    return 0;
}
/*
1
3 2
1 10 49 142
6 7
95000 100000
*/