目錄

• @0 - 參考資料@
• @1 - 多點(diǎn)求值@
  • @理論推導(dǎo)@
  • @參考代碼@
  • @例題與應(yīng)用@
• @2 - 快速插值@
  • @理論推導(dǎo)@
  • @（不建議參考的）代碼@
  • @例題與應(yīng)用@（暫無）

---

@0 - 參考資料@

Cyhlnj 的博客

@1 - 多點(diǎn)求值@

@理論推導(dǎo)@

假設(shè)已知多項(xiàng)式 \(A(x)\)，使用 FFT 可以將 \(A(w_n^0)\)，\(A(w_n^1)\)，...，\(A(w_n^{n-1})\) 的值在 \(O(n\log n)\) 的時(shí)間內(nèi)快速求出。

那么問題來了，假如我現(xiàn)在要求解任意的 \(A(x_0)\)，\(A(x_1)\)，...，\(A(x_{n-1})\) 該怎么辦呢？難道只有 \(O(n^2)\) 的算法嗎？
當(dāng)然不是。

我們依照直覺把 \(x_{0\dots {n-1}}\) 分為兩個(gè)集合：
\[S_l=\{x_0,x_1,\dots x_{\frac{n}{2}}\}\\ S_r=\{x_{\frac{n}{2}+1},x_{\frac{n}{2}+2},\dots x_{n-1}\}\\\]
然后，記多項(xiàng)式 \(P_l(x)=\prod_{x_i\in S_L}(x-x_i)\)，并用 \(A(x)\) 除以 \(P_l(x)\)，得到：
\[A(x)=P_l(x)*Q(x)+R_l(x)\]
這樣當(dāng) \(x_i\in S_l\) 時(shí)，\(P_l(x_i)=0\)，\(A(x_i) = R_l(x_i)\)。
同理我們可以對(duì) \(S_r\) 進(jìn)行相應(yīng)的處理。

問題轉(zhuǎn)換為 \(S_l\) 對(duì) \(R_l(x)\) 求值，\(S_r\) 對(duì) \(R_r(x)\) 求值。
然后就可以分治了。到底層時(shí) \(R(x)\) 就是一個(gè)常數(shù)，可以直接賦值。

\(P_l(x)\) 與 \(P_r(x)\) 可以利用類線段樹的方法處理與儲(chǔ)存下來。
一個(gè)細(xì)節(jié)：如果 \(A(x)\) 的最高次數(shù)高于 \(\prod_{x=0}^{N-1}(x-x_i)\)，則先 \(A(x)\) 模一下 \(\prod_{x=0}^{N-1}(x-x_i)\)。

復(fù)雜度為 \(O(n\log^2n)\)。

@參考代碼@

本代碼為 luoguP5050 的 AC 代碼。

#include<cstdio> #include<vector> #include<algorithm> using namespace std; const int G = 3; const int MOD = 998244353; const int MAXN = 64000*10; int pow_mod(int b, int p) {int ret = 1;while( p ) {if( p & 1 ) ret = 1LL*ret*b%MOD;b = 1LL*b*b%MOD;p >>= 1;}return ret; } struct Polynomial{void poly_copy(int *A, int *B, int n) {for(int i=0;i<n;i++)A[i] = B[i];}void poly_clear(int *A, int l, int r) {for(int i=l;i<r;i++)A[i] = 0;}void poly_revcopy(int *A, int *B, int n) {for(int i=0;i<n;i++)A[i] = B[n-i-1];}void ntt(int *A, int n, int type) {for(int i=0,j=0;i<n;i++) {if( i < j ) swap(A[i], A[j]);for(int l=(n>>1);(j^=l)<l;l>>=1);}for(int s=2;s<=n;s<<=1) {int t = (s>>1);int u = (type == -1) ? pow_mod(G, (MOD-1) - (MOD-1)/s) : pow_mod(G, (MOD-1)/s);for(int i=0;i<n;i+=s) {for(int j=0,p=1;j<t;j++,p=1LL*p*u%MOD) {int x = A[i+j], y = 1LL*p*A[i+j+t]%MOD;A[i+j] = (x + y)%MOD, A[i+j+t] = (x + MOD - y)%MOD;}}}if( type == -1 ) {int inv = pow_mod(n, MOD-2);for(int i=0;i<n;i++)A[i] = 1LL*A[i]*inv%MOD;}}int tmp1[MAXN + 5], tmp2[MAXN + 5], tmp3[MAXN + 5];void poly_mul(int *A, int *B, int *C, int n, int m) {int len; for(len = 1;len < n+m-1;len <<= 1);poly_copy(tmp1, A, n); poly_clear(tmp1, n, len);poly_copy(tmp2, B, m); poly_clear(tmp2, m, len);ntt(tmp1, len, 1); ntt(tmp2, len, 1);for(int i=0;i<len;i++) tmp3[i] = 1LL*tmp1[i]*tmp2[i]%MOD;ntt(tmp3, len, -1); poly_copy(C, tmp3, n+m-1);}void poly_inv(int *A, int *B, int n) {if( n == 1 ) {B[0] = pow_mod(A[0], MOD-2);return ;}int len; for(len = 1;len < (n<<1);len <<= 1);poly_inv(A, B, (n + 1) >> 1);poly_copy(tmp3, A, n); poly_clear(tmp3, n, len);ntt(tmp3, len, 1); ntt(B, len, 1);for(int i=0;i<len;i++) B[i] = 1LL*B[i]*(2 + MOD - 1LL*tmp3[i]*B[i]%MOD)%MOD;ntt(B, len, -1); poly_clear(B, n, len);}int tmp4[MAXN + 5], tmp5[MAXN + 5], tmp6[MAXN + 5];void poly_divide(int *A, int *B, int *C, int *R, int n, int m) {poly_revcopy(tmp4, B, m); poly_clear(tmp5, 0, 2*(n-m+1)); poly_inv(tmp4, tmp5, n-m+1);poly_revcopy(tmp4, A, n); poly_mul(tmp4, tmp5, tmp6, n-m+1, n-m+1);poly_revcopy(C, tmp6, n-m+1);poly_copy(tmp4, C, n-m+1); poly_copy(tmp5, B, m);poly_mul(tmp4, tmp5, tmp6, n-m+1, m);for(int i=0;i<m-1;i++) R[i] = (A[i] + MOD - tmp6[i])%MOD;}void poly_mod(int *A, int *B, int *R, int n, int m) {if( n < m ) {for(int i=0;i<m-1;i++) R[i] = A[i];return ;}poly_revcopy(tmp4, B, m); poly_inv(tmp4, tmp5, n-m+1);poly_revcopy(tmp4, A, n); poly_mul(tmp4, tmp5, tmp6, n-m+1, n-m+1);poly_revcopy(tmp4, tmp6, n-m+1); poly_copy(tmp5, B, m);poly_mul(tmp4, tmp5, tmp6, n-m+1, m);for(int i=0;i<m-1;i++) R[i] = (A[i] + MOD - tmp6[i])%MOD;poly_clear(tmp4, 0, n); poly_clear(tmp5, 0, n);}vector<int>P[MAXN + 5];void poly_build(int *X, int x, int l, int r) {if( l == r ) {P[x].push_back(MOD-X[l]), P[x].push_back(1);return ;}int mid = (l + r) >> 1;poly_build(X, x<<1, l, mid), poly_build(X, x<<1|1, mid+1, r);for(int i=0;i<P[x<<1].size();i++) tmp4[i] = P[x<<1][i];for(int i=0;i<P[x<<1|1].size();i++) tmp5[i] = P[x<<1|1][i];poly_mul(tmp4, tmp5, tmp6, P[x<<1].size(), P[x<<1|1].size());for(int i=0;i<P[x<<1].size()+P[x<<1|1].size()-1;i++) P[x].push_back(tmp6[i]);poly_clear(tmp4, 0, P[x<<1].size()), poly_clear(tmp5, 0, P[x<<1|1].size());}int tmp7[MAXN + 5], tmp8[25][MAXN + 5];void poly_eval(int *A, int *Y, int dep, int n, int l, int r, int x) {for(int i=0;i<P[x].size();i++) tmp7[i] = P[x][i];poly_mod(A, tmp7, tmp8[dep], n, P[x].size());if( l == r ) {Y[l] = tmp8[dep][0];return ;}int mid = (l + r) >> 1;poly_eval(tmp8[dep], Y, dep+1, P[x].size()-1, l, mid, x<<1);poly_eval(tmp8[dep], Y, dep+1, P[x].size()-1, mid+1, r, x<<1|1);poly_clear(tmp8[dep], 0, P[x].size()-1);} }oper; int f[MAXN + 5], a[MAXN + 5], g[MAXN + 5]; int main() {int n, m; scanf("%d%d", &n, &m); n++;for(int i=0;i<n;i++)scanf("%d", &f[i]);for(int i=0;i<m;i++)scanf("%d", &a[i]);oper.poly_build(a, 1, 0, m-1);oper.poly_eval(f, g, 0, n, 0, m-1, 1);for(int i=0;i<m;i++)printf("%d\n", g[i]); }

@例題與應(yīng)用@

@快速插值@

就是下面那個(gè)東西，對(duì)，它就是多點(diǎn)求值的一個(gè)應(yīng)用。

@2 - 快速插值@

@理論推導(dǎo)@

這個(gè)才是真正的毒瘤……

我們已知任意的 N+1 個(gè)點(diǎn) \((x_0,y_0)\)，\((x_1,y_1)\)，...，\((x_N,y_N)\)，想要還原一個(gè)多項(xiàng)式 \(A(x)\) 使得 \(A(x_i)=y_i\)。

首先來看一個(gè)東西，拉格朗日插值：
\[f(x) = \sum_{i=0}^{N}(\dfrac{\prod_{j=0,j\not =i}^{N}(x-x_j)}{\prod_{j=0,j\not =i}^{N}(x_i-x_j)})*y_i\]

怎么證明呢？其實(shí)比較簡單。首先證明其正確性：
記 \(P_i(x)=(\dfrac{\prod_{j=0,j\not =i}^{N}(x-x_j)}{\prod_{j=0,j\not =i}^{N}(x_i-x_j)})*y_i\)，
當(dāng) \(i = j\) 時(shí)，分子等于分母， \(P_i(x_j) = y_i\)；
當(dāng) \(i\not = j\) 時(shí)，分子等于零， \(P_i(x_j) = 0\)。
故 \(f(x_i)=\sum_{j=0}^{N}P_j(x_i)=y_i\)。

然后證明其唯一性，可以采用線性方程組的方法。

但是……我們?nèi)绻┝砘喩厦婺莻€(gè)式子，時(shí)間復(fù)雜度就高上天了。
所以我們接下來就來亂搞一下。

首先求解分母的部分，即 \(\prod_{j=0,j\not =i}^{N}(x_i-x_j)\) 的值。
記 \(g(x) = \prod_{i=0}^N(x-x_i)\)，則我們相當(dāng)于是求解：
\[\lim_{x->x_i}\dfrac{g(x)}{(x-x_i)}\]
當(dāng) \(x = x_i\) 時(shí)分母是等于 0 的，沒有意義，所以我們寫成極限的形式。
根據(jù)洛必達(dá)法則：
\[\lim_{x->x_i}\dfrac{g(x)}{(x-x_i)}=g'(x)\]
即我們相當(dāng)于要求解 \(g'(x_0)\)，\(g'(x_1)\)，...，\(g'(x_N)\)，快速插值即可。

……我相信你們應(yīng)該會(huì)一些簡單的微積分吧。
算了不重要，我繼續(xù)往下講。

【高能預(yù)警】下面的公式非常密集，請(qǐng)做好心理準(zhǔn)備。

記 \(k_i = \dfrac{y_i}{\prod_{j=0,j\not =i}^{N}(x_i-x_j)}\)，分子已知，分母可以通過我們上邊的方法求解。

則原式變?yōu)?#xff1a;
\[f(x)=\sum_{i=0}^{N}k_i*(\prod_{j=0,j\not =i}^{N}(x-x_j))\]

記 \(P_l(x)=\prod_{i=0}^{\frac{N}{2}}(x-x_i)\)，\(P_r(x)=\prod_{i=\frac{N}{2}+1}^{N}(x-x_i)\)。

再記 \(f_l(x)=\sum_{i=0}^{\frac{N}{2}}k_i*(\prod_{j=0,j\not =i}^{\frac{N}{2}}(x-x_j))\)，\(f_r(x)=\sum_{i=\frac{N}{2}+1}^{N}k_i*(\prod_{j=\frac{N}{2}+1,j\not =i}^{N}(x-x_j))\)

則 \(f(x)\) 可以寫成：
\[f(x)=P_r(x)*f_l(x)+P_l(x)*f_r(x)\]
預(yù)處理出 \(P(x)\)，分治即可。最底層直接返回 \(k_i\)。

怎么理解呢？對(duì)于 \(h_i(x) = k_i*(\prod_{j=0,j\not =i}^{N}(x-x_j))\)，它 “缺” \((x-x_i)\) 這一項(xiàng)。
對(duì)于 \(0\le i\le \frac{N}{2}\)，它們都不 “缺” \((x-x_j),\frac{N}{2}+1\le j\le N\) 這一項(xiàng)。
所以利用乘法分配律，就可以得到上面那個(gè)式子。

一個(gè)細(xì)節(jié)：多點(diǎn)求值用的 \(P(x)\) 與快速插值后面分治用的 \(P(x)\) 是同一個(gè)玩意兒，所以可以不用重復(fù)求解。

時(shí)間復(fù)雜度：\(O(n\log^2n)\)，常數(shù)極大。

@（不建議參考的）代碼@

本代碼為 luoguP5518 的 AC 代碼

#include<cstdio> #include<vector> #include<algorithm> using namespace std; const int G = 3; const int MOD = 998244353; const int MAXN = 100000*10; int pow_mod(int b, int p) {int ret = 1;while( p ) {if( p & 1 ) ret = 1LL*ret*b%MOD;b = 1LL*b*b%MOD;p >>= 1;}return ret; } int inv[MAXN + 5]; void init() {inv[1] = 1;for(int i=2;i<=MAXN;i++)inv[i] = 1LL*(MOD - MOD/i)*inv[MOD%i]%MOD; } struct Polynomial{void poly_copy(int *A, int *B, int n) {for(int i=0;i<n;i++)A[i] = B[i];}void poly_copy(int *A, vector<int> B) {for(int i=0;i<B.size();i++)A[i] = B[i];}void poly_clear(int *A, int l, int r) {for(int i=l;i<r;i++)A[i] = 0;}void poly_revcopy(int *A, int *B, int n) {for(int i=0;i<n;i++)A[i] = B[n-i-1];}void ntt(int *A, int n, int type) {for(int i=0,j=0;i<n;i++) {if( i < j ) swap(A[i], A[j]);for(int l=(n>>1);(j^=l)<l;l>>=1);}for(int s=2;s<=n;s<<=1) {int t = (s>>1);int u = (type == -1) ? pow_mod(G, (MOD-1) - (MOD-1)/s) : pow_mod(G, (MOD-1)/s);for(int i=0;i<n;i+=s) {for(int j=0,p=1;j<t;j++,p=1LL*p*u%MOD) {int x = A[i+j], y = 1LL*p*A[i+j+t]%MOD;A[i+j] = (x + y)%MOD, A[i+j+t] = (x + MOD - y)%MOD;}}}if( type == -1 ) {int inv = pow_mod(n, MOD-2);for(int i=0;i<n;i++)A[i] = 1LL*A[i]*inv%MOD;}}int tmp1[MAXN + 5], tmp2[MAXN + 5], tmp3[MAXN + 5];void poly_mul(int *A, int *B, int *C, int n, int m) {int len; for(len = 1;len < n+m-1;len <<= 1);poly_copy(tmp1, A, n); poly_clear(tmp1, n, len);poly_copy(tmp2, B, m); poly_clear(tmp2, m, len);ntt(tmp1, len, 1); ntt(tmp2, len, 1);for(int i=0;i<len;i++) tmp3[i] = 1LL*tmp1[i]*tmp2[i]%MOD;ntt(tmp3, len, -1); poly_copy(C, tmp3, n+m-1);}void poly_inv(int *A, int *B, int n) {if( n == 1 ) {B[0] = pow_mod(A[0], MOD-2);return ;}int len; for(len = 1;len < (n<<1);len <<= 1);poly_inv(A, B, (n + 1) >> 1);poly_copy(tmp3, A, n); poly_clear(tmp3, n, len);ntt(tmp3, len, 1); ntt(B, len, 1);for(int i=0;i<len;i++) B[i] = 1LL*B[i]*(2 + MOD - 1LL*tmp3[i]*B[i]%MOD)%MOD;ntt(B, len, -1); poly_clear(B, n, len);}int tmp4[MAXN + 5], tmp5[MAXN + 5], tmp6[MAXN + 5];void poly_divide(int *A, int *B, int *C, int *R, int n, int m) {poly_revcopy(tmp4, B, m); poly_clear(tmp5, 0, 2*(n-m+1)); poly_inv(tmp4, tmp5, n-m+1);poly_revcopy(tmp4, A, n); poly_mul(tmp4, tmp5, tmp6, n-m+1, n-m+1);poly_revcopy(C, tmp6, n-m+1);poly_copy(tmp4, C, n-m+1); poly_copy(tmp5, B, m);poly_mul(tmp4, tmp5, tmp6, n-m+1, m);for(int i=0;i<m-1;i++) R[i] = (A[i] + MOD - tmp6[i])%MOD;}void poly_mod(int *A, int *B, int *R, int n, int m) {if( n < m ) {for(int i=0;i<m-1;i++) R[i] = A[i];return ;}poly_revcopy(tmp4, B, m); poly_inv(tmp4, tmp5, n-m+1);poly_revcopy(tmp4, A, n); poly_mul(tmp4, tmp5, tmp6, n-m+1, n-m+1);poly_revcopy(tmp4, tmp6, n-m+1); poly_copy(tmp5, B, m);poly_mul(tmp4, tmp5, tmp6, n-m+1, m);for(int i=0;i<m-1;i++) R[i] = (A[i] + MOD - tmp6[i])%MOD;poly_clear(tmp4, 0, n); poly_clear(tmp5, 0, n);}vector<int>P[MAXN + 5];void poly_build(int *X, int x, int l, int r) {if( l == r ) {P[x].clear(); P[x].push_back(MOD-X[l]), P[x].push_back(1);return ;}int mid = (l + r) >> 1;poly_build(X, x<<1, l, mid), poly_build(X, x<<1|1, mid+1, r);poly_copy(tmp4, P[x<<1]), poly_copy(tmp5, P[x<<1|1]);poly_mul(tmp4, tmp5, tmp6, P[x<<1].size(), P[x<<1|1].size()); P[x].clear();for(int i=0;i<P[x<<1].size()+P[x<<1|1].size()-1;i++) P[x].push_back(tmp6[i]);poly_clear(tmp4, 0, P[x<<1].size()), poly_clear(tmp5, 0, P[x<<1|1].size());}int tmp7[MAXN + 5], tmp8[25][MAXN + 5];void poly_eval(int *A, int *Y, int dep, int n, int l, int r, int x) {for(int i=0;i<P[x].size();i++) tmp7[i] = P[x][i];poly_mod(A, tmp7, tmp8[dep], n, P[x].size());if( l == r ) {Y[l] = tmp8[dep][0];return ;}int mid = (l + r) >> 1;poly_eval(tmp8[dep], Y, dep+1, P[x].size()-1, l, mid, x<<1);poly_eval(tmp8[dep], Y, dep+1, P[x].size()-1, mid+1, r, x<<1|1);poly_clear(tmp8[dep], 0, P[x].size()-1);}void poly_dif(int *A, int *B, int n) {for(int i=1;i<n;i++)B[i-1] = 1LL*A[i]*i%MOD;B[n-1] = 0;}void poly_itgr(int *A, int *B, int n) {for(int i=n-1;i>=0;i--)B[i+1] = 1LL*A[i]*inv[i+1]%MOD;}int tmp9[MAXN + 5], tmp10[MAXN + 5], tmp11[25][MAXN + 5], tmp12[MAXN + 5];void poly_itplt(int *A, int dep, int x, int l, int r) {poly_clear(A, 0, P[x].size());if( l == r ) {A[0] = tmp10[l];return ;}int mid = (l + r) >> 1;poly_itplt(tmp11[dep], dep+1, x<<1, l, mid);poly_copy(tmp9, P[x<<1|1]); poly_mul(tmp9, tmp11[dep], tmp12, P[x<<1|1].size(), P[x<<1].size());for(int i=0;i<P[x].size();i++) A[i] = (A[i] + tmp12[i])%MOD;poly_itplt(tmp11[dep], dep+1, x<<1|1, mid+1, r);poly_copy(tmp9, P[x<<1]); poly_mul(tmp9, tmp11[dep], tmp12, P[x<<1].size(), P[x<<1|1].size());for(int i=0;i<P[x].size();i++) A[i] = (A[i] + tmp12[i])%MOD;}void poly_itplt(int *X, int *Y, int *A, int n) {poly_build(X, 1, 0, n-1); poly_copy(tmp9, P[1]);poly_dif(tmp9, tmp9, P[1].size());poly_eval(tmp9, tmp10, 0, P[1].size()-1, 0, n-1, 1);for(int i=0;i<n;i++) tmp10[i] = 1LL*Y[i]*pow_mod(tmp10[i], MOD-2)%MOD;poly_itplt(A, 0, 1, 0, n-1);} }oper; int x[MAXN + 5], y[MAXN + 5], f[MAXN + 5]; int main() {int n; scanf("%d", &n);for(int i=0;i<n;i++)scanf("%d%d", &x[i], &y[i]);oper.poly_itplt(x, y, f, n);for(int i=0;i<n;i++)printf("%d ", f[i]); }

@例題與應(yīng)用@（暫無）

這個(gè)東西……能夠應(yīng)用？？？

轉(zhuǎn)載于:https://www.cnblogs.com/Tiw-Air-OAO/p/10198693.html

總結(jié)

以上是生活随笔為你收集整理的@总结 - 4@ 多项式的多点求值与快速插值的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。