拉格朗日插值法学习笔记

胡说八道

啥都不会
啥都想学
大概废了

拉格朗日插值

给你$n$个点$(x_i,y_i)$，求经过这些点的$n+1$次函数及求某个$x$在该函数上的取值
由于$n$个点组成的函数是固定的
所以你可以考虑构造一种函数，使得当$x=x_i$时取值$y=y_i$，显然这种函数就是解
可以写成
$$\sum _{i=1}^n\frac{{\Pi }_{j!=i}(x-x_j)y_i}{{\Pi }_{j!=i}(x_i-x_j)}$$
用心感受一下
他就能满足上述条件
单次搞是$O(n^2)$的

优化

显然分数线上方式子是可以$O(n)$预处理$O(1)$求出的
问题在于下方式子
如果$x_i$是连续的话，下面也是可以$O(n)$的
比如$x_i=2$
下面是$(1)\ast (-1)\ast (-2)\ast ...$
$x_i=3$
下面就是$(2)\ast (1)\ast (-1)\ast (-2)....$
就这样吧.

代码

LL y[2100],inv1[2100],inv2[2100];
LL pre[2100],suf[2100];
LL polation(int K,LL x)
{
	LL ret=0;
	LL s2=1,bak=1-K,fro=0;
	pre[0]=1;suf[K+1]=1;
	for(int i=1;i<=K;i++)pre[i]=pre[i-1]*((x-i+mod)%mod)%mod;
	for(int i=K;i>=1;i--)suf[i]=suf[i+1]*((x-i+mod)%mod)%mod;
	for(int i=2;i<=K;i++)s2=s2*((1-i+mod)%mod)%mod;
	for(int i=1;i<=K;i++)
	{
		LL t1=pre[i-1]*suf[i+1]%mod*y[i]%mod;
		LL t2=s2;
		ret=(ret+t1*pow_mod((t2+mod)%mod,mod-2)%mod)%mod;
		s2=s2*(++fro)%mod;
		if(bak<0)s2=s2*inv1[-bak]%mod;
		else s2=s2*inv2[bak]%mod;
		bak++;
	}
	return ret;
}

计算系数

貌似找了半天都没有这个板子…就丢一个上来
刚好最近abc还出了这个…虽然我因为预处理挂$E$没时间做啊我的大好上分机会
考虑拉格朗日的式子
$$\sum _{i=1}^n\frac{{\Pi }_{j!=i}(x-x_j)y_i}{{\Pi }_{j!=i}(x_i-x_j)}$$
考虑对于每个$i$搞出一个多项式然后求和
我们发现对于每个$i$，下面和右边乘的那个$y_i$都是常数，需要计算的多项式实际上只是
$$\frac{\Pi (x-x_j)}{x-x_i}$$
所以预处理$\Pi (x-x_i)$的多项式，然后枚举每个把$(x-x_i)$的贡献退掉
过程类似一个退背包?
就从后往前做，每次自己先退$1$的贡献，再给前面的退掉$0$的贡献
代码中$a_i$表示$x=i$时多项式值为$a_i$

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#include<queue>
#include<vector>
#include<ctime>
#include<map>
#include<bitset>
#include<set>
#define LL long long
#define mp(x,y) make_pair(x,y)
#define pll pair<long long,long long>
#define pii pair<int,int>
using namespace std;
inline LL read()
{
	LL f=1,x=0;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
	return x*f;
}
int stack[20];
template<typename T>inline void write(T x)
{
	if(x<0){putchar('-');x=-x;}
    if(!x){putchar('0');return;}
    int top=0;
    while(x)stack[++top]=x%10,x/=10;
    while(top)putchar(stack[top--]+'0');
}
template<typename T>inline void pr1(T x){write(x);putchar(' ');}
template<typename T>inline void pr2(T x){write(x);putchar('\n');}
const int MAXN=3005;
int f[MAXN],mod,a[MAXN];
int pow_mod(int a,int b)
{
	int ret=1;
	for(;b;b>>=1,a=1LL*a*a%mod)if(b&1)ret=1LL*ret*a%mod;
	return ret;
}
int g[MAXN],h[MAXN];
void ad(int &x,int y){x+=y;if(x>=mod)x-=mod;}
void dl(int &x,int y){x-=y;if(x<0)x+=mod;}
void mul(int *x,int *y,int len)
{
	static int A[MAXN];memset(A,0,sizeof(A));
	for(int i=0;i<len;i++)for(int j=0;j<=1;j++)
		ad(A[i+j],1LL*x[i]*y[j]%mod);
	for(int i=0;i<=len;i++)x[i]=A[i];
}
void div(int *x,int *y,int len,int mu)
{
	static int A[MAXN],B[MAXN];memset(B,0,sizeof(B));
	int inv=pow_mod(y[1],mod-2);
	for(int i=0;i<len;i++)A[i]=x[i];
	for(int i=len-1;i>0;i--)
	{
		ad(B[i-1],A[i]);dl(A[i-1],1LL*A[i]*y[0]%mod);
	}
	for(int i=0;i<len;i++)ad(h[i],1LL*B[i]*mu%mod);
}
int main()
{
	mod=read();
	for(int i=0;i<mod;i++)a[i]=read();
	f[0]=1;
	for(int i=0;i<mod;i++)
	{
		g[0]=(-i+mod)%mod;g[1]=1;
		mul(f,g,i+1);
	}
	for(int i=0;i<mod;i++)
	{
		int mu=1;
		for(int j=0;j<mod;j++)if(i!=j)mu=1LL*mu*(i%mod-j%mod+mod)%mod;
		mu=1LL*pow_mod(mu,mod-2)*a[i]%mod;g[0]=(-i+mod)%mod;g[1]=1;
		div(f,g,mod+1,mu);
	}
	for(int i=0;i<mod;i++)pr1(h[i]);
	return 0;
}