[Study Notes] Lagrangian Interpolation

The first step in learning polynomials.

References:

attack's luogu blog

oi wiki Lagrangian interpolation

Apocryphal's luogu blog

1. Introduction to Lagrangian Interpolation

question:

luogu P4781 [template] Lagrangian interpolation

Solution 1: Gaussian elimination

Obviously, there are infinite polynomials of \(\deg \geqslant n\), because according to Gaussian elimination, there must be free elements.

List these \(n\) points directly to the equation system and solve it with Gaussian elimination.

After finding the polynomial, you can find \(f(k)\).

Time complexity\(O(n^3).\)

Solution 2: Lagrangian Interpolation

Give a basis function for a point \((x_i,y_i)\):

\[g(k)=\prod_{j\not=i}^{1\leqslant j\leqslant n} \dfrac{k-x_j}{x_i-x_j} \]

It is easy to find:\(\forall j\not=i,g (x_j)=0.\) Because there is always\(k=x_j\) in the cumulative multiplication, so that\(k-x_j=x_j-x_j=0,g (x_j)=0.\)

For \(j=i,g(x_j)=1.\)because each item in the cumulative product is\(\dfrac{x_i-x_j}{x_i-x_j}=1,g(x_j)=1.\ )

So our polynomial can be expressed as:

\[f(k)=\sum_{i=1}^{n}y_i\times\prod_{j\not=i}^{1\leqslant j\leqslant n} \dfrac{k-x_j}{x_i-x_j} \]

In this way\(\forall(x_i,y_i),f(x_i)=y_i\), you can also find\(f(k).\)

Probably because the inverse element is required, the time complexity is \(O(n^2+n\log n)=O(n^2).\)

\(\tt{code:}\)

#include <bits/stdc++.h>
#define ll long long
using namespace std;
int n;
const ll mod = 998244353;
ll x[2011], y[2011], k, ans;
ll ksm(ll s1, ll s2) {
	if(!s2) return 1;
	if(s2 & 1) return s1 * ksm(s1, s2 - 1) % mod;
	ll ret = ksm(s1, s2 / 2);
	return ret * ret % mod;
}
int main() {
	scanf("%d%lld", &n, &k);
	for(int i = 1; i <= n; i++) scanf("%lld%lld", &x[i], &y[i]);
	for(int i = 1; i <= n; i++) {
		ll sum = 1, mul = 1;
		for(int j = 1; j <= n; j++) {
			if(i == j) continue;
			sum *= (k - x[j]); sum %= mod; sum += mod; sum %= mod;
			mul *= (x[i] - x[j]); mul %= mod; mul += mod; mul %= mod;
		}
		sum *= ksm(mul, mod-2); sum %= mod;
		ans += (sum * y[i]) % mod; ans %= mod;
	}
	printf("%lld\n", ans);
	return 0; 
}

2. Extension of Lagrangian interpolation method

question:

Same as above, adding \(x_i\) to the limit of continuous value.

solution:

Suppose \(x_i\in[1,n],\) is formulated as:

\[f(k)=\sum_{i=1}^{n}y_i\times\prod_{j\not=i}^{1\leqslant j\leqslant n} \dfrac{k-j}{i-j} \]

It is possible to maintain \(num=\prod_{1\leqslant j\leqslant n}k-j,\) so the numerator is partially reduced to \(\dfrac{num}{k-i},\) because doing so requires an inverse element, with a log, Very bad, so you can preprocess the prefix and \(pre_i\) and the suffix and \(suf_i.\)

For the denominator, it can be transformed into the form of \(1\times 2\times...\times i-1\times (-1)\times (-2)\times...\times(i-n)\).

In essence, it is a factorial, |denominator|\(=fac_{i-1}\times fac_{n-i},\) and then deal with the positive and negative.

This plus some \(O(n)\) preprocessing \(pre,suf,fac,inv...\), we only need \(O(n)\) to complete.

3. Specific application of Lagrange interpolation & Ex amp les

[example 1] CF622F The Sum of the k-th Powers

The classic model of Lagrangian interpolation, that is to find \(\sum_{i=1}^n i^k\).

This is a polynomial of degree \(k+1\), prove it perceptually:

Consider a polynomial of degree \(n\)\(g(x)=a_nx^n+a_{n-1}x^{n-1}+...+a_0\).

Its first-order difference\(\Delta g(x)=g(x+1)-g(x)=a_n(x+1)^n-a_nx^n+...\).

Through the binomial theorem, we find that the \(n\) item will be eliminated at this time, namely:

Each time a difference is made, the degree of the polynomial is reduced by one.

Let\(f(x)=\sum_{i=1}^ni^k,\) then\(\Delta f(x)=f(x+1)-f(x)=(x+1)^ k,\) is a polynomial of degree \(k\).

So \(f(x)\) is \(k+1\) degree polynomial.

At this point, we only need to determine \(k+2\) points and do an interpolation to find \(f(n)\).

Because \(i\) is continuous, the Lagrange interpolation extension in 2.2 can be used to achieve \(O(k\log k).\)

\(i^k\) can also be obtained linearly, so I won't explain it here.

\(\tt{code:}\)

#include <bits/stdc++.h>

#define ll long long

using namespace std;

const int N = 1e6;
const ll mod = 1e9 + 7;

ll n, k, fac[N+11], pre[N+11], suf[N+11], ans;

ll ksm(ll s1, ll s2) {
	if(!s2) return 1ll;
	if(s2 % 2) return s1 * ksm(s1, s2-1) % mod;
	ll ret = ksm(s1, s2/2);
	return ret * ret % mod;
}

int main() {
	scanf("%lld%lld", &n, &k);
	ll sum = 0; 
	fac[0] = 1;
	for(int i = 1; i <= k+2; i++) 
		fac[i] = fac[i-1] * (ll)i % mod;
	pre[0] = suf[k+3] = 1;
	for(int i = 1; i <= k+2; i++)
		pre[i] = pre[i-1] * (n - i) % mod;
	for(int i = k+2; i >= 1; i--)
		suf[i] = suf[i+1] * (n - i) % mod;
	for(int i = 1; i <= k+2; i++) {
		sum += ksm(i, k); sum %= mod;
		ll num = pre[i-1] * suf[i+1] % mod;
		ll fm = fac[i-1] * fac[k+2-i] % mod;
		num *= sum; num %= mod;
		num *= ksm(fm, mod-2); num %= mod;
		if((k+2-i) % 2 == 0) ans += num;
		else ans -= num; 
		ans %= mod; ans = (ans + mod) % mod;
	}
	printf("%lld\n", ans);
	return 0;
}

[Example 2] luogu P4593 [TJOI2018] Textbook-like profanity

Problem transformation: first find the number of \(k,\) that is "blasphemy". Ignore the \(a_i\) greater than \(n\) (although I don't know if there are any), it is easy to find that one "desecration" can make all monsters with continuous HP die, so \(k=m+0?1: (\exists a_i=1).\)

And because the blood volume of the monsters is different from each other, we can convert the calculated value into \(\sum_{i=1}^n i^k,\), so it is similar to [Example 1].

The code will not be released.

Posted by slindstr on Wed, 03 Aug 2022 01:34:27 +0930