ケロシのポリ - 多项式学习笔记
ケロシのポリ - 多项式学习笔记
持续更新中……
多项式乘法
给出 \(F(x)\) 和 \(G(x)\)，求出 \(H(x)=F(x)G(x)\)。
考虑用 NTT 求出点值，而点值可以直接相乘。
然后对乘出来的点值做一遍逆变换即可。
时间复杂度 \(O(n \log n)\)。
ケロシの代码
void solve() {
cin >> n >> m;
FOR(i, 0, n) cin >> a[i];
FOR(i, 0, m) cin >> b[i];
int lim = 1, l = 0;
while(lim <= (n + m)) lim <<= 1, l++;
REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
NTT(a, lim, 1);
NTT(b, lim, 1);
REP(i, lim) Mul(a[i], b[i]);
NTT(a, lim, -1);
FOR(i, 0, n + m) cout << a[i] << " ";
cout << endl;
}
多项式乘法逆
给出 \(F(x)\)，求 \(G(x)\) 使得 \(F(x)G(x) \equiv 1 (\bmod x^n)\)。
首先 \(G(x)_0=\frac{1}{F(x)_0}\)，然后考虑倍增，用 \(\bmod x^{ \left \lceil \frac{n}{2} \right \rceil }\) 的答案 \(G_0(x)\) 推 \(\bmod x^n\) 的答案：
需要用到 \(O(n \log n)\) 的多项式乘法，计算一下复杂度：
所以复杂度 \(T(n)=O(n \log n)\)。
ケロシの代码
int c[N];
void invs(int * a, int * b, int d) {
if(d == 1) {
b[0] = fp(a[0], P - 2);
return;
}
invs(a, b, (d + 1) >> 1);
int lim = 1, l = 0;
while(lim < (d << 1)) lim <<= 1, l ++;
REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
REP(i, d) c[i] = a[i];
FOR(i, d, lim - 1) c[i] = 0;
NTT(b, lim, 1);
NTT(c, lim, 1);
REP(i, lim) Mul(b[i], sub(2, mul(b[i], c[i])));
NTT(b, lim, -1);
FOR(i, d, lim - 1) b[i] = 0;
}
多项式开根
给出 \(F(x)\)，求 \(G(x)\) 使得 \(G^2(x) \equiv F(x) (\bmod x^n)\)。
先考虑常数项，那么 \(G_0(x)\) 就是 \(F_0(x)\) 的一个二次剩余。
然后还是考虑倍增，用 \(\bmod x^{ \left \lceil \frac{n}{2} \right \rceil }\) 的答案 \(G_0(x)\) 推 \(\bmod x^n\) 的答案：
时间复杂度 \(O(n \log n)\)。
多项式在 \(\ln\) 和 \(\exp\) 上的意义
一般来说，\(\ln\) 和 \(\exp\) 都是对于实数的，那么要计算多项式等形式的 \(\ln\) 和 \(\exp\) 需要用到泰勒公式进行展开：
我们把式子中的 \(x\) 替换成多项式即可使运算有意义。
为了使其不无限循环下去，式子中的多项式的 \(x_0\) 常数项必须为 \(0\)，使得 \(x^n\) 及以后 \(\bmod x^n\) 都等于 \(0\)。
所以当多项式 \(F(x)\) 满足 \(F(x)_0=1\) 时 \(\ln(F(x))\) 有意义，
当多项式 \(F(x)\) 满足 \(F(x)_0=0\) 时 \(\exp(F(x))\) 有意义。
多项式对数函数 \(\ln\)
给出 \(F(x)\)，求 \(G(x)\) 使得 \(G(x) \equiv \ln(F(x)) (\bmod x^n)\)。
对于 \(\ln(x)\)，注意到其求导后即为 \(\frac{1}{x}\)，比较好算，所以考虑两边求导。
对于复合函数 \(f(g(x))\)，其求导需要遵循链式法则：\(f(g(x))'=f'(g(x))g'(x)\)。
需要多项式求逆和多项式乘法，时间复杂度 \(O(n \log n)\)。
ケロシの代码
int p[N], q[N];
void lns(int * a, int * b, int d) {
int lim = 1, l = 0;
while(lim < (d << 1)) lim <<= 1, l++;
REP(i, lim) q[i] = 0;
invs(a, q, d);
FOR(i, 1, d - 1) p[i - 1] = mul(i, a[i]);
FOR(i, d - 1, lim - 1) p[i] = 0;
REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
NTT(p, lim, 1);
NTT(q, lim, 1);
REP(i, lim) Mul(p[i], q[i]);
NTT(p, lim, - 1);
FOR(i, 1, d) b[i] = mul(p[i - 1], inv[i]);
b[0] = 0;
}
多项式牛顿迭代
给出多项式 \(F(x)\)，求 \(F(x)\) 的零点 \(G(x)\)，即满足 \(F(G(x))\equiv 0(\bmod x^n)\)。
考虑使用牛顿迭代来进行倍增。
先单独考虑常数项 \(F(G(x)_0)_0 = 0\)。
接下来使用牛顿迭代，用 \(\bmod x^{ \left \lceil \frac{n}{2} \right \rceil }\) 的答案 \(G_0(x)\) 推 \(\bmod x^n\) 的答案。
考虑泰勒展开公式：
考虑将泰勒展开公式中的 \(x\) 和 \(x_0\) 分别代入 \(G(x)\) 和 \(G_0(x)\)，那么 \((G(x)-G_0(x))^2\) 及以后的项 \(\bmod x^n\) 都等于 \(0\)，所以直接代入泰勒展开公式：
因为 \(F(G(x)) \equiv 0 (\bmod x^n)\)
所以可以推出 \(G(x)\)：
多项式指数函数 \(\exp\)
给出 \(F(x)\)，求 \(G(x)\) 使得 \(G(x) \equiv \exp(F(x)) (\bmod x^n)\)。
考虑将求 \(exp\) 转换成关于 \(\ln\) 的零点，
那么 \(G(x)\) 就是函数 \(H(x)=\ln(x)+F\) 的零点。
接下来使用牛顿迭代，代入公式。
但是需要求得 \(H'(x)\)，因为自变量 \(x\) 是多项式，\(F\) 也是多项式，所以对于 \(x\) 来说 \(F\) 就是常数项，所以 \(H'(x)=\ln'(x)=\frac{1}{x}\)
然后直接代入公式即可：
需要用到多项式乘法，多项式乘法逆和多项式对数函数 \(\ln\)。
时间复杂度 \(O(n \log n)\)。
ケロシの代码
int h[N];
void exps(int * a, int * b, int d) {
if(d == 1) {
b[0] = 1;
return;
}
exps(a, b, (d + 1) >> 1);
REP(i, d) h[i] = 0;
lns(b, h, d);
REP(i, d) h[i] = sub(a[i], h[i]);
Add(h[0], 1);
int lim = 1, l = 0;
while(lim < (d << 1)) lim <<= 1, l++;
REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
FOR(i, d, lim - 1) h[i] = 0;
NTT(b, lim, 1);
NTT(h, lim, 1);
REP(i, lim) Mul(b[i], h[i]);
NTT(b, lim, - 1);
FOR(i, d, lim - 1) b[i] = 0;
}
多项式快速幂
给出 \(F(x)\) 和 \(k\)，求 \(G(x)\) 使得 \(G(x) \equiv F(x)^k (\bmod x^n)\)，保证 \(F(x)_0=1\)。
考虑两边同时取对数下放指数：
时间复杂度 \(O(n \log n)\)。
多项式连续点值平移
给定不超过 \(n\) 次多项式的 \(n+1\) 个点值：\(f(0),f(1),f(2),\cdots,f(n+1)\)，并且给定 \(m\)，满足 \(m>n\)，求出 \(f(m),f(m+1),f(m+2),\cdots,f(m+n)\)。
考虑拉格朗日插值求出 \(x>n\) 的点值：
考虑卷积优化，把式子拆成关于 \(i\)，\(x-i\) 和 \(x\) 的部分：
为了方便卷积，对齐下标：
拆成两个函数，注意下标必须要正的，所以将 \(x-i\) 改成 \(x+n-i\) 即可：
然后卷起来即可：
时间复杂度 \(O(n \log n)\)。
ケロシの代码
int n, m, a[N], F[N], G[N];
int fac[N], finv[N];
void solve() {
cin >> n >> m;
FOR(i, 0, n) cin >> a[i];
fac[0] = 1;
FOR(i, 1, n) fac[i] = mul(fac[i - 1], i);
finv[n] = fp(fac[n], P - 2);
ROF(i, n - 1, 0) finv[i] = mul(finv[i + 1], i + 1);
FOR(i, 0, n) F[i] = mul({a[i], finv[i], finv[n - i], sgn(n - i)});
FOR(i, 0, n << 1) G[i] = fp(m - n + i, P - 2);
int lim = 1, l = 0;
while(lim <= (n << 1)) lim <<= 1, l ++;
REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
NTT(F, lim, 1);
NTT(G, lim, 1);
REP(i, lim) Mul(F[i], G[i]);
NTT(F, lim, - 1);
int res = 1;
FOR(i, 0, n) Mul(res, m - i);
FOR(i, 0, n) {
cout << mul(F[n + i], res) << " ";
Mul(res, fp(m + i - n, P - 2));
Mul(res, m + i + 1);
}
cout << endl;
}
Chirp Z-Transform
给定一个 \(n\) 次多项式 \(F(x)\)，给定 \(c\) 和 \(m\)，求出 \(F(c^0),F(c^1),F(c^2),\cdots,F(c^{m-1})\)。
考虑：
我们希望用卷积加速这个式子，但是有有关 \(ij\) 的项，不好处理。
注意到 \({ij}\) 是在指数上，拆成加减后会变成乘法，比较好做，所以考虑拆解 \(ij\)：
但是这个式子的指数有除以二，需要二次剩余，不好做。
考虑另外一种拆法：
这样就都是整数了，代入式子：
现在就是一部分只和 \(j\) 有关，另一部分之和 \(i+j\) 有关，变成卷积形式了，直接卷积即可。
假设 \(n\) 与 \(m\) 同阶，时间复杂度 \(O(n \log n)\)。
ケロシの代码
int n, c, m;
int a[N], pw1[B + 1], pw2[B + 1];
int F[N], G[N];
inline int pw(int x) {
return mul(pw1[x % B], pw2[x / B]);
}
inline int C2(int x) {
return (1ll * x * (x - 1) / 2) % (P - 1);
}
void solve() {
cin >> n >> c >> m;
REP(i, n) cin >> a[i];
pw1[0] = pw2[0] = 1;
FOR(i, 1, B) pw1[i] = mul(pw1[i - 1], c);
FOR(i, 1, B) pw2[i] = mul(pw2[i - 1], pw1[B]);
REP(i, n + m) F[i] = pw(C2(i));
REP(i, n) G[i] = mul(a[i], pw(P - 1 - C2(i)));
reverse(G, G + n);
int lim = 1, l = 0;
while(lim <= (n + m)) lim <<= 1, l ++;
REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
NTT(F, lim, 1);
NTT(G, lim, 1);
REP(i, lim) Mul(F[i], G[i]);
NTT(F, lim, - 1);
REP(i, m) cout << mul(F[i + n - 1], pw(P - 1 - C2(i))) << " "; cout << endl;
}
任意模数多项式乘法
能做 NTT 的模数需要大量单位根，所以不是所有的数都能做 NTT 的模数。
考虑使用三个原根都是 \(3\) 的模数：\(998244353,469762049,1004535809\)。
对于每个模数都做多项式乘法，然后使用 CRT 把答案解出来即可。
时间复杂度 \(O(n \log n)\)。
ケロシの代码
const int N = 1e6 + 5;
const int P1 = 998244353;
const int P2 = 469762049;
const int P3 = 1004535809;
const int I1 = 208783132;
const int I2 = 395249030;
inline int add1(int x, int y) { return (x + y < P1 ? x + y : x + y - P1); }
inline void Add1(int & x, int y) { x = (x + y < P1 ? x + y : x + y - P1); }
inline int sub1(int x, int y) { return (x < y ? x - y + P1 : x - y); }
inline void Sub1(int & x, int y) { x = (x < y ? x - y + P1 : x - y); }
inline int mul1(int x, int y) { return (1ll * x * y) % P1; }
inline void Mul1(int & x, int y) { x = (1ll * x * y) % P1; }
inline int add2(int x, int y) { return (x + y < P2 ? x + y : x + y - P2); }
inline void Add2(int & x, int y) { x = (x + y < P2 ? x + y : x + y - P2); }
inline int sub2(int x, int y) { return (x < y ? x - y + P2 : x - y); }
inline void Sub2(int & x, int y) { x = (x < y ? x - y + P2 : x - y); }
inline int mul2(int x, int y) { return (1ll * x * y) % P2; }
inline void Mul2(int & x, int y) { x = (1ll * x * y) % P2; }
inline int add3(int x, int y) { return (x + y < P3 ? x + y : x + y - P3); }
inline void Add3(int & x, int y) { x = (x + y < P3 ? x + y : x + y - P3); }
inline int sub3(int x, int y) { return (x < y ? x - y + P3 : x - y); }
inline void Sub3(int & x, int y) { x = (x < y ? x - y + P3 : x - y); }
inline int mul3(int x, int y) { return (1ll * x * y) % P3; }
inline void Mul3(int & x, int y) { x = (1ll * x * y) % P3; }
int fp1(int x, int y) {
int res = 1;
for(; y; y >>= 1) {
if(y & 1) Mul1(res, x);
Mul1(x, x);
}
return res;
}
int fp2(int x, int y) {
int res = 1;
for(; y; y >>= 1) {
if(y & 1) Mul2(res, x);
Mul2(x, x);
}
return res;
}
int fp3(int x, int y) {
int res = 1;
for(; y; y >>= 1) {
if(y & 1) Mul3(res, x);
Mul3(x, x);
}
return res;
}
struct cint {
int x1, x2, x3;
cint operator + (const cint & A) {
return {add1(x1, A.x1), add2(x2, A.x2), add3(x3, A.x3)};
}
cint operator - (const cint & A) {
return {sub1(x1, A.x1), sub2(x2, A.x2), sub3(x3, A.x3)};
}
cint operator * (const cint & A) {
return {mul1(x1, A.x1), mul2(x2, A.x2), mul3(x3, A.x3)};
}
cint & operator += (const cint & A) {
Add1(x1, A.x1); Add2(x2, A.x2); Add3(x3, A.x3);
return * this;
}
cint & operator -= (const cint & A) {
Sub1(x1, A.x1); Sub2(x2, A.x2); Sub3(x3, A.x3);
return * this;
}
cint & operator *= (const cint & A) {
Mul1(x1, A.x1); Mul2(x2, A.x2); Mul3(x3, A.x3);
return * this;
}
int crt(int P) {
ll x4 = 1ll * P1 * mul2(sub2(x2, x1 % P2), I1) + x1;
int k4 = mul3(sub3(x3, x4 % P3), I2);
return (x4 + 1ll * k4 * (1ll * P1 * P2 % P)) % P;
}
};
int rev[N];
void NTT(cint * a, int lim, int o) {
REP(i, lim) if(i < rev[i]) swap(a[i], a[rev[i]]);
for(int i = 1; i < lim; i <<= 1) {
cint wn = {
fp1(o == 1 ? 3 : (P1 + 1) / 3, (P1 - 1) / (i << 1)),
fp2(o == 1 ? 3 : (P2 + 1) / 3, (P2 - 1) / (i << 1)),
fp3(o == 1 ? 3 : (P3 + 1) / 3, (P3 - 1) / (i << 1))
};
for(int j = 0; j < lim; j += (i << 1)) {
cint w = {1, 1, 1};
for(int k = 0; k < i; k ++, w *= wn) {
cint x = a[j + k];
cint y = a[j + k + i] * w;
a[j + k] = x + y;
a[j + k + i] = x - y;
}
}
}
if(o == 1) return;
cint inv = {fp1(lim, P1 - 2), fp2(lim, P2 - 2), fp3(lim, P3 - 2)};
REP(i, lim) a[i] *= inv;
}
int n, m, P; cint F[N], G[N];
void solve() {
cin >> n >> m >> P;
FOR(i, 0, n) {
int x; cin >> x; x %= P;
F[i] = {x, x, x};
}
FOR(i, 0, m) {
int x; cin >> x; x %= P;
G[i] = {x, x, x};
}
int lim = 1, l = 0;
while(lim <= (n + m)) lim <<= 1, l ++;
REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
NTT(F, lim, 1);
NTT(G, lim, 1);
REP(i, lim) F[i] *= G[i];
NTT(F, lim, - 1);
FOR(i, 0, n + m) cout << F[i].crt(P) << " "; cout << endl;
}
快速阶乘
给定 \(n\) 和 \(P\)，求 \(n! \bmod P\)。
考虑分块，设快长 \(B=\sqrt{n}\)，并考虑求解一段数连乘的值，考虑设多项式 \(F_d(x)=\prod_{i=1}^{d}(x+i)\)。
那么 \(n!=\prod_{i=0}^{B-1}F_B(iB)\prod_{i=B^2+1}^{n}i\)，也就是求出 \(F_B(x)\) 即可。
考虑求解 \(F_B(x)\)，但是在 \(F_B(x)\) 需要的下标都是 \(iB\)，这启发直接计算点值。
考虑倍增求解点值，考虑 \(F_d(iB)\) 如何转移到 \(F_{2d}(iB)\)，发现 \(F_{2d}(x)=F_d(x)F_d(x+d)\)。
所以直接将 \(F_d(iB)\) 进行连续点值平移到 \(F_d(iB+d)\) 也就是 \(F_d((i+\frac{d}{B})B)\) 即可。
但这样是需要两只 \(\log\)。
考虑缩减范围，\(F_d(x)\) 就是一个 \(d+1\) 次多项式，只保留 \(F_d(0),F_d(B),F_d(2B),\cdots,F_d(dB)\) 一共 \(d+1\) 即可。
倍增到 \(F_{2d}(x)\) 的时候再做一次点值平移得出 \(F_d((d+1)B),F_d((d+1)B),\cdots,F_d(2dB)\) 即可。
这样时间复杂度就是 \(O(\sqrt{n} \log n)\)。
ケロシの代码
int n, P, B;
int a[N], b[N];
cint F[N], G[N];
int fac[N], finv[N];
void moves(int * a, int * b, int n, int m) {
int lim = 1, l = 0;
while(lim <= (n << 1)) lim <<= 1, l ++;
REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
REP(i, lim) F[i] = G[i] = Cint(0);
FOR(i, 0, n) F[i] = Cint(mul({a[i], finv[i], finv[n - i], sgn(n - i)}));
FOR(i, 0, n << 1) G[i] = Cint(inv(sub(m, n - i)));
NTT(F, lim, 1);
NTT(G, lim, 1);
REP(i, lim) F[i] *= G[i];
NTT(F, lim, - 1);
int res = 1;
FOR(i, 0, n) Mul(res, sub(m, i));
FOR(i, 0, n) {
b[i] = mul(F[n + i].crt(P), res);
Mul(res, inv(sub(m, n - i)));
Mul(res, add(m, i + 1));
}
}
void facs(int d) {
if(d == 1) {
a[0] = 1; a[1] = B + 1;
return;
}
facs(d >> 1);
moves(a, b, d >> 1, (d >> 1) + 1);
REP(i, d) a[(d >> 1) + i + 1] = b[i];
moves(a, b, (d >> 1) << 1, mul(d >> 1, inv(B)));
FOR(i, 0, (d >> 1) << 1) Mul(a[i], b[i]);
if(d & 1) {
REP(i, d) Mul(a[i], i * B + d);
a[d] = 1;
FOR(i, 1, d) Mul(a[d], d * B + i);
}
}
void solve() {
cin >> n >> P;
B = sqrt(n);
fac[0] = 1;
FOR(i, 1, B) fac[i] = mul(fac[i - 1], i);
finv[B] = fp(fac[B], P - 2);
ROF(i, B - 1, 0) finv[i] = mul(finv[i + 1], i + 1);
facs(B);
int ans = 1;
REP(i, B) Mul(ans, a[i]);
FOR(i, B * B + 1, n) Mul(ans, i);
cout << ans << endl;
}
第二类斯特林数行
给出 \(n\)，对于 \(i\in[0,n]\)，求出 \(\begin{Bmatrix} n \\ i \end{Bmatrix}\)。
考虑第二类斯特林数的通项公式：
直接设多项式 \(F_i=\frac{i^n}{i!}\)，\(G_i=\frac{(-1)^{i}}{i!}\)，直接卷积求解即可。
时间复杂度 \(O(n \log n)\)。
ケロシの代码
int n;
int fac[N], finv[N];
int F[N], G[N];
void init(int n) {
fac[0] = 1;
FOR(i, 1, n) fac[i] = mul(fac[i - 1], i);
finv[n] = fp(fac[n], P - 2);
ROF(i, n - 1, 0) finv[i] = mul(finv[i + 1], i + 1);
}
void solve() {
cin >> n;
init(n);
int lim = 1, l = 0;
while(lim <= (n << 1)) lim <<= 1, l ++;
REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
FOR(i, 0, n) F[i] = mul(fp(i, n), finv[i]);
FOR(i, 0, n) G[i] = mul(sgn(i), finv[i]);
NTT(F, lim, 1);
NTT(G, lim, 1);
REP(i, lim) Mul(F[i], G[i]);
NTT(F, lim, - 1);
FOR(i, 0, n) cout << F[i] << " "; cout << endl;
}
第二类斯特林数列
给出 \(n\) 和 \(k\)，对于 \(i\in[0,n]\)，求出 \(\begin{Bmatrix} i \\ k \end{Bmatrix}\)。
考虑组合意义，先乘上 \(k!\) 变成 \(i\) 个有标号球放进 \(k\) 个有标号非空盒子里，考虑每次加一个盒子，相当于选择放多少个球再分配标号，不难发现 EGF 的乘法可以做这个。
所以设只放一个盒子的 EGF，除了放 \(0\) 个都有 \(1\) 的方案，所以设 EGF \(F=\sum_{i=1}\frac{x^i}{i!}=e^x-1\)。
所以答案的 EGF 即为 \(\frac{F^k}{k!}\)。
使用多项式快速幂即可。
时间复杂度 \(O(n \log n)\)。
ケロシの代码
int n, k;
int F[N], G[N];
void solve() {
cin >> n >> k;
if(n < k) {
REP(_, n + 1) cout << 0 << " ";
return;
}
init(n);
int len = n - k + 1;
REP(i, len) F[i] = finv[i + 1];
lns(F, G, len);
REP(i, len) Mul(G[i], k);
exps(G, F, len);
REP(_, k) cout << 0 << " ";
REP(i, len) cout << mul({F[i], fac[i + k], finv[k]}) << " "; cout << endl;
}
第一类斯特林数行
给出 \(n\)，对于 \(i\in[0,n]\)，求出 \(\begin{bmatrix} n \\ i \end{bmatrix}\)。
所以只要求 \(x^\overline{n}\) 即可，考虑倍增：
然后考虑 \((x+n)^\overline{n}\) 怎么求：
右边很明显是一个卷积形式，设 \(F_i=i! \begin{bmatrix} n \\ i\end{bmatrix}\) 和 \(G_{n-i}=\frac{n^i}{i!}\)，直接卷即可。
时间复杂度 \(O(n \log n)\)。
ケロシの代码
int n, a[N], b[N];
int fac[N], finv[N];
int F[N], G[N];
void init(int n) {
fac[0] = 1;
FOR(i, 1, n) fac[i] = mul(fac[i - 1], i);
finv[n] = fp(fac[n], P - 2);
ROF(i, n - 1, 0) finv[i] = mul(finv[i + 1], i + 1);
}
void slv(int d) {
if(d == 1) {
a[1] = 1;
return;
}
if(d & 1) {
slv(d - 1);
FOR(i, 0, d) F[i] = 0;
REP(i, d) Add(F[i + 1], a[i]);
REP(i, d) Add(F[i], mul(a[i], d - 1));
FOR(i, 0, d) a[i] = F[i];
return;
}
slv(d >> 1);
int len = d >> 1;
int lim = 1, l = 0;
while(lim <= d) lim <<= 1, l ++;
REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
REP(i, lim) F[i] = G[i] = 0;
FOR(i, 0, len) F[i] = mul(fac[i], a[i]);
FOR(i, 0, len) G[len - i] = mul(finv[i], fp(len, i));
NTT(F, lim, 1);
NTT(G, lim, 1);
REP(i, lim) Mul(F[i], G[i]);
NTT(F, lim, - 1);
REP(i, lim) b[i] = 0;
FOR(i, 0, len) b[i] = mul(finv[i], F[len + i]);
NTT(a, lim, 1);
NTT(b, lim, 1);
REP(i, lim) Mul(a[i], b[i]);
NTT(a, lim, - 1);
}
void solve() {
cin >> n;
init(n);
slv(n);
FOR(i, 0, n) cout << a[i] << " "; cout << endl;
}
第一类斯特林数列
给出 \(n\) 和 \(k\)，对于 \(i\in[0,n]\)，求出 \(\begin{bmatrix} i \\ k \end{bmatrix}\)。
考虑使用与第二类斯特林数列相同的做法，但是注意单个环的方案是 \((i-1)!\)，所以 EGF \(F=\sum_{i=1}\frac{(i-1)!x^i}{i!}=\sum_{i=1}\frac{x^i}{i}=\ln(\frac{1}{1-x})\)。
时间复杂度 \(O(n \log n)\)。
ケロシの代码
int n, k;
int F[N], G[N];
void solve() {
cin >> n >> k;
if(n < k) {
REP(_, n + 1) cout << 0 << " ";
return;
}
init(n);
int len = n - k + 1;
REP(i, len) F[i] = inv[i + 1];
lns(F, G, len);
REP(i, len) Mul(G[i], k);
exps(G, F, len);
REP(_, k) cout << 0 << " ";
REP(i, len) cout << mul({F[i], fac[i + k], finv[k]}) << " "; cout << endl;
}
浙公网安备 33010602011771号