[BZOJ4126] 国王奇遇记 [拉格朗日插值]

最新推荐文章于 2021-07-18 10:17:16 发布

原创最新推荐文章于 2021-07-18 10:17:16 发布 · 375 阅读

0 ·

CC 4.0 BY-SA版权

拉格朗日插值专栏收录该内容

5 篇文章

订阅专栏

本文介绍了一种使用多项式差分和拉格朗日插值解决特定数学问题的算法，通过计算多项式的高阶差分找到函数间的关系，并运用拉格朗日插值公式求解未知点的值。代码实现涵盖了差分计算、拉格朗日插值及多项式幂运算。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

传送门

记 $S_n =\sum_{i=0}^{n-1} f_im^i$ ，那么有结论 $S_n = m^ng(n) - g(0)$

$S_{n}-S_{n-1}=m^ng(n)-m^{n-1}g(n-1)=f_{n-1}m^{n-1}\Rightarrow g(n)=\frac{g(n-1)+f(n-1)}{m}$

然后可以用 g(0) 表示 g(i)，最后利用差分找出 g(0), g(1)..., g(n+1) 的关系

关于差分，一个 k 次多项式 k + 1 阶差分为 0 ，也就是

$\sum_{i=0}^{k+1}(-1)^{k+1-i}\binom{k+1}{i}f(i)=0$

于是可以解出 g(0) 然后拉格朗日插出来

#include<bits/stdc++.h>
#define N 500050
using namespace std;
typedef long long ll;
const int Mod = 1000000007;
int n, m; ll f[N], g[N], a[N][2];
ll fac[N], ifac[N];
ll pre[N], suf[N];
ll add(ll a, ll b){ return (a + b) % Mod;}
ll mul(ll a, ll b){ return a * b % Mod;}
ll power(ll a, ll b){ ll ans = 1;
	for(;b;b>>=1){if(b&1) ans = mul(ans, a); a = mul(a, a);}
	return ans;
}
ll lagrange(int x){
	if(x <= m) return g[x];
	pre[0] = x; suf[m+1] = 1;
	for(int i = 1; i <= m; i++) pre[i] = mul(pre[i-1], x-i);
	for(int i = m; i >= 0; i--) suf[i] = mul(suf[i+1], x-i);
	ll ans = 0;
	for(int i = 0; i <= m; i++){
		ll v1 = mul(suf[i+1], (i==0)?1:pre[i-1]), v2 = ((m-i)&1) ? Mod-1 : 1;
		ans = add(ans, mul(v1, mul(mul(g[i], v2), mul(ifac[m-i], ifac[i]))));
	} return ans;
}
ll C(int n, int m){ return mul(fac[n], mul(ifac[n-m], ifac[m]));}
int main(){
	scanf("%d%d", &n, &m);
	if(m == 1){ cout << (n * (n+1) / 2) % Mod; return 0;}
	for(int i = 0; i <= m; i++) f[i] = power(i, m);
	fac[0] = fac[1] = ifac[0] = ifac[1] = 1;
	for(int i = 2; i <= m+1; i++) fac[i] = mul(fac[i-1], i);
	ifac[m+1] = power(fac[m+1], Mod-2);
	for(int i = m; i >= 2; i--) ifac[i] = mul(ifac[i+1], i+1);
	ll invm = power(m, Mod-2);
	a[0][0] = 0; a[0][1] = 1;
	for(int i = 1; i <= m+1; i++){
		a[i][1] = mul(a[i-1][1], invm);
		a[i][0] = mul(add(f[i-1], a[i-1][0]), invm); 
	}
	ll c = 0, d = 0;
	for(int i = 0; i <= m+1; i++){
		int v = (i&1) ? 1 : Mod-1;
		c = add(c, mul(v, mul(a[m+1-i][1], C(m+1, i))));
		d = add(d, mul(v, mul(a[m+1-i][0], C(m+1, i))));
	}
	g[0] = Mod - mul(power(c, Mod-2), d);
	for(int i = 1; i <= m+1; i++) g[i] = mul(add(g[i-1], f[i-1]), invm);
	ll ans = lagrange(n+1);
	cout << add(mul(power(m, n+1), ans), Mod - g[0]);
	return 0;
}