多项式
算法
一些常用函数导数:
\((sin(x))'\ =\ cos(x)\)
\((cos(x))'\ =\ -sin(x)\)
\((a^x)'\ =\ a^x\ ln(a)\)
\(log_a(x)\ =\ \frac{1}{x\ ln(a)}\)
\((uv)'\ =\ u'v\ +\ uv'\)
\((\frac{u}{v})'\ =\ \frac{u'v\ -\ uv'}{v^2}\)
\(f(g(x))'\ =\ f'(g(x))\ g'(x)\)
求逆:
\(A(x)B(x)\ =\ 0\ \bmod\ x^n\)
\(B(x)\ -\ B_0(x)\ =\ 0\ \bmod\ x^{\frac{n}{2}}\)
\(B(x)\ =\ 2B_0(x)\ -\ B_0^2(x)\ A(x)\ \bmod\ x^n\)
除法:
\(A(x)\ =\ B(x)D(x)\ +\ R(x)\)
\(x^{n\ -\ 1}A(x^{-1})\ =\ x^{m\ -\ 1}\ B(x^{-1})x^{n\ -\ m}\ D(x^{-1})\ +\ x^{m\ -\ 2}\ R(x^{-1})\ x^{n\ -\ m\ +\ 1}\)
\(rA(x)\ =\ rB(x)\ =\ rD(x)\ +\ rR(x)\ x^{n\ -\ m\ +\ 1},\ C(x)\ rB(x)\ =\ 1\ \bmod\ x^{n\ -\ m\ +\ 1}\)
\(rD(x)\ =\ C(x)\ rA(x)\ \bmod\ x^{n\ -\ m\ +\ 1}\)
ln:
\(ln(1\ -\ A(x))\ =\ -\sum_{i\ \geq\ 1}\ \frac{A^i(x)}{i}\)
\(B(x)\ =\ ln(A(x))\)
\(B(x)'\ =\ (ln(A(x)))'\ =\ ln'(A(x))\ A(x)\ =\ \frac{A'(x)}{A(x)}\)
牛顿迭代:
\(F(x)\ =\ F_0(x)\ -\ \frac{G(F_0(x))}{G'(F_0(x))}\ \bmod\ x^n\)
exp:
\(B(x)\ =\ e^A(x)\ \bmod\ x^n\)
\(ln(B(x))\ =\ A(x)\ \bmod\ x^n\)
\(B(x)\ =\ B_0(x)(1\ -\ ln(B_0(x))\ +\ A(x))\ \bmod\ x^n\)
应用
例1.BZOJ3456 对于一种划分集合的方案,另 \(G(x)\) 为被划分的集合,\(F(x)\) 为划分的集合,有 \(e^{F(x)}\ =\ G(x)\)。
例2.多项式乘方 \(A(x)^N\ =\ e^{N\ ln(A(x))}\),注意要使得 \([x^0]A(x)\ =\ 1\)
代码
ntt版
fft版
ntt版
const int maxn = (1 << 19) | 10; // deal with n <= 1e5
const int mod = 998244353, proot = 3;
inline int pow_mod(int x, int n) { int y = 1; while(n) { if(n & 1) y = (long long) y * x % mod; x = (long long) x * x % mod; n >>= 1; } return y; }
namespace Poly {
struct poly {
vector<int> p;
poly() { p.clear(); }
poly(const vector<int> &q): p(q) {}
poly(int n, int *a) { p.resize(n); for (int i = 0; i < n; ++i) p[i] = a[i]; }
inline int size() const { return p.size(); }
inline void resize(int n) { p.resize(n); return; }
};
void dft(int n, int *a, bool rev) {
for (int i = 0, j = 0; i < n; ++i) {
if(i < j) swap(a[i], a[j]);
for (int k = n >> 1; (j ^= k) < k; k >>= 1);
}
for (int hl = 1, l = 2; l <= n; hl = l, l <<= 1) {
int wn = pow_mod(proot, (mod - 1) / l); if(rev) wn = pow_mod(wn, mod - 2);
for (int i = 0; i < n; i += l) for (int j = 0, *x = a + i, *y = x + hl, w = 1; j < hl; ++j, ++x, ++y, w = (long long) w * wn % mod) {
int t = (long long) *y * w % mod; *y = (*x - t) % mod; *x = (*x + t) % mod;
}
}
if(rev) { int inv = pow_mod(n, mod - 2); for (int i = 0; i < n; ++i) a[i] = (long long) a[i] * inv % mod; }
return;
}
poly operator * (const poly &A, const poly &B) {
static int a[maxn], b[maxn];
int n = A.size(), m = B.size();
if(n < 10 || m < 10 || n + m - 1 < 80) {
int N = n + m - 1; static int c[maxn]; memset(c, 0, sizeof(c[0]) * N);
for (int i = 0; i < n; ++i) for (int x = A.p[i], j = 0; j < m; ++j) c[i + j] = ((long long) x * B.p[j] + c[i + j]) % mod;
return poly(N, c);
}
int N = 1; while(N < n + m - 1) N <<= 1;
for (int i = 0; i < N; ++i) a[i] = i < n ? A.p[i] : 0, b[i] = i < m ? B.p[i] : 0;
dft(N, a, 0); dft(N, b, 0); for (int i = 0; i < N; ++i) a[i] = (long long) a[i] * b[i] % mod; dft(N, a, 1);
return poly(n + m - 1, a);
}
poly operator * (const int &a, const poly &B) {
static int b[maxn]; int n = B.size();
for (int i = 0; i < n; ++i) b[i] = (long long) a * B.p[i] % mod;
return poly(n, b);
}
poly operator - (const poly &A, const poly &B) {
static int a[maxn]; int n = A.size(), m = B.size(), N = max(n, m);
for (int i = 0; i < N; ++i) a[i] = ((i < n ? A.p[i] : 0) - (i < m ? B.p[i] : 0)) % mod;
return poly(N, a);
}
poly operator + (const poly &A, const poly &B) {
static int a[maxn]; int n = A.size(), m = B.size(), N = max(n, m);
for (int i = 0; i < N; ++i) a[i] = ((i < n ? A.p[i] : 0) + (i < m ? B.p[i] : 0)) % mod;
return poly(N, a);
}
poly inv(int n, const poly &A) { // A(x)^-1 mod x^n
if(n == 1) { return poly(vector<int>{pow_mod(A.p[0], mod - 2)}); }
static poly B0, B, tA;
B0 = inv(n + 1 >> 1, A);
if(A.size() < n) tA = A, tA.resize(n);
else tA.p.clear(), tA.p.insert(tA.p.end(), A.p.begin(), A.p.begin() + n);
B = 2 * B0 - B0 * B0 * tA;
B.resize(n);
return B;
}
poly rev(const poly &A) { static poly B; B = A; reverse(B.p.begin(), B.p.end()); return B; }
poly operator / (const poly &A, const poly &B) {
static poly rA, rB, C, D; int n = A.size(), m = B.size();
rA = rev(A); rB = rev(B); C = inv(n - m + 1, rB); D = C * rA; D.resize(n - m + 1);
return rev(D);
}
poly operator % (const poly &A, const poly &B) { static poly D, ret; D = A / B; ret = A - B * D; ret.resize(B.size() - 1); return ret; }
poly Rem(const poly &A, const poly &B, const poly &D) { static poly ret; ret = A - B * D; ret.resize(B.size() - 1); return ret; }
poly dao(const poly &A) {
static int a[maxn]; int n = A.size() - 1;
for (int i = 0; i < n; ++i) a[i] = (long long) A.p[i + 1] * (i + 1) % mod;
return poly(n, a);
}
poly ji(const poly &A) {
static int a[maxn]; int n = A.size();
for (int i = 1; i <= n; ++i) a[i] = (long long) A.p[i - 1] * pow_mod(i, mod - 2) % mod; a[0] = 0;
return poly(n + 1, a);
}
poly ln(int n, const poly &A) { // ln(A(x)) mod x^n
static poly B, C;
C = dao(A); B = inv(n, A);
B = B * C; B = ji(B); B.resize(n);
return B;
}
poly exp(int n, const poly &A) { // e^A(x) mod x^n
if(n == 1) { return poly(vector<int>{1}); }
static poly B0, B, C, tA;
B0 = exp(n + 1 >> 1, A);
if(A.size() < n) tA = A, tA.resize(n);
else tA.p.clear(), tA.p.insert(tA.p.end(), A.p.begin(), A.p.begin() + n);
C = ln(n, B0);
B = B0 * (poly(vector<int>{1}) - C + tA);
B.resize(n);
return B;
}
poly pow_mod(int n, const poly &A, int N) { // A(x)^N mod x^n
static poly B;
int na = A.size(), i = 0;
while(i < na && A.p[i] == 0) ++i;
if(i == na || (long long) N * i >= n) return poly(vector<int>(n, 0));
if(i) {
B.p.clear();
for (int j = i; j < na; ++j) B.p.push_back(A.p[j]);
static poly C;
C = pow_mod(n - N * i, B, N);
B.p = vector<int>(N * i, 0);
B.p.insert(B.p.end(), C.p.begin(), C.p.end());
return B;
}
if(A.p[0] != 1) {
int t = ::pow_mod(A.p[0], mod - 2), s = ::pow_mod(A.p[0], N); B.resize(na);
for (int i = 0; i < na; ++i) B.p[i] = (long long) A.p[i] * t % mod;
B = pow_mod(n, B, N);
for (int i = 0; i < n; ++i) B.p[i] = (long long) s * B.p[i] % mod;
return B;
}
B = ln(n, A);
for (int i = 0; i < n; ++i) B.p[i] = (long long) B.p[i] * N % mod;
return exp(n, B);
}
}
fft版
const int maxn = (1 << 19) | 10, lgn = 20; // deal with n <= 1e5
const int mod = 1e9 + 7;
inline int pow_mod(int x, int n) { int y = 1; while(n) { if(n & 1) y = (long long) y * x % mod; x = (long long) x * x % mod; n >>= 1; } return y; }
namespace Poly { // don't forget to init
struct comp {
double r, i;
comp() {}
comp(double r, double i): r(r), i(i) {}
} w[lgn][maxn];
comp operator * (const comp &a, const comp &b) { return comp(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r); }
comp operator + (const comp &a, const comp &b) { return comp(a.r + b.r, a.i + b.i); }
comp operator - (const comp &a, const comp &b) { return comp(a.r - b.r, a.i - b.i); }
comp conj(const comp &a) { return comp(a.r, -a.i); }
void init() {
static const double pi = acos(-1);
for (int i = 0; i < lgn; ++i) {
int n = 1 << i;
if(n >= maxn) break;
for (int j = 0; j < n; ++j) {
if(i && (~j & 1)) w[i][j] = w[i - 1][j >> 1];
else {
double ang = 2. * pi * j / (double) n;
w[i][j] = comp(cos(ang), sin(ang));
}
}
}
return;
}
inline int toint(double x, double invn) { return llround(x * invn) % mod; }
void dft(int n, comp *a) {
for (int i = 0, j = 0; i < n; ++i) {
if(i < j) swap(a[i], a[j]);
for (int k = n >> 1; (j ^= k) < k; k >>= 1);
}
for (int hl = 1, l = 2, wi = 1; l <= n; hl = l, l <<= 1, ++wi) for (int i = 0; i < n; i += l) {
comp *x = a + i, *y = x + hl, *z = w[wi];
for (int j = 0; j < hl; ++j, ++x, ++y, ++z) {
comp t = *y * *z; *y = *x - t; *x = *x + t;
}
}
return;
}
struct poly {
vector<int> p;
poly() { p.clear(); }
poly(const vector<int> &q): p(q) {}
poly(int n, int *a) { p.resize(n); for (int i = 0; i < n; ++i) p[i] = a[i]; }
inline int size() const { return p.size(); }
inline void resize(int n) { p.resize(n); return; }
};
poly operator * (const poly &A, const poly &B) {
int n = A.size(), m = B.size();
if(n < 10 || m < 10 || n + m - 1 < 80) {
int N = n + m - 1; static int c[maxn]; memset(c, 0, sizeof(c[0]) * N);
for (int i = 0; i < n; ++i) for (int x = A.p[i], j = 0; j < m; ++j) c[i + j] = ((long long) x * B.p[j] + c[i + j]) % mod;
return poly(N, c);
}
int N = 1; while(N < n + m - 1) N <<= 1;
static comp a[maxn], b[maxn], da[maxn], db[maxn], dc[maxn], dd[maxn];
for (int i = 0; i < n; ++i) { int x = (A.p[i] + mod) % mod; a[i] = comp(x & 32767, x >> 15); } for (int i = n; i < N; ++i) a[i] = comp(0, 0);
for (int i = 0; i < m; ++i) { int x = (B.p[i] + mod) % mod; b[i] = comp(x & 32767, x >> 15); } for (int i = m; i < N; ++i) b[i] = comp(0, 0);
dft(N, a); dft(N, b);
for (int i = 0; i < N; ++i) {
int j = N - i & N - 1; static comp DA, DB, DC, DD;
DA = (a[i] + conj(a[j])) * comp(.5, 0); DB = (a[i] - conj(a[j])) * comp(0, -.5);
DC = (b[i] + conj(b[j])) * comp(.5, 0); DD = (b[i] - conj(b[j])) * comp(0, -.5);
da[j] = DA * DC; db[j] = DA * DD;
dc[j] = DB * DC; dd[j] = DB * DD;
}
for (int i = 0; i < N; ++i) a[i] = da[i] + db[i] * comp(0, 1), b[i] = dc[i] + dd[i] * comp(0, 1);
dft(N, a); dft(N, b); double iN = 1.0 / N; N = n + m - 1; static int c[maxn];
for (int i = 0; i < N; ++i) {
static int da, db, dc, dd;
da = toint(a[i].r, iN), db = toint(a[i].i, iN), dc = toint(b[i].r, iN), dd = toint(b[i].i, iN);
c[i] = (((long long) dd << 30) + ((long long) db + dc << 15) + da) % mod;
}
return poly(N, c);
}
poly operator * (const int &a, const poly &B) {
static int b[maxn]; int n = B.size();
for (int i = 0; i < n; ++i) b[i] = (long long) a * B.p[i] % mod;
return poly(n, b);
}
poly operator - (const poly &A, const poly &B) {
static int a[maxn]; int n = A.size(), m = B.size(), N = max(n, m);
for (int i = 0; i < N; ++i) a[i] = ((i < n ? A.p[i] : 0) - (i < m ? B.p[i] : 0)) % mod;
return poly(N, a);
}
poly operator + (const poly &A, const poly &B) {
static int a[maxn]; int n = A.size(), m = B.size(), N = max(n, m);
for (int i = 0; i < N; ++i) a[i] = ((i < n ? A.p[i] : 0) + (i < m ? B.p[i] : 0)) % mod;
return poly(N, a);
}
poly inv(int n, const poly &A) { // A(x)^-1 mod x^n
if(n == 1) { return poly(vector<int>{pow_mod(A.p[0], mod - 2)}); }
static poly B0, B, tA;
B0 = inv(n + 1 >> 1, A);
if(A.size() < n) tA = A, tA.resize(n);
else tA.p.clear(), tA.p.insert(tA.p.end(), A.p.begin(), A.p.begin() + n);
B = 2 * B0 - B0 * B0 * tA;
B.resize(n);
return B;
}
poly rev(const poly &A) { static poly B; B = A; reverse(B.p.begin(), B.p.end()); return B; }
poly operator / (const poly &A, const poly &B) {
static poly rA, rB, C, D; int n = A.size(), m = B.size();
rA = rev(A); rB = rev(B); C = inv(n - m + 1, rB); D = C * rA; D.resize(n - m + 1);
return rev(D);
}
poly operator % (const poly &A, const poly &B) { static poly D, ret; D = A / B; ret = A - B * D; ret.resize(B.size() - 1); return ret; }
poly Rem(const poly &A, const poly &B, const poly &D) { static poly ret; ret = A - B * D; ret.resize(B.size() - 1); return ret; }
poly dao(const poly &A) {
static int a[maxn]; int n = A.size() - 1;
for (int i = 0; i < n; ++i) a[i] = (long long) A.p[i + 1] * (i + 1) % mod;
return poly(n, a);
}
poly ji(const poly &A) {
static int a[maxn]; int n = A.size();
for (int i = 1; i <= n; ++i) a[i] = (long long) A.p[i - 1] * pow_mod(i, mod - 2) % mod; a[0] = 0;
return poly(n + 1, a);
}
poly ln(int n, const poly &A) { // ln(A(x)) mod x^n
static poly B, C;
C = dao(A); B = inv(n, A);
B = B * C; B = ji(B); B.resize(n);
return B;
}
poly exp(int n, const poly &A) { // e^A(x) mod x^n
if(n == 1) { return poly(vector<int>{1}); }
static poly B0, B, C, tA;
B0 = exp(n + 1 >> 1, A);
if(A.size() < n) tA = A, tA.resize(n);
else tA.p.clear(), tA.p.insert(tA.p.end(), A.p.begin(), A.p.begin() + n);
C = ln(n, B0);
B = B0 * (poly(vector<int>{1}) - C + tA);
B.resize(n);
return B;
}
poly pow_mod(int n, const poly &A, int N) { // A(x)^N mod x^n
static poly B;
int na = A.size(), i = 0;
while(i < na && A.p[i] == 0) ++i;
if(i == na || (long long) N * i >= n) return poly(vector<int>(n, 0));
if(i) {
B.p.clear();
for (int j = i; j < na; ++j) B.p.push_back(A.p[j]);
static poly C;
C = pow_mod(n - N * i, B, N);
B.p = vector<int>(N * i, 0);
B.p.insert(B.p.end(), C.p.begin(), C.p.end());
return B;
}
if(A.p[0] != 1) {
int t = ::pow_mod(A.p[0], mod - 2), s = ::pow_mod(A.p[0], N); B.resize(na);
for (int i = 0; i < na; ++i) B.p[i] = (long long) A.p[i] * t % mod;
B = pow_mod(n, B, N);
for (int i = 0; i < n; ++i) B.p[i] = (long long) s * B.p[i] % mod;
return B;
}
B = ln(n, A);
for (int i = 0; i < n; ++i) B.p[i] = (long long) B.p[i] * N % mod;
return exp(n, B);
}
}