多项式基本算法
多项式乘积(NTT)、多项式逆元、多项式取对数、多项式取指数、多项式求导、多项式积分：
#pragma GCC optimize("O2")
#pragma GCC optimize("O3")
#pragma GCC optimize("Ofast")
#pragma GCC optimize("unroll-loops")
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
const int mod = 998244353;
const int D = 19;
const int N = (1 << D) + 1;
const int G = 3;
inline int add(int a, int b) { a += b; if (a >= mod) a -= mod; return a; }
inline int sub(int a, int b) { a -= b; if (a < 0) a += mod; return a; }
int qp(int a, int b) {
int res = 1;
while (b) {
if (b & 1) {
res = 1LL * res * a % mod;
}
a = 1LL * a * a % mod;
b >>= 1;
}
return res;
}
int rev[N], inv[N], invfac[N], fac[N], f[N], g[N], temp[N], temp2[N], temp3[N], GG[2][D][N >> 1];
const int invG = qp(G, mod - 2);
void init_GG() {
for (int p = 0; p < D; ++p) {
int buf0 = qp(G, (mod - 1) / (1 << p + 1));
int buf1 = qp(invG, (mod - 1) / (1 << p + 1));
GG[0][p][0] = GG[1][p][0] = 1;
for (int i = 1; i < 1 << p; ++i) {
GG[0][p][i] = 1LL * GG[0][p][i - 1] * buf0 % mod;
GG[1][p][i] = 1LL * GG[1][p][i - 1] * buf1 % mod;
}
}
}
void NTT_init(int n) {
for (int i = 0; i < n; ++i)
rev[i] = ((rev[i >> 1] >> 1) | ((i & 1) ? n >> 1 : 0));
}
void init_QAQ(int n) {
inv[1] = fac[0] = invfac[0] = 1;
for (int i = 1; i <= n; ++i)
fac[i] = 1LL * fac[i - 1] * i % mod;
invfac[n] = qp(fac[n], mod - 2);
for (int i = n - 1; i >= 1; --i)
invfac[i] = 1LL * invfac[i + 1] * (i + 1) % mod;
for (int i = 2; i <= n; i++) inv[i] = 1LL * (mod - mod / i) * inv[mod % i] % mod;
}
void NTT(int *c, int n, int op) {
for (int i = 0; i < n; i++)
if (i < rev[i]) swap(c[i], c[rev[i]]);
for (int d = 1, p = 0; d < n; d <<= 1, p++) {
for (int i = 0; i < n; i += d << 1) {
int *wn = GG[op][p];
for (int k = 0; k < d; k++, ++wn) {
int x = c[i + k], y = 1LL * (*wn) * c[i + d + k] % mod;
c[i + k] = add(x, y); c[i + d + k] = sub(x, y);
}
}
}
int invn = qp(n, mod - 2);
if (op) for (int i = 0; i < n; ++i) c[i] = 1LL * c[i] * invn % mod;
}
void Poly_Mul(int *a, int *b, int n, int m) {
int len = n + m - 1;
int lim = 1; while (lim < len) lim <<= 1;
for (int i = n; i < lim; i++) a[i] = 0;
NTT_init(lim);
for (int i = 0; i < lim; i++) temp[i] = (i < m ? b[i] : 0);
NTT(a, lim, 0); NTT(temp, lim, 0);
for (int i = 0; i < lim; i++) a[i] = 1LL * a[i] * temp[i] % mod;
NTT(a, lim, 1);
}
void Poly_Inv(int *a, int *b, int n) {
int lim = 1; while (lim < n) lim <<= 1;
for (int i = n; i < lim; i++) a[i] = 0;
for (int i = 0; i < lim; i++) b[i] = 0;
assert(a[0] != 0);
b[0] = qp(a[0], mod - 2);
for (int d = 2; d <= lim; d <<= 1) {
NTT_init(d << 1);
for (int i = 0; i < d; i++) temp[i] = a[i];
for (int i = d; i < d << 1; i++) b[i] = temp[i] = 0;
NTT(b, d << 1, 0); NTT(temp, d << 1, 0);
for (int i = 0; i < d << 1; i++) b[i] = 1LL * b[i] * sub(2LL, 1LL * temp[i] * b[i] % mod) % mod;
NTT(b, d << 1, 1);
for (int i = d; i < d << 1; i++) b[i] = 0;
}
for (int i = n; i < lim; i++) b[i] = 0;
}
void Poly_Dev(int *a, int n) {
for (int i = 1; i < n; i++)
a[i - 1] = 1LL * a[i] * i % mod;
a[n - 1] = 0;
}
void Poly_InvDev(int *a, int n) {
for (int i = n - 2; i >= 0; i--) a[i + 1] = 1LL * a[i] * inv[i + 1] % mod;
a[0] = 0;
}
void Poly_Ln(int *a, int n) {
Poly_Inv(a, temp3, n);
Poly_Dev(a, n);
Poly_Mul(a, temp3, n, n);
Poly_InvDev(a, n);
}
void Poly_Exp(int *a, int *b, int n) {
int lim = 1; while (lim < n) lim <<= 1;
for (int i = n; i < lim; i++) a[i] = 0;
for (int i = 0; i < lim; i++) b[i] = 0;
b[0] = 1;
for (int d = 2; d <= lim; d <<= 1) {
for (int i = 0; i < d; i++) temp2[i] = b[i];
Poly_Ln(temp2, d);
for (int i = 0; i < d; i++) temp[i] = a[i];
for (int i = d; i < d << 1; i++) b[i] = temp[i] = temp2[i] = 0;
NTT_init(d << 1); NTT(temp, d << 1, 0); NTT(temp2, d << 1, 0); NTT(b, d << 1, 0);
for (int i = 0; i < d << 1; i++)
b[i] = 1LL * b[i] * add(sub(1LL, temp2[i]), temp[i]) % mod;
NTT(b, d << 1, 1);
for (int i = d; i < d << 1; i++)
b[i] = 0;
}
for (int i = n; i < lim; i++) b[i] = 0;
}
void Poly_LeftMove(int *a, int n, int k) {
for (int i = n - 1; i >= k; i--)
a[i] = a[i - k];
for (int i = 0; i < k; i++) a[i] = 0;
}
void Poly_RightMove(int *a, int n, int k) {
for (int i = 0; i <= n - k - 1; i++)
a[i] = a[i + k];
for (int i = max(0, n - k); i < n; i++) a[i] = 0;
}
void solve() {
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
int _ = 1;
srand((i64)time(0));
init_GG();
init_QAQ(N - 1);
// std::cin >> _;
while (_--) {
solve();
}
return 0;
}
求第二类斯特林列固定的所有行的值，即 \(\text{S}(i,j)\)，其中 \(j\) 列固定不变：
第一种方法，多项式的 \(ln/exp\) 求值：
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
const int mod = 167772161;
const int D = 19;
const int N = (1 << D) + 1;
const int G = 3;
inline int add(int a, int b) { a += b; if (a >= mod) a -= mod; return a; }
inline int sub(int a, int b) { a -= b; if (a < 0) a += mod; return a; }
int qp(int a, int b) {
int res = 1;
while (b) {
if (b & 1) {
res = 1LL * res * a % mod;
}
a = 1LL * a * a % mod;
b >>= 1;
}
return res;
}
int rev[N], inv[N], invfac[N], fac[N], f[N], g[N], temp[N], temp2[N], temp3[N], GG[2][D][N >> 1];
const int invG = qp(G, mod - 2);
void init_GG() {
for (int p = 0; p < D; ++p) {
int buf0 = qp(G, (mod - 1) / (1 << p + 1));
int buf1 = qp(invG, (mod - 1) / (1 << p + 1));
GG[0][p][0] = GG[1][p][0] = 1;
for (int i = 1; i < 1 << p; ++i) {
GG[0][p][i] = 1LL * GG[0][p][i - 1] * buf0 % mod;
GG[1][p][i] = 1LL * GG[1][p][i - 1] * buf1 % mod;
}
}
}
void NTT_init(int n) {
for (int i = 0; i < n; ++i)
rev[i] = ((rev[i >> 1] >> 1) | ((i & 1) ? n >> 1 : 0));
}
void init_QAQ(int n) {
inv[1] = fac[0] = invfac[0] = 1;
for (int i = 1; i <= n; ++i)
fac[i] = 1LL * fac[i - 1] * i % mod;
invfac[n] = qp(fac[n], mod - 2);
for (int i = n - 1; i >= 1; --i)
invfac[i] = 1LL * invfac[i + 1] * (i + 1) % mod;
for (int i = 2; i <= n; i++) inv[i] = 1LL * (mod - mod / i) * inv[mod % i] % mod;
}
void NTT(int *c, int n, int op) {
for (int i = 0; i < n; i++)
if (i < rev[i]) swap(c[i], c[rev[i]]);
for (int d = 1, p = 0; d < n; d <<= 1, p++) {
for (int i = 0; i < n; i += d << 1) {
int *wn = GG[op][p];
for (int k = 0; k < d; k++, ++wn) {
int x = c[i + k], y = 1LL * (*wn) * c[i + d + k] % mod;
c[i + k] = add(x, y); c[i + d + k] = sub(x, y);
}
}
}
int invn = qp(n, mod - 2);
if (op) for (int i = 0; i < n; ++i) c[i] = 1LL * c[i] * invn % mod;
}
void Poly_Mul(int *a, int *b, int n, int m) {
int len = n + m - 1;
int lim = 1; while (lim < len) lim <<= 1;
for (int i = n; i < lim; i++) a[i] = 0;
NTT_init(lim);
for (int i = 0; i < lim; i++) temp[i] = (i < m ? b[i] : 0);
NTT(a, lim, 0); NTT(temp, lim, 0);
for (int i = 0; i < lim; i++) a[i] = 1LL * a[i] * temp[i] % mod;
NTT(a, lim, 1);
}
void Poly_Inv(int *a, int *b, int n) {
int lim = 1; while (lim < n) lim <<= 1;
for (int i = n; i < lim; i++) a[i] = 0;
for (int i = 0; i < lim; i++) b[i] = 0;
assert(a[0] != 0);
b[0] = qp(a[0], mod - 2);
for (int d = 2; d <= lim; d <<= 1) {
NTT_init(d << 1);
for (int i = 0; i < d; i++) temp[i] = a[i];
for (int i = d; i < d << 1; i++) b[i] = temp[i] = 0;
NTT(b, d << 1, 0); NTT(temp, d << 1, 0);
for (int i = 0; i < d << 1; i++) b[i] = 1LL * b[i] * sub(2LL, 1LL * temp[i] * b[i] % mod) % mod;
NTT(b, d << 1, 1);
for (int i = d; i < d << 1; i++) b[i] = 0;
}
for (int i = n; i < lim; i++) b[i] = 0;
}
void Poly_Dev(int *a, int n) {
for (int i = 1; i < n; i++)
a[i - 1] = 1LL * a[i] * i % mod;
a[n - 1] = 0;
}
void Poly_InvDev(int *a, int n) {
for (int i = n - 2; i >= 0; i--) a[i + 1] = 1LL * a[i] * inv[i + 1] % mod;
a[0] = 0;
}
void Poly_Ln(int *a, int n) {
Poly_Inv(a, temp3, n);
Poly_Dev(a, n);
Poly_Mul(a, temp3, n, n);
Poly_InvDev(a, n);
}
void Poly_Exp(int *a, int *b, int n) {
int lim = 1; while (lim < n) lim <<= 1;
for (int i = n; i < lim; i++) a[i] = 0;
for (int i = 0; i < lim; i++) b[i] = 0;
b[0] = 1;
for (int d = 2; d <= lim; d <<= 1) {
for (int i = 0; i < d; i++) temp2[i] = b[i];
Poly_Ln(temp2, d);
for (int i = 0; i < d; i++) temp[i] = a[i];
for (int i = d; i < d << 1; i++) b[i] = temp[i] = temp2[i] = 0;
NTT_init(d << 1); NTT(temp, d << 1, 0); NTT(temp2, d << 1, 0); NTT(b, d << 1, 0);
for (int i = 0; i < d << 1; i++)
b[i] = 1LL * b[i] * add(sub(1LL, temp2[i]), temp[i]) % mod;
NTT(b, d << 1, 1);
for (int i = d; i < d << 1; i++)
b[i] = 0;
}
for (int i = n; i < lim; i++) b[i] = 0;
}
void Poly_LeftMove(int *a, int n, int k) {
for (int i = n - 1; i >= k; i--)
a[i] = a[i - k];
for (int i = 0; i < k; i++) a[i] = 0;
}
void Poly_RightMove(int *a, int n, int k) {
for (int i = 0; i <= n - k - 1; i++)
a[i] = a[i + k];
for (int i = max(0, n - k); i < n; i++) a[i] = 0;
}
void init_S(int n, int m) {
for (int i = 0; i < n; i++) g[i] = 0;
if (n < m + 1) return;
for (int i = 0; i < n; i++)
f[i] = invfac[i + 1];
Poly_Ln(f, n);
for (int i = 0; i < n; i++)
f[i] = 1LL * f[i] * m % mod;
Poly_Exp(f, g, n);
Poly_LeftMove(g, n, m);
}
int S(int n, int m) {
int res = 1LL * fac[n] * invfac[m] % mod * g[n] % mod;
return res;
}
void solve() {
int n, m; cin >> n >> m;
init_S(n + 1, m);
for (int i = 0; i <= n; i++)
cout << S(i, m) << ' ';
cout << '\n';
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
int _ = 1;
srand((i64)time(0));
init_GG();
init_QAQ(N - 1);
// std::cin >> _;
while (_--) {
solve();
}
return 0;
}
不过这个常数有点大，时间复杂度 \(O(n\log^2 n)\)，挺慢的双 \(\log\)。
第二个方法，\(cdq\) 分治+NTT：
#pragma GCC optimize("O2")
#pragma GCC optimize("O3")
#pragma GCC optimize("Ofast")
#pragma GCC optimize("unroll-loops")
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
const int mod = 167772161;
const int D = 19;
const int N = (1 << D) + 1;
const int G = 3;
inline int add(int a, int b) { a += b; if (a >= mod) a -= mod; return a; }
inline int sub(int a, int b) { a -= b; if (a < 0) a += mod; return a; }
int qp(int a, int b) {
int res = 1;
while (b) {
if (b & 1) {
res = 1LL * res * a % mod;
}
a = 1LL * a * a % mod;
b >>= 1;
}
return res;
}
int rev[N], inv[N], invfac[N], fac[N], f[N], g[N], temp[N], GG[2][D][N >> 1], v[2][D][N];
const int invG = qp(G, mod - 2);
void init_GG() {
for (int p = 0; p < D; ++p) {
int buf0 = qp(G, (mod - 1) / (1 << p + 1));
int buf1 = qp(invG, (mod - 1) / (1 << p + 1));
GG[0][p][0] = GG[1][p][0] = 1;
for (int i = 1; i < 1 << p; ++i) {
GG[0][p][i] = 1LL * GG[0][p][i - 1] * buf0 % mod;
GG[1][p][i] = 1LL * GG[1][p][i - 1] * buf1 % mod;
}
}
}
void NTT_init(int n) {
for (int i = 0; i < n; ++i)
rev[i] = ((rev[i >> 1] >> 1) | ((i & 1) ? n >> 1 : 0));
}
void init_QAQ(int n) {
inv[1] = fac[0] = invfac[0] = 1;
for (int i = 1; i <= n; ++i)
fac[i] = 1LL * fac[i - 1] * i % mod;
invfac[n] = qp(fac[n], mod - 2);
for (int i = n - 1; i >= 1; --i)
invfac[i] = 1LL * invfac[i + 1] * (i + 1) % mod;
for (int i = 2; i <= n; i++) inv[i] = 1LL * (mod - mod / i) * inv[mod % i] % mod;
}
void NTT(int *c, int n, int op) {
for (int i = 0; i < n; i++)
if (i < rev[i]) swap(c[i], c[rev[i]]);
for (int d = 1, p = 0; d < n; d <<= 1, p++) {
for (int i = 0; i < n; i += d << 1) {
int *wn = GG[op][p];
for (int k = 0; k < d; k++, ++wn) {
int x = c[i + k], y = 1LL * (*wn) * c[i + d + k] % mod;
c[i + k] = add(x, y); c[i + d + k] = sub(x, y);
}
}
}
int invn = qp(n, mod - 2);
if (op) for (int i = 0; i < n; ++i) c[i] = 1LL * c[i] * invn % mod;
}
void Poly_Mul(int *a, int *b, int n, int m) {
int len = n + m - 1;
int lim = 1; while (lim < len) lim <<= 1;
for (int i = n; i < lim; i++) a[i] = 0;
NTT_init(lim);
for (int i = 0; i < lim; i++) temp[i] = (i < m ? b[i] : 0);
NTT(a, lim, 0); NTT(temp, lim, 0);
for (int i = 0; i < lim; i++) a[i] = 1LL * a[i] * temp[i] % mod;
NTT(a, lim, 1);
}
void Poly_Inv(int *a, int *b, int n) {
int lim = 1; while (lim < n) lim <<= 1;
for (int i = n; i < lim; i++) a[i] = 0;
for (int i = 0; i < lim; i++) b[i] = 0;
assert(a[0] != 0);
b[0] = qp(a[0], mod - 2);
for (int d = 2; d <= lim; d <<= 1) {
NTT_init(d << 1);
for (int i = 0; i < d; i++) temp[i] = a[i];
for (int i = d; i < d << 1; i++) b[i] = temp[i] = 0;
NTT(b, d << 1, 0); NTT(temp, d << 1, 0);
for (int i = 0; i < d << 1; i++) b[i] = 1LL * b[i] * sub(2LL, 1LL * temp[i] * b[i] % mod) % mod;
NTT(b, d << 1, 1);
for (int i = d; i < d << 1; i++) b[i] = 0;
}
for (int i = n; i < lim; i++) b[i] = 0;
}
void Poly_LeftMove(int *a, int n, int k) {
for (int i = n - 1; i >= k; i--)
a[i] = a[i - k];
for (int i = 0; i < min(n, k); i++) a[i] = 0;
}
void Poly_RightMove(int *a, int n, int k) {
for (int i = 0; i <= n - k - 1; i++)
a[i] = a[i + k];
for (int i = max(0, n - k); i < n; i++) a[i] = 0;
}
void cdq(int l, int r, int o, int dep) {
if (l == r) {
v[o][dep][0] = 1;
v[o][dep][1] = mod - l;
return;
}
int mid = (l + r) >> 1;
cdq(l, mid, 0, dep + 1); cdq(mid + 1, r, 1, dep + 1);
Poly_Mul(v[0][dep + 1], v[1][dep + 1], mid - l + 2, r - mid + 1);
for (int i = 0; i <= r - l + 1; i++) v[o][dep][i] = v[0][dep + 1][i];
}
void init_S(int n, int m) {
for (int i = 0; i < n; i++) g[i] = 0;
if (n < m + 1) return;
cdq(1, m, 0, 0);
for (int i = m + 1; i < n; i++) v[0][0][i] = 0;
Poly_Inv(v[0][0], g, n);
Poly_LeftMove(g, n, m);
}
void solve() {
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
int _ = 1;
srand(time(0));
init_GG();
init_QAQ(N - 1);
// std::cin >> _;
while (_--) {
solve();
}
return 0;
}
这个时间复杂度也是：\(O(n\log^2 n)\)，但是常数较小，比上面那个效率高。
浙公网安备 33010602011771号