refactor: 重写 Lagrange 插值 (#5387)

* refactor: rename

* refactor: rewrite

* feat: code
pull/5334/head
Tifa 2024-02-05 17:23:06 +08:00 committed by GitHub
parent 7964b721f0
commit a4ca3d9172
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
22 changed files with 628 additions and 192 deletions

View File

@ -98,7 +98,8 @@
/math/poly/ln-exp /math/poly/elementary-func
/math/poly/tri-func /math/poly/elementary-func
/math/poly/inv-tri-func /math/poly/elementary-func
/math/poly/lagrange /math/numerical/lagrange
/math/poly/lagrange /math/numerical/interp
/math/numerical/lagrange /math/numerical/interp
/math/poly/cntt /math/poly/ntt
/math/linear-algebra/gauss /math/numerical/gauss
/math/integral /math/numerical/integral

View File

@ -0,0 +1,109 @@
#include <bits/stdc++.h>
using namespace std;
constexpr uint32_t MOD = 998244353;
struct mint {
uint32_t v_;
mint() : v_(0) {}
mint(int64_t v) {
int64_t x = v % (int64_t)MOD;
v_ = (uint32_t)(x + (x < 0 ? MOD : 0));
}
friend mint inv(mint const &x) {
int64_t a = x.v_, b = MOD;
if ((a %= b) == 0) return 0;
int64_t s = b, m0 = 0;
for (int64_t q = 0, _ = 0, m1 = 1; a;) {
_ = s - a * (q = s / a);
s = a;
a = _;
_ = m0 - m1 * q;
m0 = m1;
m1 = _;
}
return m0;
}
mint &operator+=(mint const &r) {
if ((v_ += r.v_) >= MOD) v_ -= MOD;
return *this;
}
mint &operator-=(mint const &r) {
if ((v_ -= r.v_) >= MOD) v_ += MOD;
return *this;
}
mint &operator*=(mint const &r) {
v_ = (uint32_t)((uint64_t)v_ * r.v_ % MOD);
return *this;
}
mint &operator/=(mint const &r) { return *this = *this * inv(r); }
friend mint operator+(mint l, mint const &r) { return l += r; }
friend mint operator-(mint l, mint const &r) { return l -= r; }
friend mint operator*(mint l, mint const &r) { return l *= r; }
friend mint operator/(mint l, mint const &r) { return l /= r; }
};
template <class T>
struct NewtonInterp {
// {(x_0,y_0),...,(x_{n-1},y_{n-1})}
vector<pair<T, T>> p;
// dy[r][l] = [y_l,...,y_r]
vector<vector<T>> dy;
// (x-x_0)...(x-x_{n-1})
vector<T> base;
// [y_0]+...+[y_0,y_1,...,y_n](x-x_0)...(x-x_{n-1})
vector<T> poly;
void insert(T const &x, T const &y) {
p.emplace_back(x, y);
size_t n = p.size();
if (n == 1) {
base.push_back(1);
} else {
size_t m = base.size();
base.push_back(0);
for (size_t i = m; i; --i) base[i] = base[i - 1];
base[0] = 0;
for (size_t i = 0; i < m; ++i)
base[i] = base[i] - p[n - 2].first * base[i + 1];
}
dy.emplace_back(p.size());
dy[n - 1][n - 1] = y;
if (n > 1) {
for (size_t i = n - 2; ~i; --i)
dy[n - 1][i] =
(dy[n - 2][i] - dy[n - 1][i + 1]) / (p[i].first - p[n - 1].first);
}
poly.push_back(0);
for (size_t i = 0; i < n; ++i) poly[i] = poly[i] + dy[n - 1][0] * base[i];
}
T eval(T const &x) {
T ans{};
for (auto it = poly.rbegin(); it != poly.rend(); ++it) ans = ans * x + *it;
return ans;
}
};
int main() {
NewtonInterp<mint> ip;
int n, k;
cin >> n >> k;
for (int i = 1, x, y; i <= n; ++i) {
cin >> x >> y;
ip.insert(x, y);
}
cout << ip.eval(k).v_;
return 0;
}

View File

@ -0,0 +1,10 @@
id: oi-wiki:math/interp-1
name: Interpolation
solution:
cpp: docs/math/code/numerical/interp/interp_1.cpp
target:
- 'test'
testcases:
answer: docs/math/examples/numerical/interp/interp_1.ans
input: docs/math/examples/numerical/interp/interp_1.in
version: 1000

View File

@ -0,0 +1,10 @@
id: oi-wiki:math/interp-2
name: Interpolation
solution:
cpp: docs/math/code/numerical/interp/interp_2.cpp
target:
- 'test'
testcases:
answer: docs/math/examples/numerical/interp/interp_2.ans
input: docs/math/examples/numerical/interp/interp_2.in
version: 1000

View File

@ -0,0 +1,10 @@
id: oi-wiki:math/interp-3
name: Interpolation
solution:
cpp: docs/math/code/numerical/interp/interp_3.cpp
target:
- 'test'
testcases:
answer: docs/math/examples/numerical/interp/interp_3.ans
input: docs/math/examples/numerical/interp/interp_3.in
version: 1000

View File

@ -1,10 +0,0 @@
id: oi-wiki:math/lagrange-1
name: Lagrange
solution:
cpp: docs/math/code/numerical/lagrange/lagrange_1.cpp
target:
- 'test'
testcases:
answer: docs/math/examples/numerical/lagrange/lagrange_1.ans
input: docs/math/examples/numerical/lagrange/lagrange_1.in
version: 1000

View File

@ -1,10 +0,0 @@
id: oi-wiki:math/lagrange-2
name: Lagrange
solution:
cpp: docs/math/code/numerical/lagrange/lagrange_2.cpp
target:
- 'test'
testcases:
answer: docs/math/examples/numerical/lagrange/lagrange_2.ans
input: docs/math/examples/numerical/lagrange/lagrange_2.in
version: 1000

View File

@ -0,0 +1 @@
10201

View File

@ -0,0 +1,4 @@
3 100
1 4
2 9
3 16

View File

@ -0,0 +1,82 @@
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
<svg width="600" height="480" viewBox="0 0 600 480" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">
<desc>Produced by GNUPLOT 4.2 patchlevel 0 </desc>
<defs>
<circle id="gpDot" r="1"/>
<path id="gpPt0" style="stroke-width:0.107" d="M-1,0 h2 M0,-1 v2"/>
<path id="gpPt1" style="stroke-width:0.107" d="M-1,-1 L1,1 M1,-1 L-1,1"/>
<path id="gpPt2" style="stroke-width:0.107" d="M-1,0 L1,0 M0,-1 L0,1 M-1,-1 L1,1 M-1,1 L1,-1"/>
<rect id="gpPt3" style="stroke-width:0.107" x="-1" y="-1" width="2" height="2"/>
<use xlink:href="#gpPt3" id="gpPt4" style="fill:currentColor; stroke:none"/>
<circle id="gpPt5" style="stroke-width:0.107" cx="0" cy="0" r="1"/>
<use xlink:href="#gpPt5" id="gpPt6" style="fill:currentColor; stroke:none"/>
<path id="gpPt7" style="stroke-width:0.107" d="M0,-1.33 L-1.33,0.67 L1.33,0.67 z"/>
<use xlink:href="#gpPt7" id="gpPt8" style="fill:currentColor; stroke:none"/>
<use xlink:href="#gpPt7" id="gpPt9" transform="rotate(180)"/>
<use xlink:href="#gpPt9" id="gpPt10" style="fill:currentColor; stroke:none"/>
<use xlink:href="#gpPt3" id="gpPt11" transform="rotate(45)"/>
<use xlink:href="#gpPt11" id="gpPt12" style="fill:currentColor; stroke:none"/>
</defs>
<g style="fill:none; color:white; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
</g>
<g style="fill:none; color:black; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
<path d="M88.3,445.2 L107.0,445.2 "/>
<g transform="translate(70.9,454.6)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:end">
<text>-1</text>
</g>
<path d="M88.3,240.0 L107.0,240.0 "/>
<g transform="translate(70.9,249.4)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:end">
<text> 0</text>
</g>
<path d="M88.3,34.8 L107.0,34.8 "/>
<g transform="translate(70.9,44.2)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:end">
<text> 1</text>
</g>
<path d="M164.7,240.0 L164.7,258.7 M164.7,240.0 L164.7,221.3 "/>
<g transform="translate(164.7,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 1</text>
</g>
<path d="M241.0,240.0 L241.0,258.7 M241.0,240.0 L241.0,221.3 "/>
<g transform="translate(241.0,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 2</text>
</g>
<path d="M317.4,240.0 L317.4,258.7 M317.4,240.0 L317.4,221.3 "/>
<g transform="translate(317.4,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 3</text>
</g>
<path d="M393.8,240.0 L393.8,258.7 M393.8,240.0 L393.8,221.3 "/>
<g transform="translate(393.8,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 4</text>
</g>
<path d="M470.1,240.0 L470.1,258.7 M470.1,240.0 L470.1,221.3 "/>
<g transform="translate(470.1,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 5</text>
</g>
<path d="M546.5,240.0 L546.5,258.7 M546.5,240.0 L546.5,221.3 "/>
<g transform="translate(546.5,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 6</text>
</g>
</g>
<g style="fill:none; color:gray; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
<path d="M88.3,240.0 L546.5,240.0 "/>
</g>
<g style="fill:none; color:black; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
<path d="M88.3,34.8 L88.3,445.2 M546.5,445.2 M546.5,34.8 M88.3,34.8 "/>
</g>
<g style="fill:none; color:red; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
<use xlink:href="#gpPt6" transform="translate(88.3,240.0) scale(4.67)"/>
<use xlink:href="#gpPt6" transform="translate(164.7,67.3) scale(4.67)"/>
<use xlink:href="#gpPt6" transform="translate(241.0,53.4) scale(4.67)"/>
<use xlink:href="#gpPt6" transform="translate(317.4,211.0) scale(4.67)"/>
<use xlink:href="#gpPt6" transform="translate(393.8,395.3) scale(4.67)"/>
<use xlink:href="#gpPt6" transform="translate(470.1,436.8) scale(4.67)"/>
<use xlink:href="#gpPt6" transform="translate(546.5,297.3) scale(4.67)"/>
</g>
<g style="fill:none; color:black; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
<path d="M88.3,240.0 L92.9,240.0 L97.6,240.0 L102.2,240.0 L106.8,240.0 L111.4,240.0 L116.1,240.0 L120.7,240.0 L125.3,240.0 L130.0,240.0 L134.6,240.0 L139.2,240.0 L143.8,240.0 L148.5,240.0 L153.1,240.0 L157.7,240.0 L162.4,240.0 L167.0,240.0 L171.6,240.0 L176.2,240.0 L180.9,240.0 L185.5,240.0 L190.1,240.0 L194.8,240.0 L199.4,240.0 L204.0,240.0 L208.6,240.0 L213.3,240.0 L217.9,240.0 L222.5,240.0 L227.1,240.0 L231.8,240.0 L236.4,240.0 L241.0,240.0 L245.7,240.0 L250.3,240.0 L254.9,240.0 L259.5,240.0 L264.2,240.0 L268.8,240.0 L273.4,240.0 L278.1,240.0 L282.7,240.0 L287.3,240.0 L291.9,240.0 L296.6,240.0 L301.2,240.0 L305.8,240.0 L310.5,240.0 L315.1,240.0 L319.7,240.0 L324.3,240.0 L329.0,240.0 L333.6,240.0 L338.2,240.0 L342.9,240.0 L347.5,240.0 L352.1,240.0 L356.7,240.0 L361.4,240.0 L366.0,240.0 L370.6,240.0 L375.3,240.0 L379.9,240.0 L384.5,240.0 L389.1,240.0 L393.8,240.0 L398.4,240.0 L403.0,240.0 L407.7,240.0 L412.3,240.0 L416.9,240.0 L421.5,240.0 L426.2,240.0 L430.8,240.0 L435.4,240.0 L440.0,240.0 L444.7,240.0 L449.3,240.0 L453.9,240.0 L458.6,240.0 L463.2,240.0 L467.8,240.0 L472.4,240.0 L477.1,240.0 L481.7,240.0 L486.3,240.0 L491.0,240.0 L495.6,240.0 L500.2,240.0 L504.8,240.0 L509.5,240.0 L514.1,240.0 L518.7,240.0 L523.4,240.0 L528.0,240.0 L532.6,240.0 L537.2,240.0 L541.9,240.0 L546.5,240.0 M88.3,34.8 L88.3,445.2 M546.5,445.2 M546.5,34.8 M88.3,34.8 "/>
</g>
<style xmlns="http://www.w3.org/1999/xhtml" id="ls008m4a.rnq">rt.katakana-terminator-rt::before { content: attr(data-rt); }</style></svg>

After

Width:  |  Height:  |  Size: 6.1 KiB

View File

@ -0,0 +1,89 @@
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
<svg width="600" height="480" viewBox="0 0 600 480" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">
<desc>Produced by GNUPLOT 4.2 patchlevel 0 </desc>
<defs>
<circle id="gpDot" r="1"/>
<path id="gpPt0" style="stroke-width:0.107" d="M-1,0 h2 M0,-1 v2"/>
<path id="gpPt1" style="stroke-width:0.107" d="M-1,-1 L1,1 M1,-1 L-1,1"/>
<path id="gpPt2" style="stroke-width:0.107" d="M-1,0 L1,0 M0,-1 L0,1 M-1,-1 L1,1 M-1,1 L1,-1"/>
<rect id="gpPt3" style="stroke-width:0.107" x="-1" y="-1" width="2" height="2"/>
<use xlink:href="#gpPt3" id="gpPt4" style="fill:currentColor; stroke:none"/>
<circle id="gpPt5" style="stroke-width:0.107" cx="0" cy="0" r="1"/>
<use xlink:href="#gpPt5" id="gpPt6" style="fill:currentColor; stroke:none"/>
<path id="gpPt7" style="stroke-width:0.107" d="M0,-1.33 L-1.33,0.67 L1.33,0.67 z"/>
<use xlink:href="#gpPt7" id="gpPt8" style="fill:currentColor; stroke:none"/>
<use xlink:href="#gpPt7" id="gpPt9" transform="rotate(180)"/>
<use xlink:href="#gpPt9" id="gpPt10" style="fill:currentColor; stroke:none"/>
<use xlink:href="#gpPt3" id="gpPt11" transform="rotate(45)"/>
<use xlink:href="#gpPt11" id="gpPt12" style="fill:currentColor; stroke:none"/>
</defs>
<g style="fill:none; color:white; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
</g>
<g style="fill:none; color:black; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
<path d="M88.3,445.2 L107.0,445.2 "/>
<g transform="translate(70.9,454.6)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:end">
<text>-1</text>
</g>
<path d="M88.3,240.0 L107.0,240.0 "/>
<g transform="translate(70.9,249.4)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:end">
<text> 0</text>
</g>
<path d="M88.3,34.8 L107.0,34.8 "/>
<g transform="translate(70.9,44.2)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:end">
<text> 1</text>
</g>
<path d="M164.7,240.0 L164.7,258.7 M164.7,240.0 L164.7,221.3 "/>
<g transform="translate(164.7,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 1</text>
</g>
<path d="M241.0,240.0 L241.0,258.7 M241.0,240.0 L241.0,221.3 "/>
<g transform="translate(241.0,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 2</text>
</g>
<path d="M317.4,240.0 L317.4,258.7 M317.4,240.0 L317.4,221.3 "/>
<g transform="translate(317.4,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 3</text>
</g>
<path d="M393.8,240.0 L393.8,258.7 M393.8,240.0 L393.8,221.3 "/>
<g transform="translate(393.8,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 4</text>
</g>
<path d="M470.1,240.0 L470.1,258.7 M470.1,240.0 L470.1,221.3 "/>
<g transform="translate(470.1,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 5</text>
</g>
<path d="M546.5,240.0 L546.5,258.7 M546.5,240.0 L546.5,221.3 "/>
<g transform="translate(546.5,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 6</text>
</g>
</g>
<g style="fill:none; color:gray; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
<path d="M88.3,240.0 L546.5,240.0 "/>
</g>
<g style="fill:none; color:black; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
<path d="M88.3,34.8 L88.3,445.2 M546.5,445.2 M546.5,34.8 M88.3,34.8 "/>
</g>
<g style="fill:none; color:black; stroke:currentColor; stroke-width:4.00; stroke-linecap:butt; stroke-linejoin:miter">
</g>
<g style="fill:none; color:blue; stroke:currentColor; stroke-width:4.00; stroke-linecap:butt; stroke-linejoin:miter">
<path d="M88.3,240.0 L164.7,67.3 L241.0,53.4 L317.4,211.0 L393.8,395.3 L470.1,436.8 L546.5,297.3 "/>
</g>
<g style="fill:none; color:blue; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
</g>
<g style="fill:none; color:red; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
<use xlink:href="#gpPt6" transform="translate(88.3,240.0) scale(6.54)"/>
<use xlink:href="#gpPt6" transform="translate(164.7,67.3) scale(6.54)"/>
<use xlink:href="#gpPt6" transform="translate(241.0,53.4) scale(6.54)"/>
<use xlink:href="#gpPt6" transform="translate(317.4,211.0) scale(6.54)"/>
<use xlink:href="#gpPt6" transform="translate(393.8,395.3) scale(6.54)"/>
<use xlink:href="#gpPt6" transform="translate(470.1,436.8) scale(6.54)"/>
<use xlink:href="#gpPt6" transform="translate(546.5,297.3) scale(6.54)"/>
</g>
<g style="fill:none; color:black; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
<path d="M88.3,240.0 L546.5,240.0 M88.3,34.8 L88.3,445.2 M546.5,445.2 M546.5,34.8 M88.3,34.8 "/>
</g>
<style xmlns="http://www.w3.org/1999/xhtml" id="ls07lbzl.mze">rt.katakana-terminator-rt::before { content: attr(data-rt); }</style></svg>

After

Width:  |  Height:  |  Size: 5.3 KiB

View File

@ -0,0 +1,89 @@
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
<svg width="600" height="480" viewBox="0 0 600 480" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">
<desc>Produced by GNUPLOT 4.2 patchlevel 0 </desc>
<defs>
<circle id="gpDot" r="1"/>
<path id="gpPt0" style="stroke-width:0.107" d="M-1,0 h2 M0,-1 v2"/>
<path id="gpPt1" style="stroke-width:0.107" d="M-1,-1 L1,1 M1,-1 L-1,1"/>
<path id="gpPt2" style="stroke-width:0.107" d="M-1,0 L1,0 M0,-1 L0,1 M-1,-1 L1,1 M-1,1 L1,-1"/>
<rect id="gpPt3" style="stroke-width:0.107" x="-1" y="-1" width="2" height="2"/>
<use xlink:href="#gpPt3" id="gpPt4" style="fill:currentColor; stroke:none"/>
<circle id="gpPt5" style="stroke-width:0.107" cx="0" cy="0" r="1"/>
<use xlink:href="#gpPt5" id="gpPt6" style="fill:currentColor; stroke:none"/>
<path id="gpPt7" style="stroke-width:0.107" d="M0,-1.33 L-1.33,0.67 L1.33,0.67 z"/>
<use xlink:href="#gpPt7" id="gpPt8" style="fill:currentColor; stroke:none"/>
<use xlink:href="#gpPt7" id="gpPt9" transform="rotate(180)"/>
<use xlink:href="#gpPt9" id="gpPt10" style="fill:currentColor; stroke:none"/>
<use xlink:href="#gpPt3" id="gpPt11" transform="rotate(45)"/>
<use xlink:href="#gpPt11" id="gpPt12" style="fill:currentColor; stroke:none"/>
</defs>
<g style="fill:none; color:white; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
</g>
<g style="fill:none; color:black; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
<path d="M88.3,426.5 L107.0,426.5 "/>
<g transform="translate(70.9,435.9)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:end">
<text>-1</text>
</g>
<path d="M88.3,240.0 L107.0,240.0 "/>
<g transform="translate(70.9,249.4)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:end">
<text> 0</text>
</g>
<path d="M88.3,53.5 L107.0,53.5 "/>
<g transform="translate(70.9,62.9)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:end">
<text> 1</text>
</g>
<path d="M164.7,240.0 L164.7,258.7 M164.7,240.0 L164.7,221.3 "/>
<g transform="translate(164.7,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 1</text>
</g>
<path d="M241.0,240.0 L241.0,258.7 M241.0,240.0 L241.0,221.3 "/>
<g transform="translate(241.0,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 2</text>
</g>
<path d="M317.4,240.0 L317.4,258.7 M317.4,240.0 L317.4,221.3 "/>
<g transform="translate(317.4,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 3</text>
</g>
<path d="M393.8,240.0 L393.8,258.7 M393.8,240.0 L393.8,221.3 "/>
<g transform="translate(393.8,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 4</text>
</g>
<path d="M470.1,240.0 L470.1,258.7 M470.1,240.0 L470.1,221.3 "/>
<g transform="translate(470.1,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 5</text>
</g>
<path d="M546.5,240.0 L546.5,258.7 M546.5,240.0 L546.5,221.3 "/>
<g transform="translate(546.5,305.5)" style="stroke:none; fill:black; font-family:Vera; font-size:25.00; text-anchor:middle">
<text> 6</text>
</g>
</g>
<g style="fill:none; color:gray; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
<path d="M88.3,240.0 L546.5,240.0 "/>
</g>
<g style="fill:none; color:black; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
<path d="M88.3,34.8 L88.3,445.2 M546.5,445.2 M546.5,34.8 M88.3,34.8 "/>
</g>
<g style="fill:none; color:black; stroke:currentColor; stroke-width:4.00; stroke-linecap:butt; stroke-linejoin:miter">
</g>
<g style="fill:none; color:blue; stroke:currentColor; stroke-width:4.00; stroke-linecap:butt; stroke-linejoin:miter">
<path d="M88.3,240.0 L89.8,236.6 L91.4,233.2 L92.9,229.7 L94.4,226.2 L96.0,222.7 L97.5,219.2 L99.0,215.7 L100.6,212.1 L102.1,208.6 L103.6,205.0 L105.2,201.4 L106.7,197.8 L108.2,194.3 L109.8,190.7 L111.3,187.1 L112.8,183.6 L114.4,180.1 L115.9,176.5 L117.4,173.0 L118.9,169.5 L120.5,166.1 L122.0,162.6 L123.5,159.2 L125.1,155.8 L126.6,152.4 L128.1,149.1 L129.7,145.8 L131.2,142.5 L132.7,139.3 L134.3,136.1 L135.8,132.9 L137.3,129.8 L138.9,126.7 L140.4,123.7 L141.9,120.7 L143.5,117.8 L145.0,114.9 L146.5,112.1 L148.1,109.3 L149.6,106.6 L151.1,103.9 L152.7,101.3 L154.2,98.8 L155.7,96.3 L157.3,93.9 L158.8,91.5 L160.3,89.2 L161.9,87.0 L163.4,84.8 L164.9,82.7 L166.5,80.6 L168.0,78.6 L169.5,76.7 L171.1,74.9 L172.6,73.1 L174.1,71.4 L175.6,69.8 L177.2,68.2 L178.7,66.8 L180.2,65.4 L181.8,64.0 L183.3,62.8 L184.8,61.6 L186.4,60.4 L187.9,59.4 L189.4,58.4 L191.0,57.6 L192.5,56.7 L194.0,56.0 L195.6,55.4 L197.1,54.8 L198.6,54.3 L200.2,53.8 L201.7,53.5 L203.2,53.2 L204.8,53.0 L206.3,52.9 L207.8,52.9 L209.4,52.9 L210.9,53.0 L212.4,53.2 L214.0,53.5 L215.5,53.8 L217.0,54.2 L218.6,54.7 L220.1,55.3 L221.6,55.9 L223.2,56.6 L224.7,57.4 L226.2,58.3 L227.8,59.2 L229.3,60.2 L230.8,61.3 L232.3,62.5 L233.9,63.7 L235.4,65.0 L236.9,66.4 L238.5,67.8 L240.0,69.3 L241.5,70.9 L243.1,72.5 L244.6,74.2 L246.1,76.0 L247.7,77.8 L249.2,79.7 L250.7,81.7 L252.3,83.7 L253.8,85.8 L255.3,88.0 L256.9,90.2 L258.4,92.5 L259.9,94.8 L261.5,97.2 L263.0,99.6 L264.5,102.1 L266.1,104.7 L267.6,107.3 L269.1,109.9 L270.7,112.7 L272.2,115.4 L273.7,118.2 L275.3,121.1 L276.8,124.0 L278.3,126.9 L279.9,129.9 L281.4,133.0 L282.9,136.0 L284.5,139.2 L286.0,142.3 L287.5,145.5 L289.0,148.7 L290.6,152.0 L292.1,155.3 L293.6,158.7 L295.2,162.0 L296.7,165.4 L298.2,168.8 L299.8,172.3 L301.3,175.8 L302.8,179.3 L304.4,182.8 L305.9,186.4 L307.4,189.9 L309.0,193.5 L310.5,197.2 L312.0,200.8 L313.6,204.4 L315.1,208.1 L316.6,211.8 L318.2,215.4 L319.7,219.1 L321.2,222.8 L322.8,226.5 L324.3,230.3 L325.8,234.0 L327.4,237.7 L328.9,241.4 L330.4,245.1 L332.0,248.9 L333.5,252.6 L335.0,256.3 L336.6,260.0 L338.1,263.7 L339.6,267.4 L341.2,271.1 L342.7,274.7 L344.2,278.4 L345.8,282.0 L347.3,285.7 L348.8,289.3 L350.3,292.8 L351.9,296.4 L353.4,300.0 L354.9,303.5 L356.5,307.0 L358.0,310.5 L359.5,313.9 L361.1,317.3 L362.6,320.7 L364.1,324.1 L365.7,327.4 L367.2,330.7 L368.7,333.9 L370.3,337.2 L371.8,340.3 L373.3,343.5 L374.9,346.6 L376.4,349.6 L377.9,352.7 L379.5,355.6 L381.0,358.6 L382.5,361.4 L384.1,364.3 L385.6,367.1 L387.1,369.8 L388.7,372.5 L390.2,375.1 L391.7,377.7 L393.3,380.2 L394.8,382.7 L396.3,385.1 L397.9,387.4 L399.4,389.7 L400.9,391.9 L402.5,394.1 L404.0,396.2 L405.5,398.3 L407.0,400.2 L408.6,402.2 L410.1,404.0 L411.6,405.8 L413.2,407.5 L414.7,409.2 L416.2,410.7 L417.8,412.3 L419.3,413.7 L420.8,415.1 L422.4,416.4 L423.9,417.6 L425.4,418.7 L427.0,419.8 L428.5,420.8 L430.0,421.8 L431.6,422.6 L433.1,423.4 L434.6,424.1 L436.2,424.7 L437.7,425.3 L439.2,425.7 L440.8,426.1 L442.3,426.5 L443.8,426.7 L445.4,426.9 L446.9,426.9 L448.4,427.0 L450.0,426.9 L451.5,426.7 L453.0,426.5 L454.6,426.2 L456.1,425.8 L457.6,425.3 L459.2,424.8 L460.7,424.2 L462.2,423.5 L463.7,422.7 L465.3,421.8 L466.8,420.9 L468.3,419.9 L469.9,418.8 L471.4,417.7 L472.9,416.4 L474.5,415.1 L476.0,413.7 L477.5,412.3 L479.1,410.7 L480.6,409.1 L482.1,407.5 L483.7,405.7 L485.2,403.9 L486.7,402.0 L488.3,400.1 L489.8,398.1 L491.3,396.0 L492.9,393.9 L494.4,391.6 L495.9,389.4 L497.5,387.1 L499.0,384.7 L500.5,382.2 L502.1,379.7 L503.6,377.2 L505.1,374.6 L506.7,371.9 L508.2,369.2 L509.7,366.4 L511.3,363.6 L512.8,360.8 L514.3,357.9 L515.9,355.0 L517.4,352.0 L518.9,349.0 L520.4,346.0 L522.0,342.9 L523.5,339.8 L525.0,336.7 L526.6,333.6 L528.1,330.4 L529.6,327.2 L531.2,324.0 L532.7,320.8 L534.2,317.6 L535.8,314.4 L537.3,311.1 L538.8,307.9 L540.4,304.6 L541.9,301.4 L543.4,298.2 L545.0,295.0 L546.5,291.8 "/>
</g>
<g style="fill:none; color:blue; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
</g>
<g style="fill:none; color:red; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
<use xlink:href="#gpPt6" transform="translate(88.3,240.0) scale(6.54)"/>
<use xlink:href="#gpPt6" transform="translate(164.7,83.0) scale(6.54)"/>
<use xlink:href="#gpPt6" transform="translate(241.0,70.4) scale(6.54)"/>
<use xlink:href="#gpPt6" transform="translate(317.4,213.7) scale(6.54)"/>
<use xlink:href="#gpPt6" transform="translate(393.8,381.2) scale(6.54)"/>
<use xlink:href="#gpPt6" transform="translate(470.1,418.9) scale(6.54)"/>
<use xlink:href="#gpPt6" transform="translate(546.5,292.1) scale(6.54)"/>
</g>
<g style="fill:none; color:black; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter">
<path d="M88.3,240.0 L89.8,240.0 L91.4,240.0 L92.9,240.0 L94.4,240.0 L96.0,240.0 L97.5,240.0 L99.0,240.0 L100.6,240.0 L102.1,240.0 L103.6,240.0 L105.2,240.0 L106.7,240.0 L108.2,240.0 L109.8,240.0 L111.3,240.0 L112.8,240.0 L114.4,240.0 L115.9,240.0 L117.4,240.0 L118.9,240.0 L120.5,240.0 L122.0,240.0 L123.5,240.0 L125.1,240.0 L126.6,240.0 L128.1,240.0 L129.7,240.0 L131.2,240.0 L132.7,240.0 L134.3,240.0 L135.8,240.0 L137.3,240.0 L138.9,240.0 L140.4,240.0 L141.9,240.0 L143.5,240.0 L145.0,240.0 L146.5,240.0 L148.1,240.0 L149.6,240.0 L151.1,240.0 L152.7,240.0 L154.2,240.0 L155.7,240.0 L157.3,240.0 L158.8,240.0 L160.3,240.0 L161.9,240.0 L163.4,240.0 L164.9,240.0 L166.5,240.0 L168.0,240.0 L169.5,240.0 L171.1,240.0 L172.6,240.0 L174.1,240.0 L175.6,240.0 L177.2,240.0 L178.7,240.0 L180.2,240.0 L181.8,240.0 L183.3,240.0 L184.8,240.0 L186.4,240.0 L187.9,240.0 L189.4,240.0 L191.0,240.0 L192.5,240.0 L194.0,240.0 L195.6,240.0 L197.1,240.0 L198.6,240.0 L200.2,240.0 L201.7,240.0 L203.2,240.0 L204.8,240.0 L206.3,240.0 L207.8,240.0 L209.4,240.0 L210.9,240.0 L212.4,240.0 L214.0,240.0 L215.5,240.0 L217.0,240.0 L218.6,240.0 L220.1,240.0 L221.6,240.0 L223.2,240.0 L224.7,240.0 L226.2,240.0 L227.8,240.0 L229.3,240.0 L230.8,240.0 L232.3,240.0 L233.9,240.0 L235.4,240.0 L236.9,240.0 L238.5,240.0 L240.0,240.0 L241.5,240.0 L243.1,240.0 L244.6,240.0 L246.1,240.0 L247.7,240.0 L249.2,240.0 L250.7,240.0 L252.3,240.0 L253.8,240.0 L255.3,240.0 L256.9,240.0 L258.4,240.0 L259.9,240.0 L261.5,240.0 L263.0,240.0 L264.5,240.0 L266.1,240.0 L267.6,240.0 L269.1,240.0 L270.7,240.0 L272.2,240.0 L273.7,240.0 L275.3,240.0 L276.8,240.0 L278.3,240.0 L279.9,240.0 L281.4,240.0 L282.9,240.0 L284.5,240.0 L286.0,240.0 L287.5,240.0 L289.0,240.0 L290.6,240.0 L292.1,240.0 L293.6,240.0 L295.2,240.0 L296.7,240.0 L298.2,240.0 L299.8,240.0 L301.3,240.0 L302.8,240.0 L304.4,240.0 L305.9,240.0 L307.4,240.0 L309.0,240.0 L310.5,240.0 L312.0,240.0 L313.6,240.0 L315.1,240.0 L316.6,240.0 L318.2,240.0 L319.7,240.0 L321.2,240.0 L322.8,240.0 L324.3,240.0 L325.8,240.0 L327.4,240.0 L328.9,240.0 L330.4,240.0 L332.0,240.0 L333.5,240.0 L335.0,240.0 L336.6,240.0 L338.1,240.0 L339.6,240.0 L341.2,240.0 L342.7,240.0 L344.2,240.0 L345.8,240.0 L347.3,240.0 L348.8,240.0 L350.3,240.0 L351.9,240.0 L353.4,240.0 L354.9,240.0 L356.5,240.0 L358.0,240.0 L359.5,240.0 L361.1,240.0 L362.6,240.0 L364.1,240.0 L365.7,240.0 L367.2,240.0 L368.7,240.0 L370.3,240.0 L371.8,240.0 L373.3,240.0 L374.9,240.0 L376.4,240.0 L377.9,240.0 L379.5,240.0 L381.0,240.0 L382.5,240.0 L384.1,240.0 L385.6,240.0 L387.1,240.0 L388.7,240.0 L390.2,240.0 L391.7,240.0 L393.3,240.0 L394.8,240.0 L396.3,240.0 L397.9,240.0 L399.4,240.0 L400.9,240.0 L402.5,240.0 L404.0,240.0 L405.5,240.0 L407.0,240.0 L408.6,240.0 L410.1,240.0 L411.6,240.0 L413.2,240.0 L414.7,240.0 L416.2,240.0 L417.8,240.0 L419.3,240.0 L420.8,240.0 L422.4,240.0 L423.9,240.0 L425.4,240.0 L427.0,240.0 L428.5,240.0 L430.0,240.0 L431.6,240.0 L433.1,240.0 L434.6,240.0 L436.2,240.0 L437.7,240.0 L439.2,240.0 L440.8,240.0 L442.3,240.0 L443.8,240.0 L445.4,240.0 L446.9,240.0 L448.4,240.0 L450.0,240.0 L451.5,240.0 L453.0,240.0 L454.6,240.0 L456.1,240.0 L457.6,240.0 L459.2,240.0 L460.7,240.0 L462.2,240.0 L463.7,240.0 L465.3,240.0 L466.8,240.0 L468.3,240.0 L469.9,240.0 L471.4,240.0 L472.9,240.0 L474.5,240.0 L476.0,240.0 L477.5,240.0 L479.1,240.0 L480.6,240.0 L482.1,240.0 L483.7,240.0 L485.2,240.0 L486.7,240.0 L488.3,240.0 L489.8,240.0 L491.3,240.0 L492.9,240.0 L494.4,240.0 L495.9,240.0 L497.5,240.0 L499.0,240.0 L500.5,240.0 L502.1,240.0 L503.6,240.0 L505.1,240.0 L506.7,240.0 L508.2,240.0 L509.7,240.0 L511.3,240.0 L512.8,240.0 L514.3,240.0 L515.9,240.0 L517.4,240.0 L518.9,240.0 L520.4,240.0 L522.0,240.0 L523.5,240.0 L525.0,240.0 L526.6,240.0 L528.1,240.0 L529.6,240.0 L531.2,240.0 L532.7,240.0 L534.2,240.0 L535.8,240.0 L537.3,240.0 L538.8,240.0 L540.4,240.0 L541.9,240.0 L543.4,240.0 L545.0,240.0 L546.5,240.0 M88.3,34.8 L88.3,445.2 M546.5,445.2 M546.5,34.8 M88.3,34.8 "/>
</g>
<style xmlns="http://www.w3.org/1999/xhtml" id="ls082ltu.w37">rt.katakana-terminator-rt::before { content: attr(data-rt); }</style></svg>

After

Width:  |  Height:  |  Size: 13 KiB

View File

@ -0,0 +1,220 @@
author: AtomAlpaca, billchenchina, Chrogeek, Early0v0, EndlessCheng, Enter-tainer, Henry-ZHR, hly1204, hsfzLZH1, Ir1d, Ghastlcon, kenlig, Marcythm, megakite, Peanut-Tang, qwqAutomaton, qz-cqy, StudyingFather, swift-zym, swiftqwq, Tiphereth-A, TrisolarisHD, x4Cx58x54, Xeonacid, xiaopangfeiyu, YanWQ-monad
## 引入
插值是一种通过已知的、离散的数据点推算一定范围内的新数据点的方法。插值法常用于函数拟合中。
例如对数据点:
| $x$ | $0$ | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ |
| ------ | --- | -------- | -------- | -------- | --------- | --------- | --------- |
| $f(x)$ | $0$ | $0.8415$ | $0.9093$ | $0.1411$ | $-0.7568$ | $-0.9589$ | $-0.2794$ |
![](../images/interp-1.svg)
其中 $f(x)$ 未知,插值法可以通过按一定形式拟合 $f(x)$ 的方式估算未知的数据点。
例如,我们可以用分段线性函数拟合 $f(x)$
![](../images/interp-2.svg)
这种插值方式叫做 [线性插值](https://en.wikipedia.org/wiki/Linear_interpolation)。
我们也可以用多项式拟合 $f(x)$
![](../images/interp-3.svg)
这种插值方式叫做 [多项式插值](https://en.wikipedia.org/wiki/Polynomial_interpolation)。
多项式插值的一般形式如下:
???+ note "多项式插值"
对已知的 $n+1$ 的点 $(x_0,y_0),(x_1,y_1),\dots,(x_n,y_n)$,求 $n$ 阶多项式 $f(x)$ 满足
$$
f(x_i)=y_i,\qquad\forall i=0,1,\dots,n
$$
下面介绍多项式插值中的两种方式Lagrange 插值法与 Newton 插值法。不难证明这两种方法得到的结果是相等的。
## Lagrange 插值法
由于要求构造一个函数 $f(x)$ 过点 $P_1(x_1, y_1), P_2(x_2,y_2),\cdots,P_n(x_n,y_n)$. 首先设第 $i$ 个点在 $x$ 轴上的投影为 $P_i^{\prime}(x_i,0)$.
考虑构造 $n$ 个函数 $f_1(x), f_2(x), \cdots, f_n(x)$,使得对于第 $i$ 个函数 $f_i(x)$,其图像过 $\begin{cases}P_j^{\prime}(x_j,0),(j\neq i)\\P_i(x_i,y_i)\end{cases}$,则可知题目所求的函数 $f(x)=\sum\limits_{i=1}^nf_i(x)$.
那么可以设 $f_i(x)=a\cdot\prod_{j\neq i}(x-x_j)$,将点 $P_i(x_i,y_i)$ 代入可以知道 $a=\dfrac{y_i}{\prod_{j\neq i} (x_i-x_j)}$,所以
$$
f_i(x)=y_i\cdot\dfrac{\prod_{j\neq i} (x-x_j)}{\prod_{j\neq i} (x_i-x_j)}=y_i\cdot\prod_{j\neq i}\dfrac{x-x_j}{x_i-x_j}
$$
那么我们就可以得出 Lagrange 插值的形式为:
$$
f(x)=\sum_{i=1}^ny_i\cdot\prod_{j\neq i}\dfrac{x-x_j}{x_i-x_j}
$$
朴素实现的时间复杂度为 $O(n^2)$,可以优化到 $O(n\log^2 n)$,参见 [多项式快速插值](../poly/multipoint-eval-interpolation.md#多项式的快速插值)。
???+ note "[Luogu P4781【模板】拉格朗日插值](https://www.luogu.com.cn/problem/P4781)"
给出 $n$ 个点对 $(x_i,y_i)$ 和 $k$,且 $\forall i,j$ 有 $i\neq j \iff x_i\neq x_j$ 且 $f(x_i)\equiv y_i\pmod{998244353}$ 和 $\deg(f(x))<n$(定义 $\deg(0)=-\infty$),求 $f(k)\bmod{998244353}$.
??? note "题解"
本题中只用求出 $f(k)$ 的值,所以在计算上式的过程中直接将 $k$ 代入即可。
$$
f(k)=\sum_{i=1}^{n}y_i\prod_{j\neq i }\frac{k-x_j}{x_i-x_j}
$$
本题中,还需要求解逆元。如果先分别计算出分子和分母,再将分子乘进分母的逆元,累加进最后的答案,时间复杂度的瓶颈就不会在求逆元上,时间复杂度为 $O(n^2)$.
因为在固定模 $998244353$ 意义下运算,计算乘法逆元的时间复杂度我们在这里暂且认为是常数时间。
??? note "代码实现"
```cpp
--8<-- "docs/math/code/numerical/interp/interp_1.cpp"
```
### 横坐标是连续整数的 Lagrange 插值
如果已知点的横坐标是连续整数,我们可以做到 $O(n)$ 插值。
设要求 $n$ 次多项式为 $f(x)$,我们已知 $f(1),\cdots,f(n+1)$$1\le i\le n+1$),考虑代入上面的插值公式:
$$
\begin{aligned}
f(x)&=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j\ne i}\frac{x-x_j}{x_i-x_j}\\
&=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j\ne i}\frac{x-j}{i-j}
\end{aligned}
$$
后面的累乘可以分子分母分别考虑,不难得到分子为:
$$
\dfrac{\prod\limits_{j=1}^{n+1}(x-j)}{x-i}
$$
分母的 $i-j$ 累乘可以拆成两段阶乘来算:
$$
(-1)^{n+1-i}\cdot(i-1)!\cdot(n+1-i)!
$$
于是横坐标为 $1,\cdots,n+1$ 的插值公式:
$$
f(x)=\sum\limits_{i=1}^{n+1}(-1)^{n+1-i}y_i\cdot\frac{\prod\limits_{j=1}^{n+1}(x-j)}{(i-1)!(n+1-i)!(x-i)}
$$
预处理 $(x-i)$ 前后缀积、阶乘阶乘逆,然后代入这个式子,复杂度为 $O(n)$.
???+ note " 例题 [CF622F The Sum of the k-th Powers](https://codeforces.com/contest/622/problem/F)"
给出 $n,k$,求 $\sum\limits_{i=1}^ni^k$ 对 $10^9+7$ 取模的值。
??? note "题解"
本题中,答案是一个 $k+1$ 次多项式,因此我们可以线性筛出 $1^i,\cdots,(k+2)^i$ 的值然后进行 $O(n)$ 插值。
也可以通过组合数学相关知识由差分法的公式推得下式:
$$
f(x)=\sum_{i=1}^{n+1}\binom{x-1}{i-1}\sum_{j=1}^{i}(-1)^{i+j}\binom{i-1}{j-1}y_{j}=\sum\limits_{i=1}^{n+1}y_i\cdot\frac{\prod\limits_{j=1}^{n+1}(x-j)}{(x-i)\cdot(-1)^{n+1-i}\cdot(i-1)!\cdot(n+1-i)!}
$$
??? note "代码实现"
```cpp
--8<-- "docs/math/code/numerical/interp/interp_2.cpp"
```
## Newton 插值法
Newton 插值法是基于高阶差分来插值的方法,优点是支持 $O(n)$ 插入新数据点。
为了实现 $O(n)$ 插入新数据点,我们令:
$$
f(x)=\sum_{j=0}^n a_jn_j(x)
$$
其中 $n_j(x):=\prod_{i=0}^{j-1}(x-x_i)$ 称为 **Newton 基**Newton basis
若解出 $a_j$,则可得到 $f(x)$ 的插值多项式。我们按如下方式定义 **前向差商**forward divided differences
$$
\begin{aligned}
\lbrack y_k\rbrack & := y_k, & k=0,\dots,n, \\
[y_k,\dots,y_{k+j}] & := \dfrac{[y_{k+1},\dots,y_{k+j}]-[y_k,\dots,y_{k+j-1}]}{x_{k+j}-x_k}, & k=0,\dots,n-j,~j=1,\dots,n.
\end{aligned}
$$
则:
$$
\begin{aligned}
f(x)&=[y_0]+[y_0,y_1](x-x_0)+\dots+[y_0,\dots,y_n](x-x_0)\dots(x-x_{n-1})\\
&=\sum_{j=0}^n [y_0,\dots,y_j]n_j(x)
\end{aligned}
$$
此即 Newton 插值的形式。朴素实现的时间复杂度为 $O(n^2)$.
若样本点是等距的(即 $x_i=x_0+ih$$i=1,\dots,n$),令 $x=x_0+sh$Newton 插值的公式可化为:
$$
f(x)=\sum_{j=0}^n \binom{s}{j}j!h^j[y_0,\dots,y_j]
$$
上式称为 **Newton 前向差分公式**Newton forward divided difference formula
???+ note
若样本点是等距的,我们还可以推出:
$$
[y_k,\dots,y_{k+j}]=\frac{1}{j!h^j}\Delta^{(j)}y_k
$$
其中 $\Delta^{(j)}y_k$ 为 **前向差分**forward differences定义如下
$$
\begin{aligned}
\Delta^{(0)}y_k & := y_k, & k=0,\dots,n, \\
\Delta^{(j)}y_k & := \Delta^{(j-1)} y_{k+1}-\Delta^{(j-1)} y_k, & k=0,\dots,n-j,~j=1,\dots,n.
\end{aligned}
$$
??? note " 代码实现([Luogu P4781【模板】拉格朗日插值](https://www.luogu.com.cn/problem/P4781)"
```cpp
--8<-- "docs/math/code/numerical/interp/interp_3.cpp"
```
### 横坐标是连续整数的 Newton 插值
例如:求某三次多项式 $f(x)=\sum_{i=0}^{3} a_ix^i$ 的多项式系数,已知 $f(1)$ 至 $f(6)$ 的值分别为 $1, 5, 14, 30, 55, 91$。
$$
\begin{array}{cccccccccccc}
1 & & 5 & & 14 & & 30 & & 55 & & 91 & \\
& 4 & & 9 & & 16 & & 25 & & 36 & \\
& & 5 & & 7 & & 9 & & 11 & \\
& & & 2 & & 2 & & 2 & \\
\end{array}
$$
第一行为 $f(x)$ 的连续的前 $n$ 项;之后的每一行为之前一行中对应的相邻两项之差。观察到,如果这样操作的次数足够多(前提是 $f(x)$ 为多项式),最终总会返回一个定值。
计算出第 $i-1$ 阶差分的首项为 $\sum_{j=1}^{i}(-1)^{i+j}\binom{i-1}{j-1}f(j)$,第 $i-1$ 阶差分的首项对 $f(k)$ 的贡献为 $\binom{k-1}{i-1}$ 次。
$$
f(k)=\sum_{i=1}^n\binom{k-1}{i-1}\sum_{j=1}^{i}(-1)^{i+j}\binom{i-1}{j-1}f(j)
$$
时间复杂度为 $O(n^2)$.
## C++ 中的实现
自 C++ 20 起,标准库添加了 [`std::midpoint`](https://en.cppreference.com/w/cpp/numeric/midpoint) 和 [`std::lerp`](https://en.cppreference.com/w/cpp/numeric/lerp) 函数,分别用于求中点和线性插值。
## 参考资料
1. [Interpolation - Wikipedia](https://en.wikipedia.org/wiki/Interpolation)
2. [Newton polynomial - Wikipedia](https://en.wikipedia.org/wiki/Newton_polynomial)

View File

@ -1,169 +0,0 @@
author: Ir1d, Marcythm, x4Cx58x54, YanWQ-monad, AtomAlpaca, billchenchina, Chrogeek, Early0v0, EndlessCheng, Enter-tainer, Ghastlcon, Henry-ZHR, hly1204, hsfzLZH1, kenlig, Peanut-Tang, qwqAutomaton, qz-cqy, rui\_er, StudyingFather, swift-zym, Tiphereth-A, TrisolarisHD, Xeonacid
???+ note " 例题 [Luogu P4781【模板】拉格朗日插值](https://www.luogu.com.cn/problem/P4781)"
给出 $n$ 个点对 $(x_i,y_i)$ 和 $k$,且 $\forall i,j$ 有 $i\neq j \iff x_i\neq x_j$ 且 $f(x_i)\equiv y_i\pmod{998244353}$ 和 $\deg(f(x))<n$(定义 $\deg(0)=-\infty$),求 $f(k)\bmod{998244353}$。
### 方法 1差分法
差分法适用于 $x_i=i$ 的情况。
如,用差分法求某三次多项式 $f(x)=\sum_{i=0}^{3} a_ix^i$ 的多项式形式,已知 $f(1)$ 至 $f(6)$ 的值分别为 $1, 5, 14, 30, 55, 91$。
$$
\begin{array}{cccccccccccc}
1 & & 5 & & 14 & & 30 & & 55 & & 91 & \\
& 4 & & 9 & & 16 & & 25 & & 36 & \\
& & 5 & & 7 & & 9 & & 11 & \\
& & & 2 & & 2 & & 2 & \\
\end{array}
$$
第一行为 $f(x)$ 的连续的前 $n$ 项;之后的每一行为之前一行中对应的相邻两项之差。观察到,如果这样操作的次数足够多(前提是 $f(x)$ 为多项式),最终总会返回一个定值。
计算出第 $i-1$ 阶差分的首项为 $\sum_{j=1}^{i}(-1)^{i+j}\binom{i-1}{j-1}f(j)$,第 $i-1$ 阶差分的首项对 $f(k)$ 的贡献为 $\binom{k-1}{i-1}$ 次。
$$
f(k)=\sum_{i=1}^n\binom{k-1}{i-1}\sum_{j=1}^{i}(-1)^{i+j}\binom{i-1}{j-1}f(j)
$$
时间复杂度为 $O(n^2)$。这种方法对给出的点的限制性较强。
### 方法 2待定系数法
设 $f(x)=\sum_{i=0}^{n-1} a_ix^i$ 将每个 $x_i$ 代入 $f(x)$,有 $f(x_i)=y_i$,这样就可以得到一个由 $n$ 条 $n$ 元一次方程所组成的方程组,然后使用 **高斯消元** 解该方程组求出每一项 $a_i$,即确定了 $f(x)$ 的表达式。
如果您不知道什么是高斯消元,请看 [高斯消元](./gauss.md)。
时间复杂度 $O(n^3)$,对给出点的坐标无要求。
### 方法 3拉格朗日插值法
在 [多项式部分简介](../poly/intro.md) 里我们已经定义了多项式除法。
那么我们会有:
$$
f(x)\equiv f(a)\pmod{(x-a)}
$$
因为 $f(x)-f(a)=(a_0-a_0)+a_1(x^1-a^1)+a_1(x^2-a^2)+\cdots +a_n(x^n-a^n)$,显然有 $(x-a)$ 这个因式。
这样我们就可以列一个关于 $f(x)$ 的多项式线性同余方程组:
$$
\begin{cases}
f(x)\equiv y_1\pmod{(x-x_1)}\\
f(x)\equiv y_2\pmod{(x-x_2)}\\
\vdots\\
f(x)\equiv y_n\pmod{(x-x_n)}
\end{cases}
$$
$$
\begin{aligned}
M(x)&=\prod_{i=1}^n{(x-x_i)},\\
m_i(x)&=\dfrac M{x-x_i}
\end{aligned}
$$
则 $m_i(x)$ 在模 $(x-x_i)$ 意义下的乘法逆元为
$$
m_i(x_i)^{-1}=\prod_{j\ne i}{(x_i-x_j)^{-1}}
$$
$$
\begin{aligned}
f(x)&\equiv\sum_{i=1}^n{y_i\left(m_i(x)\right)\left(m_i(x_i)^{-1}\right)}&\pmod{M(x)}\\
&\equiv\sum_{i=1}^n{y_i\prod_{j\ne i}{\dfrac {x-x_j}{x_i-x_j}}}&\pmod{M(x)}
\end{aligned}
$$
又因为 $\deg\left(f(x)\right)<n$ 所以在模 $M(x)$ 意义下 $f(x)$ 就是唯一的,即:
$$
f(x)=\sum_{i=1}^n{y_i\prod_{j\ne i}{\dfrac {x-x_j}{x_i-x_j}}}
$$
这就是拉格朗日插值的表达式。
??? note "通常意义下拉格朗日插值的一种推导"
由于要求构造一个函数 $f(x)$ 过点 $P_1(x_1, y_1), P_2(x_2,y_2),\cdots,P_n(x_n,y_n)$。首先设第 $i$ 个点在 $x$ 轴上的投影为 $P_i^{\prime}(x_i,0)$。
考虑构造 $n$ 个函数 $f_1(x), f_2(x), \cdots, f_n(x)$,使得对于第 $i$ 个函数 $f_i(x)$,其图像过 $\begin{cases}P_j^{\prime}(x_j,0),(j\neq i)\\P_i(x_i,y_i)\end{cases}$,则可知题目所求的函数 $f(x)=\sum\limits_{i=1}^nf_i(x)$。
那么可以设 $f_i(x)=a\cdot\prod_{j\neq i}(x-x_j)$,将点 $P_i(x_i,y_i)$ 代入可以知道 $a=\dfrac{y_i}{\prod_{j\neq i} (x_i-x_j)}$,所以
$f_i(x)=y_i\cdot\dfrac{\prod_{j\neq i} (x-x_j)}{\prod_{j\neq i} (x_i-x_j)}=y_i\cdot\prod_{j\neq i}\dfrac{x-x_j}{x_i-x_j}$。
那么我们就可以从另一个角度推导出通常意义下(而非模意义下)拉格朗日插值的式子为:
$f(x)=\sum_{i=1}^ny_i\cdot\prod_{j\neq i}\dfrac{x-x_j}{x_i-x_j}$。
??? note "代码实现"
因为在固定模 $998244353$ 意义下运算,计算乘法逆元的时间复杂度我们在这里暂且认为是常数时间。
```cpp
--8<-- "docs/math/code/numerical/lagrange/lagrange_1.cpp"
```
本题中只用求出 $f(k)$ 的值,所以在计算上式的过程中直接将 $k$ 代入即可。
$$
f(k)=\sum_{i=1}^{n}y_i\prod_{j\neq i }\frac{k-x_j}{x_i-x_j}
$$
本题中,还需要求解逆元。如果先分别计算出分子和分母,再将分子乘进分母的逆元,累加进最后的答案,时间复杂度的瓶颈就不会在求逆元上,时间复杂度为 $O(n^2)$。
### 横坐标是连续整数的拉格朗日插值
如果已知点的横坐标是连续整数,我们可以做到 $O(n)$ 插值。
设要求 $n$ 次多项式为 $f(x)$,我们已知 $f(1),\cdots,f(n+1)$$1\le i\le n+1$),考虑代入上面的插值公式:
$$
\begin{aligned}
f(x)&=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j\ne i}\frac{x-x_j}{x_i-x_j}\\
&=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j\ne i}\frac{x-j}{i-j}
\end{aligned}
$$
后面的累乘可以分子分母分别考虑,不难得到分子为:
$$
\dfrac{\prod\limits_{j=1}^{n+1}(x-j)}{x-i}
$$
分母的 $i-j$ 累乘可以拆成两段阶乘来算:
$$
(-1)^{n+1-i}\cdot(i-1)!\cdot(n+1-i)!
$$
于是横坐标为 $1,\cdots,n+1$ 的插值公式:
$$
f(x)=\sum\limits_{i=1}^{n+1}y_i\cdot\frac{\prod\limits_{j=1}^{n+1}(x-j)}{(x-i)\cdot(-1)^{n+1-i}\cdot(i-1)!\cdot(n+1-i)!}
$$
预处理 $(x-i)$ 前后缀积、阶乘阶乘逆,然后代入这个式子,复杂度为 $O(n)$。
???+ note " 例题 [CF622F The Sum of the k-th Powers](https://codeforces.com/contest/622/problem/F)"
给出 $n,k$,求 $\sum\limits_{i=1}^ni^k$ 对 $10^9+7$ 取模的值。
本题中,答案是一个 $k+1$ 次多项式,因此我们可以线性筛出 $1^i,\cdots,(k+2)^i$ 的值然后进行 $O(n)$ 插值。
也可以通过组合数学相关知识由差分法的公式推得下式:
$$
f(x)=\sum_{i=1}^{n+1}\binom{x-1}{i-1}\sum_{j=1}^{i}(-1)^{i+j}\binom{i-1}{j-1}y_{j}=\sum\limits_{i=1}^{n+1}y_i\cdot\frac{\prod\limits_{j=1}^{n+1}(x-j)}{(x-i)\cdot(-1)^{n+1-i}\cdot(i-1)!\cdot(n+1-i)!}
$$
??? note "代码实现"
```cpp
--8<-- "docs/math/code/numerical/lagrange/lagrange_2.cpp"
```

View File

@ -89,7 +89,7 @@ $$
### Lagrange 插值公式法
考虑 Lagrange 插值公式
考虑 [Lagrange 插值公式](../numerical/interp.md#lagrange-插值法)
$$
\begin{aligned}

View File

@ -322,7 +322,7 @@ nav:
- 非公平组合游戏: math/game-theory/partizan-game.md
- 反常游戏: math/game-theory/misere-game.md
- 数值算法:
- 拉格朗日插值: math/numerical/lagrange.md
- 插值: math/numerical/interp.md
- 数值积分: math/numerical/integral.md
- 高斯消元: math/numerical/gauss.md
- 牛顿迭代法: math/numerical/newton.md