拉格朗日插值法

传送门

我们遇到的问题是,给定一个多项式的点值表示和一个数,求出这个数带入多项式后的值。

这个问题如果用待定系数法,可以使用高斯消元,但是复杂度是(O(n^3))的,无法通过本题。

所以我们来引入拉格朗日插值法。

它的关键在于,有一个拉格朗日基本公式:

$$f(k) = sum_{i=0}^ny_iprod_{i eq j} frac{k - x_j}{x_i - x_j}$$

这个公式使得(f(x_i) = y_i) ,也就是同时满足这些点的表示。

原因在于,对于一个(x_i),只有在第(i)项,他后面的值是1,剩下的时候全是0,相当于全部被消掉了。(感觉这个地方和(CRT)的构造好像……),所以这个公式自然就满足所有点的表示了。

之后,我们就可以直接把值带入,模拟计算即可。时间复杂度(O(n^2))

特殊的,一般来讲我们遇到的题目,所有的(x_i)都是连续的。这样我们可以进一步优化这个算法,使之在(O(n))内计算出值。

具体的,我们再观察一下上面的式子。我们发现分母其实就是一个阶乘的形式。可以先预处理出来。然后对于分子,对于每个(k),我们可以先维护其前缀积和后缀积,即(pre_i = prod_{j = 0}^i k - j),(suf_i = prod_{j = i}^n j - k),式子就变成了

$$f(k) = sum_{i=0}^ny_i frac{pre_{i-1}suf_{i+1}}{fac_i fac_n-i}$$

之后就可以(O(n))计算了。 注意分母,在(n-i)为奇数的时候,分母是负数。

先这些吧……剩下的坑以后再填。

看一下代码。

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<iostream>
#include<cmath>
#include<set>
#include<vector>
#include<map>
#include<queue>
#define rep(i,a,n) for(int i = a;i <= n;i++)
#define per(i,n,a) for(int i = n;i >= a;i--)
#define enter putchar('
')
#define fr friend inline
#define y1 poj
#define mp make_pair
#define pr pair<int,int>
#define fi first
#define sc second
#define pb push_back

using namespace std;
typedef long long ll;
const int M = 100005;
const int INF = 1000000009;
const double eps = 1e-7;
const ll mod = 998244353;

ll read()
{
    ll ans = 0,op = 1;char ch = getchar();
    while(ch < '0' || ch > '9') {if(ch == '-') op = -1;ch = getchar();}
    while(ch >= '0' && ch <= '9') ans = ans * 10 + ch - '0',ch = getchar();
    return ans * op;
}

ll n,k,x[M],y[M],ans;

ll qpow(ll a,ll b)
{
   ll p = 1;
   while(b)
   {
      if(b & 1) p *= a,p %= mod;
      a *= a, a %= mod;
      b >>= 1;
   }
   return p;
}

int main()
{
   n = read(),k = read();
   rep(i,1,n) x[i] = read(),y[i] = read();
   rep(i,1,n)
   {
      ll cur = 1;
      rep(j,1,n)
      {
	 if(i == j) continue;
	 ll now = (x[i] - x[j] + mod) % mod;
	 cur = cur * (k - x[j] + mod) % mod * qpow(now,mod-2) % mod;
      }
      ans = ans + (y[i] * cur % mod),ans %= mod;
   }
   printf("%lld
",ans);
   return 0;
}

原文地址:https://www.cnblogs.com/captain1/p/10107031.html