多项式板子
dft
可以在dft是用 \(unsigned\ long\ long\) 存下每个数，注意到每次一个数只会变最多 \(mod^2\)，故可以减少取模次数。
for (int hl = 1, l = 2, _ = 0; l <= n; hl = l, l <<= 1, ++_) {
if(!rev) {
for (int i = 0; i < n; i += l) {
unsigned long long *x = pool + i, *y = x + hl;
for (int j = 0, *w = w_pool[_]; j < hl; ++j, ++x, ++y, ++w) {
unsigned long long t = *y % mod * (unsigned long long) *w;
*y = *x + mod2 - t;
*x = *x + t;
if(_ == 17 && l < n) *x %= mod, *y %= mod;
}
}
}
}
求逆
令 \(N\ =\ 2n\)。
注意到 \(g_0\ =\ g\ \bmod\ x^n\)，故
考虑优化 \(f\cdot g_0\ \bmod\ x^N\)，\(f\) 是 \(N\) 次多项式，\(g_0\) 是 \(n\) 次多项式，所得结果为 \(3n\) 次多项式。由于 \(f\cdot g_0\ =\ 1\ \bmod\ x^n\)，故，只用求出 \([x^i]f\cdot g_0\ (n\leq i<N)\) 即可。做长度为 \(N\) 的 \(dft\) 即可，即求 \(f\cdot g_0\ \bmod\ x^N-1\)。
求 \((2\ -\ f\cdot g_0)\ g_0\) 可以用类似的方法进行优化。
t = inv(n, a);
for (int i = 0; i < n; ++i) foo[i] = t.p[i];
dft(N, foo, 0);
for (int i = 0; i < N; ++i) bar[i] = a.p[i];
dft(N, bar, 0);
for (int i = 0; i < N; ++i) bar[i] = (long long) bar[i] * foo[i] % mod;
dft(N, bar, 1);
for (int i = 0; i < n; ++i) bar[i] = 0;
for (int i = n; i < N; ++i) bar[i] = mod - bar[i];
dft(N, bar, 0);
for (int i = 0; i < N; ++i) foo[i] = (long long) foo[i] * bar[i] % mod;
dft(N, foo, 1);
for (int i = 0; i < n; ++i) foo[i] = t.p[i];
return poly(N, foo);
除法
令 \(A^R(x)\) 表示 \(A\) 系数逆序的多项式。
则有 $$AR(x) = BR(x)\cdot QR(x) + RR(x)\cdot x^{n-m+1}$$
令 \(a(x)\ =\ A^R(x),\ b(x)\ =\ B^R(x),\ N\ =\ n-m+1\)。
令 $n\ =\ \frac{N}{2},\ b_0(x)\ =\ b{-1}(x) \bmod xn,\ b_1(x)\ =\ \frac{b{-1}(x) - b_0(x)}{xn},\ $$a_0(x)\ =\ a(x)\ \bmod\ xn, a_1(x) = \frac{a(x) - a_0(x)}{xn}$
先求逆得出 \(b^{-1}(x)\)。
\(b_0\) 的 \(dft\) 结果在算 \(b^{-1}(x)\) 时已经计算过，故只需要算 \(a_0,\ a_1,\ b_1\) 的 \(dft\)。
b0 = b;
b0.resize(n);
b0 = inv(b0);
static int foo[maxn], bar[maxn];
memset(foo, 0, sizeof(foo[0]) * N);
for (int i = 0; i < n; ++i) foo[i] = b0.p[i];
dft(N, foo, 0); // b0
for (int i = 0; i < N; ++i) bar[i] = b.p[i];
dft(N, bar, 0);
for (int i = 0; i < N; ++i) bar[i] = (long long) bar[i] * foo[i] % mod;
dft(N, bar, 1);
for (int i = 0; i < n; ++i) bar[i] = 0;
for (int i = n; i < N; ++i) bar[i] = mod - bar[i];
dft(N, bar, 0);
for (int i = 0; i < N; ++i) bar[i] = (long long) foo[i] * bar[i] % mod;
dft(N, bar, 1);
for (int i = 0; i < n; ++i) bar[i] = bar[i + n], bar[i + n] = 0;
dft(N, bar, 0); // b1
static int Foo[maxn], Bar[maxn];
memset(Foo + n, 0, sizeof(Foo[0]) * n);
memset(Bar + n, 0, sizeof(Bar[0]) * n);
for (int i = 0; i < n; ++i) Foo[i] = a.p[i], Bar[i] = a.p[i + n];
dft(N, Foo, 0); // a0
dft(N, Bar, 0); // a1
for (int i = 0; i < N; ++i) bar[i] = (long long) bar[i] * Foo[i] % mod;
dft(N, bar, 1);
for (int i = 0; i < N; ++i) Bar[i] = (long long) Bar[i] * foo[i] % mod;
dft(N, Bar, 1);
for (int i = 0; i < N; ++i) foo[i] = (long long) foo[i] * Foo[i] % mod;
dft(N, foo, 1);
q.resize(N);
for (int i = 0; i < n; ++i) q.p[i] = foo[i], q.p[i + n] = ((bar[i] + Bar[i]) % mod + foo[i + n]) % mod;
q.resize(A.size() - B.size() + 1);
reverse(q.p.begin(), q.p.end());
return q;
ln
令 \(f_0\ =\ f^{-1}\ mod\ x^n,\ f_1\ =\ \frac{f\cdot f_0\ -\ 1}{x^n},\ f'_0\ =\ f'\ mod\ x^n,\ f'_1\ =\ \frac{f'\ -\ f'_0}{x^n}\)。
\(g\ =\ f'\cdot(1\ -\ (f\cdot f_0\ -\ 1))\cdot f_0\)
\(\ =\ (f'_0\ +\ f'_1\cdot x^n)(1\ -\ f_1\cdot x^n)\cdot f_0\)
\(\ =\ f'_0\cdot f_0\ +\ (f'_1\ -\ f'_0\cdot f_1)\cdot f_0\ x^n\bmod x^N\)
f.resize(N);
f0 = f;
f0.resize(n = N >> 1);
f0 = inverse(f0);
for (int i = 0; i < n; ++i) F0[i] = f0.p[i];
for (int i = 0; i < N; ++i) F1[i] = f.p[i];
dft(N, F0, 0);
dft(N, F1, 0);
for (int i = 0; i < N; ++i) F1[i] = (long long) F1[i] * F0[i] % mod;
dft(N, F1, 1);
for (int i = 0; i < n; ++i) F1[i] = F1[i + n], F1[i + n] = 0;
f = derivative(f);
for (int i = 0; i < n; ++i) DF0[i] = f.p[i], DF1[i] = (i == n - 1 ? 0 : f.p[i + n]);
dft(N, DF0, 0);
dft(N, F1, 0);
for (int i = 0; i < N; ++i) F1[i] = (long long) F1[i] * DF0[i] % mod;
dft(N, F1, 1);
for (int i = 0; i < n; ++i) F1[i] = (DF1[i] - F1[i]) % mod, F1[i + n] = 0;
dft(N, F1, 0);
for (int i = 0; i < N; ++i) F1[i] = (long long) F1[i] * F0[i] % mod;
dft(N, F1, 1);
for (int i = 0; i < N; ++i) F0[i] = (long long) F0[i] * DF0[i] % mod;
dft(N, F0, 1);
for (int i = 0; i < n; ++i) F0[i + n] = (F0[i + n] + F1[i]) % mod;
f = poly(N, F0);
f = integration(f);
f.resize(ff.size());
return f;
开方
\(g\ =\ f^{\frac{1}{2}}\bmod x^N,\ h\ =\ g^{-1}\bmod x^N,\ g_0\ =\ g\bmod x^n\)。
\(g\ =\ g_0\ -\ \frac{(g_0^2\ -\ f)\cdot h_0}{2}\)
计算 \(g_0^2\ -\ f\) 时只用求出 \(g_0\) 长度为 \(n\) 的dft，然后循环卷积一下，再减去 \(g_0\) 就行了。
算 \((g_0^2\ -\ f)\cdot h_0\) 同理。
每次求 \(g\) 时保留 \(h_0\) 长度为 \(N\) 的 \(dft\) 结果，优化求 \(h\)。
每次求 \(h\) 时需要计算 \(g\) 长度为 \(N\) 的 \(dft\) 结果，可以在算 \(2N\) 时的 \(g_0^2\) 时用到。
g0 = sqrt1(n, f, h, lastN);
for (int i = 0; i < n; ++i) dftg0[i] = (long long) dftg0[i] * dftg0[i] % mod;
dft(n, dftg0, 1);
for (int i = 0; i < n; ++i) {
dftg[i] = 0;
dftg[i + n] = (dftg0[i] - f.p[i]) % mod;
}
for (int i = n; i < N; ++i) dftg[i] = (dftg[i] - f.p[i]) % mod;
dft(N, dftg, 0);
for (int i = 0; i < n; ++i) dfth[i] = h.p[i], dfth[i + n] = 0;
dft(N, dfth, 0);
for (int i = 0; i < N; ++i) {
dftg[i] = (long long) dftg[i] * dfth[i] % mod;
dftg[i] = (dftg[i] & 1 ? dftg[i] + mod : dftg[i]) >> 1;
}
dft(N, dftg, 1);
g = g0;
g.resize(N);
for (int i = n; i < N; ++i) g.p[i] = -dftg[i];
if(N < lastN) {
for (int i = 0; i < N; ++i) dftg0[i] = g.p[i];
dft(N, dftg0, 0);
for (int i = 0; i < N; ++i) dftg[i] = (long long) dfth[i] * dftg0[i] % mod;
dft(N, dftg, 1);
for (int i = 0; i < n; ++i) dftg[i] = 0, dftg[i + n] = -dftg[i + n];
dft(N, dftg, 0);
for (int i = 0; i < N; ++i) dfth[i] = (long long) dfth[i] * dftg[i] % mod;
dft(N, dfth, 1);
h.resize(N);
for (int i = n; i < N; ++i) h.p[i] = dfth[i];
}
