拉格朗日插值

So far 没啥好说的。终于有东西了。

原本的插值

给定 \(n\) 个坐标系上的点（\(x\) 坐标两两不同），让你求一个最高不超过 \(n-1\) 次的多项式函数经过这些点。

设 \(n\) 个点的坐标为 \((x_1,y_1),\dots,(x_n,y_n)\)，则函数

\[\large{ f(x)=\sum_{i=1}^ny_i\prod_{j\ne i}\frac{x-x_j}{x_i-x_j} }\]

即为所求。

证明：

只要将 \(x=x_1,x_2,\dots,x_n\) 代入原式一切就明白了，每次代入都只有一项为 \(y_i\)，其他为 \(0\)。

若 \(x_i\) 没有规律的话复杂度为 \(O(n^2)\)，也比高斯消元优秀。

P4781 【模板】拉格朗日插值

远古的提交

P5667 拉格朗日插值2

配合多项式。

提交

CF622F The Sum of the k-th Powers

求 \(\sum_{i=1}^ni^k \bmod (10^9+7)\)。

\(n\) 很大，要求 \(O(k\log k)\)。

不妨设答案为 \(S_k(n)\)。

其实，\(S_k(n)\) 是关于 \(n\) 的 \(k+1\) 次多项式，给出证明。

证明 1

做一下差分，发现 \(\Delta S_k(n)=n^k\)，是一个 \(k\) 次多项式，而 \(k\) 次多项式的前缀和为 \(k+1\) 次多项式。

证明 2（仅供参考）

参考 我写的

\[S_{k}(n)=\frac{1}{k+1}\sum_{r=1}^{k+1}\binom{k+1}{r}B_{k+1-r}(n+1)^r \]

其中 \(B_i\) 为伯努利数。

---

所以说我们任意找不同的 \(k+2\) 个点进行插值即可求出 \(S_{k}(n)\)。

我们就选取好算的 \(1\sim k+2\) 来插值，发现每个点的值和插值时的系数可以快速求出，就解决了。

点击查看代码

//https://www.luogu.com.cn/problem/CF622F
#include<bits/stdc++.h>
#define mod 1000000007
#define N 1000010
using namespace std;

int pl[N], pr[N], fac[N];
int qpow(int a, int b) {
	int ans = 1;
	for(; b >= 1; b >>= 1, a = 1ll * a * a % mod)
		if(b & 1) ans = 1ll * ans * a % mod;
	return ans;
}
signed main() {
	int n, k, y = 0, ans = 0;
	scanf("%d%d", &n, &k);
	pl[0] = pr[k + 3] = fac[0] = 1;
	for(int i = 1; i <= k + 2; i ++)
		pl[i] = 1ll * pl[i - 1] * (n - i) % mod;
	for(int i = k + 2; i >= 1; i --)
		pr[i] = 1ll * pr[i + 1] * (n - i) % mod;
	for(int i = 1; i <= k + 2; i ++)
		fac[i] = 1ll * fac[i - 1] * i % mod;
	for(int i = 1; i <= k + 2; i ++) {
		y = (y + qpow(i, k)) % mod;
		int a = pl[i - 1] * 1ll * pr[i + 1] % mod;
		int b = fac[i - 1] * ((k - i) & 1 ? -1ll : 1ll) * fac[k + 2 - i] % mod;
		ans = (ans + 1ll * y * a % mod * qpow(b, mod - 2) % mod) % mod;
	}
	printf("%d\n", (ans + mod) % mod);
	return 0;
}

BZOJ3453. tyvj 1858 XLkxc

给定 \(k,a,n,d,p\)。

\(f(x)=\sum_{i=1}^{x}i^k\)

\(g(x)=\sum_{i=1}^{x}f(i)\)

\(h(x)=\sum_{i=0}^{x}g(a+id)\)

求 \(h(n)\bmod p\)。

\(1\le k\le 123,0\le a,n,d\le 123456789\)

\(p=1234567891\) 一个素数。

和上一道差不多，就是多套了几层，插值 \(k+4\) 个点即可。

时间 \(O(k^3)\)。

细节有一些。

点击查看代码

#pragma GCC optimize("Ofast")
#include<bits/stdc++.h>
using namespace std;
#define For(i,j,k) for(int i=(j),i##_=(k);i<=i##_;i++)
#define Rof(i,j,k) for(int i=(j),i##_=(k);i>=i##_;i--)
#define ll long long
const ll P=1234567891;//prime
const int N=130;
ll pw(ll x,ll y){
	ll res=1;
	while(y){
		if(y&1){
			(res*=x)%=P;
		}
		y>>=1;
		(x*=x)%=P;
	}
	return res;
}
int k;
ll A,n,d,v[N],f[N],g[N],h[N];
ll calc(int lim,ll x){//[0,lim) 内整数已知点值，求 f(x)（x 对 P 取模）
	ll ans=0,res;
	For(i,0,lim-1){
		res=1;
		For(j,0,lim-1) if(i!=j) (res*=i-j+P)%=P;
		res=pw(res,P-2)*v[i]%P;
		For(j,0,lim-1) if(i!=j) (res*=x-j+P)%=P;
		(ans+=res)%=P;
	}
	return ans;
}
void work(){
	cin>>k>>A>>n>>d;
	For(i,1,k+2) f[i]=(f[i-1]+pw(i,k))%P;
	For(i,1,k+2) g[i]=(g[i-1]+f[i])%P;
	copy(g,g+k+3,v);
	h[0]=calc(k+3,A);
	For(i,1,k+3) h[i]=(h[i-1]+calc(k+3,(A+i*d)%P))%P;
	copy(h,h+k+4,v);
	cout<<calc(k+4,n)<<"\n";
}
signed main(){ios::sync_with_stdio(false),cin.tie(nullptr);
	int T;cin>>T;
	while(T--)work();
return 0;}
//插值式子上面是 x-x_j 不是 x_i-x，公式记牢
//这里模数两倍爆 int！！

[AGC061F] Perfect Strings

看题解区吧，我也不是非常懂，属于拉插高端运用了。