拉格朗日插值法

传送门

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

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

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

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

$$f(k) = \sum_{i=0}^ny_i\prod_{i \neq 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('\n')
#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\n",ans);
   return 0;
}

posted @ 2018-12-12 10:37  CaptainLi  阅读(384)  评论(0编辑  收藏  举报