求 \(b(x)\) 滿足 \(a(x)b(x) \equiv 1 \pmod n\)。牛頓迭代得
\[b_(x) \equiv b_t(x)\left(2-a(x)b_t(x)\right)\pmod}}
\]樸素的實現需要做3次長度為 \(2^\) 的fft，把多餘的部分捨去，常數較大。
發現 \(b_(x)\) 的前 \(2^t\) 項和 \(b_t(x)\) 一樣，所以只需要求後 \(2^t\) 項，即求 \(a(x)b_t^2(x)\) 的 \(x^\dots x^-1}\) 項係數。
設\[a(x)b_t(x) \equiv 1 + b_0x^+b_1x^+\dots+b_x^-1} \pmod}}
\]這個結果的前 \(2^t\) 項確定是1了，所以只有後半部分是有用的。
由於 \(\deg b_t = 2^t\)，如果做長度為 \(2^\) 的卷積，多餘的部分會迴圈到前半部分，不會影響後半部分的結果。
同樣的， \(a(x)b_t^2(x) = a(x)b_t(x) \times b_t(x)\) ，卷積多餘的部分會迴圈到前 \(2^t\) 項，後半部分不會受到影響。
所以只需要做5次長度為 \(2^\) 的fft，在實際測試中（用100000的資料測試）常數約為正常寫法的三分之二。
下面這個是遞迴寫法：
polynom inverse(polynom a) ;
int m = n >> 1;
polynom b = inverse(polynom(a.begin(), a.begin() + m)), c = b;
b.resize(n);
dft(a), dft(b);
for (int i = 0; i < n; i++) a[i] = 1ll * a[i] * b[i] % p;
idft(a);
for (int i = 0; i < m; i++) a[i] = 0;
for (int i = m; i < n; i++) a[i] = p - a[i];
dft(a);
for (int i = 0; i < n; i++) a[i] = 1ll * a[i] * b[i] % p;
idft(a);
for (int i = 0; i < m; i++) a[i] = c[i];
return a;
}
求 \(b(x)\) 使 $a(x)\equiv b^2(x) \pmod $。牛頓迭代得
\[b_(x) \equiv b_t(x) - \frac \pmod}}
\]只要求 \((b_t^2(x) - a(x))\cdot(x)}\) 的後半部分即可。
容易發現 \(b_t^2(x) - a(x)\) 前 \(2^t\) 項是 \(0\)，所以 \(b_t^(x)\) 的後 \(2^t\) 項並不會對答案的後半部分產生貢獻，只求在模 \(x^\) 意義下的 \(b_t^(x)\) 即可。
對於 \(b_t^2(x)\)，由於前 \(2^t\) 項已知，所以只需要求長度為 \(2^t\) 的迴圈卷積，
實現時，不用每一層都呼叫多項式求逆，可以和開根一起迭代，具體看**。
實際測試常數優化效果顯著，可以在洛谷 p5205 的前幾個看見我的提交。
polynom sqrt(polynom a)
int len = a.size();
assert((len & len - 1) == 0);
assert(a[0] == 1); // warning: sqrtmod is needed if a[0] > 1.
polynom b(len), binv, bsqr; // sqrt, sqrt_inv, sqrt_sqr
polynom foo, bar; // temp
b[0] = 1;
auto shift = (int x) ; // quick div 2
for (int m = 1, n = 2; n <= len; m <<= 1, n <<= 1)
binv.resize(n);
dft(foo), dft(binv);
for (int i = 0; i < n; i++) foo[i] = 1ll * foo[i] * binv[i] % p;
idft(foo);
for (int i = m; i < n; i++) b[i] = shift(foo[i]);
// inv
if (n == len) break;
for (int i = 0; i < n; i++) foo[i] = b[i];
bar.resize(n), binv = bar;
dft(foo), dft(bar);
bsqr.resize(n);
for (int i = 0; i < n; i++) bsqr[i] = 1ll * foo[i] * foo[i] % p;
idft(bsqr);
for (int i = 0; i < n; i++) foo[i] = 1ll * foo[i] * bar[i] % p;
idft(foo);
for (int i = 0; i < m; i++) foo[i] = 0;
for (int i = m; i < n; i++) foo[i] = p - foo[i];
dft(foo);
for (int i = 0; i < n; i++) foo[i] = 1ll * foo[i] * bar[i] % p;
idft(foo);
for (int i = m; i < n; i++) binv[i] = foo[i];
} return b;
}
迴圈卷積與線性卷積
迴圈卷積 針對的是兩個長度都為n的序列，對兩個序列做fft，然後再做ifft得到的結果就是迴圈卷積，結果的長度也是n。直接計算步驟 序列a與序列b，長度都是n，新的序列c 1 把b倒過來。翻轉 2 把b向右平移乙個元素。最右側的元素補到左邊。3 計算此時a和b對應元素的積的和。將其加到c的末尾。4 ...
迴圈卷積 總結
一般的線性卷積 f i sum i a j b i j 如果將 b 陣列迴圈複製得到 b n 就能得到週期卷積 f i sum a j b n i j 而一般比較常見的迴圈卷積其實就是週期卷積的主值序列 0,n 1 項 f i sum a j b i j n 其中下標 n 表示其主值序列限定在 0,...
迴圈卷積和線性卷積的關係
迴圈卷積和線性卷積的關係 一般訊號處理濾波器 時域資料 與 濾波器係數的線性卷積，但卷積的運算量比較大，所以用頻域的相乘來替代時域卷積，而頻域的相乘等於時域迴圈卷積，所以要有乙個迴圈卷積和線性卷積轉換的過程，如下matlab例子 x 3 2 1 2 5 y 7 1 8 5 1 n length x ...