1. 程式人生 > 其它 >OI 中的多項式(填坑中)

OI 中的多項式(填坑中)

OI 中的多項式

目前已有知識點:

  • 拉格朗日插值
  • 快速傅立葉變換 FFT

拉格朗日插值

模板

給定 \(n\) 個點 \((x_i,y_i)\),要求求出一個多項式函式 \(f(x)\),滿足 \(f(x)\) 的影象經過這 \(n\) 個點,並且 \(f(x)\) 的度為 \(n-1\)。給出一個 \(k\),求 \(f(k)(\bmod 998\ 244\ 353)\)。(from Luogu P4781

一個很樸素的辦法就是直接列出方程組,然後用高斯消元求解。這樣的時間複雜度是 \(\mathcal O(n^3)\) 的,並不算一個高效的演算法。而使用拉格朗日插值可以在 \(\mathcal O(n^2)\)

的時間內解決這個問題,是很高效的一種演算法。

設拉格朗日的基本多項式為

\[l_j(t) = \prod\limits_{i\neq j}\dfrac{t-x_i}{x_j-x_i} \]

那麼會發現對於一個 \(l_j\),會有 \(l_j(x_j)=1\) 並且 \(l_j(x_i)=0,\forall i \neq j\)

有了這個基本多項式,那麼可以就構造出函式 \(f(x)\)

\[f(x) = \sum\limits_{i} y_i \cdot l_i(x) \]

容易發現,這個函式一定是滿足題目中給出的 \(n\) 個點的。將兩個式子合併到一起,就可以得到答案:

\[f(x)=\sum\limits_i (y_i\cdot\prod\limits_{j\neq i} \dfrac{x-x_j}{x_i-x_j}) \]

由於題目讓我們求 \(f(k)\)

的值,所以可以直接將上面式子中的 \(x\) 替換成 \(k\) 進行計算。

程式碼

快速傅立葉變換 FFT

前置知識

複數

複數 \(z\) 是形如 \(a+bi\) 的數,其中 \(i^2=-1\)\(a\) 稱為實部,\(b\) 稱為虛部。

複平面是一個笛卡爾平面,有兩條座標軸,橫軸是實部,縱軸是虛部。那麼根據複數的形式可知,複數可以用複平面上的一個向量 \((a,b)\) 表示。以橫軸正方向為始邊,向量為終邊,逆時針轉過的角度叫做複數 \(z\) 的輻角 \(\theta\),並且 \(\theta\) 必須為正數。

複數的模 \(|z|=\sqrt{a^2+b^2}\),與向量的模相同。

複數運算
加減法

加減法很好理解,實部和虛部對應相加即可。

\[z_1\pm z_2=(a_1+b_1i)\pm(a_2+b_2i)=(a_1\pm a_2)+(b_1\pm b_2)i \]
乘法

類似二項式乘法的計算方法。

\[\begin{array}{lllll} z_0&=&z_1\times z_2 &=& (a_1+b_1i) \times (a_2 + b_2i) \\ &&&=& a_1a_2+a_1b_2i+a_2b_1i+b_1b_2i^2\\ &&&=& (a_1a_2-b_1b_2)+(a_1b_2+a_2b_1)i \end{array} \]

看一下在複平面上的意義。發現 \(\theta_0=\theta_1+\theta_2\),並且 \(|z_0|=|z_1|\times|z_2|\)

尤拉公式

\[e^{i\pi}=1 \]

單位根

複數域下,滿足 \(x^n=1\)\(x\) 被稱作 \(n\) 次單位根。

顯然這樣的 \(x\) 會有 \(n\) 個。按照輻角大小排序後第 \(k\) 大的根是 \(e^{i\frac{2k\pi}{n}}\)

可以發現,這 \(n\) 個單位根在複平面上剛好 \(n\) 等分單位圓。

本原單位根

代數上的不好理解,就不提了。複平面上輻角最小的單位根就叫做本原單位根。

\(n\) 次本原單位根記作 \(\omega_n\),那麼有 \(\omega_n=e^{i\frac{2k}n}\)。根據虛數乘法在複平面上的表示方法,可以發現第 \(k\) 大的 \(n\) 次單位根可以表示成 \(\omega_n^k\)

關於本原單位根,有一個很重要的性質:

\(n = 2m\),則有:

\[\begin{array}{lll}\omega_n^{2k}&=&\omega_m^k\\\omega_n^{m+k}&=&-\omega_n^k\end{array} \]

關於這兩個式子的證明,可以畫出複平面過後然後用旋轉的想法去證明,這裡不再給出(其實是懶)

離散傅立葉變換 DFT

考慮兩個多項式做乘法,如果直接用係數表示法的話顯然時間複雜度是 \(\mathcal O(n^2)\) 的,並且很難優化。而使用點值表示法的話,計算乘法的時間複雜度僅是 \(\mathcal O(n)\),但是將係數表示法轉化為點值表示法的時間複雜度是 \(\mathcal O(n^2)\),並且將點值表示法轉化稱為係數表示法需要高斯消元,時間複雜度是 \(\mathcal O(n^3)\) 的,比直接係數表示法算慢多了。那麼如果可以找到一種能夠快速將係數表示法轉化成點值表示法並且同樣快速的轉化回來的演算法,那麼用點值表示法計算將會快很多。

將係數表示法轉化為點值表示法的變換叫做離散傅立葉變換(DFT),將點值表示法轉化為係數表示法的變換叫做離散傅立葉逆變換(IDFT)

將一個 \(n-1\) 次多項式表示成為 \(A(x)\),方便下面的書寫。

考慮構造出一個長度為 \(n\) 的數列 \(\{b_i\}\),其中 \(b_k=A(\omega_n^k)\),即(將 \(\sum\) 內的變數寫作 \(j\) 的原因是避免與虛數的 \(i\) 混淆):

\[b_k=\sum\limits_{j=0}^{n-1}a_j\times\omega_n^{kj} \]

這裡的 \(b_k\) 就表示了多項式在 \(\omega_n^k\) 處的點值。如果直接求這一過程,時間複雜度是 \(\mathcal O(n^2)\) 的,如果使用快速傅立葉變換(FFT)就可以在 \(\mathcal O(n\log n)\) 的時間內求出這一點值。

快速傅立葉變換 FFT

\(A(x)\) 按照指數的奇偶性分成兩組(\(n=2m\)):

\[\begin{array}{lll}A(x)&=&\displaystyle\sum\limits_{j=0}^{n-1} a_j\cdot x^j\\&=& \displaystyle\sum\limits_{j=0}^m(a_{2j}\cdot x^{2j}) + \sum\limits_{j=0}^m(a_{2j+1}\cdot x^{2j+1})\end{array} \]

將後面那個求和式中的最後一個 \(x\) 拿出來:

\[\begin{array}{lll}A(x) &=& \displaystyle \sum\limits_{j=0}^{m-1}(a_{2j}\cdot x^{2j}) + x\sum\limits_{j=0}^{m-1}(a_{2j+1}\cdot x^{2j}) \\&=& \displaystyle \sum\limits_{j=0}^{m-1}(a_{2j}\cdot(x^2)^j) + x\sum\limits_{j=0}^{m-1}(a_{2j+1}\cdot (x^2)^j)\end{array} \]

假設:

\[\begin{array}{lll}A_0(x) &=& \displaystyle \sum\limits_{j=0}^{m-1} a_{2j}\cdot x^j\\A_1(x) &=& \displaystyle \sum\limits_{j=0}^{m-1} a_{2j+1}\cdot x^j\end{array} \]

那麼有:

\[A(x) = A_0(x^2)+x\times A_1(x^2) \]

如果知道了 \(A_0\)\(A_1\) 的點值,那麼就可以 \(\mathcal O(n)\) 去計算 \(A(x)\) 在各點的點值了,並且會發現求 \(A_0,A_1\) 的過程與求解 \(A\) 的方式是一樣的,所以計算 \(A_0,A_1\) 可以直接遞迴處理。但是計算一下時間複雜度,會發現這樣仍然是 \(\mathcal O(n^2)\) 的,因此考慮去優化這個做法。

在前置知識中提到了關於本原單位根的兩個公式:

\[\begin{array}{rll}\omega_n^{2k}&=&\omega_m^k\\\omega_n^{m+k}&=&-\omega_n^k\end{array} \]

那麼根據第二個公式,啟發我們只計算 \(A(x)\)\(m\) 個位置的點值,然後嘗試發掘前 \(m\) 個位置和後 \(m\) 個位置的關係。

\(0\le k<m\),那麼上面推出來的式子可以化成:

\[A(\omega_n^k)=A_0(\omega_n^{2k})+\omega_n^k\cdot A_1(\omega_n^{2k}) \]

那麼根據第一個式子可以將這個式子化成:

\[A(\omega_n^k)=A_0(\omega_m^k)+\omega_n^k\cdot A_1(\omega_m^k) \]

然後對於後 \(m\) 個位置,可以寫成:

\[A(\omega_n^{m+k}) = A_0(\omega_n^{2(m+k)}) + \omega_n^{m+k} \cdot A_1 (\omega_n^{2(m+k)}) \]

用第二個式子可以化成:

\[A(\omega_n^{m+k}) = A_0(\omega_n^{2k}) - \omega_n^k \cdot A_1(\omega_n^{2k}) \]

然後效仿前一個的推法:

\[A(\omega_n^{m+k}) = A_0(\omega_m^k) - \omega_n^k \cdot A_1(\omega_m^k) \]

會發現求 \(A(x)\) 的時候只需要求到 \(A_0,A_1\)\(m\) 項的值即可,相比於之前,每次遞迴計算的值少了一半,所以時間複雜度變成了 \(\mathcal O(n\log n)\)

離散傅立葉逆變換 IDFT

目前已知了一個多項式經過 DFT 後變成的點值表示式,即上面所說的 \(\{b_k\}\) 數列,也就是:

\[b_k = \sum\limits_{j=0}^{n-1} a_j\cdot \omega_n^{kj} \]

要求原係數表示式數列 \(\{a_k\}\)

這裡直接給出結論(不會證)

\[a_k = \dfrac 1 n \sum\limits_{j=0}^{n-1} b_j\cdot \omega^{-kj} \]

如果直接求解的話仍然是 \(\mathcal O(n^2)\) 的,所以考慮用快速傅立葉逆變換(IFFT)來進行這一操作。

快速傅立葉逆變換 IFFT

會發現求解的這個式子和 FFT 是高度一致的。因為 \(\omega_n^{-k} = \omega_n^{n+k}\),所以可以直接用 FFT,然後將計算出來的 \(\{b_k\}\) 翻轉過來就是 IFFT 得到的式子。


這樣,我們得到了一個可以 \(\mathcal O(n\log n)\) 計算多項式乘法的演算法。用洛谷的模板題給出程式碼。

給定多項式 \(F(x),G(x)\),求 \(F(x)\times G(x)\)。(from Luogu P3803

值得一說的是,這道題在之前遞迴做法因為常數過大是被卡了的,但是現在可以比較輕鬆跑過去了。

程式碼