前言
數學在\(\text{OI}\)中十分重要。其中大多都是數論。
什么是數論？
本文包含所有側邊目錄中呈現的內容。絕對豐富！！！
下面直奔主題。
整除、同余、裴蜀定理、輾轉相除法、擴展歐幾里得
群
\((G,\cdot)\)表示\(G\)對定義在集合\(G\)上的運算\(\cdot\)構成群，它需要滿足：
（1）封閉性：若\(a,b \in G\)，則存在唯一確定的\(c \in G\)，使得\(a\cdot b=c\)。
（2）結合律：對於\(a,b,c \in G : (a\cdot b)\cdot c = a\cdot (b\cdot c)\)
（3）存在單位元：存在\(e \in G\)，使得\(\forall a \in G : a\cdot e=e\cdot a\)，\(e\)被稱作單位元。
（4）存在逆元：\(\forall a \in G, \exist b \in G : a\cdot b=b\cdot a=e\)，\(b\)又寫作\(a^{-1}\)。
逆元
歐拉函數
歐拉定理及拓展形式
威爾遜定理
同余方程（組）、中國剩余定理（CRT、EXCRT）
階和原根
BSGS、exBSGS
剩余問題
解決形如\(x^k \equiv a \pmod p\)的問題。暴力找似乎代價太高了。我們先從簡單的入手。
二次剩余
若\(x^2 \equiv a \pmod p\)有解，則稱\(a\)是\(p\)的二次剩余。否則稱其為非二次剩余。
定義\(\text{Legendre}\)符號：
（長相奇怪，其實只是記號）
先討論\(p\)是奇素數的情況。此時有：
\(Proof\)：
我們證明了所有的二次剩余得到的結果都是\(1\)，所有的非二次剩余得到的結果都是\(0\)，而又不存在第三者，所以命題能成立。
然后有\(p\)是奇質數下，二次剩余和非二次剩余的數量一樣多。這個很顯然，讓\(1\)到\(p-1\)的平方在模\(p\)意義下構成集合，而我們有：
這個東西展開下就證完了，所以我們只要考慮\(1\)到\(\frac{p-1}{2}\)之間的平方即可。
若有\(x^2 \equiv y^2 \pmod p\)，\(x \neq y\)，則\(x^2-y^2 \equiv 0 \pmod p \Rightarrow p \mid (x-y)(x+y)\)，由於\(p\)是質數，所以有\(p \mid x-y\)或\(p \mid x+y\)而得到\(x \equiv y\)或\(-y\)，第一種情況不存在，第二種情況發現兩數一定不在同一個部分（即在\(\frac{p-1}{2}\)兩邊），不合假設。所以得到\(1\)到\(\frac{p-1}{2}\)之間，兩兩的平方不同，而這樣得到的二次剩余只有\(\frac{p-1}{2}\)個，非二次剩余也只有\(\frac{p-1}{2}\)個。所以得證。
Cipolla
該算法是用來解決二次剩余問題的，思想基於構造。注意下文解決的是\(x^2 \equiv n \pmod p\)，右邊的是\(n\)而不是\(a\)了。
如果\(n\)是二次非剩余，顯然無解；
如果\(n\)是二次剩余，我們隨機一個\(a\)，使得\(\left( \frac{a^2-n}{2} \right) = -1\)。由非二次剩余和二次剩余數量一樣多，所以有\(\frac{1}{2}\)的概率取到，而期望下\(2\)次即能得到。證明不重要，放在后面。定義\(\omega = \sqrt{a^2-n}\)，事實上\(\omega\)在整數域中沒有數存在，所以類似於定義虛數單位\(i\)一樣，存數也先存成\(a+b\omega\)這種形式，之后會消掉這個東西。
我們要\(x^2 \equiv a \pmod p\)，先要知道這幾個引理。
引理1：\((x+\omega)^p \equiv x^p + \omega^p \pmod p\)。
\(Proof\)：
引理2：\(\omega^p \equiv -\omega\)。
\(Proof\)：
結論：\(x=(a+\omega)^\frac{p+1}{2}\)
\(Proof\)：
那這個式子咋求？直接上復數。根據拉格朗日定理，最后虛部一定為\(0\)。不會拉格朗日定理？其實我也不會那我們反證下吧。
\(Proof\)：
還有一種情況就是\(p=2\)，此時它是質數但不是奇質數，所以特判掉即可。
如果\(x\)有解，其實存在兩組解，另一組解為\(p-x\)。
int T, N, P;
namespace Cipo {
#define rd() (1ll * rand() * RAND_MAX + rand()) // 隨機數可能范圍不夠
int n, P; // x^2=n (%p)
int a, t; // t = w^2 = a^2-n
struct comp { // 自定義復數類
int x, y;
comp operator * (const comp &rhs) const { // (x1,y1)(x2,y2)=(x1x2+y1y2t,x1y2+x2y1)，這是復數乘，注意取模
return (comp){(1ll*x*rhs.x + 1ll*y*rhs.y%P*t)%P, (1ll*x*rhs.y+1ll*y*rhs.x)%P};
}
};
comp qpow(comp a, int b) { // 復數類快速乘
comp res = (comp){1, 0};
for (comp k = a; b; k = k*k, b >>= 1)
if (b & 1) res = res*k;
return res;
}
int qpow(int a, int b) { // 普通快速乘
int res = 1;
for (register int k = a; b; k = (ll)k*k%P, b >>= 1)
if (b & 1) res = (ll)res*k%P;
return res;
}
int sqrt(int n, int P) { // 開根號
Cipo::n = n, Cipo::P = P;
if (qpow(n, (P-1)>>1) == P-1) return -1; // n不是二次剩余
while (a = rd() % P, qpow(t = (1ll*a*a-n+P)%P, (P-1)>>1) == 1); // 判斷a^2-n得到的是不是二次剩余，如果是重新隨機a
return qpow((comp){a, 1}, (P+1)>>1).x; // 返回結果
}
}
int main() {
srand(time(0));
T = read();
while (T--) {
N = read(), P = read(); N %= P;
if (!N) { printf("0\n"); continue; } // 如果n=0表明x=0(%p)
if (P == 2) { printf("1\n"); continue; } // p=2特判
int ans = Cipo::sqrt(N, P); // 開根
if (ans == -1) printf("Hola!\n"); // 無解
else printf("%d %d\n", min(ans, P-ans), max(ans, P-ans)); // 否則有兩解
}
return 0;
}
該算法的復雜度：\(\text{O}(\log p)\)。
附期望是\(2\)的證明：
K次剩余
定義同上，只不過把\(2\)換成了\(k\)，但無法套用上面的解法求解。這里仍然假定\(p\)是奇質數。
發現構造走不通了！沒關系，我們學習了原根！
眾所周知原根遍歷了所有與\(p\)互質的數，所以我們用原根替換它們。
比如說\(x^k \equiv a \pmod p\)，設\(g\)為原根，\(a=g^s\)，\(x=g^m\)，則原式變成了\((g^m)^k \equiv g^s \pmod p\)，也就是\(g^{km} \equiv g^s \pmod p\)，\(s\)可以通過\(\text{BSGS}\)解得；\(m\)未知。然后由階的性質得到\(km \equiv s \pmod{\varphi(p)}\)，也就是\(km \equiv s \pmod{p-1}\)，解\(m\)即可。
想要所有解？剛剛的同余方程的解在模\(p-1\)的意義下的所有取值即為所有解。事實上題目一般要求找出任意解，正常解同余方程即可。如果上面的任意一步都出了問題（比如說\(s\)無解，或者同余方程沒有解），那就沒有解了。
// 碼量驚人，博主寫了140+行。這里只放最核心部分
// 這里用K次剩余解決二次剩余問題
int solve(int k, int a, int P) { // x^k=a (%p)
int s = BSGS(g, a, P); // 找出s使得g^s=a
int d = gcd(k, P-1);
if (s == -1 || s % d) return -1; // BSGS無解或同余方程無解，說明該問題無解
return qpow(g, 1ll*(s/d)*inv(k/d, (P-1)/d)%P, P); // 解同余方程km=s(%phi(p))
}
int main() {
init(); // 初始化質數表
T = read();
while (T--) {
N = read(), P = read(); N %= P;
if (!N) { printf("0\n"); continue; }
if (P == 2) { printf("1\n"); continue; }
getG(P); // 得到模P的原根
int ans = solve(2, N, P); // 解決K次剩余就寫成solve(K, N, P)
if (ans == -1) printf("Hola!\n");
else printf("%d %d\n", min(ans, P-ans), max(ans, P-ans));
}
return 0;
}
Miller Rabin
如果我們想要快速判斷一個數是不是質數，你可能會想到\(\text{O}(\sqrt n)\)復雜度的算法，但這個算法太劣了；或者想到運用費馬小定理，因為當\(p\)為質數時有\(a^{p-1} \equiv 1 \pmod p\)，所以猜測逆命題也成立。我們通過隨機\(a\)判斷式子是否成立來判斷\(p\)是不是質數。如果不准確，一次中了不代表一定是質數，那么我們多測幾次，可以將其為合數的概率降到很低。但有一類特殊的合數能讓式子恆成立，無論你怎么選擇\(a\)，該同余式子一直成立，這類合數叫做\(\text{Carmichael}\)數，例如\(561\)。
那怎么辦呢？先介紹一下二次探測原理。
若\(p\)是質數，且有\(a^2 \equiv 1 \pmod p\)，那么\(a \equiv \pm 1 \pmod p\)。前文有類似的證明，也就是移項化成整除式，這里就不證了。
我們通過這個方法來驗證。
把\(p-1\)分解成\(2^k \times t\)。然后隨機一個\(a\)然后讓其自乘成\(a^2\)，結合二次探測原理來判斷，如果在某一時刻自乘出\(1\)，則前一次必須是\(1\)或\(p-1\)，否則違背了原理，這個一定是合數。測試\(t\)輪后就得到了\(a^{p-1}=a^{2^k \times t}\)。如果最后得到的數不是\(1\)，那么就違背了費馬小定理，\(p\)是合數。
據說測試\(n\)次的正確率為\(1-\left( \frac{1}{4} \right)^n\)。
int test[10] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29}; // 據說用這幾個數在OI的ll范圍內出錯率為0%
int MR(ll P) {
if (P == 1) return 0; // 毒瘤的1:他不是質數，但會卡后面
ll t = P-1; int k = 0;
while (!(t & 1)) k++, t >>= 1; // 分解出2
rep(i, 0, 9) {
if (P == test[i]) return 1; // test中的數都是質數。如果這個數出現過那它一定是質數
else if (test[i] > P) return 0; // 測試的數比模數大，說明測試失敗了
ll a = pow(test[i], t, P), nxt = a; // 初始化a=test[i]^t%P
rep(j, 1, k) {
nxt = Mult(a, a, P); // 計算nxt=a^2%P。由於模數可能很大，這里的Mult用的是gui速乘
if (nxt == 1 && a != P-1 && a != 1) return 0; // 如果違背了二次探測定理，測試失敗
a = nxt;
}
if (a != 1) return 0; // 違背了費馬小定理，測試失敗
}
return 1; // 通過了所有的測試
}
忽略快（gui）速乘的復雜度，該復雜度為\(\text{O}(\log P)\)。
Pollard Rho
\(\text{Pollard-Rho}\)是一個用來快速分解質因數的算法，它巧妙的運用了“生日悖論”這一現象將\(N\)分解，配合着\(\text{Miller-Rabin}\)可以將復雜度降為\(\text{O}(N^{\frac{1}{4}})\)。
描述這個算法，我們先從“生日悖論”說起。
班級里有\(50\)名學生，保證生日是隨機的，不存在暗箱操作。假定一年為\(365\)天，問出現兩個人生日相同的概率是多少？你可能脫口而出：\(\frac{50}{365}\)。
但事實上概率非常之高以至於幾乎肯定有兩個人的生日在同一天，這個概率\(P=1-\frac{365 \choose 50}{365^{50}} \approx 97.142\%\)。為啥直覺跟計算出來的差異如此之大？原因在於要保證兩兩生日都不等的概率占樣本空間的比重非常之小。換個問題：如果你走進一個\(50\)人的班級里，你與班級里同學生日相等的概率又是多少呢？這時就是\(\frac{50}{365}\)了。
然后有一個結論：假設一年有\(N\)天，出現兩個人生日在同一天的期望人數是\(\text{O}(\sqrt N)\)。這個證明博主不會，直接跳過。
進入到\(\text{Pollard-Rho}\)前，還要介紹一個算法中的一個隨機函數。
設隨機函數\(f(x)\)，在模\(P\)的意義下，滿足\(f(x)=x^2+c\)，\(c\)為隨機參數，那么構成了\([0,P)\rightarrow[0, P)\)下的映射，且有\(f(x)=f(x+P)\)。借助初始參數\(x_0\)構造隨機數列\(\{ x_i\}\)，滿足：\(\forall i>0 : x_i=f(x_{i-1})\)。顯然其在模\(N\)意義下必定會走向一個循環，因為定義域和值域大小有限。這個函數生成的數列還是很隨機的，然后根據根據生日悖論有期望環大小為\(\text{O}(\sqrt P)\)，原因是如果出現了之前出現過的數字，顯然就會走入循環，根據上面的結論，出現兩個相同數字的期望次數就為\(\text{O}(\sqrt P)\)。把\(x\)的一步一步映射畫出來圖，發現長得十分像\(\rho\)，算法由此得名。
對於將要被分解的數\(N\)，我們選取兩個數字\(s\)和\(t\)，如果有\((\mid s-t \mid, N) \neq 1,N\)，很顯然我們找到了一個因子，可以遞歸分解兩部分了。但是畢竟數字太多，如何合理地隨機\(s\)和\(t\)呢？就運用隨機函數。
如果我們隨機選取一些數字兩兩判斷發現時間復雜度和空間復雜度又會彈回去。所以我們采用“\(\text{Floyd}\)判圈法”：有兩個人在起點\(x_0\)處，在環上一個人一次走一步，另一個一次走兩步。當他們再次相遇的時候表明一個人至少走完了一圈。這時一個人的數字記為\(s\)，另一個記為\(t\)，再繼續如上判斷。運用此方法能更好判環。
期望下我們需要\(\text{O}(N^{\frac{1}{4}})\)找出來一個因子。假設\(P\)可被分解，有\(P=a_1a_2,a_1 \leq a_2\)，在模\(P\)意義下，生成出來的序列出現環的期望為\(\text{O}(\sqrt N)\)。對於環上的任意兩個數\(x'\)和\(x''\)，\((\mid x'-x'' \mid, P)=a_1\)，當且僅當：\(x' \equiv x'' \pmod {a_1}\)成立。而生成的數列\(\{x_i\}\)在模\(N\)下是隨機的，在模\(a_1\)下也是隨機的，再次運用結論得到出現\(x'\equiv x''\)的期望為\(\text{O}(\sqrt{a_1})\)，又因為有\(a_1 \leq \sqrt P\)，所以得到期望為\(\text{O}(N^{\frac{1}{4}})\)。
但是還有種情況：\(x'=x''\)，會導致\((\mid x' - x'' \mid, P)=P\)。這表明兩個人相遇，選取的是環上的同一個數。遇到這種情況表明我們隨機數選取的不好，要更換隨機參數\(c\)重新操作。至此最原始的\(\text{Pollard-Rho}\)就出來了。
分析復雜度時我們只關心最大的數分解的復雜度，分解后的數字繼續分解對總復雜度的貢獻不大。但要不斷求\(\gcd\)，總復雜度為\(\text{O}(N^{\frac{1}{4}}\log N)\)。有\(\log\)的復雜度不好看，考慮削掉。
為了減少\(\gcd\)的次數，我們將若干次得到的\(\mid s-t \mid\)的結果乘起來。\((a,b)=d\)時，\((ac \bmod b,b)\geq d\)一定成立。記乘積模\(N\)的值為\(prod\)，如果\((prod,N)=1\)，顯然每一次\(\mid s-t \mid\)都與\(N\)互質；如果\(=N\)，表明其中必定有因子或\(0\)；否則就是\(N\)的因子。此時我們退回來一步一步判斷，找出因子或者是發現相遇了。調整適當的次數連乘，均攤下來復雜度就變成了\(\text{O}(N^{\frac{1}{4}})\)。
// 下列程序段實現求N的最大質因子
inline ll f(ll x, ll c, ll p) { return inc(Mult(x, x, p), c, p); } // 隨機函數
ll Floyd(ll N) {
ll c = rand()&15, s = 1, t = f(s, c, N); // 隨機c，s=1表示x0=1，t表示從起點跳一步
for (;;) {
ll _s = s, _t = t, prod = 1;
for (int i = 1; i <= 128; i++)
prod = Mult(prod, abs(_s - _t), N), _s = f(_s, c, N), _t = f(f(_t, c, N), c, N); // 連續處理128步
ll d = gcd(prod, N); // d=(prod, N)
if (d == 1) { s = _s, t = _t; continue; } // 如果d=1，表示沒找到因子，再來一組
for (int i = 1; i <= 128; i++) {
if ((d = gcd(abs(s - t), N)) != 1) return d; // !=1返回
s = f(s, c, N), t = f(f(t, c, N), c, N); // 否則繼續
}
}
}
ll PR(ll N) { // Pollard-Rho主過程
if (MR(N)) return N; // 如果N是質數，直接返回
ll ans;
while ((ans = Floyd(N)) == N); // 找到的是同一個點，重新來
return max(PR(ans), PR(N / ans)); // 遞歸，返回最大值
}
int main() {
srand(time(0));
scanf("%d", &T);
while (T--) {
scanf("%lld", &N);
ll ans = PR(N); // 找最大的質因子
if (ans == N) printf("Prime\n"); else printf("%lld\n", ans);
}
return 0;
}
Möbius反演相關
博主曾經寫過一篇，可能不夠完整，在這里就放個鏈接了：Mobius反演學習
（准備填天坑）
前置知識
前文談了很多關於數論的知識，這里加上一些嚴謹的定義。
數論函數：定義域為\(\Z\)，值域\(\subseteq \Complex\)的函數。多數情況下值域也為\(\Z\)。
積性函數：設\(f\)為數論函數，\(\forall a,b:(a,b)=1\)，有\(f(ab)=f(a)f(b)\)恆成立，則稱\(f\)為積性函數。特殊的，\(\forall a,b\)，有\(f(ab)=f(a)f(b)\)，則稱\(f\)為完全積性函數，注意此時沒有互質的要求。
質因數分解：將任意正整數化成若干個質數的冪的乘積，也就是將\(n\)化成這樣的形式：\(n=p_1^{\alpha_1}p_2^{\alpha_a}\dotsb p_s^{\alpha_s}\)。
根據以上三條，由於\(f(n)=f(p_1^{\alpha_1})f(p_2^{\alpha_2})\dotsb f(p_s^{\alpha_s})\)，我們可以用線性篩在\(\text{O}(n)\)的時間內求解前\(n\)個積性函數的值。
提到這里，我們引入幾個常用的數論函數。
單位函數：\(\epsilon (n)=[n=1]=\begin{cases}1&n=1\\ 0&n\neq 1\end{cases}\)。\([\text{state}]\)表示\(\text{state}\)為真時結果為\(1\)，否則為\(0\)。是完全積性函數。
除數函數：\(\sigma_k(n)=\sum\limits_{d \mid n}d^k\)。\(k=0\)表示求其約數個數此時簡寫成\(d(n)\)；\(k=1\)表示求約數和，此時簡寫成\(\sigma(n)\)。不難證明，除數函數也是積性函數。
歐拉函數：就是前文講的\(\varphi(n)\)，也是積性函數。歐拉函數有一個重要的性質：\(n=\sum\limits_{d \mid n}\varphi(d)\)。
\(Proof\)：
冪函數：\(\text{Id}_k(n)=n^k\)。這里是數論中的冪函數。\(k=1\)時可以省略。
理論基礎貌似差不多了。
Dirichlet卷積
設\(f\)，\(g\)是數論函數，對於數論函數\(h\)如果\(\forall n\)，\(h(n)=\sum\limits_{d \mid n}f(d)g\left(\dfrac{n}{d}\right)\)或者\(h(n)=\sum\limits_{i\times j=n}f(i)g(j)\)（兩式意義一樣），稱\(h\)為\(f\)和\(g\)的\(\text{Dirichlet}\)卷積，簡記\(h=f*g\)。
\(\text{Dirichlet}\)卷積有以下性質。
交換律和結合律：\(f*g=g*f\)，\(f*(g*h)=(f*g)*h\)。這些直接拿定義證明即可。
\(f*\epsilon=f\)，即任何數論函數與單位函數的卷積為本身。
以及有非常重要的應用：\(\sigma_k(n)=\text{Id}_k * 1\)和\(Id=\varphi * 1\)。前面的是定義得來的，后面的是剛剛講的證明，這些都是卷積形式。
Möbius反演
提到\(\text{Möbius}\)反演之前，需要了解\(\bold{Möbius}\)函數。
\(\bold{Möbius}\)函數：\(\mu(n)=\begin{cases} (-1)^s &n=p_1p_2\dotsb p_s,s可以為0，此時為1\\0 &\exist p_i: p_i^2 \mid n \end{cases}\)。這里\(p_i\)都是質數。\(\mu\)是積性函數，同樣也可以用線性篩求解。
\(\mu\)有什么用？一個很重要的性質：\(\mu * 1 = \epsilon\)。
\(Proof\)：
然后定義\(\text{Möbius}\)變換：如果\(g=f*1\)，則稱\(g\)為\(f\)的\(\text{Möbius}\)變換，而\(f\)為\(g\)的\(\text{Möbius}\)逆變換。往往我們想運用逆變換，此時稱作\(\text{Möbius}\)反演。此時有：
\(Proof\)：
該定理稱為\(\bold{Möbius}\)反演定理。
數論函數篩法
泰勒公式
給一個函數\(f(x)\)和其在\(x_0\)處的取值，用下面的式子無限逼近\(f(x)\)：
相當於通過數學的手段不斷擬合。怎么擬合？
下文的定義可能不那么嚴謹，不影響后文應用。
拉格朗日中值定理：
這個定理很顯然，隨便證。（原命題好像更嚴謹，這里不會影響下文）
於是又有了有限增量定理：
所以我們有：
但是當\(\Delta x\)很大時往往\(a\)的影響會變大，我們需要更進一步。
考慮構造一個多項式：
我們不斷向函數\(f(x)\)擬合，也就是讓：
所以有：
得到表達式：
然后還有個誤差分析關於\(f(x)-P(x)\)。感覺有點超其實我不會所以就不管了。最后得到的誤差很小小於我們定的階數\(n\)，所以在\(n\)趨於\(+\infin\)下，命題成立。
運用這個公式的例子：
多項式
多項式\(A(x)=a_0+a_1x+a_2x^2+\dotsb\)在平常十分常見，其應用十分廣泛，如多項式乘法、生成函數等等。多項式與多項式之間可以運算。多項式的度數即為最高項的次數，記作\(\text{deg}_A\)。對於\(A(x)=\sum_{i=0}^n a_ix^i,a_n\neq 0\)有時直接稱\(A(x)\)為\(n\)次多項式
有兩個多項式\(A,B\)，定義多項式加法：\(A(x)+B(x)=\sum_{i=0}^{\max(\deg_A,\deg_B)}(A_i+B_i)x^i\)；
定義多項式乘法：\(A(x)B(x)=\sum_{i=0}^{\deg_A} \sum_{j=0}^{\deg_B} A_iB_jx^{i+j}=\sum_{i=0}^{\deg_A+\deg_B} \sum_{j+k=i}A_jB_kx^{j+k}\)。后面的變形在之后應用上會比較方便。
觀察發現多項式加法的復雜度為\(\text{O}(n)\)，乘法為\(\text{O}(n^2)\)。乘法的復雜度過高了，有沒有辦法讓復雜度降下來呢？
我們平常所采用的多項式表示方法為系數表示法。此外，還有點值表示法，由於有一個定理：\(n\)個點唯一確定最高次數為\(n-1\)的多項式。我們就用\(n+1\)個點表示一個\(n\)次多項式。也就是說對於一個n次多項式\(A(x)\)，找一組數\(\{x_0, x_1,\dotsb, x_n\}\)，用\(n+1\)個點\((x_0,A(x_0))\)、\((x_1,A(x_1))\dotsb(x_n,A(x_n))\)。采用這個方式有什么好處呢？對於多項式加法\(A(x)+B(x)\)，假定次數都為\(n\)，那么得到的\(n+1\)個點即為\((x_i,A(x_i)+B(x_i)),i\in [0,n]\)；同理對於多項式乘法\(A(x)B(x)\)，得到的\(n+1\)個點即為\((x_i,A(x_i)B(x_i))\)。這樣操作一次的復雜度為\(\text{O}(n)\)了，優點就在此。
但是我們還面臨着問題：就是系數表示和點值表示之間如何轉換？系數向點值的轉化很簡單，用一下秦九韶公式（自行翻閱數學必修三），這需要\(\text{O}(n^2)\)的復雜度；而點值向系數需要用到拉格朗日插值，同樣也需要\(\text{O}(n^2)\)。可是這樣的復雜度並不優秀，如果我們采用轉換表示方法的角度來加快運算，雖然運算的復雜度降低，但是轉換的復雜度仍然很高，即使我們采用后文的快速變換改進公式也不夠優秀，這樣似乎變得不可行，其實之后我們會使用一組特殊的點使得轉換的速度加快。對於一般的點而言只能運用拉格朗日插值還原成系數表示，下面先介紹這個方法。
拉格朗日插值
給你\(n\)個點\((x_i,y_i)\)（保證\(x\)互不相同），這樣能確定一個多項式\(f(x)\)。給你\(x\)，求出\(f(x)\)。該算法解決這樣的問題。
用中學的知識考慮：對於\(1\)個點，確定一條水平的線；對於\(2\)個點，可以解出一條直線。對於\(3\)個點，可以解出一條拋物線。但是隨着點數的增多，我們需要求解一個\(n+1\)元一次方程，而解方程要用到高斯消元（后文會提到），這個的復雜度是\(\text{O}(n^3)\)，不夠優秀。下面我們有巧妙的方法來構造使得復雜度為\(\text{O}(n^2)\)，加上多項式快速變換可以在\(\text{O}(n \log^2 n)\)的時間內還原出多項式，但這個的基礎就是要用到拉格朗日插值。
構造函數\(l_i(x)=\prod_{j \neq i}\frac{x-x_j}{x_i-x_j}\)。這樣有什么好處呢？當\(x=x_j,j=1,2 \dotsb,j \neq i\)時，\(l_i(x)=0\)；反之如果\(x=x_i\)，\(l_i(x)=1\)。
構造多項式\(f(x)=\sum_{i=1}^n y_i l_i(x)\)，則這個多項式就是我們要求的多項式。往往只要求其中一個點時的值，就把\(x\)代入求解，畢竟直接還原多項式復雜度比較高。
這為什么是對的？由於構造的函數一定經過這\(n\)個點，而且保證過\(n\)個點一定確定小於\(n\)次的多項式，而這個就是我們要的多項式。
最終式子：
求解這個式子的復雜度為\(\text{O}(n^2)\)。
for (int i = 1; i <= n; i++) {
int mul = y[i], div = 1; // mul表示y[i]和式子分母的乘積
for (int j = 1; j <= n; j++)
if (i != j) mul = 1ll*mul*dec(k, x[j])%P, div = 1ll*div*dec(x[i], x[j])%P; // 分母分子分別乘上對應的數
ans = (ans + 1ll*mul*qpow(div, P-2)) % P; // 算到ans中
}
優化在后面會見到。
快速傅里葉變換(FFT)
快速在哪呢？轉換的過程只需要\(\text{O}(n \log n)\)，也就是說多項式乘法的復雜度就會降至\(\text{O}(n \log n)\)。\(\text{FFT}\)分為\(\text{DFT}\)和\(\text{IDFT}\)，分別稱作離散傅里葉變換和逆離散傅里葉變換，這些變換是時域與頻域之間的轉換，不過我們並不關心。學習該算法前我們需要學習復數的相關知識。
復數
復數是為了完善數學體系而新增加的一類數，用\(\Complex\)表示這類數的集合。具體來說，引入了\(i=\sqrt{-1}\)，於是所有復數都可以寫成\(z=a+bi\)這種形式，其中\(a,b\in \R\)，\(a\)被稱作實部，\(b\)被稱作虛部。特殊的，當\(b=0\)時這類數稱為實數；當\(a=0\)時這類數稱為純虛數。
復數可以對應成向量。
首先復數有一些基本定義和運算：
對於復數\(z=a+bi\)，
定義復數的模\(\mid z \mid=\sqrt{a^2+b^2}\)。這個運用勾股定理即可得出。
定義復數的幅角：\(z\)在平面上對應的點與原點的連線上，如果\(x\)的正半軸逆時針旋轉\(\theta\)角度與連線重合，則模角為\(\theta\)。注意是逆時針方向，不能從順時針方向。
定義兩個復數\(x=a_x+b_xi\)，\(y=a_y+b_yi\)，復數的加法\(x+y=(a_x+a_y)+(b_x+b_y)i\)；減法同理。幾何意義就是向量的加減，遵循平行四邊形法則。
定義復數的乘法：\(xy=(a_xa_y-b_xb_y)+(a_xb_y+a_yb_x)i\)。推導運用上文對\(i\)的定義即可（\(i^2=-1\)）。除法類似，對下文不重要便略去。復數乘法在幾何上的意義非常有趣且很特別：對於乘出來的復數，其模為兩個乘數的模相乘，其幅角等於兩個乘數的幅角相加。
證明可以運用三角變換。同時我們得到了\(\mid z \mid = \mid x \mid \mid y \mid\)。
定義共軛復數\(\overline z\)表示與\(z\)實部相等，虛部相反，即\(\overline z=a-bi\)。有一個推導公式\(\overline{ab}=\overline a \cdot \overline b\)，據說這個公式能夠優化掉\(\text{FFT}\)的\(\frac{1}{2}\)的常數，博主不會。幾何意義是兩個復數乘法前后在平面上關於\(x\)軸對稱。
歐拉公式：\(e^{i\theta}=\cos \theta + i\sin \theta\)。證明要用到泰勒公式，直接展開即可。這個代數寫法能很大程度代替幾何上的描述。比如復數乘法就可以轉化成指數上角的相加了。這個公式很重要，后文將會用到。
由以上的知識我們可以進入\(\text{FFT}\)的學習了。
正文：DFT、IDFT和優化
DFT
\(\text{DFT}\)實現的是將系數表達轉變成點值表達，對於一個次數小於\(2^n\)的多項式，我們選擇特殊的\(2^n\)個點使得轉換的速度加快。
首先定義單位根\(\omega_n=e^{\frac{2\pi i}{n}}\)，其中\(n\)是\(2\)的冪。后文中默認\(n\)都為\(2\)的冪次，包括實現。相當於將\(1\)沿逆時針方向旋轉\(\frac{2\pi}{n}\)的弧度。假設\(n=8\)，如下圖。
其模為\(1\)，幅角為\(45^\circ\)。對於其的任意次冪如下。
實際上這個是\(x^n=1\)的所有解。
發現\(\omega_n\)到\(\omega_n^n\)八個復數正好平分了半徑為\(1\)的單位圓，且有\(\omega_n^k=\omega_n^{k \bmod n}\)。特別的，\(\omega_n^{kn}=\omega_n^0=1\)。
對於一個多項式\(A(x)\)，我們希望快速求出\(A(\omega_n^0)\)、\(A(\omega_n^1)\)、\(A(\omega_n^2)\dotsb A(\omega_n^{n-1})\)的值，這些共能組成\(n\)個點。
接着介紹幾個有用的定理。這些可以直接運用定義證明。
定理\(1\)：\(\omega_{dn}^{dk}=\omega_n^k\)。（消去引理）
\(Proof\)：
定理\(2\)：\((\omega_{n}^k)^2=(\omega_n^{k+n/2})^2\)。（折半引理）
\(Proof\)：
這個引理非常重要。
定理\(3\)：\(\sum_{i=0}^{n-1}\omega_n^i=0\)。（求和引理）
\(Proof\)：
這個在\(\text{IDFT}\)中會用到。
有了上面的東西，我們就可以考慮對多項式進行變形而運用分治求解。
對於一個多項式\(A(x)=\sum_{i=0}^{n-1}a_ix^i\)，按照\(i\)的奇偶分成兩組得到：
然后假設：
所以我們有：
這樣我們可以運用以上引理了。
\(A(x)\)在\(\omega_n^k\)和\(\omega_n^{k+n/2}\)處的取值：
因為折半引理和消去引理，\(\omega_n^k = -\omega_n^{k+n/2}\)，所以有：
這樣我們可以通過求出\(A^{[0]}(\omega_{n/2}^k)\)和\(A^{[1]}(\omega_{n/2}^k)\)來求出\(A(\omega_n^k)\)和\(A(\omega_n^{k+n/2})\)。而對於每一個長度為\(n\)的多項式，我們在遞歸后都需要\(\text{O}(n)\)的時間合並這一層的信息，所以總復雜度就是\(\text{T}(n)=2\text{T}(\frac{n}{2})+\text{O}(n)=\text{O}(n \log n)\)了。
// 代碼有些難懂，需要細細體會
void DFT(comp *a, int l, int r) { // 長度一定要為2的冪
if (l == r) return; // 區間只有一個元素退出
int mid = l+r>>1, len = mid-l+1; // mid表示區間中點，len表示每組的長度，對應了n/2
static comp t[maxn << 2]; // 用來存交換的東西
for (int i = 0; i < len<<1; i += 2) t[l+(i>>1)] = a[l+i]; // 按照奇偶分組。先是偶數放前面
for (int i = 1; i < len<<1; i += 2) t[l+len+(i>>1)] = a[l+i]; // 再是奇數
for (int i = l; i <= r; i++) a[i] = t[i]; // 拷貝回去
DFT(a, l, mid), DFT(a, mid+1, r); // 先遞歸算出左右區間的插值結果
comp Wn = (comp){cos(PI/len), sin(PI/len)}, W = (comp){1, 0}; // 計算單位根Wn，單位根的冪用W表示
for (int i = 0; i < len; i++)
t[l+i] = a[l+i] + W*a[l+len+i], t[l+len+i] = a[l+i] - W*a[l+len+i], W = W * Wn; // 算出l+i，l+len+i處的插值結果，對應的冪為i。
for (int i = l; i <= r; i++) a[i] = t[i]; // 將算出的插值結果復制替換系數
}
IDFT
\(\text{IDFT}\)是將點值表達轉成系數表達。先給出如下的定義。
定義\(M=\left[\begin{matrix}\omega_n^0 & \omega_n^0 & \omega_n^0 & \dotsb & \omega_n^0 \\ \omega_n^0 & \omega_n^1 & \omega_n^2 & \dotsb & \omega_n^{n-1} \\ \omega_n^0 & \omega_n^2 & \omega_n^4 & \dotsb & \omega_n^{n-2} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \omega_n^0 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \dotsb & \omega_n^1 \end{matrix}\right]\)，定義向量矩陣\(x=\left[\begin{matrix}a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_n \end{matrix}\right]\)，則\(\text{DFT}\)的作用就是求出\(y=Mx\)，而\(\text{IDFT}\)的作用就是利用\(y\)反求\(x\)。
仍然定義\(\overline{M}\)，\(\overline M_{ij}=\overline{M_{ij}}\)，表示\(\overline M_{ij}\)為對應的\(M_{ij}\)的共軛復數。
定理：\(\frac{1}{n}\overline My=x\)。
\(Proof\)：
void IDFT(comp *a, int l, int r) { // 跟DFT差不多
if (l == r) return;
int mid = l+r>>1, len = mid-l+1;
static comp t[maxn << 2];
for (int i = 0; i < len<<1; i += 2) t[l+(i>>1)] = a[l+i];
for (int i = 1; i < len<<1; i += 2) t[l+len+(i>>1)] = a[l+i];
for (int i = l; i <= r; i++) a[i] = t[i];
IDFT(a, l, mid), IDFT(a, mid+1, r);
comp Wn = (comp){cos(PI/len), -sin(PI/len)}, W = (comp){1, 0}; // 與上面的不同之處在於虛部有負號
for (int i = 0; i < len; i++)
t[l+i] = a[l+i] + W*a[l+len+i], t[l+len+i] = a[l+i] - W*a[l+len+i], W = W * Wn;
for (int i = l; i <= r; i++) a[i] = t[i];
}
...
for (int i = 0; i < N; i++) a[i].x /= N, a[i].y /= N; // 這句話把n除掉，要寫在外面
...
FFT的改進
\(\text{FFT}\)的核心上文已經闡述完了，但是發現這樣遞歸下來空間消耗很大，常數很大，因此我們需要優化。
首先發現\(\text{DFT}\)和\(\text{IDFT}\)幾乎沒有什么代碼上的差別，我們可以合並。
然后將遞歸寫法改成非遞歸寫法。可是我們發現輸入按單位根的冪是奇數還是偶數來分開算，從二進制上看，發現是最后一位是\(0\)的分一組，最后是\(1\)的分一組。接着考慮倒數第二位，繼續如上操作。當\(n=8\)，分組如下。
最下面一層順序是\(\{0,4,2,6,1,5,3,7\}\)，做\(\text{FFT}\)的過程中每層連續兩組合並起來，每組越合並越大，合並是將兩組的數從左往右每次取出一個排成的。第一層的單位根為\(\omega_8\)，第二層是\(\omega_4\)，第三層是\(\omega_2\)，而每組對應的都是\(A(x)\)在\(\omega_n^0\)到\(\omega_n^{n-1}\)處的取值，這樣一來計算就可以從底往上合並了。
但是如何得到最底下一層的順序呢？
把最底下一層二進制寫出來：\(\{000,100,010,110,001,101,011,111\}\)，再把下標\(0\)到\(7\)的二進制列出來：\(\{000,001,010,011,100,101,110,111\}\)。發現正好是將下標的二進制翻過來得到的序列！這個證明很簡單，顯然會先根據二進制末尾分組，然后再根據倒數第二位分組，以此類推，相當於下標的反串排序然后再翻過來，其實等於直接把下標翻過來（下標的反串構成的集合等於下標的集合）。
這樣就能改進\(\text{FFT}\)了。具體細節看代碼。
void FFT(comp *a, int n, int opt) { // n為長度，opt表示是DFT(1)還是IDFT(-1)
static int rev[maxn << 2];
for (int i = 0; i < n; i++)
if ((rev[i] = rev[i>>1]>>1|(i&1?n>>1:0)) > i) swap(a[i], a[rev[i]]); // 細細品味線性遞推出下標i對應的位置，通過二進制翻轉求得
for (int q = 1; q < n; q <<= 1) { // 考慮組的大小，從下往上合並
comp Wn = (comp){cos(PI/q), opt*sin(PI/q)}; // 處理單位根，opt=1和-1時不一樣
for (int p = 0; p < n; p += q<<1) {
comp W = (comp){1, 0}; // 單位根的冪
for (int i = 0; i < q; i++) {
comp tmp = a[p+q+i]*W; // 用tmp記錄可以無需原來的t數組，這個是蝴蝶操作，用來優化。
a[p+q+i] = a[p+i] - tmp;
a[p+i] = a[p+i] + tmp;
W = W * Wn;
}
}
}
if (opt == -1) for (int i = 0; i < n; i++) a[i].x /= n, a[i].y /= n; // 如果是IDFT要除n
}
// FFT(a, n, 1) 表示DFT
// FFT(a, n, -1) 表示IDFT
總結
\(\text{FFT}\)不僅僅能做多項式乘法，事實上，它以及后面的一系列變換都有一套模型：處理卷積問題。
定義運算\(\circ\)，則卷積為\(F_n=\sum_{i\circ j=n}a_ib_j\)，這些運算可以在\(\text{O}(n\log n)\)的時間求出\(F_0\dotsb F_{n-1}\)。在上面的問題中，運算符為\(+\)，多項式乘法就變成了求解\(F\)數組了。\(\text{diricklet}\)卷積對應運算符是\(\times\)。
拓展：分治FFT
（占坑）
// 洛谷P4721 【模板】分治 FFT
// 采用cdq分治解決，亦可采用多項式求逆
// 這個的核心思想就是求f[l,r]，在完成f[l,mid]的基礎上，將f[l,mid]對f[mid+1,r]的貢獻用fft掃出來，然后遞歸求解[mid+1,r]，此時一定有f[mid+1]解決。
// 前置知識：NTT
void divide(int l, int r) {
if (l == r) return;
int mid = l+r>>1;
divide(l, mid); // 先完成f[l,mid]
int len = r-l+1, N = 1;
while (N < len) N <<= 1;
memcpy(x, f+l, sizeof(int)*len); // 將f[l,r]拷到x
memcpy(y, g, sizeof(int)*len); // 將g[0,len]拷到y
for (int i = mid-l+1; i < N; i++) x[i] = 0; // 防止不必要的影響，將一些位置設成0
fft(x, N, 1), fft(y, N, 1);
for (int i = 0; i < N; i++) x[i] = 1ll*x[i]*y[i]%P; // 卷積
fft(x, N, -1);
for (int i = mid-l+1; i <= r-l; i++) f[l+i] = inc(f[l+i], x[i]); // 讓f[mid+1,r]加上f[l,mid]對它們的貢獻
divide(mid+1, r); // 再完成f[mid+1,r]
}
復雜度：\(\text{O}(n \log^2 n)\)。運用生成函數加上多項式求逆可以做到\(\text{O}(n \log n)\)。
快速數論變換(NTT)
就是將多項式乘法限制在模意義下，通常模數都是\(\text{NTT}\)模數，即形如\(P=s\times 2^k+1\)的數，首先這個數必須是質數，其次要滿足\(k\)盡量大，為什么呢？下面細細說來。
\(\text{NTT}\)與\(\text{FFT}\)十分相像：它也有單位根。單位根要滿足其\(n\)次冪為\(1\)，小於\(n\)次都不為\(1\)，而原根能很好地勝任，因為有\(g^{\varphi(P)}=1\)且冪為任意\(\varphi(P)\)的因子都不為\(1\)。\(P\)是質數，所以\(\varphi(P)=P-1\)。我們可以構造\(\omega_n=g^{\frac{P-1}{n}}\)，這樣\(\omega_n^n\equiv 1\)而且\(m<n,\omega_n^m \not\equiv 1\)。這樣就可以完全替換掉\(\text{FFT}\)里面的東西了。由於\(n\)是\(2\)的冪且可能很大，所以對\(P-1\)的要求就是有很多的因子\(2\)。
\(\text{NTT}\)的優點在於：能在特定模數意義下運算、沒有任何的精度丟失。但無法直接運用到任意模數\(\text{MTT}\)上。這樣的情況后面再討論。
代碼如下。
// 這里NTT模數P取的是998244353，g=3。當然也可以取別的合適的
void NTT(int *a, int n, int opt) {
static int rev[maxn << 2];
for (int i = 0; i < n; i++)
if ((rev[i] = rev[i>>1]>>1|(i&1?n>>1:0)) > i) swap(a[i], a[rev[i]]);
for (int q = 1; q < n; q <<= 1) {
int Wn = qpow(g, (P-1)/q>>1); // 這里用原根替換了復數
if (opt == -1) Wn = qpow(Wn, P-2); // 如果是DFT，要求g^(-(P-1)/2)，實際上就是求上面那個數的逆元
for (int p = 0; p < n; p += q<<1) {
int W = 1;
for (int i = 0; i < q; i++) {
int tmp = 1ll*a[p+q+i]*W%P;
a[p+q+i] = inc(a[p+i], P-tmp);
a[p+i] = inc(a[p+i], tmp);
W = 1ll*W*Wn%P;
}
}
}
if (opt == -1) for (int i = 0, inv = qpow(n, P-2); i < n; i++) a[i] = 1ll*a[i]*inv%P; // 這里除n變成乘逆元
}
常數也不小（主要是取模），但比\(\text{FFT}\)小。
快速沃爾什變換(FWT)
任意模數NTT(MTT)
具體原理詳見myy的集訓隊論文。（先咕了~）
#include <bits/stdc++.h>
#define rep(i, a, b) for (int i = a, i##end = b; i <= i##end; ++i)
#define per(i, a, b) for (int i = a, i##end = b; i >= i##end; --i)
#define rep0(i, a) for (int i = 0, i##end = a; i < i##end; ++i)
#define per0(i, a) for (int i = a-1; ~i; --i)
#define chkmin(a, b) a = std::min(a, b)
#define chkmax(a, b) a = std::max(a, b)
typedef long long ll;
const int maxn = 100000 + 5;
const double PI = acos(-1);
int P;
inline int read() {
int w = 0, f = 1; char c;
while (!isdigit(c = getchar())) c == '-' && (f = -1);
while (isdigit(c)) w = w*10+(c^48), c = getchar();
return w * f;
}
struct comp {
double x, y;
comp(double x=0, double y=0) : x(x), y(y) {}
};
comp operator + (comp a, comp b) { return comp(a.x+b.x, a.y+b.y); }
comp operator - (comp a, comp b) { return comp(a.x-b.x, a.y-b.y); }
comp operator * (comp a, comp b) { return comp(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x); }
comp conj(comp a) { return comp(a.x, -a.y); }
comp W[maxn << 2];
void prework(int n) {
for (int i = 1; i < n; i <<= 1)
for (int j = 0; j < i; j++) W[i+j] = comp(cos(j*PI/i), sin(j*PI/i));
}
void fft(comp *a, int n, int opt) {
static int rev[maxn << 2] = {0};
for (int i = 1; i < n; i++)
if ((rev[i] = rev[i>>1]>>1|(i&1?n>>1:0)) > i) std::swap(a[i], a[rev[i]]);
for (int q = 1; q < n; q <<= 1)
for (int p = 0; p < n; p += q<<1)
for (int i = 0; i < q; i++) {
comp t = a[p+q+i]*W[q+i]; a[p+q+i] = a[p+i]-t, a[p+i] = a[p+i]+t;
}
if (~opt) return;
for (int i = 0; i < n; i++) a[i].x /= n, a[i].y /= n;
std::reverse(a+1, a+n);
}
void conv(int *x, int *y, int *z, int n) {
for (int i = 0; i < n; i++) x[i] %= P, y[i] %= P;
static comp a[maxn << 2], b[maxn << 2];
static comp dfta[maxn << 2], dftb[maxn << 2], dftc[maxn << 2], dftd[maxn << 2];
for (int i = 0; i < n; i++)
a[i] = comp(x[i] >> 15, x[i] & 32767), b[i] = comp(y[i] >> 15, y[i] & 32767);
fft(a, n, 1), fft(b, n, 1);
for (int i = 0; i < n; i++) {
int j = (n-i) & (n-1);
static comp a1, a2, b1, b2;
a1 = (a[i] + conj(a[j])) * comp(0.5, 0);
a2 = (a[i] - conj(a[j])) * comp(0, -0.5);
b1 = (b[i] + conj(b[j])) * comp(0.5, 0);
b2 = (b[i] - conj(b[j])) * comp(0, -0.5);
dfta[i] = a1 * b1, dftb[i] = a1 * b2, dftc[i] = a2 * b1, dftd[i] = a2 * b2;
}
for (int i = 0; i < n; i++)
a[i] = dfta[i] + dftb[i] * comp(0, 1), b[i] = dftc[i] + dftd[i] * comp(0, 1);
fft(a, n, -1), fft(b, n, -1);
for (int i = 0; i < n; i++) {
int ax = (ll)(a[i].x + 0.5) % P, ay = (ll)(a[i].y + 0.5) % P, bx = (ll)(b[i].x + 0.5) % P, by = (ll)(b[i].y + 0.5) % P;
z[i] = (((ll)ax << 30) + ((ll)(ay + bx) << 15) + by) % P;
}
}
int getsz(int n) { int x = 1; while (x < n) x <<= 1; return x; }
int n, m, N, a[maxn << 2], b[maxn << 2], c[maxn << 2];
int main() {
n = read() + 1, m = read() + 1, P = read();
rep0(i, n) a[i] = read();
rep0(i, m) b[i] = read();
N = getsz(n+m-1); prework(N);
conv(a, b, c, N);
rep0(i, n+m-1) printf("%d ", c[i]);
return 0;
}
多項式基本操作
多項式求逆
給一個長度為\(n\)的多項式\(A\)（無關時\((x)\)略去），求\(B\)，滿足\(AB\equiv 1 \pmod {x^n}\)。\(n\)不一定是2的冪。
求解需要用倍增的方法。注意下文模數的變化。
首先我們有：
接着考慮構造一個解\(B'\)滿足：
所以有：
兩邊同時平方得：
這里模數能平方是因為同余\(0\)，否則不可以直接平方。
展開：
所以我們就有：
這樣不斷遞歸下去，當\(n=1\)時就是\(A\)常數項的逆元。在此基礎上倍增構造即可。復雜度\(\text{O}(n \log n)\)。
// 多項式以及一些支持的操作
struct poly {
int a[maxn << 1], len; // 各項系數和長度
poly(int a0=0) : len(1) { memset(a, 0, sizeof a); a[0] = a0; } // 初始化常數項，或者為0
int &operator [] (int i) { return a[i]; } // 定義下標訪問
void load(int *from, int n) { // 傳入
len = n;
memcpy(a, from, sizeof(int)*n);
}
void cpyto(int *to, int begin, int len) { // 傳出
memcpy(to, &a[begin], sizeof(int)*len);
}
void resize(int n=-1) { // 擴展或收縮長度
if (n == -1) { while (len > 1 && !a[len-1]) len--; return; } // 不傳入參數表明去掉最高位的0的真實長度
while (len < n) a[len++] = 0; // 保證新項為0
while (len > n) len--;
}
void print() { // 輸出
for (int i = 0; i < len; i++) printf("%d ", a[i]);
}
} A;
// 多項式求逆。自底向上實現
poly poly_inv(poly A) {
poly B = poly(qpow(A[0], P-2)); // 邊界
for (int len = 1; len < A.len; len <<= 1) { // 倍增
static int x[maxn << 2], y[maxn << 2];
for (int i = 0; i < len << 2; i++) x[i] = y[i] = 0; // 清空數組防止錯誤
A.cpyto(x, 0, len << 1); B.cpyto(y, 0, len); // 將信息拷到數組中。注意長度不能太長也不能過短
ntt(x, len << 2, 1), ntt(y, len << 2, 1);
for (int i = 0; i < len << 2; i++) x[i] = inc(y[i], inc(y[i], P - 1ll*x[i]*y[i]%P*y[i]%P)); // 求2B'-AB'^2
ntt(x, len << 2, -1);
B.load(x, len << 1); // 更新len*2位的值
}
B.resize(A.len); // 保證A和B長度相等
return B;
}
多項式開根
給定長度為\(n\)的多項式\(A\)，求項數在\(n\)以內的多項式\(B\)，使得\(B^2\equiv A \pmod{x^n}\)。
仍然采用倍增的方法。由於與求逆相似，步驟上的簡略參考上文。構造\(B'\)使得\(B'\equiv B \pmod{x^{\lceil \frac{n}{2} \rceil}}\)。然后就有：
而\(B^2 \equiv A \pmod{x^n}\)，所以有：
需要多項式求逆。然后直接套。當\(n=1\)時常數項為模意義下的\(\sqrt{A_0}\)。模板題保證\(A_0=1\)，省去了求解二次剩余的麻煩。復雜度\(\text{O}(n \log n)\)。
// 多項式開根
poly poly_sqrt(poly A) {
poly B = poly(1); int inv2 = P+1>>1; // 邊界和2的逆元
for (int len = 1; len < A.len; len <<= 1) {
static int x[maxn << 2], y[maxn << 2], z[maxn << 2];
for (int i = 0; i < len << 2; i++) x[i] = y[i] = z[i] = 0;
A.cpyto(x, 0, len << 1), B.cpyto(y, 0, len); // 拷貝信息
B.resize(len << 1); // 重定B多項式的大小。這里為了使得求解出的逆元長度為2*len
poly_inv(B).cpyto(z, 0, len << 1); // 拷貝逆的信息
ntt(x, len << 2, 1), ntt(y, len << 2, 1), ntt(z, len << 2, 1);
for (int i = 0; i < len << 2; i++) x[i] = 1ll*(x[i]+1ll*y[i]*y[i]%P)*z[i]%P*inv2%P; // 求解(A+B'^2)/(2B')
ntt(x, len << 2, -1);
B.load(x, len << 1); // 更新
}
B.resize(A.len);
return B;
}
多項式除法
給一個\(n\)次多項式\(F(x)\)和\(m\)次的多項式\(G(x)\)，求\(Q(x)\)和\(R(x)\)，使得滿足\(Q(x)\)的次數為\(n-m\)，\(R(x)\)的次數小於\(m\)，且\(F(x)=Q(x)G(x)+R(x)\)。
首先定義\(n\)次多項式\(A(x)\)的翻轉多項式\(A_{rev}(x)=x^nA(\frac{1}{x})\)，直觀上相當於將系數調換。然后我們來轉化：
根據定義替換得到上式：
於是\(Q\)通過求逆就能求出來了，然后\(R\)再帶回去減一下就出來了。復雜度\(\text{O}(n \log n)\)。
int getsize(int n) { // 獲取剛好大於n的2的冪
int N = 1;
while (N < n) N <<= 1;
return N;
}
poly operator - (poly A, poly B) { // 定義減法
int len = max(A.len, B.len);
A.resize(len), B.resize(len);
for (int i = 0; i < len; i++) A[i] = inc(A[i], P-B[i]);
A.resize(); return A;
}
poly operator * (poly A, poly B) { // 定義乘法
int len = getsize(A.len + B.len); poly C;
static int x[maxn << 2], y[maxn << 2];
for (int i = 0; i < len; i++) x[i] = y[i] = 0;
A.cpyto(x, 0, A.len), B.cpyto(y, 0, B.len);
ntt(x, len, 1), ntt(y, len, 1);
for (int i = 0; i < len; i++) x[i] = 1ll*x[i]*y[i]%P;
ntt(x, len, -1);
C.load(x, len), C.resize(); return C;
}
poly operator / (poly A, poly B) { // 定義除法
A.rev(), B.rev(); // .rev()用來翻轉多項式
int n = A.len - B.len + 1; // 縮短多項式
A.resize(n), B.resize(n);
poly C = A * poly_inv(B); C.resize(n); C.rev(); // 求完之后要翻轉回來
return C;
}
...
{
Q = F / G; // 求Q
R = F - Q * G; // 求R
}
多項式求ln
給一個多項式\(A\)，求\(B\)滿足\(B\equiv \ln A \pmod{x^n}\)。
前置知識：微積分。
然后多項式求導中\(C\)一般為\(0\)。復雜度\(\text{O}(n \log n)\)。
poly poly_deri(poly A) { // 多項式求導
for (int i = 0; i < A.len - 1; i++) A[i] = 1ll*A[i+1]*(i+1)%P;
A.resize(A.len - 1); return A;
}
poly poly_int(poly A) { // 多項式積分
for (int i = A.len-1; i; i--) A[i] = 1ll*A[i-1]*inv[i]%P; // 為了減小常數，預處理逆元
A[0] = 0; return A;
}
poly poly_ln(poly A) { // 多項式求ln
poly B = poly_inv(A) * poly_deri(A); B.resize(A.len);
return poly_int(B);
}
多項式求exp
給一個多項式\(A\)，求\(B\)滿足\(B\equiv e^A \pmod{x^n}\)。
這個有點復雜，需要用到泰勒公式和牛頓迭代。泰勒公式前文有講，牛頓迭代運用到了泰勒展開的東西，這里只提多項式的牛頓迭代。
牛頓迭代用來求解滿足\(G(F(x)) \equiv 0 \pmod{x^n}\)的多項式\(F(x)\)。怎么迭代？仍然是倍增。
當\(n=1\)時，\(G(F(x))\equiv 0\pmod x\)，這個單獨解。
當\(n>1\)時，假設我們已經求解出\(F_0(x)\)使得\(G(F_0(x)) \equiv 0 \pmod{x^{\lceil \frac{n}{2} \rceil}}\)，根據泰勒展開我們有：
因為\(F(x) \equiv F_0(x) \pmod{x^{\lceil \frac{n}{2} \rceil}}\)，所以\((F(x)-F_0(x))^2 \equiv 0 \pmod{x^n}\)，從二階導之后的全部略去，就有了：
又因為\(G(F(x)) \equiv 0 \pmod{x^n}\)，就有了：
於是我們就可以求解任意多項式函數的零點了。
首先變下形：
令\(G(B(x))=\ln B(x)-A(x)\)，就可以用牛頓迭代求出\(B(x)\)了。中間我們要求\(G'(F_0(x))\)。注意這里自變量是多項式\(F_0(x)\)。
所以\(G'(F_0(x)) \equiv \ln F_0(x) - A(x) \equiv \frac{1}{F_0(x)} \pmod{x^n}\)，由於自變量，這里\(A(x)\)求導后為\(0\)。
這樣就可以解決了。復雜度\(\text{O}(n \log n)\)。
// 這里保證了A[0]=0，可以得到B[0]=1
poly poly_exp(poly A) {
poly B = poly(1), C;
for (int len = 1; len < A.len; len <<= 1) { // 倍增求解
B.resize(len << 1);
C = poly(1) - poly_ln(B) + A; C.resize(len << 1); // 分布算防止長度溢出
B = B * C; B.resize(len << 1);
}
B.resize(A.len); return B;
}
如果\(A_0 \neq 0\)呢？做不了。
多項式快速冪
給一個多項式\(A\)和指數\(k\)，求\(B \equiv A^k \pmod{x^n}\)。
直接快速冪？\(\text{O}(n \log n \log k)\)！這個不夠優秀，我們有\(\text{O}(n\log n)\)的做法！
往下看：
這樣一來我們只需要：多項式求\(\ln\)，多項式求\(\exp\)即可。
// 沒什么東西，就是各種板子來一遍
poly poly_pow(poly A, int k) {
return poly_exp(k * poly_ln(A));
}
多項式多點求值
就是給你一個\(n\)次多項式\(A(x)\)和\(m\)個數，每個數\(x_i\)，求出\(A(x_1)\)、\(A(x_2)\dotsb A(x_m)\)處的取值。
直接秦九韶做是\(\text{O}(nm)\)的，不夠優秀（但是常數小），我們用構造和分治來解決。
構造多項式\(F(x)=\prod_{i=1}^{\lfloor \frac{m}{2} \rfloor} (x-x_i)\)和\(G(x)=\prod_{i=\lfloor \frac{m}{2} + 1\rfloor}^m(x-x_i)\)，這樣對於\(\forall i\in [1,\lfloor \frac{m}{2} \rfloor],F(x_i)=0\)、\(\forall i\in [\lfloor \frac{m}{2} \rfloor + 1, m],G(x_i)=0\)，然后求出\(A(x) \bmod F(x)\)的值，也就是多項式除法中剩余的\(R(x)\)。這樣有什么用呢？代入我們要求的值，發現\(R(x_i)=A(x_i)\)，因為\(F(x)\)與除法商的乘積為\(0\)。\(A(x)\)是\(n\)次多項式，而\(R(x)\)是次數小於\(\lfloor \frac{m}{2} \rfloor\)多項式，然后將左邊一半遞歸求解，以\(R(x)\)作為新的\(A(x)\)，這樣下去多項式次數自然會越來越小。
但仍然有一個問題：如何快速求出\(F(x)\)？\(F(x)=F_{left}(x)F_{right}(x)\)。提前分治求解然后存起來即可。注意這里要用\(\text{vector}\)存，否則會爆。
時間復雜度：\(\text{O}(n \log^2 n)\)。常數很大？小范圍用一用秦九韶有益身心健康。
// 這里poly內用vector代替數組
poly Q[maxn << 4]; // 分治預處理
void divide(int o, int l, int r, int *x) {
if (l == r) { Q[o].resize(2), Q[o][0] = P-x[l], Q[o][1] = 1; return; } // 邊界
int mid = l+r>>1;
divide(o<<1, l, mid, x); divide(o<<1|1, mid+1, r, x);
Q[o] = Q[o<<1] * Q[o<<1|1];
}
void calc(poly A, int o, int l, int r, int *x) {
if (r-l <= 512) { // 一定范圍內用秦九韶展開，我的代碼最適宜在512
for (int i = l; i <= r; i++) {
int ans = 0; // 運用了循環展開進一步減小常數
static ull pw[16];
pw[0] = x[i];
for (int j = 1; j < 16; j++) pw[j] = 1ll*pw[j-1]*x[i]%P;
#define I(a,b) A[a]*pw[b]
for (int j = A.len-1; j >= 16; j -= 16)
ans = (ans*pw[15]+I(j,14)+I(j-1,13)+I(j-2,12)+I(j-3,11)+I(j-4,10)+I(j-5,9)+I(j-6,8)+I(j-7,7)+I(j-8,6)+I(j-9,5)+I(j-10,4)+I(j-11,3)+I(j-12,2)+I(j-13,1)+I(j-14,0)+A[j-15]) % P;
#undef I
for (int j = (A.len-1)&15; ~j; j--) ans = (1ll*ans*x[i]+A[j])%P;
x[i] = ans;
}
return;
}
int mid = l+r>>1;
calc(A % Q[o<<1], o<<1, l, mid, x); calc(A % Q[o<<1|1], o<<1|1, mid+1, r, x); // 否則分治
}
void poly_calc(poly A, int *x, int n) {
divide(1, 0, n-1, x); calc(A, 1, 0, n-1, x); // 主過程
}
多項式快速插值
前文我們討論過拉格朗日插值。這里的問題變成了求出\(n\)個點插值出來的多項式。
看一看拉格朗日插值：
記\(g(x)=\prod_i (x-x_i)\)，那么式子變成：
這個\(g(x_i)/(x_i-x_j)\)分母分子都是\(0\)，怎么求呢？
引理：洛必達法則。
證明直接用柯西中值定理。感性理解就是在點無限靠近\(0\)時，\(f(x) \approx f'(x)\)，\(g(x) \approx g'(x)\)，但此時兩者都不為\(0\)。
於是\(g(x_i)/(x_i-x_j)\)可以直接用\(g'(x_i)/(x_i-x_j)'\)代替，即為\(g'(x_i)\)。
式子又變成了：
令\(f_{l,r}(x)\)為第\(l\)個點到第\(r\)個點的插值結果（注意這里並不是真正的插值結果），有：
然后分治求解即可。對於\(\prod\)的求解方法與上一篇相同。復雜度：\(\text{O}(n \log^2 n)\)。
void poly_calc(poly A, int *x, int n) {
divide(1, 0, n-1, x); calc(A, 1, 0, n-1, x);
}
poly inter(int o, int l, int r, int *x, int *y) {
if (l == r) return poly(1ll*y[l]*qpow(x[l], P-2)%P);
int mid = l+r>>1;
return inter(o<<1, l, mid, x, y)*Q[o<<1|1] + inter(o<<1|1, mid+1, r, x, y)*Q[o<<1]; // 如上
}
poly poly_inter(int *x, int *y, int n) {
return divide(1, 0, n-1, x), calc(poly_deri(Q[1]), 1, 0, n-1, x), inter(1, 0, n-1, x, y); // 先求出若干連續(x-xi)的乘積，然后根據洛必達法則求值，最后遞歸算出多項式
} // 傳入點，傳出插值后的多項式
組合數學
主要用來數數求方案數，然后就是一堆球的有(luan)趣(lun)問題。
組合：\(n\)個人選出\(m\)個人的方案數，記作\(\begin{aligned}n \choose m\end{aligned}\)。\(\begin{aligned}{n \choose m}=\frac{n!}{m!(n-m)!}\end{aligned}\)。
排列：\(n\)個人選出\(m\)個人有順序排隊的不同方案數，記作\(A_n^m\)。顯然\(\begin{aligned}A_n^m={n \choose m}m!\end{aligned}\)
然后有大名鼎鼎的楊輝三角。
第\(i\)行第\(j\)個的值：\(\begin{aligned}i-1 \choose j-1\end{aligned}\)。所以我們還有\(\begin{aligned}{n \choose m}={n-1 \choose m-1}+{n-1 \choose m}\end{aligned}\)。組合意義上就是新加入一個人然后看這個人選或不選來分類討論。
然后如果\(\forall x_i:0<x_1<x_2<x_3<\dotsb<x_m\leq n\)，則不等式的整數解個數有\(\begin{aligned}n \choose m\end{aligned}\)種。
插板法：將長度為\(n\)的區間划分成\(m\)端，每段不空。問方案數。很顯然要在\(n-1\)個端中插入\(m-1\)個板，方案數為\(\begin{aligned}n-1 \choose m-1\end{aligned}\)。如果段的長度可以為\(0\)呢？先往每段塞一個數，就變成上面的問題了，之后取回來即可，方案數為\(\begin{aligned}n+m-1 \choose m-1\end{aligned}\)。
有時數學模型是\(\forall x_i:0\leq x_1\leq x_2\leq x_3\leq \dotsb\leq x_m\leq n\)。沒關系，把所有的\(x_i\)加上\(i\)就變成了\(0<x_1<x_2<x_3<\dotsb<x_m\leq n+m\)，然后組合數求解即可。
組合數性質公式
1、\(\begin{aligned}{n \choose m}={n \choose n-m}\end{aligned}\)。很直觀看楊輝三角是對稱的。也可以從補集上來理解。
2、\(\begin{aligned}{n \choose m}={n-1 \choose m-1}+{n-1 \choose m}\end{aligned}\)。
3、（二項式定理）\(\begin{aligned}(a+b)^n=\sum_{i=0}^n{n \choose i}a^ib^{n-i}\end{aligned}\)。根據多項式乘多項式，\(a^ib^{n-i}\)的系數可以看成從\(n\)個數中挑出來\(i\)個\(a\)相乘，方案數即為\(\begin{aligned}n \choose i\end{aligned}\)。還有廣義二項式定理：\(\begin{aligned}(a+b)^n=\sum_{i=0}^{\infty}{n \choose i}a^ib^{n-i}\end{aligned}\)，其中\(\begin{aligned}{n \choose i}\end{aligned}\)的\(n\)可以為實數，此時\(\begin{aligned}{n \choose m}=\frac{n(n-1)(n-2)\dotsb (n-m+1)}{i!}\end{aligned}\)。
4、（頂變底不變）\(\begin{aligned}\sum_{i=0}^m{n+i \choose r}={n+m+1 \choose r+1}\end{aligned}\)。可以用性質\(2\)歸納證明，也可以從組合意義上下手：假設最右邊選的人是第\(n+i+1\)位，還剩\(r\)個人沒選，就在前面的\(n+i\)個中選。顯然最右邊位置范圍在\(n+1\)到\(n+m+1\)，所以式子顯然成立。
5、（底變頂不變）\(\begin{aligned}\sum_{i=0}^n{n \choose i}\end{aligned}=2^n\)。用二項式定理看一看\((1+1)^n\)的展開吧。
6、\(\begin{aligned}\sum_{i=0}^n(-1)^n{n \choose i}=0\end{aligned}\)。看一看\((1-1)^n\)？
7、\(\begin{aligned}\sum_{i=0}^ni\times {n \choose i}=n2^{n-1}\end{aligned}\)。發現是一個組合數構成的三角形，首項\(+\)末項搞一搞就行了。
8、（范德蒙恆等式）\(\begin{aligned}\sum_{i=0}^k{n \choose i}{m \choose k-i}={n+m \choose k}\end{aligned}\)。直接從組合意義上看，相當於枚舉左邊選\(i\)個，右邊選\(k-i\)個。最后所有可能加起來等價於\(n+m\)中選\(k\)個。還有推論：\(\begin{aligned}\sum_{i=0}^n{n \choose i}^2={2n \choose n}\end{aligned}\)。用一用性質\(1\)就能發現玄機了。如果想要嚴謹的證明，通過觀察\((1+x)^{m+n}\)的\(x^k\)的系數和\((1+x)^m\)和\((1+x)^n\)乘起來\(x^k\)的系數也能得到。
9、\(\begin{aligned}{n \choose m}{m \choose r}={n \choose r}{n-r \choose m-r}\end{aligned}\)。直接從組合意義入手，也就是相當於先划一部分人，再從這部分人中選出來\(r\)個人，反過來考慮如果選中的就是這\(r\)個人，有多少種可能先前划人的可能？實際上就是\(\begin{aligned}{n-r \choose m-r}\end{aligned}\)種。最后選\(r\)個人再乘上方案數\(\begin{aligned}{n \choose r}\end{aligned}\)即可。
Lucas定理
定理：當\(p\)是質數時，\(\begin{aligned}{n \choose m} \equiv {\lfloor \dfrac{n}{p} \rfloor \choose \lfloor \dfrac{m}{p} \rfloor}{n \bmod p \choose m \bmod p} \pmod p\end{aligned}\)
引理：\((1+x)^p \equiv 1+x^p \pmod p\)
\(Proof\)：
下面開始證明。
\(Proof\)：
int fac[maxp], ifac[maxp], inv[maxp]; // 階乘、階乘逆元、逆元
void prework(int P) {
fac[0] = ifac[0] = inv[1] = 1; // 邊界
for (int i = 2; i < P; i++) inv[i] = 1ll*(P-P/i)*inv[P%i]%P; // 求逆元
for (int i = 1; i < P; i++) fac[i] = 1ll*fac[i-1]*i%P, ifac[i] = 1ll*ifac[i-1]*inv[i]%P; // 求階乘及其逆元
}
int C(int n, int m, int p) {
return 1ll*fac[n]*ifac[m]%p*ifac[n-m]%p; // 組合數公式
}
int lucas(int n, int m, int p) {
return !m ? 1 : 1ll*lucas(n/p, m/p, p)*C(n%p, m%p, p)%p; // lucas定理，注意邊界m=0
}
復雜度\(\text{O}(p+\log_p n)\)。主要的復雜度在預處理階乘及其逆元。
如果\(p\)不是質數呢？上面的方法變得不太可行了。怎么辦呢？
發現問題出在求組合數上，逆元不一定存在！原因很簡單，因為不互質。
設\(p=p_1^{\alpha_1}p_2^{\alpha_2}\dotsb\)，我們只需求
然后再按照中國剩余定理合並即可。
對於單個的情況
我們把與\(p_i^{\alpha_i}\)互質的因子提出來（即\(p_i\)），得到：
考慮如何求\(\dfrac{n!}{P^x}\bmod P^k\)，接下來的很復雜。
發現\(P^{\lfloor \frac{n}{P} \rfloor}\)是其中的因子，\(\lfloor \dfrac{n}{P} \rfloor!\)中可能還含有因子\(P\)，所以我們要遞歸求解因子個數，這樣可以得到\(x,y,z\)。
對於剩下部分的值，我們也要求出來。發現\(\begin{aligned}\prod\limits_{i=1,P\nmid i}^n i\end{aligned}\)有循環，可以寫成：
之后就可以求出來了。細節看代碼。
void prework(int P, int Pk) {
pdt[0] = 1; // 預處理與Pk互質的數前綴積
for (int i = 1; i <= Pk; i++)
pdt[i] = i%P ? 1ll*pdt[i-1]*i%Pk : pdt[i-1];
}
int fac(ll n, int &x, int P, int Pk) {
return !n ? 1 : (x += n/P, 1ll*fac(n/P, x, P, Pk)*qpow(pdt[Pk], n/Pk, Pk)%Pk*pdt[n%Pk]%Pk); // 處理出n!/P^x，順道求出來x。公式詳細見上文
}
int a[20], b[20];
// p[]是提前篩出來的素數
int exlucas(ll n, ll m, int P) {
int sz = 0;
for (int i = 1; P > 1; i++) {
if (P % p[i]) continue; // p[i]不是P的因子
b[++sz] = 1;
int k = 0;
while (!(P % p[i])) P /= p[i], b[sz] *= p[i]; // 分解因數
prework(p[i], b[sz]); // 預處理
int x = 0, y = 0, z = 0;
a[sz] = 1ll*fac(n, x, p[i], b[sz])*inv(fac(m, y, p[i], b[sz]), b[sz])%b[sz]*inv(fac(n-m, z, p[i], b[sz]), b[sz])%b[sz]*qpow(p[i], x-y-z, b[sz])%b[sz]; // 表達式比較長，但其實就是公式的翻譯
}
return excrt(a, b, sz); // crt合並出來
}
復雜度相差不大，主要在預處理和分解因數上。
此被稱作擴展\(\bold{Lucas}\)定理。
生成函數
(learned from _rqy)
生成函數又稱母函數，是計數的一大有力工具。生成函數分為普通生成函數（\(\text{OFG}\)）和指數生成函數（\(\text{EFG}\)）。
下面先定義下降冪和上升冪。
定義下降冪：\(x^{\underline n}=x\times (x-1) \times\dotsb\times (x-n+1)=\prod\limits_{i=0}^{n-1} (x-i)\)
定義上升冪：\(x^{\overline n}=x\times (x+1) \times\dotsb\times (x+n-1)=\prod\limits_{i=0}^{n-1} (x+i)\)
前置芝士：多項式。
普通生成函數
定義普通生成函數
生成函數看着好像很簡單，也不知道有什么威力。先從最簡單的看起。
數列\(\{1,1,1,1,\dotsb\}\)的生成函數為\(F(x)=\sum\limits_{i=0}^{\infty}x^i\)，我們發現它是個等比數列，即\(F(x)=\dfrac{1}{1-x}\)，這里我們假設\(\mid x \mid<1\)。
就像這個生成函數一樣，我們可以將一些生成函數化成一個簡單的式子，以有限表示無限。有趣吧！再來幾個！
數列\(\{0,0,1,1,1,\dotsb\}\)的生成函數？相當於最普通的向右平移了\(x^2\)，即\(\dfrac{x^2}{1-x}\)。
數列\(\{1,-1,1,-1,\dotsb\}\)？\(\dfrac{1}{1+x}\)。
數列\(\{1,c,c^2,c^3,\dotsb\}\)？\(\dfrac{1}{1-cx}\)。
數列\(\begin{aligned}\{{n \choose 0},{n \choose 1},{n \choose 2},{n \choose 3},\dotsb\}\end{aligned}\)？\(=(1+x)^n\)。
數列\(\{0,1,2,3,\dotsb\}\)？這個有點難，其實是\(\dfrac{1}{(1-x)^2}\)，相當於\(F(x)\)與自身卷一次積得到的結果，也可以從微分的角度理解；如果暫時沒看明白，下面會將運算法則。
數列\(\{0,1,\dfrac{1}{2},\dfrac{1}{3},\dotsb\}\)？這個更難了，為\(\ln\dfrac{1}{1-x}\)。兩邊積分即可得到。
同理數列\(\{0,1,-\dfrac{1}{2},\dfrac{1}{3},-\dfrac{1}{4},\dotsb\}\)對應的生成函數為\(\ln(1+x)\)，是由\(F(x)=-\dfrac{1}{1+x}\)兩邊積分得來的。
還有很多就不說了，下面談談生成函數的運算法則，這樣我們可以用來推導更多的生成函數。
發現最后一個等式：左邊乘積對應右邊卷積。
指數生成函數
定義指數生成函數
我們發現\(\begin{aligned}{n \choose m}=\frac{n!}{m!(n-m)!}\end{aligned}\)，如果指數生成函數的第\(m\)項和第\(n-m\)項的值相乘，會得到\(\dfrac{f_mf_{n-m}}{m!(n-m)!}\)，如果再乘上\(n!\)就可以得到\(\begin{aligned}{n \choose m}f_mf_{n-m}\end{aligned}\)。
所以可以得到
這個適用於排列的計算。
同樣數列對應指數生成函數。
數列\(\{1,1,1,1,\dotsb\}\)對應的生成函數為\(e^x\)。看一看泰勒展開？
數列\(\{1,c,c^2,c^3,\dotsb\}\)的生成函數為\(e^{cx}\)。
數列\(\{0,1,2,3,\dotsb\}\)的生成函數\(xe^x\)。求個導再平移一下就出來了。
指數生成函數的運算法則
求導積分注意有的因數被約掉了。再加上上面的乘法，基本上構成了指數生成函數的運算法則，這些不難推導。
生成函數的運用
求通項公式
利用遞歸式的特性，構造生成函數，然后通過解方程的形式得到式子后再展開得到信息。
Fibonacci數列
數列\(\{0,1,1,2,3,5,8,\dotsb\}\)被稱為\(\text{Fibonacci}\)數列，它有以下遞推式
定義\(\text{Fibonacci}\)生成函數為\(\sum\limits_{i\geq 0}f_ix^i\)，則根據遞推式，我們有：
加上一個單位函數是為了保證邊界\(i=1\)時為\(1\)。這樣我們可以解出來\(F(x)=\dfrac{x}{1-x-x^2}\)，進一步展開得到：
其中\(\phi=\dfrac{1+\sqrt 5}{2}\)，\(\hat \phi=\dfrac{1-\sqrt 5}{2}\)。
然后因為\(\dfrac{1}{1-cx}=1+cx+c^2x^2+\dotsb\)，得到
得知第\(n\)項的系數為\(\dfrac{\phi^i-\hat\phi^i}{\sqrt 5}\)，即為\(\text{Fibonacci}\)的第\(n\)項的值。
Catalan數
\(\text{Catalan}\)數的遞推式：
定義卡特蘭數的生成函數：\(F(x)=\sum\limits_{i \geq 0}f_ix^i\)，由遞推式得到：
解得\(F(x)=\dfrac{1\pm \sqrt{1-4x}}{2x}\)，有一個\(\pm\)怎么辦？
（坑坑坑）先看用母函數推卡特蘭數通項公式 知乎吧（博主瘋了）
斯特林數
組合數學中有兩類斯特林數：第一類斯特林數、第二類斯特林數。
第一類斯特林數
第一類斯特林數表示將\(n\)個不同元素構成\(m\)個圓排列的數目，記作\(s(n,m)\)。能推出遞推式
特別的，\(s(0,0)=1\)，\(s(n,0)=0\)，\(s(n,n)=1\)。
斯特林數一般求一行或一列。
求一行
有一個式子：
\(Proof\)：
如何求解呢？可以用分治\(\text{FFT}\)，\(\text{O}(n\log^2 n)\)。但不夠優秀，我們還有\(\text{O}(n\log n)\)的方法。
考慮倍增。假設我們求出了\(x^{\overline n}\)對應的多項式\(\sum\limits_{i=0}^n a_ix^i\)，我們來求一下\((x+n)^{\overline n}\)。
后面式子翻轉一下卷積即可求出。復雜度為\(\text{T}(n)=\text{T}(n/2)+\text{O}(n\log n)=\text{O}(n\log n)\)。
遞歸求解。對於余項直接乘。
(咕咕咕)
求一列
定義環排列指數生成函數
因為\(n\)個數的環排列有\((n-1)!\)種。用指數生成函數滿足划分環大小時，在同一個環內的方案不被重復計算。
所以有：
最后是由於划分\(m\)個環大小的過程中有\(m!\)種相同的情況，所以要除去。
多項式快速冪即可。注意要偏移下標，因為常數項為\(0\)。
(咕咕咕)
第二類斯特林數
第二類斯特林數表示將\(n\)個不同的元素放到\(m\)個相同集合里的方案數，記作\(S(n,m)\)。也有遞推式
二項式反演
如果：
那么有：
此被稱作二項式反演。
\(Proof\)：
最后一步是根據\(\begin{aligned}\sum\limits_{i=0}^n(-1)^i{n \choose i}=\begin{cases}0&(n>0)\\1&(n=0)\end{cases}\end{aligned}\)得來的。
還有一個比較常用的式子：
由：
得到：
證明類似。
拉格朗日反演
定義兩個多項式\(f(x)\)和\(g(x)\)，如果\(f(g(x))=x\)成立，此時顯然\(f(g(x))=g(f(x))\)，且\([x^0]f(x)=0\)，\([x^1]f(x)\times [x^1]g(x)=1\)，那么有
定理\(1\)：
\(Proof\)
單位根反演
定理：\(\dfrac{\sum\limits_{k=1}^n \omega_n^{dk}}{n}=[n\mid d]\)。\([state]\)是單位函數。
不證，在\(\text{FFT}\)的求和引理中有。