進軍多項式。
1. 拉格朗日插值
1.1. 普通插值
首先給出公式：
解釋：對於每對點值 \((x_k,y_k)\)，我們需要構造出一個函數 \(G(x)\)，使得其在 \(x=x_k\) 處的取值為 \(y_k\)，其余處取值為 \(0\)。
首先構造函數 \(D(x)=\prod_{i=1,i\neq k}^n x-x_i\)。顯然當 \(i\neq k\) 時，有 \(D(x_i)=0\)。但是現在我們不能保證 \(D(x_k)=y_k\)。為了使 \(D(x_k)=y_k\)，我們只需要先將其除以 \(D(x_k)\)，再乘以 \(y_k\) 即可，這就有了上面的拉格朗日插值公式。
通常情況下，題目會要求我們求出 \(F(x)\) 在給定某個 \(x\) 處的取值，此時我們不把 \(x\) 看做函數的一個元，而是直接將 \(x\) 帶入上式即可。時間復雜度為 \(\mathcal{O}(n^2)\)，代碼見例題 I。
1.2. 連續取值插值
很多情況下，我們求出的點值 \(x_i\) 滿足 \(x_i=i\)，即 \(x_i\) 是連續的。此時我們重新寫一下公式：
記 \(p_i=\prod_{j=1}^ix-i\)，\(s_i=\prod_{j=i+1}^n x-i\)，這些可以線性預處理，那么上述柿子右邊就變成了
預處理階乘就可以線性插值，應用見例題 II。
1.3. 求 \(F(x)\) 各項系數
構造函數 \(D(x)=\prod_{i=1}^n(x-x_i)\)，設 \(d_k=y_i\prod_{i=1,i\neq k}^n\dfrac{1}{x_k-x_i}\)，則 \(F(x)=\sum_{k=1}^n \dfrac{d_kD(x)}{x-x_k}\)。注意到 \(D(x)\) 的各項系數可以在 \(n^2\) 的時間內暴力處理出來，而對於每個 \(k\)，我們可以線性將 \(D(x)\) 除以一個一次多項式。最后加和即可。時間復雜度 \(\mathcal{O}(n^2)\)。
1.4. 例題
I. P4781 【模板】拉格朗日插值
板子題。
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const ll mod=998244353;
const int N=2e3+5;
ll ksm(ll a,ll b){
ll s=1;
while(b){
if(b&1)s=s*a%mod;
a=a*a%mod,b>>=1;
} return s;
}
int n,k,x[N],y[N],ans;
int main(){
cin>>n>>k;
for(int i=1;i<=n;i++)cin>>x[i]>>y[i];
for(ll i=1,s1=1,s2=1;i<=n;i++,s1=s2=1){
for(int j=1;j<=n;j++)if(i!=j)s1=s1*(k-x[j])%mod,s2=s2*(x[i]-x[j])%mod;
ans=(ans+y[i]*s1%mod*ksm(s2,mod-2))%mod;
} cout<<(ans%mod+mod)%mod<<endl;
return 0;
}
II. CF622F The Sum of the k-th Powers
經典題。一個結論是自然數 \(k\) 次方和是 \(k+1\) 次多項式，那么只需要帶 \((i,i^k)\ (i\in [0,k+1])\) 插值即可。注意到當 \(i=0\) 時對答案無貢獻，所以在插值的時候可以跳過（但是預處理 \(p,s\) 的時候仍應考慮 \(i=0\) 時 \(n-i\) 的這個 \(n\)）。
時間復雜度 \(\mathcal{O}(k\log k)\)，使用線性篩篩 \(i^k\) 可以做到線性。
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const ll mod=1e9+7;
const int N=1e6+5;
ll ksm(ll a,ll b){
ll s=1;
while(b){
if(b&1)s=s*a%mod;
a=a*a%mod,b>>=1;
} return s;
} ll inv(ll x){return ksm(x,mod-2);}
ll n,k,ans;
ll p[N],s[N],fc[N];
int main(){
cin>>n>>k,s[k+2]=1;
for(int i=0;i<=k+1;i++)fc[i]=i?fc[i-1]*i%mod:1,p[i]=i?p[i-1]*(n-i)%mod:n;
for(int i=k+1;~i;i--)s[i]=i==k+1?n-i:s[i+1]*(n-i)%mod;
for(ll i=1,res=0;i<=k+1;i++){
res=(res+ksm(i,k))%mod;
ans=(ans+p[i-1]*s[i+1]%mod*res%mod*ksm(fc[i]*((k-i)&1?1:-1)*fc[k+1-i]%mod,mod-2))%mod;
} cout<<(ans%mod+mod)%mod<<endl;
return 0;
}
III. P4463 [集訓隊互測 2012] calc
不妨設 \(a_i<a_{i+1}\ (1\leq i<n)\)，那么只需將答案乘 \(n!\) 即可。
首先把轉移方程寫出來：\(f_{i,j}\) 表示值域落在 \([1,i]\) 且長度為 \(j\) 時的答案（即 \(k=i\)，\(n=j\)）。顯然有 \(f_{i,j}=f_{i-1,j}+[j>0]f_{i-1,j-1}\times i\)。
-
接下來我們分析一下 \(f_{i,j}\) 的次數：
首先我們有 \(f_{i,0}\) 是關於 \(i\) 的 \(0\) 次多項式。接下來使用一個小 Trick。
Trick：差分分析次數。
根據轉移方程 \(f_{i,j}=f_{i-1,j}+f_{i-1,j-1}\times i\)，有 \(f_{i,j}-f_{i-1,j}=f_{i-1,j-1}\times i\)。由於將兩個相差 \(1\) 的數帶入一個 \(c\) 次多項式得到的差是 \(c-1\) 次多項式，我們有 \(f_{i,j}\) 的次數 \(-1\) 等於 \(f_{i,j-1}\) 的次數 \(+1\)，因此 \(f_{i,j}\) 是關於 \(i\) 的 \(2j\) 次多項式。
所以拉格朗日插值即可。
const int N=1.5e3+5;
ll f[N][N],k,n,p,ans;
ll ksm(ll a,ll b){
ll s=1;
while(b){
if(b&1)s=s*a%p;
a=a*a%p,b>>=1;
}
return s;
}
int main(){
cin>>k>>n>>p,f[0][0]=1;
for(ll i=1;i<=n*3+4;i++)
for(ll j=0;j<=min(i,n);j++)
f[i][j]=(f[i-1][j]+(j?f[i-1][j-1]*i:0))%p;
for(ll i=n;i<=n*3+3;i++){
ll s1=f[i][n],s2=1;
for(ll j=n;j<=n*3+3;j++)if(i!=j)
s1=s1*(k+p-j)%p,
s2=s2*(i+p-j)%p;
ans=(ans+s1*ksm(s2,p-2))%p;
}
for(ll i=1;i<=n;i++)ans=ans*i%p;
cout<<ans<<endl;
return 0;
}
*IV. [BZOJ2137]submultiple
題意簡述：對於 \(m\) 的所有約數 \(d\)，求 \(\sum \mathrm{d}^k(d)\)。\(m\) 由唯一分解給出。
對於 \(45\%\) 的數據，\(k\leq 2^{31}-1\)，\(p\leq 10^5\)；對於剩下 \(55\%\) 的數據，\(k\leq 12\)，\(p\leq 2^{63}-1\)。
\(n\leq 10^5\)，\(n\) 表示 \(m\) 由前 \(n\) 小的質數組成。
一個顯然的結論是我們只關心唯一分解后每個質因子的次數 \(p_i\)，而底數是什么根本不重要：約數個數僅與質因子次數有關。考慮枚舉 \(d\) 每個質因子的次數 \(q_i\in [0,p_i]\)，得到如下柿子：
根據加法對乘法的分配律，稍作化簡：
設 \(S_k(n)\) 表示 \(\sum_{i=1}^n i^k\)，則題目就是在求 \(\prod_{i=1}^nS_k(p_i+1)\)。一個經典結論是自然數的 \(k\) 次方前綴和是 \(k+1\) 次多項式，因此對於 \(k\) 很大，\(p\) 很小的部分直接 \(\mathcal{O}(p\log k)\) 預處理所有 \(S_k(i)\)，而 \(k\) 很小，\(p\) 很大的部分直接單次 \(\mathcal{O}(k^2)\) 或 \(\mathcal{O}(k)\)（連續取值）插值即可。
綜上，時間復雜度為 \(\mathcal{O}(\min(n+p\log k,nk^2))\)。
template <class T> void cmin(T &a, T b){a = a < b ? a : b;}
template <class T> void cmax(T &a, T b){a = a > b ? a : b;}
bool Mbe;
const int N = 1e5 + 5;
const ll mod = 1e9 + 7;
ll ksm(ll a, ll b = mod - 2) {
ll s = 1;
while(b) {
if(b & 1) s = s * a % mod;
a = a * a % mod, b >>= 1;
} return s;
}
ll n, k, type = 1, pw[N];
void solve1() {
static ll f[N], ans = 1; f[0] = 0;
for(int i = 1; i < N; i++)
f[i] = (f[i - 1] + ksm(i, k)) % mod;
for(int i = 1; i <= n; i++)
ans = ans * f[pw[i] + 1] % mod;
cout << ans << endl;
}
void solve2() {
static ll f[N], ans = 1; f[0] = 0;
for(int i = 1; i <= k + 1; i++) f[i] = (f[i - 1] + ksm(i, k)) % mod;
for(int i = 1; i <= n; i++) pw[i] %= mod;
for(int i = 1; i <= n; i++) {
ll val = 0;
for(int j = 0; j <= k + 1; j++) {
ll nume = f[j], deno = 1;
for(int p = 0; p <= k + 1; p++) if(p != j)
deno = deno * (mod + j - p) % mod,
nume = nume * (pw[i] + 1 + mod - p) % mod;
val = (val + nume * ksm(deno)) % mod;
}
ans = ans * val % mod;
} cout << ans << endl;
}
bool Med;
int main(){
fprintf(stderr, "%.2lf\n", (&Mbe - &Med) / 1048576.0);
cin >> n >> k;
for(int i = 1; i <= n; i++) pw[i] = read(), type &= pw[i] <= 1e5;
if(type) solve1();
else solve2();
return 0;
}
*V. [BZOJ3453]tyvj 1858 XLkxc
題意簡述：設 \(f_k(n)\) 表示 \(\sum_{i=1}^ni^k\)，\(g_k(n)\) 表示 \(\sum_{i=1}^n f_k(i)\)。給定 \(k,a,n,d\)，求 \(\sum_{i=0}^{n}g_{k}(a+id)\bmod 1234567891\ (p)\)。
\(k\leq 123\)，\(a,n,d\leq 123456789\)。
感覺數據范圍像是出題人隨手打的（大霧。
經典結論：\(n^k\) 是關於 \(n\) 的 \(k\) 次多項式，其前綴和（積分）\(f_k(n)\) 是關於 \(n\) 的 \(k+1\) 次多項式，\(f_{k}(n)\) 的積分 \(g_{k}(n)\) 是關於 \(n\) 的 \(k+2\) 次多項式，次數隨着前綴和次數是線性增長的，所以再套幾層娃也沒關系。
考慮只算一個 \(g_k(a)\) 怎么做：直接拉格朗日插值即可，把插值公式寫出來：
如果我們把后面那個 \(\prod\) 看做一個關於 \(p\) 的 \(k+1\) 次多項式 \(h_i(p)\)，那么柿子可以寫作：
不慌，對 \(h_i\) 做一遍前綴和，記作 \(H_i\)，這是一個關於 \(n\) 的 \(k+2\) 次多項式，可以繼續插值插出來（取 \(0\sim k+2\) 處的點值）。即求：
其中 \(h_i(j)\) 是一個插值形式的柿子。這叫 插 值 套 差 值。注意到兩次插值都是連續取值插值，因此可以做到 \(\mathcal{O}(k)\) 插值，總時間復雜度 \(\mathcal{O}(k^2)\)，比 tzc 直接求系數 + 二項式展開不知道簡單到哪里去了（大霧。
代碼偷懶使用了 \(k^2\) 插值。
template <class T> void cmin(T &a, T b){a = a < b ? a : b;}
template <class T> void cmax(T &a, T b){a = a > b ? a : b;}
bool Mbe;
const int N = 1e3 + 5;
const ll mod = 1234567891;
ll ksm(ll a, ll b = mod - 2) {
ll s = 1;
while(b) {
if(b & 1) s = s * a % mod;
a = a * a % mod, b >>= 1;
} return s;
}
ll k, a, n, d, c, y[N];
void solve() {
ll ans = 0;
cin >> k >> a >> n >> d, c = k + 3;
for(int i = 1; i <= c; i++) y[i] = (y[i - 1] + ksm(i, k)) % mod;
for(int i = 1; i <= c; i++) y[i] = (y[i - 1] + y[i]) % mod;
for(int i = 1; i <= c; i++) {
ll coef = 0;
static ll y2[N];
for(int j = 0; j <= c + 2; j++) {
y2[j] = j ? y2[j - 1] : 0;
ll deno = 1, nume = 1;
for(int p = 1; p <= c; p++) if(i != p)
nume = nume * ((a + j * d) % mod + mod - p) % mod,
deno = deno * (i + mod - p) % mod;
y2[j] = (y2[j] + nume * ksm(deno)) % mod;
} // 第一遍插值插出 h_i
for(int j = 1; j <= c + 2; j++) {
ll deno = 1, nume = y2[j];
for(int p = 1; p <= c + 2; p++) if(j != p)
nume = nume * (n + mod - p) % mod,
deno = deno * (j + mod - p) % mod;
coef = (coef + nume * ksm(deno)) % mod;
} // 第二遍插值插出 H_i
ans = (ans + coef * y[i]) % mod;
} cout << ans << endl;
}
bool Med;
signed main(){
fprintf(stderr, "%.2lf\n", (&Mbe - &Med) / 1048576.0);
int T; cin >> T;
while(T--) solve();
return 0;
}
2. 多項式乘法：加法卷積
2.1. FFT
對於一個 \(n\) 次多項式 \(F(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n\)，我們可以用 \(n+1\) 個點值 \((x_k,y_k)\) 唯一確定該多項式。即 \(y_k=\sum_{i=0}^na_ix_k^i\)（注意下文的 \(n\) 表示不小於多項式次數 \(+1\) 的最小的 \(2\) 的冪）。設 \(A=F\times G\)，注意到 \(A(x)=F(x)G(x)\)，其中 \(x\) 是一個確定的值。因此，我們只需要將一個多項式快速（\(n\log n\)）轉點值，再快速轉成系數表示，就可以做到時間復雜度 \(n\log n\) 的多項式乘法。
點值的取法很有講究，高明的方法能夠極大化地減小時間復雜度。這里我們采用 \(n\) 次單位根 \(\omega_n\)，並主要利用以下性質：
- \(\omega_n^k=\omega_{2n}^{2k}\)。
- \(\omega_n^k=-\omega_n^{k+n/2}\)。
- \(\omega_n^n=1\)。
- \(\omega_n=\cos(\dfrac{2\pi}{n})+i\sin(\dfrac{2\pi}{n})\)，從而計算 \(n\) 次單位根。
當 \(n=4\) 時，\(F(x)=a_0+a_1x+a_2x^2+a_3x^3=(a_0+a_2x^2)+x(a_1+a_3x^2)\)。設 \(L(x)=a_0+a_2x\)，\(R(x)=a_1+a_3x\)，那么 \(F(x)=L(x^2)+xR(x^2)\)。將單位根 \(\omega_n\) 帶入，那么當 \(k<\dfrac{n}{2}\) 時，\(F(\omega_n^k)=L(\omega_{n/2}^k)+\omega_{n}^kR(\omega_{n/2}^k)\)，\(F(\omega_n^{k+n/2})=L(\omega_{n/2}^k)-\omega_n^kR(\omega_{n/2}^k)\)。注意到兩式只有一個正負號不同。因此，如果我們已經知道了 \(L,R\) 在 \(\omega_{n/2}^i,\ i\in[0,\dfrac{n}{2})\) 處的取值，那么我們就可以在線性時間內求出 \(F\) 在 \(\omega_n^i, \ i\in[0,n)\) 處的取值。考慮遞歸樹的每一層都是線性的，因此總復雜度為 \(n\log n\)。
遞歸處理很慢，於是我們使用迭代：考慮系數 \(a_i\) 向 \(L\) 分治為 \(0\)，向 \(R\) 分治為 \(1\)，那么顯然 \(a_i\) 最終形成的 “分治序列” 形成的二進制數倒過來就是 \(i\) 的二進制表示。考慮求出 \(r_i\) 表示 \(i\) 二進制翻轉后得到的數。假設 \(r_0,r_1,\cdots,r_{i-1}\) 都已經求出 ，那么有 \(r_i=\lfloor\dfrac{r_{\lfloor\dfrac{i}{2}\rfloor}}{2}\rfloor+\dfrac{n}{2}\times(i\bmod 2)\)。左邊是 \(i\) 不考慮最低位（即假設其為 \(0\)）時二進制翻轉得到的數，然后再考慮 \(i\) 的最低位的影響即可。然后對於每一對無序對 \((i,r_i)\)，將 \(a_i\) 與 \(a_{r_i}\) 交換，那么最終的分治樹的形態類似線段樹，直接從最底層向上迭代即可。該操作被稱為蝴蝶迭代。
卡常技巧：
- 結構體不要寫構造函數。
- 重載加減乘運算符放到
struct里面。
2.2. NTT
由於復數運算很慢，而通常情況下我們是在模意義下進行多項式運算，所以當模數取一些特殊值時，我們可以用 \(\mathbb{Z}\) 中的數 \(g\) 代替單位根的復數運算。具體地，\(g\) 需要是模 \(p\) 的原根。
鴿着。
2.3. FFT 優化字符串匹配
Trick 1：翻轉一個序列常常可以使關於它的某些計算變成卷積的形式。
對於一個文本串 \(s\) 與匹配串 \(t\)（下標從 \(0\) 開始），設它們的長度分別為 \(n\)，\(m\)。稱它們在位置 \(p\) 匹配，當且僅當對於任意 \(i\in[0,m)\)，有 \(s_{p-m+i+1}=t_i\)。不難發現它的充分條件為 \(\sum_{i=0}^{m-1}(s_{p-m+i+1}-t_i)^2=0\)。展開，得到 \(\sum_{i=0}^{m-1}(s^2_{p-m+i+1}+t^2_{i}-2s_{p-m+i+1}t_i)=0\)。注意到前面兩項容易預處理得到，但后面一項同時關於兩個字符串，比較麻煩。又發現它的形式類似卷積，但又不是卷積：翻轉字符串 \(t\) 即可。因此柿子變為 \(\sum_{i=0}^{m-1}(s^2_{p-m+i+1}+t_{m-i-1}^2-2s_{p-m+i+1}t_{m-i-1})\)，后面一項即 \(2\sum_{0\leq i<m,i+j=p}s_jt_i\)，FFT 計算即可。
對於有通配符的字符串，我們不妨將該位置上的值設為 \(0\)，然后乘到上面的柿子里去，即 \(\sum_{i=0}^{m-1}(s_{p-m+i+1}-t_i)^2s_{p-m+i+1}t_i\)。化簡得到 \(\sum_{0\leq i<m,i+j=p}(s^3_jt_i-2s^2_jt^2_i+s_jt^3_i)\)，做 6 次 DFT + 1 次 IDFT 即可。
Trick 2：一般情況下，上述方法足夠應付多數題目。但是如果萬惡的出題人卡了 FFT 精度就涼涼了。因此，為了保險，應盡量使用 NTT。但是 NTT 也有一個致命問題：如果計算出來的值剛好是 \(998244353\) 的倍數，那么就會在不匹配的地方判定匹配。看起來好像沒有解決辦法了？非也。在 P4173 殘缺的字符串 這題的討論區 https://www.luogu.com.cn/discuss/show/303076，我找到了一個高明的手段解決這個問題：注意到上述柿子最大值可達到 \(3\times m\times (|\mathbb{\Sigma}|-1)^4\)，如果 \(m\) 取 \(5\times 10^5\)，字符集取大小 \(26\)，那么數量級為 \(1.5\times 10^6\times 25^4\approx 6\times 10^{11}\)，顯然不太行。但是實際上我們沒必要將整個 \(s\) 和 \(t\) 乘進去，因為我們關心的只是某一位是否是通配符，而具體這一位是什么並不重要。因此，我們令 \(S_i=[s_i\neq\texttt{*}]\)，\(T_i=[t_i\neq\texttt{*}]\)，只需要將 \(S,T\) 而非 \(s,t\) 乘入即可。 這時最大值僅為 \(3\times 5\times 10^5\times 25^2=9.375\times 10^8<998244353\)，有驚無險地保證了正確性。當然，不同題目還應根據 \(m\) 和字符集大小的不同取值具體分析正確性。
也許看了上述分析的你認為：既然將 \(S,T\) 乘進去，那么不如直接將是通配符的位置當做 \(0\) 來計算 \(\sum_{i=0}^{m-1}(s_{p-m+i+1}-t_i)^2\) 不就好了？非也。因為這樣在計算只關於某一個字符串的項時，考慮不到另外一個字符串的對應位置是否是通配符。因此，6 次 DFT 是逃不掉的。最終的柿子即為 \(\sum_{0\leq i<m,i+j=p}((s_j^2S_j)\times T_i-2(s_jS_j)\times (t_iT_i)+S_j\times (t^2_iT_i))\)。
例題 & 代碼可以看 II.
2.4. 任意模數：MTT
MTT 又稱為任意模數 FFT。
2.4.1. 拆系數：7 次 FFT
對於值域為 \(V\) 的兩個 \(n\) 次多項式，它們相乘后值域為 \(V^2n\)。對於一般的題目，\(V\approx 10^9\)，\(n\approx 10^5\)，所以 \(V^2n\approx 10^{23}\)，long double 都無法承受。
對此，我們可以將系數拆分成 \(xp+y\) 的形式，其中 \(p\) 一般取平方大於值域的最小的 \(2\) 的冪，這里可以使用 \(2^{15}\)。這樣，我們設 \(a=a_0p+a_1\)，\(b=b_0p+b_1\)，那么 \(c=a\times b=(a_0p+a_1)(b_0p+b_1)=a_0b_0p^2+(a_0b_1+a_1b_0)p+a_1b_1\)。值域約為 \((2^{15})^2\times n\approx 10^{14}\)，可以接受。
綜上，我們有了一個顯然的 4 DFT + 3 IDFT = 7 FFT 的做法，但是不夠快。
2.4.2. 優化 Trick：合並 DFT 兩個實系數多項式
構造多項式 \(p=a+bi\)，\(q=a-bi\)，由於 \(a,b\) 的系數都是實數，故 \(p,q\) 對應項系數互為共軛。為了讓它們對應位置點值表示的結果也是共軛的，由於積的共軛等於共軛的積，我們需要保證帶入的單位根也互相共軛。例如，設 \(p\) 第 \(j\) 項前的系數為 \(c+di\)，則 \(q\) 的對應項系數為 \(c-di\)。為了讓 \((c+di)\omega_p^j\) 與 \((c-di)\omega_q^j\) 互為共軛，則 \(\omega_p\) 與 \(\omega_q\) 也需互為共軛。
如果我們先對 \(p\) 進行 DFT，也就是算出 \(p(\omega_n^i)\) 的點值，那么根據上面的分析，它與 \(q(\overline{\omega_n^i})\) 互為共軛。由於 \(\omega_{n}^i\) 當 \(i=0\) 時的共軛為它本身，當 \(i\neq 0\) 時的共軛為 \(\omega_n^{n-i}\)，並且 \(p\) 在位置 \(i\) 處求得的點值為 \(p(\omega_n^i)\)，因此為了得到 \(q\) 在位置 \(i\) 處的點值 \(q(\omega_n^i)\)，我們只需要將 \(p\) 的 \(1\sim n-1\) 項的點值進行序列翻轉，再求共軛即可。\(a=\dfrac{p+q}{2}\)，\(b=\dfrac{p-q}{2i}=\dfrac{(q-p)i}{2}\)。
同時 DFT 兩個實系數多項式，可以壓縮成一次 DFT。大大減小了常數。
2.4.3. 4 次 FFT
根據 2.4.2. 的 Trick，原來的 4 DFT 就被優化成了 2 DFT。四維變二維。
有了 \(a_0,a_1,b_0,b_1\) 的點值，接下來可以計算 \(a_0b_0,a_0b_1+a_1b_0,a_1b_1\) 的點值。可是 IDFT 回去仍要做三遍 FFT，並不是很優秀。
我們借鑒上述思想，構造多項式 \(p=a_0b_0+a_0b_1i\)。將其進行 IDFT 之后，由於 \(a_0b_0\) 與 \(a_0b_1\) 的系數表示法的每一項系數都為實數，因此 \(p\) 的實部就是 \(a_0b_0\) 的系數表示，\(p\) 的虛部就是 \(a_0b_1\) 的系數表示。一個直觀的解釋就是：給出 \(3+4i\)，你知道是它是由一個實數和一個純虛數相加得到，那么很容易就能得出原來的兩個數分別是 \(3\) 和 \(4i\)。相反，如果你知道它是由一個復數和另一個復數相加得到，那么顯然無法還原原來的兩個數是什么。而我們所做的就是前者。
對於 \(a_1b_0\) 與 \(a_1b_1\)，同理。
綜上，我們只需要進行 2 DFT + 2 IDFT = 4 FFT，常數非常優秀。代碼見例題 III.
2.5. 任意模數：三模數 NTT
考慮對 \(f(x)\) 在可 NTT 模數下進行 NTT，然后使用中國剩余定理 CRT 合並即可。這個還是比較好理解的，就不細說了。
一般選模數分別為 \(469762049,998244353,1004535809\)，這樣可以處理到值域為大約 \(4\times 10^8\times 10^9\times 10^9\approx 4\times 10^{26}\) 的多項式乘法，對付一般的題目應該是綽綽有余了。
不過常數稍大，可以手寫結構體將對於這三個模數分別取模后的值打包，這樣常數會小一些。
2.6. 任意長度：Bluestein 算法
2.7. 例題
I. P3803 【模板】多項式乘法（FFT）
FFT：
#include <bits/stdc++.h>
using namespace std;
#define ld double
const int N=1<<21;
const ld Pi=acos(-1);
struct com{
ld x,y;
com operator + (com b){return (com){x+b.x,y+b.y};}
com operator - (com b){return (com){x-b.x,y-b.y};}
com operator * (com b){return (com){x*b.x-y*b.y,x*b.y+y*b.x};}
}a[N],b[N];
int n,m,lim=1,bit,r[N];
void FFT(com *a,int tp){
for(int i=0;i<lim;i++)if(i<r[i])swap(a[i],a[r[i]]);
for(int l=1;l<lim;l<<=1){
com wn={cos(Pi/l),tp*sin(Pi/l)};
for(int j=0;j<lim;j+=(l<<1)){
com w={1,0},x,y;
for(int k=0;k<l;k++,w=w*wn)
x=a[j+k],y=w*a[j+k+l],a[j+k]=x+y,a[j+k+l]=x-y;
}
}
}
int main(){
cin>>n>>m;
for(int i=0;i<=n;i++)scanf("%lf",&a[i].x);
for(int i=0;i<=m;i++)scanf("%lf",&b[i].x);
while(lim<=n+m)lim<<=1,bit++;
for(int i=1;i<lim;i++)r[i]=(r[i>>1]>>1)|((i&1)<<bit-1);
FFT(a,1),FFT(b,1);
for(int i=0;i<lim;i++)a[i]=a[i]*b[i];
FFT(a,-1);
for(int i=0;i<=n+m;i++)cout<<(int)(a[i].x/lim+0.5)<<" ";
return 0;
}
NTT：
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned long long
#define gc getchar()
inline int read(){
int x=0; char s=gc;
while(!isdigit(s))s=gc;
while(isdigit(s))x=s-'0',s=gc;
return x;
}
const ll mod=998244353;
const int N=1<<21;
ll ksm(ll a,ll b){
ll s=1;
while(b){
if(b&1)s=s*a%mod;
a=a*a%mod,b>>=1;
} return s;
} ll inv(ll x){return ksm(x,mod-2);}
const int G=3;
const int ivG=inv(3);
ll tr[N],lim=1,l;
ll n,m,f[N],g[N];
void NTT(ll *a,bool tp){
static ull f[N],w[N]; w[0]=1;
for(int i=0;i<lim;i++)f[i]=a[tr[i]];
for(int l=1;l<lim;l<<=1){
ll wn=ksm(tp?G:ivG,(mod-1)/(l+l));
for(int i=1;i<l;i++)w[i]=w[i-1]*wn%mod;
for(int i=0;i<lim;i+=l<<1){
for(int j=0;j<l;j++){
int y=w[j]*f[i|j|l]%mod;
f[i|j|l]=f[i|j]+mod-y,f[i|j]+=y;
}
} if(l==(1<<17))for(int i=0;i<lim;i++)f[i]%=mod;
}
if(!tp){
ll iv=inv(lim);
for(int i=0;i<lim;i++)a[i]=f[i]%mod*iv%mod;
} else for(int i=0;i<lim;i++)a[i]=f[i]%mod;
}
int main(){
cin>>n>>m;
for(int i=0;i<=n;i++)f[i]=read();
for(int i=0;i<=m;i++)g[i]=read();
while(lim<=n+m)lim<<=1,l++;
for(int i=1;i<lim;i++)tr[i]=(tr[i>>1]>>1)|((i&1)<<l-1);
NTT(f,1),NTT(g,1);
for(int i=0;i<lim;i++)f[i]=f[i]*g[i]%mod;
NTT(f,0);
for(int i=0;i<=n+m;i++)printf("%lld ",f[i]);
}
II. P4173 殘缺的字符串
經典字符串匹配題。
#include <bits/stdc++.h>
using namespace std;
typedef double db;
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
#define gc getchar()
#define pb push_back
#define mem(x,v,n) memset(x,v,sizeof(int)*n)
#define cpy(x,y,n) memcpy(x,y,sizeof(int)*n)
const ld Pi=acos(-1);
const ll mod=998244353;
inline int read(){
int x=0; char s=gc;
while(!isdigit(s))s=gc;
while(isdigit(s))x=x*10+s-'0',s=gc;
return x;
}
ll ksm(ll a,ll b){
ll s=1;
while(b){
if(b&1)s=s*a%mod;
a=a*a%mod,b>>=1;
}
return s;
}
ll inv(ll x){return ksm(x,mod-2);}
const int N=1<<19;
const ll G=3;
const ll ivG=inv(3);
int r[N],pren;
void pre(int n){
if(n==pren)return;
for(int i=1;i<n;i++)r[i]=(r[i>>1]>>1)|(i&1?n>>1:0);
}
void NTT(int *g,int n,bool op){
pre(n);
static ull f[N],w[N]; w[0]=1;
for(int i=0;i<n;i++)f[i]=g[r[i]];
for(int l=1;l<n;l<<=1){
ull wn=ksm(op?G:ivG,(mod-1)/(l+l));
for(int i=1;i<l;i++)w[i]=w[i-1]*wn%mod;
for(int i=0;i<n;i+=l<<1)
for(int j=0;j<l;j++){
int t=w[j]*f[i|j|l]%mod;
f[i|j|l]=f[i|j]+mod-t,f[i|j]+=t;
}
if(l==(1<<16))for(int i=0;i<n;i++)f[i]%=mod;
}
if(op)for(int i=0;i<n;i++)g[i]=f[i]%mod;
else{
ll iv=inv(n);
for(int i=0;i<n;i++)g[i]=f[i]%mod*iv%mod;
}
}
int n,m,lim=1,a[N],b[N],res[N];
int ans[N],cnt;
string A,B;
int main(){
cin>>n>>m>>A>>B,reverse(A.begin(),A.end());
while(lim<m)lim<<=1;
for(int i=0;i<n;i++)if(A[i]!='*')a[i]=(A[i]-'a')*(A[i]-'a');
for(int i=0;i<m;i++)if(B[i]!='*')b[i]=1;
NTT(a,lim,1),NTT(b,lim,1);
for(int i=0;i<lim;i++)res[i]=1ll*a[i]*b[i]%mod;
mem(a,0,lim),mem(b,0,lim);
for(int i=0;i<n;i++)if(A[i]!='*')a[i]=A[i]-'a';
for(int i=0;i<m;i++)if(B[i]!='*')b[i]=B[i]-'a';
NTT(a,lim,1),NTT(b,lim,1);
for(int i=0;i<lim;i++)res[i]=(res[i]-2ll*a[i]*b[i]%mod+mod)%mod;
mem(a,0,lim),mem(b,0,lim);
for(int i=0;i<n;i++)if(A[i]!='*')a[i]=1;
for(int i=0;i<m;i++)if(B[i]!='*')b[i]=(B[i]-'a')*(B[i]-'a');
NTT(a,lim,1),NTT(b,lim,1);
for(int i=0;i<lim;i++)res[i]=(res[i]+1ll*a[i]*b[i])%mod;
NTT(res,lim,0);
for(int i=n-1;i<m;i++)if(!res[i])ans[++cnt]=i-n+2;
cout<<cnt<<endl;
for(int i=1;i<=cnt;i++)cout<<ans[i]<<" ";
return 0;
}
III. P4245 【模板】任意模數多項式乘法
const int N=(1<<18)+1;
const int B=32767;
struct cp{
double re,im;
cp operator + (cp x){return (cp){re+x.re,im+x.im};}
cp operator - (cp x){return (cp){re-x.re,im-x.im};}
cp operator * (cp x){return (cp){re*x.re-im*x.im,re*x.im+im*x.re};}
cp operator * (double x){return (cp){re*x,im*x};}
cp operator / (double x){return (cp){re/x,im/x};}
cp conj (){return (cp){re,-im};}
}I,w[N],a0[N],b0[N],a1[N],b1[N],P[N],Q[N];
int lim=1,r[N];
void init(){
for(int i=1;i<lim;i++)r[i]=(r[i>>1]>>1)|(i&1?lim>>1:0);
for(int i=0;i<=lim;i++)w[i]=(cp){cos(2*i*Pi/lim),sin(2*i*Pi/lim)};
}
void FFT(cp *a,int op){
static cp f[N];
for(int i=0;i<lim;i++)f[i]=a[r[i]];
for(int l=1;l<lim;l<<=1)
for(int i=0,b=lim/l>>1;i<lim;i+=l<<1)
for(int j=0,c=0;j<l;j++,c+=b){
cp k=f[i|j|l]*w[~op?c:lim-c];
f[i|j|l]=f[i|j]-k,f[i|j]=f[i|j]+k;
}
for(int i=0;i<lim;i++)a[i]=op==1?f[i]:f[i]/lim;
}
void FFFT(cp *a,cp *b){
for(int i=0;i<lim;i++)a[i].im=b[i].re; FFT(a,1);
for(int i=0;i<lim;i++)b[i]=a[i?lim-i:0].conj();
for(int i=0;i<lim;i++){
cp p=a[i],q=b[i];
a[i]=(p+q)*0.5,b[i]=(q-p)*0.5*I;
}
}
int n,m,p;
int main(){
cin>>n>>m>>p,I={0,1};
for(int i=0,v;i<=n;i++)v=read()%p,a0[i].re=v>>15,a1[i].re=v&B;
for(int i=0,v;i<=m;i++)v=read()%p,b0[i].re=v>>15,b1[i].re=v&B;
while(lim<=n+m)lim<<=1; init();
FFFT(a0,a1),FFFT(b0,b1);
for(int i=0;i<lim;i++){
P[i]=a0[i]*b0[i]+a0[i]*b1[i]*I,
Q[i]=a1[i]*b0[i]+a1[i]*b1[i]*I;
}
FFT(P,-1),FFT(Q,-1);
for(int i=0;i<=n+m;i++){
ll p1=(ll)(P[i].re+0.5)%p+p;
ll p2=(ll)(P[i].im+Q[i].re+0.5)%p+p;
ll p3=(ll)(Q[i].im+0.5)%p+p;
printf("%lld ",((p1<<30)+(p2<<15)+p3)%p);
}
return 0;
}
IV. P3723 [AH2017/HNOI2017]禮物
設第 \(1\) 個手環與第 \(2\) 個手環增加亮度的差為 \(c\)，那么有 \((x_i+c-y_i)^2=x_i^2+y_i^2+c^2+2x_ic-2y_ic-2x_iy_i\)，注意到如果我們枚舉 \(c\)，那么就要求 \(\sum x_iy_i\) 的最小值。進行翻轉后變成卷積形式，再破環成鏈即倍長數組，注意到值域 \(\leq 5\times 10^8\)，使用 NTT 即可。
const ll G=3;
const ll ivG=inv(3);
const int N=1<<18;
int n,m,len,r[N];
void init(int n){for(int i=1;i<n;i++)r[i]=(r[i>>1]>>1)|(i&1?n>>1:0);}
void NTT(ll *a,int n,bool op){
static ull f[N],w[N]; w[0]=1;
for(int i=0;i<n;i++)f[i]=a[r[i]];
for(int l=1;l<n;l<<=1){
ll wn=ksm(op?G:ivG,(mod-1)/(l+l));
for(int i=1;i<l;i++)w[i]=w[i-1]*wn%mod;
for(int i=0;i<n;i+=l+l)
for(int j=0;j<l;j++){
int y=f[i|j|l]*w[j]%mod;
f[i|j|l]=f[i|j]+mod-y,f[i|j]+=y;
} if(l==(1<<16))for(int i=0;i<n;i++)f[i]%=mod;
} if(!op){
ll iv=inv(n);
for(int i=0;i<n;i++)a[i]=f[i]%mod*iv%mod;
} else for(int i=0;i<n;i++)a[i]=f[i]%mod;
}
ll a[N],b[N];
ll ssq,ss,res,ans=1e10;
int main(){
cin>>n>>m;
for(int i=0;i<n;i++)scanf("%d",&a[i]),ssq+=a[i]*a[i],ss+=a[i];
for(int i=0;i<n;i++)scanf("%d",&b[i]),ssq+=b[i]*b[i],ss-=b[i],b[i+n]=b[i];
reverse(a,a+n);
len=1; while(len<3*n)len<<=1;
init(len),NTT(a,len,1),NTT(b,len,1);
for(int i=0;i<len;i++)a[i]=a[i]*b[i]%mod;
NTT(a,len,0);
for(int i=n-1;i<2*n;i++)res=max(res,a[i]);
for(int i=-m;i<=m;i++)ans=min(ans,ssq+2*ss*i+n*i*i-2*res);
cout<<ans<<endl;
return 0;
}