[學習筆記&教程] 信號, 集合, 多項式, 以及各類卷積性變換 (FFT,NTT,FWT,FMT)

目錄

• 信號, 集合, 多項式, 以及卷積性變換
  • 卷積
    • 卷積性變換
  • 傅里葉變換與信號
    • 引入: 信號分析
    • 變換的基礎: 複數
    • 傅里葉變換
    • 離散傅里葉變換
  • FFT 與多項式
    • \(n\) 次單位復根
    • 消去引理, 折半引理與求和引理
    • 從新定義
    • 多項式的表示
    • 快速傅里葉變換FFT
    • 經過 FFT 在單位複數根處插值
    • FFT的速度優化與迭代實現
  • 炸精現場與 NTT
    • 原根
    • NTT
    • 任意模數 NTT
  • 卷積狀物體與分治 FFT
  • FWT 與位運算卷積
    • FWT 與 \(\text{or}\) 卷積
    • FWT 與 \(\text{and}\) 卷積
    • FWT 與 \(\text{xor}\) 卷積
  • FMT 與子集卷積
  • 卷積與字符串
    • 卷積與迴文序列
    • 卷積與通配符
    • 一點優化?
  • 寫在最後

信號, 集合, 多項式, 以及卷積性變換

寫在前面的須知內容: 這篇博客主要講的是板子以及一些理解變換本質的內容, 並無幾道例題.

卷積

理論上你們對卷積這個東西應該並不陌生了吧? 不過你們對於卷積都有什麼理解?

• 多項式乘法?
• 狄利克雷卷積?

其實卷積能夠統一成同一個形式:

給定兩個 \(n\) 維向量 \(\vec a,\vec b\) 和一個下標運算 \(\circ\), 那麼若 \(\vec c\) 知足下式:
\[ c_k=\sum_{i\circ j=k}a_i b_j \]
那麼稱 \(\vec c\) 爲 \(\vec a\) 和 $\vec b $ 的離散卷積, 簡稱卷積(題外話: 連續卷積是個積分的形式), 記做 \(\vec c = \vec a *_{\circ} \vec b\). 在不引發混淆的狀況下可記做 \(\vec c = \vec a * \vec b\), 其中 \(*\) 稱爲卷積算子. (因此不要在 \(\LaTeX\) 裏亂用 \(*\) ...老老實實打 \(\times\) ...)

咱們將下標集合設爲 \(S\), 則 \((S,\circ)\) 組成了一個代數系統, 咱們稱爲卷積的下標系統.

常見的卷積:

• 多項式的乘積: 等價於下標系統爲 \(\mathbb{Z}\) 上 (OI中其實通常是 \(\mathbb{N}\) ) 的加法運算 (這是FFT/NTT要作的事情).
• 循環卷積: 下標系統爲模 \(n\) 意義下的加法. (這(大概)也是FFT/NTT要作的事情)
• 位運算卷積: 下標系統爲 \(n\) 位二進制數的運算 (與, 或, 異或...) (這是FWT要作的事情)
• 子集卷積: 下標是集合, 下標運算對於兩個交集非空的集合無定義, 對於兩個不相交集合定義爲兩者的並. (這是FMT要作的事情)
• 狄利克雷卷積: 下標系統即爲正整數乘法. (這是Chaiki_Cage要作的事情)

卷積性變換

對於變換 \(T\) 以及向量(數列) \(\vec a,\vec b\) , 若是知足:
\[ T(\vec a)\cdot T(\vec b)=T(\vec a * \vec b) \]
那麼稱變換 \(T\) 知足卷積定理, 這種性質被稱爲卷積性.

咱們接下來要講的 DFT, FFT, NTT, FWT, FMT, BFT, FST 都是卷積性變換, 用途就是根據上式將卷積轉化爲點積來加速卷積的計算.

傅里葉變換與信號

引入: 信號分析

講FFT怎麼能直接從多項式開始講呢

如下內容截取自 3Blue1Brown 的視頻 形象化展現傅立葉變換. 原視頻超棒的! 若是有條件的話必定要看看!

首先咱們先從一個看起來和OI絕不相關的(然而其實炒雞重要的)研究課題開始.

給定一個波, 如何將它分解爲若干正弦波的疊加?

或者說如今給定了一個信號, 咱們要把它分解成幾個頻率的疊加.

(其實這就是FFT本來要作的事情)

好比下面這兩個信號:

它們疊加起來就是下面這個鬼東西:

其實這兩個波的狀況最後疊加出來的仍是能看出些許規律的. 然而若是是長這樣的:

當場去世.png

按照常理, 咱們確定只能探測到最上面的那個波. 那麼咱們要怎麼根據它來反推出它的組成成分呢?

這個時候咱們須要想辦法引入另外一個頻率, 而且讓這個波在特定組分頻率下能體現出不一樣.

怎麼體現頻率?

要頻率, 就得有周期.

怎麼要週期? 循環?

對!

讓這個信號轉圈圈!

好比對於下面這個 3bps 的純正弦波, 咱們把它以 0.5rps 的速度卷在原點附近:

由於 3 正好是 0.5 的整倍數因此看起來好像變少了...實際上:

沒錯, 它真的是捲起來的...

因而咱們如今完成了第一步: 引入一個頻率.

第二步是想辦法讓這個信號在對應組分頻率上表現得特別一些(否則咱們怎麼知道當前頻率是否是這個信號的一個組分)

咱們首先看看若是轉圈的頻率和信號組分頻率同樣會發生什麼:

額密密麻麻...

它如今有啥特徵?

以前纏繞頻率和組分頻率不一樣的時候信號在相對原點的各個方向上分佈得比較均勻, 而如今這個圖像明顯要偏向右側.

爲啥?

你看當纏繞頻率和組分頻率一致的時候, 咱們每次轉完一圈的, 波峯和波谷都在同一條穿過原點的直線上. 固然他們要整個偏向右側啦

那麼咱們怎麼表示這個特徵?

咱們先用通俗一點的概念來表示它: 重心!

想象這個波是些有質量的東西 (好比鐵絲), 而後咱們就能夠用重心距離原點的距離來度量整個波捲起來以後偏離中心的程度.

當纏繞頻率和組分頻率不一樣的時候:

能夠看到重心距離原點很近.

當纏繞頻率和組分頻率相同 (或者接近) 的時候:

重心就距離原點比較遠.

或者說, 這個纏繞的過程能夠想辦法把波峯和波谷的影響在一條固定的直線上疊加 (波峯最大, 放在正方向上; 波谷最小, 放在負方向上, 因而影響疊加).

因而咱們能夠嘗試用這個原理來對這個純正弦波分解一下:

左下角的圖中橫座標是纏繞頻率, 縱座標是重心的 \(x\) 座標值.

能夠看到在組分頻率 \(3\) 那裏產生了一個峯.

等等爲啥 \(0\) 那裏也有個大峯?

其實上面那個波並非 "純正弦波". 真正的純正弦波應該是長這樣的:

這樣的話纏一纏就會變成這樣:

(纏出來其實就是個圓, 由於波谷是負的)

顯然咱們只用單一組分的純正弦波是沒法體現出 "分解" 頻率的過程的. 因而咱們拿下面兩個波的疊加來試驗一下

咱們把它卷卷卷...

組分 \(2\texttt{Hz}\) 的波峯和波谷產生的貢獻疊加了, 在最終的圖像中產生了一個峯. 咱們繼續纏:

實際上咱們發現, 把波疊加的同時, 最終產生的變換圖像也會疊加 (顯然它們在對應的點上都是算術相加的關係對吧)

仔細想想這個地方.

那麼咱們怎麼在數學上實現這個頻率分解過程呢?

首先咱們講一下怎麼讓信號轉圈圈

變換的基礎: 複數

如下是高考內容, 請掉線的同窗當即嘗試重連(怎麼可能掉線)

最開始引入複數是爲了解決這種東西:
\[ x^2=-1 \]
顯然它在實數範圍內無解對吧

那麼咱們強行讓它有解

怎麼強行讓它有解呢?

定義一個唄

因而 \(i\) 的定義就是這樣的:
\[ i^2=-1 \]
\(i\) 被稱爲虛數單位. 由於在創造它的時候認爲它是個在天然界裏沒卵用的假數, 因此叫虛數.

(實際上物理裏交流電分析裏要用虛數的. 並且有一部分原本能作功的能量還會被虛部吃掉(變焦耳熱了))

然而只有虛數還不行, 咱們可能還得有實數部分, 因而咱們把它們複合起來的複雜的東西叫作複數(complex number).

任意一個複數均可以表示爲 \(a+bi\) 的形式 (\(a,b\in \mathbb{R}\)). (全體複數集定義爲 \(\mathbb{C}\)).

其中 \(a\) 爲複數 \(a+bi\) 的實部, \(b\) 爲複數 \(a+bi\) 的虛部.

複數的普通運算

首先是簡單一點的, 複數加減法:
\[ (a+bi)+(c+di)=(a+c)+(b+d)i \\ (a+bi)-(c+di)=(a-c)+(b-d)i \]
複數乘法按照定義的話應該是長這樣的:
\[ \begin{align} (a+bi)(c+di)&=ac+a(di)+(bi)c+(bi)(di) \\ &=ac+adi+bci+bdi^2\\ &=ac+adi+bci-bd\\ &=(ac-bd)+(ad+bc)i \end{align} \]
然而實際上它還有一種解釋方式.

咱們引入一個座標系, 橫座標爲複數的實部, 縱座標爲複數的虛部, 咱們稱其爲複平面. 那麼一個複數均可以對應複平面上的一個點(或者說一個向量).

既然咱們有了一個座標系, 咱們就能夠接着定義一個複數的另外兩個參數: 模長表示當前複數對應點到原點的距離, 幅角表示原點到當前複數所在點與 \(x\) 軸正方向的有向夾角( \([0,2\pi)\) ).

因而一個複數也能夠表示爲 \(r\angle\theta\) 的形式(吧)(反正卡西歐這麼表示的(逃))

那麼複數乘法能夠定義爲將兩個複數的模長相乘, 幅角相加, 或者說:
\[ (r\angle\theta)(s\angle\varphi)=(rs)\angle(\theta+\varphi) \]
爲啥?

咱們來證實一下:

設 \(\alpha\) 爲複數 \(a+bi\) 的幅角, \(\beta\) 爲複數 \(c+bi\) 的幅角, 則:
\[ \tan\alpha=\frac b a, \tan\beta=\frac d c \]
那麼根據 \(\tan\) 的和角公式, 咱們有:
\[ \begin{align} \tan(\alpha+\beta)&=\frac{\tan\alpha + \tan\beta}{1-\tan\alpha\tan\beta} \\ &=\frac{\frac b a+\frac d c}{1-\frac{bd}{ac}} \\ &=\frac{\frac {bc} {ac}+\frac {ad} {ac}}{1-\frac{bd}{ac}} \\ &=\frac{bc+ad}{ac-bd} \end{align} \]
而 \((a+bi)(c+di)\) 就等於 \((ac-bd)+(ad+bc)i\).

(其實用 \(\tan\) 來證實是不夠嚴謹的, 由於它的週期是 \(\pi\) 而不是 \(2\pi\) ...然而 \(\sin\) 和 \(\cos\) 的式子太大懶得搞了...有興趣本身試試吧)

模長相乘很是好搞, 直接就出來了
\[ \begin{aligned} \sqrt{a^2+b^2}\times\sqrt{c^2+d^2}&=\sqrt{(a^2+b^2)(c^2+d^2)} \\ &=\sqrt{a^2c^2+a^2d^2+b^2c^2+b^2d^2}\\ \sqrt{(ac-bd)^2+(ad+bc)^2}&=\sqrt{(a^2c^2-2abcd+b^2d^2)+(a^2d^2+2abcd+b^2c^2)}\\ &=\sqrt{a^2c^2+b^2d^2+a^2d^2+b^2c^2}\\ \\ \sqrt{a^2+b^2}\times\sqrt{c^2+d^2}&=\sqrt{(ac-bd)^2+(ad+bc)^2} \end{aligned} \]

複數的高端運算

下面有請宇宙第一裝B公式:
\[ e^{i\pi}+1=0 \]
這個公式實際上是Taylor展開搞出來的, 實際上更廣義的公式是長這樣的:
\[ e^{i\varphi}=\cos\varphi+i\sin\varphi=1\angle\varphi \]
當 \(\varphi\) 取 \(\pi\) 的時候就是歐拉公式辣 (815班旗上的元素之一)

也就是下面這個圖:

泰勒展開的具體式子就不放了...放了大概也沒卵用...

有了上面的運算法則, 咱們就能夠定義傅里葉變換了

傅里葉變換

咱們剛剛搞出了 \(e^{i\varphi}=\cos\varphi+i\sin\varphi=1\angle\varphi\) 這個式子, 而後又知道根據複數運算法則, 複數相乘時模長相乘幅角相加, 那麼一個複數乘以 \(e^{i\varphi}\) 其實就至關於在複平面上旋轉了 \(\varphi\) 的角度.

這就是轉圈圈的數學實現了!

假設這個波形是函數 \(g\), 咱們只要計算 \(F(x)=g(x)e^{ikx}\) 的圖像的重心就能夠了!

其中 \(k\) 做爲係數決定纏繞速度. 爲了讓它和頻率直接相關, 咱們改造一下就星了:
\[ F(t)=g(t)e^{2\pi ift} \]
其中 \(f\) 是當前頻率.

不過這個鬼東西是逆時針纏繞的, 咱們認爲這不是很清真(實際上是由於相位的關係), 因而咱們把它換成順時針:
\[ F(t)=g(t)e^{-2\pi ift} \]
至於重心, 咱們能夠考慮"採樣". 在纏繞的圖像上選取若干個點代數相加

其中 \(N\) 就是採樣點的個數, 咱們將它取個平均值就是重心位置了.

當咱們採樣點愈來愈多的時候:

咱們即可以用微分的思想, 直接把它積出來

因而咱們就獲得了一個假傅里葉變換的表達式
\[ \frac 1 {t_2-t_1}\int_{t_1}^{t_2}g(t)e^{-2\pi ift}\text{d}t \]
爲啥是假傅里葉變換呢?

真傅里葉變換是長這樣的
\[ \mathcal{F}(f)=\int_{-\infty}^{\infty}g(t)e^{-2\pi ift}\text{d}t \]
從物理意義上講, 當一個信號持續的時間越長, 你便越能肯定這個信號的頻率. 不知道你們有沒有這樣的經歷: 當你在等紅燈的時候, 碰巧看見兩輛車的轉向燈好像是同步的. 然而等了一會以後, 它們又變得不一樣步了. 因而咱們將前面的求平均值部分捨棄, 來表達咱們對這個成分頻率的肯定程度.

好比當信號很長的時候:

當咱們纏繞的頻率稍有誤差, 就會產生這樣的效果:

能夠看到, 波形很是迅速地比較均勻地分佈在了原點周圍, 這會讓咱們的變換波形上出現一個尖銳的峯.

然而這個東西的缺陷在於什麼呢?

它處理的是連續信號. 也就是說你得有個平滑的可積的函數表達式才能讓計算機處理.

可天然界哪來那麼多平滑函數? 尤爲是這種隨機性極大的信號?

因而咱們就有了...

離散傅里葉變換

顯然在計算機中, 咱們只能經過原波形上的一系列的點來搞事情 (由於原波形通常複雜的一匹). 這些點稱爲採樣.

實際上WC2018課上wys說採樣應該是單位衝激信號和原信號的連續卷積

因而咱們 ROLLBACK 到採樣點那一部分

而後把前面的 \(\frac 1 N\) 去掉, 因而就獲得了離散傅里葉變換(DFT)的表達式
\[ \mathcal{F}(f)=\sum_{k=1}^Ng(t_k)e^{-2\pi i f t_k} \]

FFT 與多項式

在上一節中咱們獲得了DFT的表達式, 而後咱們來分(魔)析(改)一下這個式子.

首先對於一個正經的採樣來講, 它的時間間隔必定是相等的. 其次對於離散傅里葉變換來講, 這個時間具體是多少其實已經沒有任何卵用. 並且對於離散的採樣來講, 它的旋轉系數(或者每個採樣的旋轉角)只有最小時間間隔的整數倍纔有意義. 再加上咱們能夠把相位共軛一下而去掉這個式子中 \(e\) 的冪上的負號, 最後再把下標換成從 \(0\) 開始, 因而咱們就能夠獲得下面這個式子:
\[ f(x)=\sum_{k=0}^{N-1}g(k)\omega_n^{kx} \]
等等這個 \(\omega\) 是啥玩意?

\(n\) 次單位復根

\(n\) 次單位復根是知足 \(\omega^n=1\) 的複數 \(\omega\). 根據複數乘法的集合意義以及同餘, 咱們能夠得出: \(n\) 次單位根剛好有 \(n\) 個, 均勻分佈在以原點爲圓心, 半徑爲 \(1\) 的圓周上. 以下圖所示:

因而咱們能夠獲得對於 \(k=0,1,\dots,n-1\), \(n\) 次單位根爲 \(e^{2\pi i k / n}\).

其中 \(\omega_n=e^{2\pi i/n}\) 爲主 \(n\) 次單位根.

爲啥呢它是這麼分佈的呢?

首先咱們有這樣一個顯然的結論:
\[ \omega_n^n=\omega_n^0=1 \]
因而咱們就有了:
\[ a\equiv b \pmod n \Rightarrow \omega_n^a=\omega_n^b \]
因而咱們能夠把乘法轉化指數上的模 \(n\) 意義下加法, 顯然任意整數 \(k\) 乘上一個 \(n\) 在模 \(n\) 意義下都是 \(0\).

其實 \(n\) 次單位根集合 \(\varOmega_n\) 和複數乘法運算組成了一個有限羣 \((\varOmega_n,\times)\).

還有一個顯然的事情, 就是對於一個 \(n\) 次單位根 \(\omega_n^k\), 它在極座標系中的幅角就是 \(\frac k n\) 個周角. 由於它們均勻分佈在單位圓上. 而對於分數, 咱們是能夠進行約分的, 這也就引出了對後續FFT十分重要的結論:

消去引理, 折半引理與求和引理

直接引用算法導論中給出的證實:

引理(消去引理): 對任意整數 $n \geq 0, k\geq 0 $, 以及 \(d\geq 0\), 有:
\[ \omega_{dn}^{dk}=\omega_n^k \]
證實: 由單位複數根的定義式可直接推出引理, 由於
\[ \omega_{dn}^{dk}=(e^{2\pi i / dn})^{dk}=(e^{2\pi i /n})^k= \omega_n^k \]

十分優雅而精簡的證實, 正是運用了單位複數根的特殊性質. 證實過程當中使用了冪的冪中指數相乘的運算法則, 將 \(d\) 乘進括號內, 與括號內在指數分母上的另外一個 \(d\) 相抵消.

根據消去引理, 咱們還能夠推出另外一個引理:

引理(折半引理): 若是 \(n>0\) 爲偶數, 那麼 \(n\) 個 \(n\) 次單位複數根的平方的集合就是 \(n/2\) 個\(n/2\) 次單位複數根的集合.

證實: 根據消去引理, 對任意非負整數 \(k\) , 咱們有 \((\omega_n^k)^2=\omega_{n/2}^k\) . 注意, 若是對全部 \(n\) 次單位複數根進行平方, 那麼得到每一個 \(n/2\) 次單位根正好兩次. 由於:
\[ (\omega_n^{k+n/2})^2=\omega_n^{2k+n}=\omega_n^{2k}\omega_n^n=\omega_n^{2k}=(\omega_n^k)^2 \]
所以, \(\omega_n^{k+n/2}\) 與 \(\omega_n^k\) 平方相同.

\(\blacksquare\)

從消去引理推出的折半引理是後面FFT中使用的分治策略的正確性的關鍵, 由於折半引理保證了遞歸子問題的規模只是遞歸調用前的一半, 從而保證 \(O(n\log(n))\) 時間複雜度.

另外一個重要的引理是求和引理, 是後文中使用FFT快速插值的基礎.

引理(求和引理): 對於任意整數 \(n\geq 1\) 與不能被 \(n\) 整除的非負整數 \(k\) , 有:
\[ \sum_{j=0}^{n-1}(\omega_n^k)^j=0 \]
證實: 根據幾何級數的求和公式:
\[ \sum_{k=0}^nx^k=\frac{x^{n+1}-1}{x-1} \]

咱們能夠獲得以下推導:

\[ \sum_{j=0}^{n-1}(\omega_n^k)^j=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}=\frac{(\omega_n^n)^k-1}{\omega_n^k-1}=\frac{(1)^k-1}{\omega_n^k-1}=0 \]

定理中要求 \(k\) 不可被 \(n\) 整除, 而僅當 \(k\) 被 \(n\) 整除時有 \(\omega_n^k=1\), 因此保證分母不爲 \(0\) .

\(\blacksquare\)

從新定義

咱們把 \(g(k)\) 看作一個數列, 那麼它就是個數論函數它就是個向量. 那麼其實後面的和式能夠看作另外一個函數:
\[ F(x)=\sum_{k=0}^{N-1}g(k)x^k \]
而 \(f(x)\) 就成了:
\[ f(x)=F(\omega_n^x) \]
而這個 \(F\) 其實就是個多項式!

而咱們所作的, 其實就是求出了多項式 \(F\) 在全部的 \(n\) 次單位根上的取值!

然而求這個有啥用?

多項式的表示

多項式的係數表達

相信你們都知道對於一個次數界爲 \(n\) 的多項式 \(A(x)\) 能夠表示爲以下形式:

\[ A(x) = \sum _{i=0} ^{n-1} a_i x^i \]

這時, 多項式 \(A(x)\) 的係數所組成的向量 \(\vec{a}=(a_0,a_1,...a_{n-1})\) 是這個多項式的係數表達. (其實是列向量不是行向量OwO)

使用係數表達來表示多項式能夠進行一些方便的運算, 好比對其進行求值, 這時咱們能夠採用秦九昭算法來在 \(O(n)\) 時間複雜度內對多項式進行求值.

\[ A(x_0) = a_0+x_0(a_1+x_0(a_2+...+x_0(a_{n-2}+x_0 a_{n-1})...)) \]

可是當咱們想把兩個採用係數表達的多項式乘起來的時候, 恭喜你, 問題來了: 一個多項式中的每一個係數都要與另外一個多項式中的每一個係數相乘. 而後時間複雜度就變成了鬼畜的 \(\Theta(n^2)\) ...

也就是對於兩個多項式的係數表示 \(\vec{a}\) 和 \(\vec{b}\) 來講, 兩個多項式乘積的係數表達 \(\vec{c}\) 中的值的求值公式以下:

\[ c_i=\sum_{j=0}^i a_jb_{i-j} \]

然而多項式乘法和卷積在不少地方都須要快速的實現, 而對於多項式乘法來講, 另外一種表示方法彷佛更爲合適:

多項式的點值表達

首先咱們對點值表達進行定義:

一個次數界爲 \(n\) 的多項式 \(A(x)\) 的點值表達爲一個由 \(n\) 個點值對所組成的集合:

\[ \{(x_0,y_0),(x_1,y_1),(x_2,y_2),...,(x_{n-1},y_{n-1})\} \]

使得對於 \(k=0,1,2,...,n-1\) , 全部的 \(x_k\) 各不相同

(咱們能夠先暫且把 \(A(x)\) 看做一個函數)

其中 \(y_k=A(x_k)\) , 也就是說選取 \(n\) 個不一樣的值分別代入求值以後產生 \(n\) 個值, 就會獲得 \(n\) 個未知數的值與多項式值的點對. 這 \(n\) 個點對就是多項式的一個點值表達.

不難看出點值表達並不像多項式表達同樣具備惟一性, 由於 \(x_k\) 是沒有限制條件的.

而咱們剛剛DFT中所作的, 就是求出了多項式 \(F\) 的一個點值表達. 因而咱們就能夠把它拿來搞事情了. 而按照定義計算出這些點值須要 \(O(n^2)\) 的時間.

....

等等...? 怎麼仍是 \(O(n^2)\) ? 別急, 很快咱們就會發現, 若是咱們像某些喪病出題人同樣精心選擇數據, 咱們就能夠在 \(O(n\log(n))\) 的時間複雜度內完成這個運算.

而從多項式的點值表達轉換爲惟一的係數表達就沒有那麼顯然了. 首先咱們介紹兩個概念:

求值與插值

從一個多項式的係數表達肯定其點值表達的過程稱爲求值(彷佛很顯然的樣子? 畢竟咱們求點值表達的過程就是取了 \(n\) 個 \(x\) 的值而後扔進了多項式求了 \(n\) 個值出來233), 而求值運算的逆運算(也就是從一個多項式的點值表達肯定多項式的係數表達)被稱爲插值. 而插值多項式的惟一性定理告訴咱們只有多項式的次數界等於已知的點值對的個數, 插值過程纔是明確的(能獲得一個肯定的多項式表達) , 聯想以前的矩陣與線性方程組和增廣矩陣能夠獲得以下證實:

定理(插值多項式的惟一性): 對於任意 \(n\) 個點值對組成的集合 \(\{(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1}\}\) ,且其中全部 \(x_k\) 不一樣, 則存在惟一的次數界爲 \(n\) 的多項式 \(A(x)\) , 知足 \(y_k=A(x_k),k=0,1,...,n-1\) .

證實 證實主要是根據某個矩陣存在逆矩陣. 多項式係數表達等價於下面的矩陣方程:
\[ \begin{bmatrix} 1 & x_0 & x_0^2 & \dots & x_0^{n-1} \\ 1 & x_1 & x_1^2 & \dots & x_1^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n-1} & x_{n-1}^2 & \dots & x_{n-1}^{n-1} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_{n-1} \end{bmatrix} \]

左邊的矩陣表示爲 \(V(x_0,x_1,...,x_{n-1})\) , 稱爲範德蒙矩陣. 而這個矩陣的行列式爲:

\[ \prod _{0\leq j < k \leq n-1} (x_k-x_j) \]

因此根據定理 "\(n\times n\) 矩陣是奇異的, 當且僅當該矩陣的行列式爲 \(0\)", 若是 \(x_k\) 皆不相同, 則這個矩陣是可逆的. 所以給定點值表達 \(\vec y\) , 則能肯定惟一的係數表達 \(\vec{a}\) 使:

\[ \vec{a}=V(x_0,x_1,...,x_{n-1})^{-1}\vec{y} \]

\(\blacksquare\)

上面的證實過程實際上也給出了插值的一種算法, 即求出範德蒙矩陣的逆矩陣. 咱們能夠在 \(O(n^3)\) 時間複雜度內求出逆矩陣因此相應地能夠在 \(O(n^3)\) 時間複雜度內計算出點值表達所對應的係數表達.

然而這特麼比求值還慢...=_=

所幸咱們還有一種基於 \(n\) 個點的插值算法, 基於拉格朗日公式:

\[ A(x)=\sum_{k=0}^{n-1} y_k \frac{\prod _{j\neq k} (x-x_j)} {\prod_{j\neq k} (x_k-x_j)} \]

基於這個公式咱們能夠在 \(O(n^2)\) 時間複雜度內計算出點值表達所對應的係數表達.

而後插值勉強追上了求值的時間複雜度.

那麼問題來了, 計算這個點值表達有啥用?

點值表達的運算

點值表達在不少多項式相關操做上比係數表達要方便不少, 好比對於加法, 若是 \(C(x)=A(x)+B(x)\) ,則對於任意點值 \(x_k\) , 知足 \(C(x_k)=A(x_k)+b(x_k)\) . 而正確性是顯而易見的.因此對於兩個點值表達的多項式進行加法操做, 時間複雜度爲 \(O(n)\).

與加法相似, 多項式乘法在點值表達的狀況下也變得更加方便: 若是 \(C(x)=A(x)B(x)\), 則對於任意點 \(x_k\) , \(C(x_k)=A(x_k)B(x_k)\) . 也就是說點值表達經常能夠大大簡化多項式間的運算!

因而這就是咱們拿它來卷積的理由了.

然而如今它是 \(O(n^2)\) 的, 咱們如今來優化它.

快速傅里葉變換FFT

FFT (Fast Fourier Transform)利用了 \(n\) 次單位根的特殊性質, 讓咱們能夠在 \(O(n\log n)\) 的時間內計算出DFT.

​ 咱們下面統一這樣定義DFT:

​ 對於 \(A\) 的係數表達 \(\vec{a}= (a_0,a_1,...,a_{n-1})\) 和 \(k=0,1,...,n-1\) , 定義結果 \(y_k\):
\[ y_k=A(\omega_n^k)=\sum_{i=0}^{n-1}a_i\omega^{ki}_n \]
​ 則 \(\vec{y}=(y_0,y_1,...,y_{n-1})\) 就是係數向量 \(\vec{a}\) 的離散傅里葉變換(DFT). 簡記爲 \(\vec{y}=DFT_n(\vec{a})\)

在下文中咱們爲了分治方便, 統一認爲 \(n\) 是 \(2\) 的整次冪. 雖然對於非整次冪的 \(n\) 的算法已經出現, 可是日常 OI 中使用的時候通常採起補一大坨值爲 \(0\) 的高次項係數強行補到 \(2\) 的整次冪23333

FFT使用的是分治策略. 根據係數下標的奇偶來拆分子問題, 將這些係數分配到兩個次數界爲 \(n/2\) 的多項式 \(A^{[0]}(x)\) 和 \(A^{[1]}(x)\) :

\[ A^{[0]}(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{n/2-1} \]

\[ A^{[1]}(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{n/2-1} \]

能夠發現 \(A^{[0]}(x)\) 中包含了 \(A\) 中下標爲偶數的係數, 而 \(A^{[1]}(x)\) 則包含了 \(A\) 中下標爲奇數的係數. 因此實際上咱們根據下標能夠獲得以下結論:

\[ A(x)=A^{[0]}(x^2)+A^{[1]}(x^2)x \]

由於次數界折半了因此代入的是 \(x^2\) , 而後對於奇數下標係數還要乘上一個 \(x\) 做爲補償.

因而咱們就把求 \(A(x)\) 在 \(\omega_n^0,\omega_n^1,...,\omega_n^{n-1}\) 處的值轉化成了這樣:

1.求多項式 \(A^{[0]}(x)\) 和 多項式 \(A^{[1]}(x)\) 在下列點上的取值:

\[ (\omega_n^0)^2,(\omega_n^1)^2,...(\omega_n^{n-1})^2 \]

2.根據剛剛推導出的將兩個按奇偶拆開的多項式的值合併的公式合併結果.

可是根據折半引理, 這 \(n\) 個要代入子多項式中求值的點顯然並非 \(n\) 個不一樣值, 而是 \(n/2\) 個 \(n/2\) 次單位複數根組成. 因此咱們遞歸來對次數界爲 \(n/2\) 的多項式 \(A^{[0]}(x)\) \(A^{[1]}(x)\) 在 \(n/2\) 個 \(n/2\) 次單位複數根處求值. 然咱們成功地縮小了了子問題的規模. 也就是說咱們將 \(DFT_n\) 的計算劃分爲兩個 \(DFT_{n/2}\) 來計算. 這樣咱們就能夠計算出 \(n\) 個元素組成的向量 \(\vec a\) 的DFT了.

咱們結合遞歸FFT的C++實現來理解一下

void FFT(Complex* a,int len){
    if(len==1)
        return;
    Complex* a0=new Complex[len/2];
    Complex* a1=new Complex[len/2];
    for(int i=0;i<len;i+=2){
        a0[i/2]=a[i];
        a1[i/2]=a[i+1];
    }
    FFT(a0,len/2);
    FFT(a1,len/2);
    Complex wn(cos(2*Pi/len),sin(2*Pi/len));
    Complex w(1,0);
    for(int i=0;i<(len/2);i++){
        a[i]=a0[i]+w*a1[i];
        a[i+len/2]=a0[i]-w*a1[i];
        w=w*wn;
    }
    delete[] a0;
    delete[] a1;
}

讓咱們逐行解讀一下這個板子:

第2,3行是遞歸邊界, 顯然一個元素的 \(DFT\) 就是它自己, 由於 \(\omega_1^0=1\), 因此

\[ y_0=a_0\omega_1^0=a_0 \]

第4~9行將多項式的係數向量按奇偶拆成兩部分.

第12,13和第17行結合起來用來更新 \(\omega\) 的值, 12行計算主 \(n\) 次單位複數根, 13行將變量 \(w\) 初始化爲 \(\omega_n^0\) 也就是 \(1\) , 17行使在組合兩個子多項式的答案時能夠避免從新計算 \(\omega_n\) 的冪, 而是採用遞推方式最大程度利用計算中間結果(OI中的經常使用技巧, 對於主 \(n\) 次單位複數根還能夠打個表來節約計算三角函數時的巨大耗時).

第10,11行將拆出的兩個多項式遞歸處理. 咱們設多項式 \(A^{[0]}(x)\) 的 \(DFT\) 結果爲 \(\vec{y^{[0]}}\) , \(A^{[1]}(x)\) 的 \(DFT\) 結果爲 \(\vec{y^{[1]}}\), 則根據定義, 對於 \(k=0,1,...,n/2-1\), 有:

\[ y_k^{[0]}=A^{[0]}(\omega_{n/2}^k) \]

\[ y_k^{[1]}=A^{[1]}(\omega_{n/2}^k) \]

而後根據消去引理, 咱們能夠獲得 \(\omega_{n/2}^k=\omega_n^{2k}\), 因此

\[ y_k^{[0]}=A^{[0]}(\omega_n^{2k}) \]

\[ y_k^{[1]}=A^{[1]}(\omega_n^{2k}) \]

設整個多項式的 \(DFT\) 結果爲 \(\vec y\), 則第15,16行用以下推導來將遞歸計算兩個多項式的計算結果綜合爲一個多項式的結果:

第15行, 對於 \(y_0,y_1,...,y_{n/2-1}\), 咱們設 \(k=0,1,...,n/2-1\) :
\[ \begin{aligned} y_k&=y^{[0]}_k+\omega_n^ky_k^{[1]} \\ &=A^{[0]}(\omega_n^{2k})+\omega_n^kA^{[1]}(\omega_n^{2k}) \\ &= A(\omega_n^k) \end{aligned} \]
上面的推導中, 第一行是咱們在代碼中進行的運算, 第二行將 \(y_k^{[0]}\) 和 \(y_k^{[1]}\) 按定義代入, 第三行基於咱們剛剛推導出的組合兩個子多項式 \(DFT\) 結果的式子.

第16行, 對於 \(y_{n/2}, y_{n/2+1}, ..., y_{n-1}\) , 咱們一樣設 \(k=0,1,...,n/2-1\) :
\[ \begin{aligned} y_{k+(n/2)} &= y^{[0]}_{k+(n/2)} - \omega_n^ky_{k+(n/2)}^{[1]} \\ &= y_k^{[0]} + \omega_n^{k+(n/2)} y_k^{[1]} \\ &= A^{[0]}(\omega_n^{2k}) +\omega_n^{k+(n/2)} A^{[1]}(\omega_n^{2k}) \\ &= A^{[0]}(\omega_n^{2k+n}) + \omega_n^{k+(n/2)} A^{[1]}(\omega_n^{2k+n}) \\ &= A(\omega_n^{k+(n/2)}) \end{aligned} \]
上面的推導中, 第一行一樣是咱們在代碼中進行的運算, 第二行從 \(\omega_n^{k+(n/2)}=-\omega_n^k\) 推導出, 第三行將 \(y_k^{[0]}\) 和 \(y_k^{[1]}\) 按定義代入, 第四行從 \(\omega_n^n=1\) 推導出, 最後一行基於推導出的組合遞歸結果用的等式.

根據這兩段推導, 咱們在代碼中所做的計算可以獲得正確的\(\vec y\).

其中咱們把 \(\omega_n^k\) 稱爲旋轉因子.

上述的算法實際上對應於下面的遞歸式:
\[ T(n)=2T(n/2)+\Theta(n)=\Theta(n\log(n)) \]
因此咱們講完了...個鬼

別忘了咱們只是求出了 DFT, 咱們算卷積的時候把它們乘起來還得再 IDFT 回去彩星. 這種時候咱們就得...

經過 FFT 在單位複數根處插值

根據咱們的四步快速卷積計算方案(倍次/求值/點積/插值), 咱們還剩下最後一步就能夠完成這一切了.

首先, \(DFT\) 的計算過程按照定義能夠表示成一個矩陣, 根據上文關於範德蒙德矩陣的等式 , 咱們能夠將 \(DFT\) 寫成矩陣乘積 \(\vec y=V_n \vec a\), 其中 \(V_n\) 是用 \(\omega_n\) 的必定冪次填充的矩陣, 這個矩陣大概長這樣:

\[ \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{bmatrix}= \begin{bmatrix} 1 & 1 & 1 & 1 & \dots & 1 \\ 1 & \omega_n & \omega_n^2 & \omega_n^ 3 & \dots & \omega_n^{n-1} \\ 1 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \dots & \omega_n^{2(n-1)} \\ 1 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \dots & \omega_n^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \dots & \omega_n^{(n-1)(n-1)} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{bmatrix} \]

觀察矩陣中 \(\omega_n\) 的冪次, 咱們彷佛發現 \(\omega_n\) 的指數彷佛組成了一張....乘法表???

由於咱們將 \(DFT\) 表示爲 \(\vec y=V_n \vec a\) , 則對於它的逆運算 \(\vec a=DFT^{-1}(\vec y)\) , 咱們能夠選擇直接求出 \(V_n\) 的逆矩陣 \(V^{-1}_n\) 並乘上 \(\vec y\) 來計算.

而後就會有一個玄學的定理來告訴咱們逆矩陣的形式, 算導上的證實以下:

定理: 對於 \(j,k=0,1,...,n-1\), \(V_n^{-1}\) 的 \((j,k)\) 處元素爲 \(\omega_n^{-kj}/n\).

證實: 咱們證實 \(V_n^{-1}V_n=I_n\), 其中 \(I_n\) 爲一個 \(n\times n\) 的單位矩陣. 考慮 \(V_n^{-1}V_n\) 中 \((j,j')\) 處的元素:

\[ [V_n^{-1}V_n]_{jj'}=\sum_{k=0}^{n-1}(\omega_n^{-kj}/n)(\omega_n^{kj'})=\sum_{k=0}^{n-1}\omega_n^{k(j'-j)}/n \]
當 \(j'=j\) 時, 此值爲 \(1\) , 不然根據求和引理, 此和爲 \(0\).

注意, 只有 \(-(n-1)\leq j'-j \leq n-1\) , 從而使得 \(j'-j\) 不能被 \(n\) 整除,才能應用求和引理.

\(\blacksquare\)

上面證實的思路是, 求出 \(V_n^{-1}V_n\) 的各元素的值, 並和單位矩陣 \(I_n\) 比較, 發現求出的值與單位矩陣相符, 原命題得證.

獲得這個定理以後咱們就能夠推導出 \(DFT^{-1}(\vec y)\) 了, 而咱們發現它是長這樣的:

\[ a_j=\frac{1}{n} \sum_{k=0}^{n-1} y_k\omega_n^{-kj} \]

咱們發現它和 \(DFT\) 的定義極其類似! 惟一的區別僅僅在於前面的多出的 \(\frac{1}{n}\) 和 \(\omega_n\) 指數上多出來的負號. 也就是說, 咱們能夠經過稍微修改一下 \(FFT\) 就能夠用來計算 \(DFT^{-1}\) 了!

最終咱們能夠獲得像這樣的代碼, 經過第三個參數指定計算 \(DFT\) 仍是 \(DFT^{-1}\) . 參數爲 \(1\) 則計算 \(DFT\) , 參數爲 \(-1\) 則計算 \(DFT^{-1}\). 除以 \(n\) 的過程在函數執行後在函數外進行便可.

void FFT(Complex* a,int len,int opt){
    if(len==1)
        return;
    Complex* a0=new Complex[len/2];
    Complex* a1=new Complex[len/2];
    for(int i=0;i<len;i+=2){
        a0[i/2]=a[i];
        a1[i/2]=a[i+1];
    }
    FFT(a0,len/2);
    FFT(a1,len/2);
    Complex wn(cos(2*Pi/len),opt*sin(2*Pi/len));
    Complex w(1,0);
    for(int i=0;i<(len/2);i++){
        a[i]=a0[i]+w*a1[i];
        a[i+len/2]=a0[i]-w*a1[i];
        w=w*wn;
    }
    delete[] a0;
    delete[] a1;
}

FFT的速度優化與迭代實現

看到上面的實現, 根據OIer的常識, 咱們顯然能夠意識到: 遞歸常數比較大. 並且上述實現中拆分多項式的時候還要從新動態分配一段空間顯然慢的一匹. 原本FFT各類浮點運算常數大如狗而後你還要遞歸那麼, 咱們是否能夠經過模擬遞歸時對係數向量的重排與操做來取得更好的執行效率呢?

襠燃能夠辣!

首先咱們能夠注意到, 在上面的實現代碼中計算了兩次 \(\omega_n^k y^{[1]}_k\), 咱們能夠將它只計算一次並將結果放在一個臨時變量中. 而後咱們的循環大概會變成這樣:

for(int i=0;i<(len/2);i++){
    int t=w*a1[i];
    a[i]=a0[i]+t;
    a[i+len/2]=a0[i]-t;
    w=w*wn;
}

咱們定義將旋轉因子與 \(y^{[1]}_k\) 相乘並存入臨時變量 \(t\) 中, 並從 \(y^{[0]}_k\) 中增長與減去 \(t\) 的操做爲一個蝴蝶操做.

接着咱們考慮如何將遞歸調用轉化爲迭代. 首先咱們觀察遞歸時的調用樹:

(沒錯這圖就是算導上的)

咱們能夠發現若是咱們自底向上觀察這棵樹, 若是將初始向量按照葉子的位置預先重排好的話, 徹底能夠自底向上一步一步合併結果. 具體的過程大概像這樣:

• 首先成對取出元素, 對於每對元素進行 \(1\) 次蝴蝶操做計算出它們的 \(DFT\) 並用它們的 \(DFT\) 替換原來的兩個元素, 這樣 \(\vec a\) 中就會存儲有 \(n/2\) 個二元素 \(DFT\) ;

• 接着繼續成對取出這 \(n/2\) 個 \(DFT\) , 對於每對 \(DFT\) 進行 \(2\) 次蝴蝶操做計算出包含 \(4\) 個元素的 \(DFT\) , 並用計算出的 \(DFT\) 替換掉原來的值.

• 重複上面的步驟, 直到計算出 \(2\) 個長度爲 \(n/2\) 的 \(DFT\) , 最後使用 \(n/2\) 次蝴蝶操做便可計算出整個向量的 \(DFT\).

實際實現時, 咱們發現計算該輪要進行蝴蝶操做的兩個元素的下標是比較方便的, 最外層循環維護當前已經計算的 \(DFT\) 長度, 每次循環完成後翻倍; 次級循環維護當前所在的 \(DFT\) 的開始位置的下標, 每次循環完成後加上 \(DFT\) 長度值; 最內層循環枚舉 \(DFT\) 內的偏移量. 若是三個循環的循環變量分別爲 \(i,j,k\), 則 \(j+k\) 指向當前要操做的第一個值, \(j+k+i\) 指向下一個 \(DFT\) 中須要計算的值.

而後關鍵的問題就在於如何快速重排獲得遞歸計算時的葉子順序.

咱們以長度爲 \(8\) 的 \(DFT\) 爲例, 讓咱們打表找規律看看能不能發現什麼:

原位置	0	1	2	3	4	5	6	7
原位置的二進制表示	000	001	010	011	100	101	110	111
重排後下標的二進制表示	000	100	010	110	001	101	011	111
重排結果	0	4	2	6	1	5	3	7

看起來彷佛...是把二進制表示給翻過來了?

這個過程在算導上彷佛叫作位逆序置換, 這個過程能夠直接遍歷一遍並交換二進制表示互爲逆序的各個下標所對應的值, 也能夠預處理出來保存在一個數組中來在屢次定長 \(FFT\) 中減小重複計算.

代碼實現相似這樣(我在這裏將位逆序置換預處理掉存在rev數組裏了)

void FFT(Complex* a,int len,int opt){
    for(int i=0;i<len;i++)
        if(i<rev[i])
            std::swap(a[i],a[rev[i]]);
    for(int i=1;i<len;i<<=1){
        Complex wn=Complex(cos(PI/i),opt*sin(PI/i));
        int step=i<<1;
        for(int j=0;j<len;j+=step){
            Complex w=Complex(1,0);
            for(int k=0;k<i;k++,w=w*wn){
                Complex x=a[j+k];
                Complex y=w*a[j+k+i];
                a[j+k]=x+y;
                a[j+k+i]=x-y;
            }
        }
    }
}

預處理部分大概這樣:

for(int i=0;i<bln;i++)
    rev[i]=(rev[i>>1]>>1)|((i&1)<<(bct-1));

其中bct是DFT長度的冪次.

實現中, 第一維枚舉當前分治區間的長度(由於遞歸的時候遞歸到底纔開始計算, 因此從小到大枚舉分治長度), 第二維枚舉分治區間的位置, 第三維在當前分治區間中進行蝴蝶操做.

上面的板子不包含最後除以DFT長度 \(n\) 的部分

而後咱們就成功地給 \(FFT\) 加了個速OwO

如今的實現已經能夠拿來作題了, 好比板子題 UOJ #34 多項式乘法

完整代碼實現:

#include <bits/stdc++.h>

const int MAXN=270000;
const int DFT=1;
const int IDFT=-1;
const double PI=acos(-1);

struct Complex{
    double real;
    double imag;
    Complex(double r=0,double i=0){
        this->real=r;
        this->imag=i;
    }
};
Complex operator+(const Complex& a,const Complex& b){
    return Complex(a.real+b.real,a.imag+b.imag);
}
Complex operator-(const Complex& a,const Complex& b){
    return Complex(a.real-b.real,a.imag-b.imag);
}
Complex operator*(const Complex& a,const Complex& b){
    return Complex(a.real*b.real-a.imag*b.imag,a.real*b.imag+a.imag*b.real);
}

Complex a[MAXN],b[MAXN],c[MAXN];

int n;             //Length of a
int m;             //Length of b
int bln=1;         //Binary Length
int bct;           //Bit Count
int len;           //Length Sum
int rev[MAXN];     //Binary Reverse Sort

void Initialize();
void FFT(Complex*,int,int);

int main(){
    Initialize();
    FFT(a,bln,DFT);
    FFT(b,bln,DFT);
    for(int i=0;i<=bln;i++){
        c[i]=a[i]*b[i];
    }
    FFT(c,bln,IDFT);
    for(int i=0;i<=len;i++){
        printf("%d ",int(c[i].real/bln+0.5));
    }
    putchar('\n');
    return 0;
}

void FFT(Complex* a,int len,int opt){
    for(int i=0;i<len;i++)
        if(i<rev[i])
            std::swap(a[i],a[rev[i]]);
    for(int i=1;i<len;i<<=1){
        Complex wn=Complex(cos(PI/i),opt*sin(PI/i));
        int step=i<<1;
        for(int j=0;j<len;j+=step){
            Complex w=Complex(1,0);
            for(int k=0;k<i;k++,w=w*wn){
                Complex x=a[j+k];
                Complex y=w*a[j+k+i];
                a[j+k]=x+y;
                a[j+k+i]=x-y;
            }
        }
    }
}

void Initialize(){
    scanf("%d%d",&n,&m);
    len=n+m;
    while(bln<=len){
        bct++;
        bln<<=1;
    }
    for(int i=0;i<bln;i++){
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bct-1));
    }
    for(int i=0;i<=n;i++){
        scanf("%lf",&a[i].real);
    }
    for(int i=0;i<=m;i++){
        scanf("%lf",&b[i].real);
    }
}

炸精現場與 NTT

有時候咱們會遇到這樣的東西:

• 求一個東西的方案數量, 計算過程當中發現須要快速卷積
• 求一個東西的機率, 輸出模一個大質數意義下的值, 計算過程當中發現須要快速卷積

這個時候咱們發現一個很是刺激的問題: FFT炸精了!

其次 FFT 全程對複數進行各類騷操做, 看着就慢的一匹, 能不能用整數代替啊?

固然能夠啦!

可是單位復根的性質並非很好在整數系統中體現...

咱們在發現 \(\omega_n\) 能夠簡化運算的時候主要用了啥性質?

別忘了它的定義是啥.

\(\omega_n^n=1\) 啊!

因此咱們要在整數中找知足 \(x^n=1\) 的非 \(1\) 值 \(x\).

等等好像只剩下 \(-1\) 了?

首先普通的乘法運算確定不會越乘絕對值越小, 因此咱們考慮怎麼能容許它越乘越小呢?

取模!

因而咱們就有了...

原根

階

定義: 若 \(m>1\) 且 \(\gcd(x,m)=1\), 咱們定義一個數 \(x\) 在剩餘系 \(m\) 下的階 \(\delta_m(x)\) 爲使得 \(x^r \equiv 1 \pmod m\) 的最小的整數 \(r\).

因而咱們有定理: 若 \(m>1\) 且 \(\gcd(x,m)=1\), \(x^n\equiv1\pmod m\), 則 \(\delta_m(x)|n\).

由於在 \(x\) 的全部冪次上, 第一次達到 \(1\) 即爲 \(\delta_m(x)\), 因而 \(n\) 顯然是 \(\delta_m(x)\) 的倍數.

這個階其實就是指 \(x^k\) 在模 \(m\) 意義下的全部取值的個數. 也就是說當 \(k\) 取 \([0,\delta_m(x))\) 中的全部整數的時候 \(x^k\) 在模 \(m\) 意義下互不相同. 大概能夠這樣證實:

• 設 \(a,b\in[0,\delta_m(x))\), \(a<b\), \(x^a\equiv x^b\pmod m\), 則咱們必定有 \(x^{b-a}\equiv1\pmod m\), 這與 \(\delta_m(x)\) 的定義矛盾.

原根的定義

設 \(m\) 是正整數, \(g\) 是整數, 若 \(\delta_m(g)=\varphi(m)\), 則稱 \(g\) 爲模 \(m\) 的一個原根.

原根的性質

• 對於正整數 \(m\), 當且僅當 \(m=2,4,p^k,2p^k\) 時 \(m\) 有原根. (\(p\) 爲奇質數)
• 對於模 \(m\) 的原根 \(g\), 知足 \(g^k\pmod m(0\le k<\varphi(m))\) 與全部與 \(m\) 互質的數一一對應

這兩個結論都挺顯然的(吧)

因而原根也構成了一個有限羣 \((G_m,\times_{\mathbb{Z}/m})\) . 其中 \(|G_m|=\varphi(m)\) .

原根的意義在於, 它能把 \([0,p-1)\) 這 \(p-1\) 個指數取不一樣的取值, 因而就成了一個最小週期爲 \(p-1\) 的環.

回想 \(\omega_n\) 的那張圖:

其實它的本質就是一個乘法環上的若干個不一樣值.

可是有一點不一樣!

\(g\) 能產生一個最小週期爲 \(p-1\) 的環, 那麼它最多隻能產生 \(p-1\) 以及 \(p-1\) 的約數次單位根. 想想, 爲啥?

咱們在求複數根的時候至關於在一個連續的圓周上均勻選了 \(n\) 個點, 但如今這個大圓不是連續的了. 它變成了 \(p-1\)個點, 而咱們要均勻選 \(n\) 個點的話一定要知足 \(n|(p-1)\).

說了這麼多, 那單位復根知足的三個引理它都知足嘛?

固然知足啦!

按照上面的定義, 咱們在乘法環上這 \(n\) 個點的方法就是令 \(n\) 次單位根爲 \(g^{(p-1)/n}\). 它的性質證實直接把 \(g^{(p-1)/n}\) 代進 \(\omega_n\) 裏就證完了.

原根的求法

其實某個數的最小原根通常都很小, 再加上模數通常是固定的因而本機暴力枚舉或者直接記住就能夠了

實際上, 對於一個質數, \(a^{p-1}\equiv 1\pmod p\) , 因此 \(\delta_p(a)|(p-1)\), 因而咱們只要枚舉 \(p-1\) 的全部質因子 \(p_k\) 看看 \(a^{\frac{p-1}{p_k}}\) 是否都不是 \(1\) 就能夠了.

擴展到模數爲合數的狀況只要用歐拉定理把 \(p-1\) 換成 \(\varphi(m)\) 就能夠了.

NTT

其實會了 FFT 以後基本就會 NTT ((Fast) Number Theory Transform) 了.

由於根據上面的各類推♂倒, 咱們只要用 \(g^{(p-1)/n}\) 代替 \(\omega_n\) 就能夠了.

可是由於咱們必需要令 \(n|(p-1)\), 而 \(n\) 由於後續分治的須要必須是 \(2\) 的整次冪, 因而咱們須要令 \(p-1\) 是 \(r2^k\) 的形式, 其中必須有 \(2^k\ge n\) .

然而在OI題目中能實現的範圍內, NTT長度不會超過 \(2^{20}\), 因而咱們只要記住幾個質數直接拿來用就星了

題目若是給定模數的話扔進計算器裏看看二進制表示也就知道它是否是NTT模數了.

好比最經常使用的 UOJ 模數 \(998244353\), 它的一個原根爲 \(3\).

它實際上等於 \(7\times17\times2^{23}+1\), 能夠作長達 \(8388608\) 的NTT(然而實際上你真的跑這麼長的 NTT 在 OI 時限內是跑不完的)

因而咱們就能夠獲得 NTT 的板子辣~

void NTT(int* a,int len,int opt){
    for(int i=0;i<len;i++){
        if(i<rev[i]){
            std::swap(a[i],a[rev[i]]);
        }
    }
    for(int i=1;i<len;i<<=1){
        int wn=Pow(G,(opt*((MOD-1)/(i<<1))+MOD-1)%(MOD-1),MOD);
        int step=i<<1;
        for(int j=0;j<len;j+=step){
            int w=1;
            for(int k=0;k<i;k++,w=(1ll*w*wn)%MOD){
                int x=a[j+k];
                int y=1ll*w*a[j+k+i]%MOD;
                a[j+k]=(x+y)%MOD;
                a[j+k+i]=(x-y+MOD)%MOD;
            }
        }
    }
    if(opt==-1){
        int r=Pow(len,MOD-2,MOD);
        for(int i=0;i<len;i++)
            a[i]=1ll*a[i]*r%MOD;
    }
}

理論上由於都是整數運算因此跑得比 FFT 要快

實際上上面這句理論是就是在扯淡.

由於瓦們有 SSE!

(SSE是一個現代CPU指令集, 全稱Streaming SIMD Extensions, 單指令多數據流擴展. 包含各類浮點運算指令/MMX整數運算強化指令/連續內存塊傳輸優化指令, 能夠理解爲這個指令集就是專門爲各類現代化圖形/科學計算設計的)

因而通常咱們有這樣的規律:

• 若是你用 std::complex, 那麼 FFT 通常都要慢於 NTT
• 若是你手寫 Complex 並且沒有開 -O2, 那麼 NTT 通常快於 FFT
• 若是你手寫 Complex 並且開了 -O2, 那麼 FFT 通常快於 NTT

除了 \(998244353\) 還有另外的一些模數, 好比:

• \(p=1004535809=479\times2^{21}+1\), \(g=3\)
• \(p=469762049=7\times2^{26}+1\), \(g=3\) (這個NTT長度更長, 可是數值有些偏小)

附一個(大概是miskcoo的) long long 之內的 NTT 模數表:

素數	\(r\)	\(k\)	\(g\)
\(3\)	\(1\)	\(1\)	\(2\)
\(5\)	\(1\)	\(2\)	\(2\)
\(17\)	\(1\)	\(4\)	\(3\)
\(97\)	\(3\)	\(5\)	\(5\)
\(193\)	\(3\)	\(6\)	\(5\)
\(257\)	\(1\)	\(8\)	\(3\)
\(7681\)	\(15\)	\(9\)	\(17\)
\(12289\)	\(3\)	\(12\)	\(11\)
\(40961\)	\(5\)	\(13\)	\(3\)
\(65537\)	\(1\)	\(16\)	\(3\)
\(786433\)	\(3\)	\(18\)	\(10\)
\(5767169\)	\(11\)	\(19\)	\(3\)
\(7340033\)	\(7\)	\(20\)	\(3\)
\(23068673\)	\(11\)	\(21\)	\(3\)
\(104857601\)	\(25\)	\(22\)	\(3\)
\(167772161\)	\(5\)	\(25\)	\(3\)
\(469762049\)	\(7\)	\(26\)	\(3\)
\(1004535809\)	\(479\)	\(21\)	\(3\)
\(2013265921\)	\(15\)	\(27\)	\(31\)
\(2281701377\)	\(17\)	\(27\)	\(3\)
\(3221225473\)	\(3\)	\(30\)	\(5\)
\(75161927681\)	\(35\)	\(31\)	\(3\)
\(77309411329\)	\(9\)	\(33\)	\(7\)
\(206158430209\)	\(3\)	\(36\)	\(22\)
\(2061584302081\)	\(15\)	\(37\)	\(7\)
\(2748779069441\)	\(5\)	\(39\)	\(3\)
\(6597069766657\)	\(3\)	\(41\)	\(5\)
\(39582418599937\)	\(9\)	\(42\)	\(5\)
\(79164837199873\)	\(9\)	\(43\)	\(5\)
\(263882790666241\)	\(15\)	\(44\)	\(7\)
\(1231453023109121\)	\(35\)	\(45\)	\(3\)
\(1337006139375617\)	\(19\)	\(46\)	\(3\)
\(3799912185593857\)	\(27\)	\(47\)	\(5\)
\(4222124650659841\)	\(15\)	\(48\)	\(19\)
\(7881299347898369\)	\(7\)	\(50\)	\(6\)
\(31525197391593473\)	\(7\)	\(52\)	\(3\)
\(180143985094819841\)	\(5\)	\(55\)	\(6\)
\(1945555039024054273\)	\(27\)	\(56\)	\(5\)
\(4179340454199820289\)	\(29\)	\(57\)	\(3\)

任意模數 NTT

剛剛咱們說到, NTT 有一個侷限就是模數必須能寫成 \(r2^k+1\) 的形式並且 \(2^k\ge n\).

那若是有毒瘤出題人搞了個題目模數爲 \(998244853\) (貓錕模數) 或者 \(99824353\) (meaty模數) 怎麼辦呢?

其實咱們能夠發現一個事情: 若是卷積的長度爲 \(n\), 模數是 \(m\), 那麼最終不取模的準確答案不會超過 \(nm^2\).

因此咱們找三個 NTT 模數分別作一遍 NTT 而後 CRT 一發就好啦!

如今泥萌知道我剛剛爲啥要多說幾個模數了吧

然而這樣搞出來的結果會爆 long long. 因此實現上仍是有點技巧的.

咱們要合併的式子是:
\[ \begin{cases} x\equiv a_1\pmod {p_1} \\ x\equiv a_2\pmod {p_2} \\ x\equiv a_3\pmod {p_3} \end{cases} \]

咱們首先在 long long 範圍內把前兩項合併, 獲得下面的東西:
\[ \left \{ \begin{aligned} x &\equiv A \pmod M \\ x &\equiv a_3 \pmod {p_3} \\ \end{aligned} \right. \]
而後咱們設最後的答案是:
\[ x=kM+A \]
且 \(k\) 要知足:
\[ kM+A\equiv a_3 \pmod{p_3} \]
那麼這個 \(k\) 必定能夠在模 \(p_3\) 意義下求出來:
\[ k \equiv (a_3 - A)M^{-1} \pmod {p_3} \]
其中 \(M^{-1}\) 是模 \(p_3\) 意義下的逆元.

因而咱們就能夠根據這個式子直接計算在對給定模數 \(p\) 取模的意義下的值了:
\[ x=kM+A \pmod{p} \]
好像很是清真?

然而並不

這個過程其實侷限也蠻多的.

首先模數 \(p\) 必須不能超過 \(10^9\) 數量級 (再大的話須要加模數, 加模數就得 NTT 更多遍, 並且最後 CRT 的時候可能沒法避免地要計算高精度取模...常數當場爆炸)

其次第一步合併兩個數的時候由於 \(M\) 是個 long long 範圍內的數因而須要快速乘...然而倍增快速乘常數好像很捉雞...因而咱們就須要這個板子了:

LL Mul(LL a,LL b,LL P){
    a=(a%P+P)%P,b=(b%P+P)%P;
    return ((a*b-(LL)((long double)a/P*b+1e-6)*P)%P+P)%P;
}

主要原理是取模以後的答案用 \(A\% B=A-\left\lfloor\frac A B\right\rfloor B\) 來算.

注意其中 \(A\) 和 \(\left\lfloor\frac A B\right\rfloor B\) 都是可能爆 long long 的.

然而由於是減法並且溢出以後通常是循環溢出因此做差以後(大概)不會鍋.

好像 \(2009\) 年集訓隊論文裏有這玩意?

相信這個板子就用它吧, 不信就倍增...

若是堅信有符號溢出是 UB 的話你能夠試試 unsigned

卷積狀物體與分治 FFT

有時候咱們會發現形如這樣的式子
\[ f(x)=\sum_{i=1}^xf(x-i)g(i) \]
它是卷積麼?

顯然它不是卷積. 可是它看上去確實長得很像卷積, 惟一不一樣的是 \(f\) 同時出如今了等式兩邊.

(注意求和下指標是從 \(1\) 開始的)

這式子怎麼求?

\(O(n^2)\) 的算法很是顯然, 根據定義式從小到大枚舉 \(x\) 就能夠計算了.

然而出題人才不會這麼放過你呢哼!

教練我會 FFT!

然而直接算由於會有本身轉移到本身這種操做, 因此咱們若是對於每一個 \(x\) 都算一遍 FFT 的話就變成 \(O(n^2\log n)\) 了...還不如 BFT 呢...

咱們考慮優化.

由於這是個求和式子, 其中每一項都是前面某一個已經求過的 \(f(x)\) 與另外一個相對比較固定的函數值的乘積. 咱們考慮將這些貢獻分別計算.

假設咱們已經知道了 \(f(1)\) 到 $ f(\frac n 2)$ 的值, 咱們要計算它對未知部分 \(f(\frac n 2+1)\) 到 \(f(n)\) 的貢獻. 咱們發現對於一個已知的 \(f(x)\) , 它對任意一個 \(i>x\) 的 \(f(i)\) 的貢獻都只有 \(f(x)g(i-x)\). 設未知部分的 \(f\) 值爲 \(0\) , 因而貢獻函數 \(F\) 就長這樣了:
\[ F(x)=\sum_{i=1}^xf(i)g(x-i) \]
好辣如今它是個卷積辣~

實際應用的時候, \(f\) 只須要填充 \([l,mid)\) 的部分, \(g\) 只須要填充 \([1,r-l]\) 的部分就能夠了 (由於其餘地方的 \(f\) 值都是 \(0\), 而貢獻只要算 \([mid,r)\) 的部分).

因而咱們CDQ分治一下, 先遞歸跑 \([l,mid)\), 而後算貢獻, 而後跑 \([mid,r)\)就完成了.

板子沒時間寫了, 竊了一份.

void cdq(int l,int r){
    if(l==r) return;
    int mid=l+r>>1;
    cdq(l,mid);

    int up=r-l-1;
//  for(lg=0,sze=1;sze<=(up<<1);sze<<=1)++lg;--lg;
    for(lg=0,sze=1;sze<=up;sze<<=1)++lg;--lg;//其實只用的到區間長度，實際卻有兩倍，因此保險用兩倍。
    for(int i=0;i<sze;i++) R[i]=(R[i>>1]>>1)|((i&1)<<lg),A[i]=B[i]=0;
//  memset(A,0,sizeof(A));memset(B,0,sizeof(B));
//  fill(A,A+sze,0);fill(B,B+sze,0);//這兩種方法賦值fill會比memset快一些（由於fill肯定了大小）
    for(int i=l;i<=mid;i++) A[i-l]=f[i];
    for(int i=1;i<=r-l;i++) B[i-1]=g[i];
    calc(A,B,sze);
    
    for(int i=mid+1;i<=r;i++) f[i]=(f[i]+A[i-l-1]%mod)%mod;
    cdq(mid+1,r);
}

calc 函數是計算 A 和 B 的卷積並把結果放在 A 裏.

FWT 與位運算卷積

咱們以前講到卷積, 說離散卷積能夠表示爲以下形式:
\[ c_k=\sum_{i\circ j=k}a_i b_j \]
咱們迷思一下: 當運算 \(\circ\) 爲位運算的時候是什麼玩意?

這即是位運算卷積了, 至關於下標系統是二進制數以及位運算的卷積.

回想咱們求多項式乘法卷積的時候, 搞了個大新聞構造了個卷積性變換 DFT, 而後把 \(A\) 和 \(B\) 都 DFT 掉而後求點積再 IDFT 回來就獲得了卷積.

因而咱們如今要構造一些新的卷積性變換 \(T\) , 知足 \(T(A)\cdot T(B)=T(A*_\oplus B)\). 並且還得有搭配的逆變換 \(T^{-1}\).

這個變換即是 Fast Walsh-Hadamard Transform.

FWT 與 \(\text{or}\) 卷積

構造變換 \(T\)

若 \(x\ \text{or}\ y=z\) , 那麼 \(x\) 和 \(y\) 必定知足它們的二進制表明的集合都是 \(z\) 的子集.

因而咱們把下標當作集合, 那麼咱們發現這個卷積對並集一致的下標一視同仁, 因而咱們能夠把全部子集都求和來讓它們能快速捲起來(乘法分配律你不會嘛?), 因而咱們嘗試如此定義變換 \(T\):
\[ T(A)_i=\sum_{j\subseteq i}A_j \]
那麼咱們就有了 \(T(A*B)=T(A)\cdot T(B)\).

爲啥?

咱們展開它:
\[ T(A*B)_i=\sum_{j\subseteq i}\sum_{a\cup b=j}A_aB_b \]

\[ (T(A)\cdot T(B))_i=\left (\sum_{a\subseteq i}A_a \right )\left(\sum_{b\subseteq i}B_b\right) \]

顯然它們都等價於:
\[ \sum_{(a\cup b)\subseteq i}A_aB_b \]
因而咱們便完成了第一步: 構造一個卷積性變換.

計算變換 \(T\)

直接按照定義式計算顯然又雙㕛叒叕變成 \(O(n^2)\) 的了...

咱們發現, 位運算有個特色就是每一位相互獨立. 因而咱們很是自信地選擇按位分治.

設咱們要變換的數列爲 \(A\), 則 \(A_0\) 爲 \(A\) 中下標的最高的二進制位爲 \(0\) 的值構成的數列, \(A_1\) 爲下標的最高二進制位爲 \(1\) 的值構成的數列. 實際上 \(A_0\) 就是前一半, \(A_1\) 就是後一半. 顯然當 \(A\) 長度爲 \(1\) 的時候變換就是它本身.

若是咱們知道了 \(T(A_0)\) 和 \(T(A_1)\) , 那是否是就能夠求 \(T(A)\) 了呢?

襠燃能夠辣!

假設當前 \(A\) 的長度爲 \(2^k\):

• \(A_0\) 至關因而第 \(k\) 位填 \(0\) (不選), 那麼子集確定仍是隻有 \(T(A_0)\) 裏的值.
• \(A_1\) 至關因而第 \(k\) 位填 \(1\) (選中), 那麼子集除了 \(T(A_1)\) 裏的值, 還有 \(T(A_0)\) 裏的值.

因而咱們就有了合併方法了.

設 \(A=(B,C)\) 表示把 \(B\) 和 \(C\) 首尾相接獲得的數列賦給 \(A\), 那麼:
\[ T(A)=(T(A_0),T(A_1)+T(A_0)) \]
其中 \(+\) 表示直接對位算術相加.

爲啥直接算術相加? 顯然分治的過程能夠保證對位的時候除了第 \(k\) 位外其餘的更低的位的選中狀況徹底一致.

其實這好像是個CDQ分治的過程, 遞歸分治而後算左半部分對右半部分的貢獻.

因而咱們就在 \(O(n\log n)\) 的時間複雜度內計算出 \(T(A)\) 辣~

計算逆變換 \(T^{-1}\)

其實有了變換的分治算法, 逆變換算法就直接解出來了.

咱們只要遞歸地把左邊對右邊的貢獻減掉就行了啊

因而式子就長這樣了:
\[ T^{-1}(A)=(T^{-1}(A_0),T^{-1}(A_1)-T^{-1}(A_0)) \]
應該沒啥很差理解的(吧)

迭代實現

咱們發現這一整個分治過程都和FFT極其類似, 那麼咱們能不能像FFT同樣迭代實現它呢?

固然能夠啦!

並且此次由於 \(A_0\) 就是前半部分 \(A_1\) 就是後半部分因此甚至根本就沒有重排過程.

因而咱們就能夠寫出這樣的板子:

void FWT(LL *a,int opt){
    for(int i=2;i<=N;i<<=1)
        for(int p=i>>1,j=0;j<N;j+=i)
            for(int k=j;k<j+p;++k)
                a[k+p]+=a[k]*opt;
}

第一維枚舉當前分治長度, 第二維枚舉當前分治區間的位置, 第三維計算貢獻. 跟 FFT 的迭代實現極爲類似.

而後就沒了.

FWT 與 \(\text{and}\) 卷積

構造變換 \(T\)

若 \(x\ \text{and}\ y=z\), 那麼 \(x\) 和 \(y\) 必定都知足它們的二進制表明的集合都是 \(z\) 的超集.

因而咱們依葫蘆畫瓢, 按照 \(\text{or}\) 卷積的形式定義變換 \(T\):
\[ T(A)_i=\sum_{j\supseteq i}A_j \]
那麼它確定也知足 \(T(A*B)=T(A)\cdot T(B)\) 了.

計算變換 \(T\)

按照 \(\text{or}\) 卷積的基本思路, 咱們一樣把 \(A\) 分爲兩部分 \(A_0\) 和 \(A_1\) , 而後計算它們之間的影響.

顯然咱們能夠獲得下面的計算式:
\[ T(A)=(T(A_0)+T(A_1),T(A_1)) \]
爲啥這樣應該不用我多說了吧

• \(A_0\) 至關因而第 \(k\) 位填 \(0\) (不選), 那麼超集除了 \(T(A_0)\) 裏的值, 還有 \(T(A_1)\) 裏的值.
• \(A_1\) 至關因而第 \(k\) 位填 \(1\) (選中), 那麼超集確定仍是隻有 \(T(A_1)\) 裏的值.

計算逆變換 \(T^{-1}\)

其實就是把左邊向右貢獻改爲右邊向左貢獻了吧
\[ T^{-1}(A)=(T^{-1}(A_0)-T^{-1}(A_1),T^{-1}(A_1)) \]
不理解的話大概須要 ROLLBACK 到 \(\text{or}\) 卷積的地方...

迭代實現

void FWT(LL *a,int opt){
    for(int i=2;i<=N;i<<=1)
        for(int p=i>>1,j=0;j<N;j+=i)
            for(int k=j;k<j+p;++k)
                a[k]+=a[k+p]*opt;
}

是的, 就是把 \(\text{or}\) 卷積的 a[k] 和 a[k+p] 給 std::swap 了一下.

FWT 與 \(\text{xor}\) 卷積

前面 \(\text{or}\) 卷積和 \(\text{and}\) 卷積都是小意思並且本質類似.

然而 \(\text{xor}\) 就開始鬼畜起來了... (雖然最後寫出來也沒長到哪去)

構造變換 \(T\)

在 \(\text{and}\) 和 \(\text{or}\) 卷積中, 咱們利用了一些集合操做成功地構造出了一個卷積性變換, 然而 \(\text{xor}\) 可就沒有這樣優秀的性質了...

咱們須要尋找一些像上面兩個卷積中那樣的性質.

這個時候咱們引入一個定理:

設 \(d(x)\) 表示 \(x\) 的二進制表示中 \(1\) 的個數的奇偶性, 則咱們有:
\[ d(i\,\text{and}\,k)\,\text{xor}\,d(j\,\text{and}\,k)=d((i\,\text{xor}\,j)\,\text{and}\,k) \]
理解一下?

其實首先裏面那個 \(\text{and}\,k\) 能夠直接扔掉無論, 咱們直接看其餘的.

這個定理的正確性咱們能夠理解爲, \(j\) 中每多包含一個 \(1\), 一定都會讓 \(d(i\,\text{xor}\,j)\) 的值翻轉一次.

因而咱們就能夠嘗試構造出長這樣的變換 \(T\):
\[ T(A)_i=\sum_{d(i\,\text{and}\,j)=0}A_j-\sum_{d(i\,\text{and}\,j)=1}A_j \]
因而咱們又有了 \(T(A*B)=T(A)\cdot T(B)\).

咱們展開來看看發生了什麼 (前方鬼畜預警):
\[ T(A*B)_i=\sum_{d(i\,\text{and}\,j)=0}\sum_{a\,\text{xor}\,b=j}A_aB_b-\sum_{d(i\,\text{and}\,j)=1}\sum_{a\,\text{xor}\,b=j}A_aB_b \\ \]

\[ \begin{aligned} (T(A)\cdot T(B))_i&= \left( \sum_{d(i\,\text{and}\,j)=0}A_j-\sum_{d(i\,\text{and}\,j)=1}A_j \right) \left( \sum_{d(i\,\text{and}\,j)=0}B_j-\sum_{d(i\,\text{and}\,j)=1}B_j \right) \\ &= \left( \left( \sum_{d(i\,\text{and}\,j)=0}A_j \right) \left( \sum_{d(i\,\text{and}\,j)=0}B_j \right) \right)+ \left( \left( \sum_{d(i\,\text{and}\,j)=1}A_j \right) \left( \sum_{d(i\,\text{and}\,j)=1}B_j \right) \right)- \left( \left( \sum_{d(i\,\text{and}\,j)=0}A_j \right) \left( \sum_{d(i\,\text{and}\,j)=1}B_j \right) \right)- \left( \left( \sum_{d(i\,\text{and}\,j)=1}A_j \right) \left( \sum_{d(i\,\text{and}\,j)=0}B_j \right) \right)\\ &= \left( \left( \sum_{d(i\,\text{and}\,j)=0}A_j \right) \left( \sum_{d(i\,\text{and}\,j)=0}B_j \right)+ \left( \sum_{d(i\,\text{and}\,j)=1}A_j \right) \left( \sum_{d(i\,\text{and}\,j)=1}B_j \right) \right)- \left( \left( \sum_{d(i\,\text{and}\,j)=0}A_j \right) \left( \sum_{d(i\,\text{and}\,j)=1}B_j \right)+ \left( \sum_{d(i\,\text{and}\,j)=1}A_j \right) \left( \sum_{d(i\,\text{and}\,j)=0}B_j \right) \right)\\ &= \left( \sum_{d(i\,\text{and}\,a)\,\text{xor}\,d(i\,\text{and}\,b)=0}A_aB_b \right)- \left( \sum_{d(i\,\text{and}\,a)\,\text{xor}\,d(i\,\text{and}\,b)=1}A_aB_b \right) \end{aligned} \]

而實際上根據上面的定理, 咱們有:
\[ \sum_{d(i\,\text{and}\,a)\,\text{xor}\,d(i\,\text{and}\,b)=0}A_aB_b =\sum_{d(i\,\text{and}\,(a\,\text{xor}\,b))=0} A_aB_b =\sum_{d(i\,\text{and}\,j)=0}\sum_{a\,\text{xor}\,b=j}A_aB_b \\ \sum_{d(i\,\text{and}\,a)\,\text{xor}\,d(i\,\text{and}\,b)=1}A_aB_b =\sum_{d(i\,\text{and}\,(a\,\text{xor}\,b))=1} A_aB_b =\sum_{d(i\,\text{and}\,j)=1}\sum_{a\,\text{xor}\,b=j}A_aB_b \]
因而咱們莫名其妙地又完成了第一步.

計算變換 \(T\)

咱們仍是照樣按位分治, 分別考慮最高位取 \(0\) 或 \(1\) 的狀況.

但此次與 \(\text{or}\) 和 \(\text{and}\) 不一樣的地方在於, 此次咱們分治兩側的區間會相互產生貢獻了.

(不過還好它們沒有變成環)

首先咱們假設咱們不考慮第 \(k\) 位的貢獻會發生什麼.

由於不考慮最高位, \(T(A_0)\) 和 \(T(A_1)\) 的對應位置的求和指標在對應的分治區間內的相對位置都是徹底同樣的.

• 當最高位是 \(0\) 的時候, \(\text{and}\) 運算的結果並不產生影響, 因此咱們直接把 \(T(A_0)\) 和 \(T(A_1)\) 的相對位置加起來就能夠了.

• 當最高位是 \(1\) 的時候, \(\text{and}\) 運算的結果可能會產生區別. 左半部分由於下標最高位是 \(0\) 因而直接產生貢獻, 右半部分下標最最高位是 \(1\), 因而多了一個翻轉奇偶性的貢獻, 翻轉奇偶性以後根據定義式, 兩個求和號的內容互換, 至關於交換符號, 因而咱們把它取反加在對應位置.

因而總的計算式就是:
\[ T(A)=(T(A_0)+T(A_1),T(A_0)-T(A_1)) \]

計算逆變換 \(T^{-1}\)

話說泥萌再看看這個變換計算式, 它是否是有點眼熟?
\[ T(A)=(T(A_0)+T(A_1),T(A_0)-T(A_1)) \]
這個 \(+/-\) 的過程是否是有點像...Manhattan距離和Chebyshev距離互化的座標系變換過程?

那個變換的逆變換是啥來着?

其實咱們直接解就能夠了, \(A_0'\) 和 \(A_1'\) 分別求和或者求差就能夠消掉另外一項, 而後再除以 \(2\) 就能夠了.

因而咱們就有了 \(T^{-1}\):
\[ T^{-1}(A)=\left(\frac{T^{-1}(A_0)+T^{-1}(A_1)}2,\frac{T^{-1}(A_0)-T^{-1}(A_1)} 2\right) \]

迭代實現

跟前面一個思路, 不過逆變換要記得除以 \(2\).

void FWT(int *a,int opt){
    for(int i=2;i<=N;i<<=1)
        for(int p=i>>1,j=0;j<N;j+=i)
            for(int k=j;k<j+p;++k)
            {
                int x=a[k],y=a[k+p];
                a[k]=(x+y)%MOD;a[k+p]=(x-y+MOD)%MOD;
                if(opt==-1)a[k]=1ll*a[k]*inv2%MOD,a[k+p]=1ll*a[k+p]*inv2%MOD; 
            }
}

inv2 是預處理的 \(2\) 的逆元.

FMT 與子集卷積

有時候咱們發如今圖論問題中有時候會搞出這種東西:
\[ F(S)=\sum_{T\subseteq S}f(T)g(S-T) \]
其中 \(S\) 和 \(T\) 都是集合. 或者另外一種(大概)等價的表達式:
\[ F(S)=\sum_{A\cup B=S,A\cap B=\varnothing}f(A)g(B) \]
第一眼咱們大概會用 \(\text{or}\) 卷積來搞, 然而咱們會發現一個很是滑稽的事情: \(A\) 和 \(B\) 不能有交集, 而 \(\text{or}\) 卷積則沒有限制這一點. 咋辦?

首先你們可能會天真地覺得 FMT 就是直接幹這個的

然而實際上 FMT 和 \(\text{or}\) 卷積一毛同樣, 由於它的定義是這樣的:
\[ F(S)=\sum_{T\subseteq S}f(T) \]

網上的 FMT 計算方法是這樣的:

inline void FMT(int *a,int len,int opt) {
    for (int i = 1; i < len; i <<= 1)
        for (int j = 0; j < len; ++j)
            if (j & i) 
                a[j]+=a[j^i]*opt;
}

其實這個計算過程也和上面的 FWT 計算 \(\text{or}\) 卷積本質徹底同樣. 想想爲啥 (注意看 i 枚舉的是什麼)

那咱們要怎麼知足交集爲 \(\varnothing\) 的限制呢?

咱們加一維表示集合大小再作 FMT. 這樣的話式子就變成了:
\[ F_k(S)=\sum_{T\subseteq S,|T|=k}f(T) \]
而後咱們對位相乘的時候把全部集合大小與當前位置相等的都手動捲起來, 獲得:
\[ \begin{align} P_k(S)&=\sum_{i+j=k}F_i(S)G_j(S) \\ &=\sum_{i+j=k} \left( \sum_{T\subseteq S,|T|=i}f(T) \right) \left( \sum_{T\subseteq S,|T|=j}g(T) \right) \\ &=\sum_{A\subseteq S,B\subseteq S,|A|+|B|=k}f(A)g(B) \\ &=\sum_{A\subseteq S,B\subseteq S}f(A)g(B)[|A|+|B|=k] \end{align} \]
而後咱們再逆變換回去, 便獲得:
\[ \begin{align} p_k(S)&=\sum_{A\cup B=S} f(A)g(B)[|A|+|B|=k] \\ &=\sum_{A\cup B=S,|A|+|B|=k} f(A)g(B) \end{align} \]
爲啥逆變換回去以後是對的呢?

咱們在變換與逆變換的過程當中, 求和式中的每一項都是獨立的, 因而逆變換的時候 \(f(A)g(B)[|A|+|B|=k]\) 就會被看作一個總體對待(由於它們之間是相乘的, 是同一項).

衆所周知, 若是 \(|A\cup B|=|A|+|B|\) 的話, \(A\cap B=\varnothing\). 因此咱們只取 \(p_k(S)\) 中 \(|S|=k\) 的部分就是答案了.

總時間複雜度由於多了一維要再乘上一個集合大小, 轉換成 FMT 長度的話就是 \(O(n\log^2n)\).

卷積與字符串

首先講個笑話

去年打NEERC的時候發現的

從前, 有三個神仙組了個隊註冊了, 他們叫...叫啥來着...(忽然忘了)

後來, 又有三我的組了個隊註冊了, 他們叫作:

因而前面三位神仙就改了個隊名:

然而剛剛講了那麼多東西, FFT 和bitset怎麼扯上的關係呢?

其實它們功能的交集是一些字符串問題.

卷積與迴文序列

記得BZOJ裏有這麼一個題: 萬徑人蹤滅

首先咱們能夠獲得一個很是顯然的思路: 對於每一個迴文中心, 計算以它爲中心對稱的字符個數 \(k\), 因而它會對答案產生 \(2^k\) 的貢獻. 同時題目要求不容許連續的迴文串, 可是咱們發現連續的迴文串必定都在剛剛的過程當中算過了, 因而咱們跑一遍 Manacher 求出迴文子串個數在答案中減掉就行了.

這個過程的複雜度瓶頸在對於每一個迴文中心都計算對稱字符個數, 直接枚舉計算顯然是 \(O(n^2)\) 的. 考慮優化這個過程.

首先咱們發現字符集大小 \(\Sigma\) 只有 \(2\) , 咱們能夠考慮把每種字符產生的貢獻分開計算. 咱們考慮相對某個迴文中心對稱的下標有什麼共同點?

它們的下標值之和是固定的!

咱們能夠理解爲, 兩個位置的下標值的平均值就是對稱中心的位置. 而平均值中除以 \(2\) 的部分咱們能夠同時忽略.

那麼咱們設原字符串爲 \(s\) , \(f(x)\) 爲對稱中心 \(\frac x2\) 產生的貢獻, 那麼咱們就能夠獲得這樣的式子:
\[ f(x)=\sum_{i+j=x}[s_i=k][s_j=k] \]
通過 \(0.5\texttt{s}\) 的思考以後咱們很是滑稽地發現: 它是一個卷積的形式!

咱們把原串 \(s\) 轉化爲數列 \(S_i=[s_i=k]\), 而後用 FFT/NTT 計算出 \(S*S\) 的結果就能夠了!

因而問題完美解決!

卷積與通配符

考慮這樣一個問題:

給定一個字符串和一個模板串, 兩個串中有一些字符是通配符, 求模板串與字符串的全部匹配位置.

你們先順着上一節的套路撕烤一下?

首先咱們得想辦法匹配字符串.

KMP/ST/SAM等等狀態轉移自動機確定是不行的, 單次匹配卻是沒問題. 可是單次匹配你爲啥不直接暴力呢?

咱們考慮什麼初等計算能夠起到 "只要進行這個運算, 無論什麼值都會獲得一個一樣的值" 這樣的通配效果?

乘 \(0\) !

咱們能夠構造一個表達式, 讓它爲 \(0\) 的時候表示該位置匹配, 這樣咱們只要讓通配符出現的時候讓它與 \(0\) 相乘就可讓整個項表明匹配了, 也就是實現了通配.

比較顯然的操做是把對應位置字符值作差, 可是差來差去萬一有正有負結果消成 \(0\) 了不就假了?

瓦們能夠套個絕對值!

因而咱們設原串爲 \(S\), 模板串爲 \(P\), 能夠這樣表示一個位置是否匹配:
\[ M(x)=\sum_{i=0}^{\text{len}(P)-1}\left|S_{x+i}-P_i\right|\times S_{x+i}\times P_i \]
若某位置是通配符則字符值爲 \(0\) , 那麼顯然當且僅當 \(M(x)=0\) 的時候位置 \(x\) 匹配.

然而絕對值可不怎麼好算...

注意到咱們只是想看看最終答案是否是 \(0\), 那麼顯然咱們能夠經過平方的方式把絕對值去掉, 因而獲得:
\[ \begin{align} M(x)&=\sum_{i=0}^{\text{len}(P)-1}(S_{x+i}-P_i)^2S_{x+i}P_i \\ &=\sum_{i=0}^{\text{len}(P)-1}(S_{x+i}^2-2S_{x+i}P_i+P_i^2)S_{x+i}P_i \\ &=\sum_{i=0}^{\text{len}(P)-1}S_{x+i}^3P_i-2S^2_{x+i}P^2_i+S_{x+i}P_i^3 \\ &= \left( \sum_{i=0}^{\text{len}(P)-1}S_{x+i}^3P_i \right) - 2\left( \sum_{i=0}^{\text{len}(P)-1}S^2_{x+i}P^2_i \right) + \left( \sum_{i=0}^{\text{len}(P)-1}S_{x+i}P_i^3 \right) \end{align} \]
如今它變成了三個看上去像卷積的東西...

然而咱們很是滑稽地發現: 這個卷積中下標運算不是加法而是減法!

咋辦?

咱們只須要用下面這個小學生都知道的式子:
\[ a-b=a+(-b) \]
咱們想辦法把其中一個串的下標變成負的!

然而直接變成負的的話好像很是不清真. 咱們再加個 \(\text{len}\) 把它補成正的吧!

因而咱們就要搞一個這樣的變換來預處理:
\[ P'_i=P_{\text{len}(P)-i} \]
emmmmm...

這不就是™ std::reverse 了一下麼

(也就是說當你發現減法卷積的時候翻轉其中一個就能夠了)

而後就能夠搞事情了! 咱們預處理出 \(S_i^2,S_i^3,P_i^2,P_i^3\) 而後再 std::reverse 一下 \(P\), 捲起來看哪位是 \(0\) 就能夠了!

由於都是整數, NTT/FFT都可.

一點優化?

剛剛咱們的作法要算三個卷積作⑨遍FFT/NTT, 常數直接起飛.

若是原串和模板串都帶通配符的話並無別的辦法, 只能像個⑨同樣作⑨遍FFT. (講道理⑨怎麼會FFT呢)

(實際上咱們本質上是多項式乘法再加法的形式, 因此按照點值表示的運算法則算就能夠了, 總計 7 遍 FFT)

可是若是隻有模板串或只有原串帶通配符呢?

咱們能夠考慮複數!

剛剛咱們爲啥不能直接作差呢? 由於正負相消, 兩個不匹配位置可能一個 \(+k\) 一個 \(-k\) 形成假匹配.

若是是複數呢?

咱們按照慣例先說個結論: \(k\) 個模長爲 \(1\) 的複數相加, 它們的和爲 \(k\) 當且僅當這些複數都是 \(1\).

因而咱們把這個差搬到 \(e^i\) 的指數上就好辣!
\[ \begin{align} M(x)&=\sum_{k=0}^{\text{len}(P)-1}e^{\frac{i\pi(S_{x+k}-P_k)}\Sigma} \\ &=\sum_{k=0}^{\text{len}(P)-1}e^{\frac{i\pi S_{x+k}}\Sigma}e^{\frac{-i\pi P_k}\Sigma} \end{align} \]
爲啥?

假設出現了假匹配, 那麼咱們顯然會獲得這種東西:

那個圓就是單位圓.

而實際上只有下面這種東西咱們纔會被認爲是匹配:

顯然假匹配的實部是不夠長的.

因而咱們對於原串作變換 \(S_i=\exp\left(\frac{s_i\pi i}\Sigma\right)\), 模板串作變換 \(P_i=\exp\left(\frac{-p_i\pi i}\Sigma\right)\), 翻轉模板串後直接卷積就能夠了.

判斷某個位置匹配只須要看卷出來以後的實部大小是否和模板長度減去通配符數量一致就能夠了.

因而咱們只須要一次卷積三次 FFT 就解決了通配符匹配問題.

此次由於動用了複數因此並不能用 NTT 來搞了...

實際上常數依然不小, 因而...

(大霧

寫在最後

卷積性變換在 OI 圈裏的普及度愈來愈高了...並且確實頗有用也很好用. 畢竟卷積這種形式在天然界也蠻常見的.

尤爲是在以 998244353 出名的UOJ(霧

這篇博客原本是一份知識分享講稿, 結果因爲時間關係最後部分寫得比較倉促, 還請如今在看這篇博客的您諒解!

若是你以爲這篇博客寫得比較全面, 還請點擊右下角推薦一下qwq~ 若是你以爲博客中有錯誤或寫得晦澀的地方也歡迎在評論區指正或詢問~