多项式相关
慢慢更，想学的时候就学一点。给自己看的，随便写写。
拉格朗日插值法
给定 \(n\) 个点值 \(f(x_i)=y_i\)，我们断言，能够唯一确定一个 \(n-1\) 项的多项式 \(f\)。
基本形式
对于任意 \(x_0\)，
这是很好证明的：显然 \(f\) 为 \(n-1\) 项的多项式；注意到 \(\mathscr F_i(x_0)=\prod\limits_{i\ne j}\dfrac{x_0-x_j}{x_i-x_j}\) 满足 \(\mathscr F_i(x_i)=1\)、\(\mathscr F_i(x_j)=0\)（\(i\ne j\)），所以构造出的函数 \(f(x_0)=\sum y_i\mathscr F_i(x_0)\) 必定符合 \(n\) 个点值。显然多项式是唯一的，所以 \(f\) 当然也是确定的啦。
显然是 \(O(n^2)\) 的，至少比高斯消元快 /bx
\(\text{ }\)
点值连续的快速做法
......好吧挺唐氏的，就是在 \(x_i\) 有一定规律时，快速预处理出上文的 \(\mathscr F_i\) 即可，做到线性。
一个经典应用：CF622F。题意：求 \(\sum_{i=1}^n i^k\)。注意到这是一个关于 \(n\) 的 \(k+1\) 次多项式，所以只要随便求出 \(k+2\) 个点值（即 \(n\) 对应的 \(\sum_{i=1}^n i^k\)），然后插值一下就搞定。我们显然是取 \(0\le n\le k+1\) 的位置，不难得到 \(O(k\log k)\) 解法。
\(\text{ }\)
\(\text{ }\)
FMT&FWT
给定序列 \(f,g\)，求解 \(h\) 满足
其中 \(\otimes\) 符号可以是 OR，AND（FMT）或 XOR（FWT）。
核心思想：构造序列变换 \(a\to T(a)\)，满足 \(T(a)_i\times T(b)_i=T(c)_i\)，这样只需要对 \(f,g\) 进行变换后直接求出 \(T(h)\)，然后快速求出 \(T(h)\)，再进行逆变换 \(T(h)\to h\) 得到序列 \(h\)。
变换方式统一为：类似高维前缀和，枚举每一位 \(b\)，看 \(i,i\oplus 2^b\) 互相造成的贡献，然后递推。
三个操作复杂度均可以做到 \(O(n2^n)\)。构造方法如下（证明略，可以看 OI-Wiki）
OR
\(T(a)_i=\sum\limits_{j|i=i}a_j\)
正变换 \(op=1\)，逆变换 \(op=-1\)。
void OR(int f[], int op) {
for (int k = 1; k < (1 << n); k <<= 1) for (int i = 0; i < (1 << n); i += (k << 1)) for (int j = 0; j < k; j++)
f[i + j + k] = ((ll)f[i + j + k] + op * f[i + j] + mod) % mod;
}
AND
\(T(a)_i=\sum\limits_{j\&i=i}a_j\)
正变换 \(op=1\)，逆变换 \(op=-1\)。
void AND(int f[], int op) {
for (int k = 1; k < (1 << n); k <<= 1) for (int i = 0; i < (1 << n); i += (k << 1)) for (int j = 0; j < k; j++)
f[i + j] = ((ll)f[i + j] + op * f[i + j + k] + mod) % mod;
}
XOR
\(T(a)_i=\left(\sum_{\operatorname{popc}(i\& j)\equiv0\pmod2} a_j\right)-\left(\sum_{\operatorname{popc}(i\& j)\equiv1\pmod2} a_j\right)\)
正变换 \(op=1\)，逆变换 \(op=\frac12\)。
void XOR(int f[], int op) {
for (int k = 1; k < (1 << n); k <<= 1) for (int i = 0; i < (1 << n); i += (k << 1)) for (int j = 0; j < k; j++) {
int x = f[i + j], y = f[i + j + k];
f[i + j] = 1ll * op * (x + y) % mod, f[i + j + k] = 1ll * op * (x - y + mod) % mod;
}
}
\(\text{ }\)
\(\text{ }\)
子集卷积
本质只是一道 FMT 的练习题而已。。
给定序列 \(f,g\)，求解 \(h\) 满足
核心思想：\(j|k=i,j\&k=0\) 的一个等价描述是 \(j|k=i,\operatorname{popc}(j)+\operatorname{popc}(k)=\operatorname{pop}(j|k)\)。
popcount 是 \(O(n)\) 的，可以直接枚举 popcount 的值消除第二个限制，然后第一个限制 OR-FMT 解决。
复杂度 \(O(n^22^n)\)。
for (int i = 0; i < (1 << n); i++) scanf("%d", &a[popc[i]][i]);
for (int i = 0; i < (1 << n); i++) scanf("%d", &b[popc[i]][i]);
for (int i = 0; i <= n; i++) fmt(a[i], 1), fmt(b[i], 1);
for (int i = 0; i <= n; i++) for (int S = 0; S < (1 << n); S++) for (int j = 0; j <= i; j++)
ans[i][S] = (ans[i][S] + 1ll * a[j][S] * b[i - j][S]) % mod;
for (int i = 0; i <= n; i++) fmt(ans[i], -1);
for (int i = 0; i < (1 << n); i++) printf("%d ", ans[popc[i]][i]);
\(\text{ }\)
\(\text{ }\)
FFT
给定两个多项式 \(F,G\)，求它们的乘积，具体即求 \(H_i=\sum\limits_{j+k=i}F_jG_k\)。
核心思想：有点类似 fwt，考虑找到一个媒介做正逆变换：对 \(F\) 取出 \(n\) 个位置的点值从而唯一确定多项式，\(G\) 同理，点值可以直接乘得到 \(H\)，然后对 \(H\) 做点值到多项式的逆变换得到 \(H\) 多项式。
核心思想2：取出 \(1,\omega_n,\omega_n^2,\cdots,\omega_n^{n-1}\) 作为点值。
多项式->点值
现在考虑对多项式 \(F\) 转为点值。设其系数为 \(a_{0\sim n-1}\)，我们有
不妨认为 \(n\) 是偶数，考虑按下标奇偶性分类：
从中取出两个函数，我们有
现在我们要引入 \(\omega_n\) 了。这个是定义域在复数上的 \(n\) 次单位根，即 \(\omega_n^n=1\)。它的性质很多，首先将它作为复数放在二维坐标系上，相当于把圆等分为 \(2n\) 段取，所以 \(\omega_n^k=\cos k\cdot\frac{2\pi}n+\sin k\cdot\frac{2\pi}ni\)。其次，很难不注意到 \(\omega_n^{k+n}=\omega_n^k\)；在 \(n\) 为偶数时，\(\omega_n^k=\omega_{n/2}^{k/2}\)；\(\omega_n^{k+n/2}=-\omega_n^k\)。这些在数学中其实都是很基础的，但是我们已经可以完成 FFT 的步骤了。
我们代入 \(F(x=\omega_n^k)\)，得到
我们代入 \(F(x=\omega_n^{k+n/2})\)，得到
我去，有没有很神奇！！\(A,B\) 的对应位置竟然不变，也就是我们只要求出 \(A(\omega_{n/2}^k)\) 与 \(B(\omega_{n/2}^k)\) 的值，就能算出 \(F(\omega_n^k)\) 与 \(F(\omega_n^{k+n/2})\) 的值。这是一个子问题结构，直接递归进去，每次 \(n\) 会除以 \(2\)，yy 一下复杂度就是对的！\(O(n\log n)\)。
\(\text{ }\)
点值->多项式
草这个我真懒得看了，直接记吧，就是搞 \(op=-1\) 直接做 ntt 逆变换，然后最后给 \(a_i\) 都除以 \(n\) 就是多项式了。
\(\text{ }\)
实现技巧
递归比较慢，无法通过洛谷模板题。于是一个实现技巧是去除递归，就是注意到递归每次会把奇偶位置的值分别存到两个位置去，打表发现二进制位会有一个反转的操作，就像这样
（图片从洛谷题解盗的）
反正就是随便怎么搞搞？嘟嘟嘟嘟嘟嘟嘟嘟，我现在怎么这么魔怔了
void fft(int n, CPX A[], int op) {
for (int i = 0; i < n; i++) if (i < r[i]) swap(A[i], A[r[i]]);
for (int k = 1; k < n; k <<= 1) {
CPX wn = {cos(M_PI / k), op * sin(M_PI / k)};
for (int i = 0; i < n; i += (k << 1)) {
CPX w(1, 0); for (int j = 0; j < k; j++, w *= wn) {CPX x = A[i + j], y = w * A[i + j + k]; A[i + j] = x + y, A[i + j + k] = x - y;}
}
}
if (op == -1) {for (int i = 0; i < n; i++) A[i] /= n;}
}
其中，正变换 \(op=1\)，逆变换 \(op=-1\)。
有没有发现 fft 的代码有点类似 fwt 的！！太神奇了多项式。
哦对有个细节，由于有大量 \(n\to \frac n2\) 的操作，最好在开始时把多项式长度补全为 \(2\) 的幂次。
哦对再给一个求 \(r\) 的代码，\(r_i\) 即 \(i\) 的二进制反转：
//O(nlogn) 直接做
for (int i = 0; i < lim; i++) {r[i] = 0; for (int j = 0; j < w; j++) if (i >> j & 1) r[i] ^= (1 << (w - j - 1));}
//O(n), 应该也很好理解，我不写原因啦
for (int i = 0; i < lim; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (w - 1));
\(\text{ }\)
\(\text{ }\)
NTT
可以把 FFT 的模板认为模数为 \(998244353\) 然后做 NTT。
核心思想：算 \(\omega_n\) 要进复数太唐了，考虑在模意义下一个类似的东西叫原根。直接把单位根替换为原根就行，再稍微改改，就能在整数域下完成操作啦，无论是常数还是 int 与 double 的精度差距，NTT 都吊打 FFT。
void ntt(int n, int A[], int op) {
for (int i = 0; i < n; i++) if (i < r[i]) swap(A[i], A[r[i]]);
for (int k = 1; k < n; k <<= 1) {
int wn = qpow(op, (mod - 1) / (k << 1));
for (int i = 0; i < n; i += (k << 1)) for (int j = 0, w = 1; j < k; j++, w = 1ll * w * wn % mod) {
int x = A[i + j], y = 1ll * w * A[i + j + k] % mod; A[i + j] = (x + y) % mod, A[i + j + k] = (x - y + mod) % mod;
}
}
if (op != G) {int inv = qpow(n, mod - 2); for (int i = 0; i < n; i++) A[i] = 1ll * A[i] * inv % mod;}
}
其中，正变换 \(op=g\)，逆变换 \(op=g^{-1}\)，\(g\) 为原根。
注意模数 \(p=k2^n+1\) 时，需保证多项式长度 \(\le2^n\)。
- \(p=998244353\)，\(g=3\)；
- \(p=1004535809\)，\(g=3\)；
- 猫がかわいい
\(\text{ }\)
\(\text{ }\)
分治 FFT / NTT
就是半在线卷积，给定 \(g\)，求 \(f_i=\sum_{j=1}^i f_{i-j}g_j\)。
直接给 FFT / NTT 套一个 cdq 就行，没啥难度。2log
void solve(int l, int r) {
if (l == r) {if (!l) f[0] = 1; return;} int mid = (l + r) >> 1, _n = 0, _m = 0;
solve(l, mid);
for (int i = l, j = 0; i <= mid; i++, j++) a[j] = f[i], _n = j;
for (int i = 0, j = 0; i <= r - l; i++, j++) b[j] = g[i], _m = j;
int lim = 1, w = 0; while (lim <= r - l + 1) lim <<= 1, w++;
for (int i = 0; i < lim; i++) rv[i] = (rv[i >> 1] >> 1) | ((i & 1) << (w - 1));
for (int i = _n + 1; i < lim; i++) a[i] = 0; //remind: clear
for (int i = _m + 1; i < lim; i++) b[i] = 0;
ntt(lim, a, G), ntt(lim, b, G);
for (int i = 0; i < lim; i++) a[i] = 1ll * a[i] * b[i] % mod;
ntt(lim, a, iG);
for (int i = mid + 1; i <= r; i++) (f[i] += a[i - l]) %= mod;
solve(mid + 1, r);
}
\(\text{ }\)
\(\text{ }\)
多项式 inv
给定多项式 \(F(x)\)，求一个多项式 \(G(x)\) 使得 \(F(x)G(x)\equiv1\pmod {x^n}\)。
核心思想：倍增。
假设我们获取了 \(F(x)\bmod x^{n/2}\) 的逆元 \(G(x)\)，要求 \(x^n\) 的逆元 \(H(x)\)：
所以 \(H(x)=2G(x)-F(x)G^2(x)\)，可以递归求解。
\(T(n)=T(n/2)+\Theta(n\log n)\)，得算法复杂度 \(O(n\log n)\)。
void solve(int n) {
if (!n) return g[0] = qpow(f[0], mod - 2), void();
solve(n >> 1);
int lim = 1, w = 0; while (lim <= n * 2) lim <<= 1, w++;
for (int i = 0; i < lim; i++) rv[i] = (rv[i >> 1] >> 1) | ((i & 1) << (w - 1));
for (int i = 0; i <= n; i++) tf[i] = f[i], tg[i] = g[i];
for (int i = n + 1; i < lim; i++) tf[i] = tg[i] = 0;
ntt(lim, tf, G), ntt(lim, tg, G);
for (int i = 0; i < lim; i++) tg[i] = ((2ll * tg[i] - 1ll * tf[i] * tg[i] % mod * tg[i]) % mod + mod) % mod;
ntt(lim, tg, iG);
for (int i = 0; i <= n; i++) g[i] = tg[i];
}
这个是递归函数，但是递归层数太少了所以对常数没影响。
\(\text{ }\)
\(\text{ }\)
多项式 ln
求 \(G(x)\equiv\ln F(x)\pmod{x^n}\)。
核心思想：导完再积回去。
多项式的求导与积分都是平凡的；除法就求逆完成即可。
void dev(int n, int A[]) {for (int i = 1; i <= n; i++) a[i - 1] = 1ll * i * a[i] % mod; a[n] = 0;}
void ved(int n, int A[]) {for (int i = n - 1; ~i; i--) a[i + 1] = 1ll * inv[i + 1] * a[i] % mod; a[0] = 0;}
//in main():
get_inv(a, ia, n); dev(n, a);
int lim = 1, w = 0; while (lim <= n * 2) lim <<= 1, w++;
for (int i = 0; i < lim; i++) rv[i] = (rv[i >> 1] >> 1) | ((i & 1) << (w - 1));
ntt(lim, a, G), ntt(lim, ia, G);
for (int i = 0; i < lim; i++) a[i] = 1ll * a[i] * ia[i] % mod;
ntt(lim, a, iG); ved(n, a);
\(\text{ }\)
\(\text{ }\)
多项式 exp
求 \(G(x)\equiv \text e^{F(x)}\pmod{x^n}\)。
【警告⚠⚠⚠ 这里是垃圾 2log 做法 ⚠⚠⚠警告】
即
就是半在线卷积的形式，直接分治 NTT，\(O(n\log^2n)\)。代码不放了。
其实很难不发现，套用这个做法前面两个 inv 与 ln 都能做，但是毕竟复杂度更劣，还是太嘟嘟了。
\(\text{ }\)
\(\text{ }\)
多项式 qpow
给定多项式 \(F(x)\)，求 \(G(x)\) 满足 \(G(x)\equiv F(x)^k\pmod{x^n}\)。
特殊情况：[x^0]F(x)=1
即 \(a_0=1\)。
一个唐氏做法就是直接 qpow，是 \(O(n\log n\log k)\) 的。模板题 \(k\le10^{10^5}\)，所以要稍微处理一下指数，直接对 \(p\) 取模即可，注意不是 \(\varphi(p)\)。然后再稍微减少一下 NTT 调用次数，就可以 1.5s 内通过了。
为什么可以令 \(k\gets k\bmod p\)？
考虑证明 \(f(x)^p\equiv f(x^p)\pmod p\)。考虑归纳法，首先 \(a_0=1\)，故零次多项式肯定成立；然后设当前的 \(f(x)\) 为 \(n\) 次多项式，且其低 \(n-1\) 次多项式 \(g(x)\) 已经满足上述命题。提取出最高位 \(f(x)=g(x)+a_kx^k\)，简单推推式子\[f(x)^p\equiv\left(g(x)+a_kx^k\right)^p\equiv g(x)^p+a_k^px^{kp}\equiv g(x^p)+a_kx^{kp}=f(x^p) \]有点跳步：1->2 定义，2->3 二项式定理拆开有用项，3->4 对应一下项显然。
注意这是 2log 的，比较菜。感觉不如：
核心思想：ln 完再 exp 回去。
直接做，1log。我 exp 板子太烂是 2log 的
ln(n, lim);
for (int i = 0; i < n; i++) f[i] = 1ll * f[i] * k % mod;
for (int i = n; i <= lim; i++) f[i] = 0;
exp(n, lim);
\(\text{ }\)
任意情况
去除了 \(a_0=1\)。
这个看起来没B用的限制，实则使得我们上面那个 \(k\gets k\bmod p\) 的证明失效了！！
而且你 ln 的一个保证也是 \(a_0=1\)，也寄寄了！！
核心思想：强行使 \(a_0=1\)。
首先若 \(a_0=0\)，那就先一直找到最低的位置 \(a_t\ne0\)，然后降次一下。即 \(F(x)^k=\left(\frac{F(x)}{x^t}\right)^k\times x^{tk}\)；
那么此时默认 \(a_0\ne0\)。直接对多项式除以 \(a_0\) 即可。即 \(F(x)^k=\left(\frac{F(x)}{a_0}\right)^k\times a_0^k\)。
然后就可以套用上面做法啦。注意一堆细节：
- 若找不到 \(a_t\ne0\)，直接 return 全零即可。
- 若 \(tk\ge n\)，同样 return 全零。
- 多项式指数 mod \(p\)；数值的指数 mod \(p-1\)（比如此处的 \(a_0^k\)）。
- 可能要特判 \(k=0\)。
code：
int n, k_m, k_m1; ll real_k;
void fuck() {
char op = getchar(); while (op < 48 || op > 57) op = getchar();
while (48 <= op && op <= 57) {
k_m = (10ll * k_m + (op ^ 48)) % mod;
k_m1 = (10ll * k_m1 + (op ^ 48)) % (mod - 1);
real_k = min(10ll * real_k + (op ^ 48), (ll)1e10);
op = getchar();
}
}
int temp[N];
int main() {
scanf("%d", &n), fuck();
if (!real_k) {for (int i = 0; i < n; i++) printf("%d ", !i ? 1 : 0); return 0;}
for (int i = 0; i < n; i++) scanf("%d", &f[i]);
int t = 0; while (t < n && !f[t]) t++;
if (t == n || (__int128)t * real_k >= n) {while (n--) printf("0 "); return 0;}
for (int i = t; i < n; i++) temp[i - t] = f[i]; n -= t;
for (int i = 0; i < n; i++) f[i] = temp[i];
int val = qpow(f[0], k_m1), iv = qpow(f[0], mod - 2);
for (int i = 0; i < n; i++) f[i] = 1ll * f[i] * iv % mod;
inv[1] = 1; for (int i = 2; i <= n; i++) inv[i] = 1ll * inv[mod % i] * (mod - mod / i) % mod;
int lim = 1, w = 0; while (lim <= 2 * n) lim <<= 1, w++;
for (int i = 0; i < lim; i++) rv[i] = (rv[i >> 1] >> 1) | ((i & 1) << (w - 1));
ln(n, lim);
for (int i = 0; i < n; i++) f[i] = 1ll * f[i] * k_m % mod;
for (int i = n; i <= lim; i++) f[i] = 0;
exp(n, lim);
for (int _ = 1; _ <= t * real_k; _++) printf("0 ");
for (int i = 0; i < n + t - t * real_k; i++) printf("%d ", 1ll * ans[i] * val % mod);
return 0;
}
这个是傻逼洛谷模板题代码。因为太唐了我就不截取了，全部丢上来（
\(\text{ }\)
\(\text{ }\)
多项式 sqrt
给定 \(F(x)\)，求 \(G(x)\) 满足 \(G(x)^2\equiv F(x)\pmod{x^n}\)。
特殊情况：[x^0]F(x)=1
保证了 \(a_0=1\)，直接当作 \(\frac12\) 次数按上面 qpow 来做即可。\(O(n\log n)\)
ln(n, lim);
for (int i = 0; i < n; i++) f[i] = 1ll * f[i] * i2 % mod;
for (int i = n; i <= lim; i++) f[i] = 0;
exp(n, lim);
\(\text{ }\)
任意情况
我不会二次剩余啊 /kk
\(\text{ }\)
\(\text{ }\)
多项式(带余)除法
给定多项式 \(F(x),G(x)\)，求 \(Q(x),R(x)\) 满足 \(F(x)=G(x)Q(x)+R(x)\)。
\(F,G,Q,R\) 多项式次数分别为 \(n,m,n-m+1,m-1\)。（洛谷模板题中 \(n,m\) 加了一）
核心思想：转置。
考虑一个 "多项式系数反转" 的操作，即 \([x^{k}]\mathcal A(x)=[x^{n-k}]A(x)\)。容易编出等价描述
考虑算 \(Q\)。就把带余除法的式子对着这个转置做就行：
对 \(F,G\) 转置，\(G\) 求逆，算出 \(Q\) 的转置，再转回去得到 \(Q\)，再直接算得到 \(R\)。\(O(n\log n)\)
for (int i = 0; i < n; i++) fR[i] = f[n - i - 1];
for (int i = 0; i < m; i++) gR[i] = g[m - i - 1];
get_inv(n - m + 1, gR, inv_gR);
int lim = 1, w = 0; while (lim <= (n << 1)) lim <<= 1, w++;
for (int i = 0; i < lim; i++) rv[i] = (rv[i >> 1] >> 1) | ((i & 1) << (w - 1));
ntt(lim, fR, G), ntt(lim, inv_gR, G);
for (int i = 0; i < lim; i++) qR[i] = 1ll * fR[i] * inv_gR[i] % mod;
ntt(lim, qR, iG);
for (int i = n - m + 1; i < lim; i++) qR[i] = 0;
for (int i = 0; i < n - m + 1; i++) q[i] = qR[(n - m + 1) - i - 1];
ntt(lim, q, G), ntt(lim, g, G);
for (int i = 0; i < lim; i++) r[i] = 1ll * q[i] * g[i] % mod;
ntt(lim, q, iG), ntt(lim, r, iG);
for (int i = 0; i < m - 1; i++) r[i] = (f[i] - r[i] + mod) % mod;
\(\text{ }\)
\(\text{ }\)
下降幂多项式乘法
下降幂多项式形如 \(F(x)=\sum a_ix^{\underline i}\)。求 \(F(x)G(x)\) 得到的下降幂多项式。
对于多项式 \(F(x)\)，定义其 点值 EGF 为
我们在此代入 \(F(x)\) 的表达式，
也就是说，\(\hat F(x)\) 只需要用下降幂系数 \(a_i\) 对应的 OGF 乘上 \(\text e^x\) 即可得出。
注意 \(\hat F(x)\) 是点值，所以分别对 \(F,G\) 算出点值 EGF \(\hat F(x),\hat G(x)\)，直接对应位相乘就能得到 \(H\) 的点值 EGF \(\hat H(x)\)，然后逆变换（乘上 \(\text e^{-x}\)）得到 \(H\)。\(O(n\log n)\)
void o2e(int A[], int n) {
for (int i = 0; i < n; i++) t[i] = ifac[i];
mul(A, t, n, n);
for (int i = 0; i < n; i++) A[i] = 1ll * A[i] * fac[i] % mod;
for (int i = n; i < N; i++) A[i] = 0;
}
void e2o(int A[], int n) {
for (int i = 0; i < n; i++) A[i] = 1ll * A[i] * ifac[i] % mod;
for (int i = 0; i < n; i++) t[i] = (i % 2 == 0 ? ifac[i] : mod - ifac[i]);
mul(A, t, n, n);
for (int i = n; i < N; i++) A[i] = 0;
}
//in main()
o2e(a, n + m - 1), o2e(b, n + m - 1);
for (int i = 0; i < n + m - 1; i++) a[i] = 1ll * a[i] * b[i] % mod;
e2o(a, n + m - 1);
关于点值相乘：考虑 \(\hat F(x)\) 搞成 OGF 的形式的系数实际是 \(\frac{F(i)}{i!}\)，我们 EGF 肯定是转回去基本形式 OGF 再算。直接对每一项乘上 \(i!\) 就得到了真正点值 \(F(i)\)，现在可以直接对应相乘了；\(F(i)G(i)=H(i)\)。然后再给 \(H(i)\) 除以 \(i!\) 就得到了 \(\hat H(x)\) 的系数。
\(\text{ }\)
\(\text{ }\)
封装多项式
待测试，估计还会有问题的，别急。
Code
const int mod = 998244353, G = 3, iG = 332748118;
//todo: add "assert"
//todo: optimize space
namespace Shit {
inline int qpow(int x, int y) {int ans = 1; while (y) {if (y & 1) ans = 1ll * ans * x % mod; x = 1ll * x * x % mod; y >>= 1;} return ans;}
int lim, w; vector <int> rv;
void init(int n) {
lim = 1, w = 0; while (lim <= n) lim <<= 1, w++; rv.resize(lim);
for (int i = 0; i < lim; i++) rv[i] = (rv[i >> 1] >> 1) | ((i & 1) << (w - 1));
}
struct poly {
vector <int> A; int &operator[] (int x) {return A[x];}
poly(int n = 0) {A.shrink_to_fit(), A = vector<int>(n, 0);} poly(vector <int> T): A(T) {}
inline void maintain() {while (!A.empty() && A.back() == 0) A.pop_back();}
inline int size() {return A.size();}
inline void rev() {reverse(A.begin(), A.end());}
inline void dev() {int n = A.size(); for (int i = 1; i < n; i++) A[i - 1] = 1ll * i * A[i] % mod; A.pop_back();}
inline void ved() {int n = A.size(); A.push_back(0); for (int i = n - 1; ~i; i--) A[i + 1] = 1ll * qpow(i + 1, mod - 2) * A[i] % mod; A[0] = 0;}
void ntt(int op) {
A.resize(lim); for (int i = 0; i < lim; i++) if (i < rv[i]) swap(A[i], A[rv[i]]);
for (int k = 1; k < lim; k <<= 1) for (int i = 0, wn = qpow(op, (mod - 1) / (k << 1)); i < lim; i += (k << 1)) for (int j = 0, w = 1; j < k; j++, w = 1ll * w * wn % mod) {int x = A[i + j], y = 1ll * w * A[i + j + k] % mod; A[i + j] = (x + y) % mod, A[i + j + k] = (x - y + mod) % mod;}
if (op != G) {int inv = qpow(lim, mod - 2); for (int i = 0; i < lim; i++) A[i] = 1ll * A[i] * inv % mod;}
}
friend poly operator + (poly x, poly y) {int len = max(x.size(), y.size()); poly z(len); for (int i = 0; i < len; i++) z[i] = (x[i] + y[i]) % mod; return z;}
friend poly operator - (poly x, poly y) {int len = max(x.size(), y.size()); poly z(len); for (int i = 0; i < len; i++) z[i] = (x[i] - y[i] + mod) % mod; return z;}
friend poly operator * (poly x, poly y) {int len = x.size() + y.size() - 1; init(len); x.ntt(G), y.ntt(G); for (int i = 0; i < lim; i++) x[i] = 1ll * x[i] * y[i] % mod; x.ntt(iG), x.A.resize(len); return x;}
friend poly operator / (poly f, poly g) { //O((n+m)log(n+m))
f.maintain(), g.maintain();
int n = f.size() - g.size() + 1; poly fR = f, gR = g; fR.rev(), gR.rev();
gR.A.resize(n), gR.get_inv(); fR *= gR; fR.A.resize(n); fR.rev(); return fR;
}
// friend poly operator / (poly f, poly g) { //O(nm), use when m=O(1)
// poly ans; f.maintain(), g.maintain();
// for (int i = f.size() - 1, gsiz = g.size(); i >= gsiz - 1; i--) {
// int k = 1ll * f[i] * qpow(g[gsiz - 1], mod - 2) % mod; ans.A.push_back(k);
// for (int j = gsiz - 1, mi = i; ~j; j--, mi--) f[mi] = ((f[mi] - 1ll * k * g[j]) % mod + mod) % mod;
// }
// ans.rev(); return ans;
// }
friend poly operator % (poly f, poly g) {poly R = f - g * (f / g); R.A.resize(g.size() - 1); return R;} //O((n+m)log(n+m))
void ksm(int M1, int M2, ll M3) { //M1=y%(mod-1), M2=y%mod, M3=min(y,10^10)
int n = A.size();
if (!M3) {for (int i = 0; i < n; i++) A[i] = !i ? 1 : 0; return;}
int t = 0; while (t < n && !A[t]) t++; if (t == n || (__int128)t * M3 >= n) {A = vector<int>(n, 0); return;}
for (int i = t; i < n; i++) A[i - t] = A[i]; A.resize(n -= t);
int val = qpow(A[0], M1); poly F; *this /= A[0]; ln(), *this *= M2, exp();
for (int i = 0; i < t * M3; i++) F.A.push_back(0);
for (int i = 0; i < n + t - t * M3; i++) F.A.push_back(1ll * A[i] * val % mod);
*this = F;
}
friend poly operator ^ (poly x, ll y) {x.ksm(y % (mod - 1), y % mod, y); return x;}
friend poly& operator += (poly &x, const poly &y) {return x = x + y;}
friend poly& operator -= (poly &x, const poly &y) {return x = x - y;}
friend poly& operator *= (poly &x, const poly &y) {return x = x * y;}
friend poly& operator /= (poly &x, const poly &y) {return x = x / y;}
friend poly& operator %= (poly &x, const poly &y) {return x = x % y;}
friend poly& operator ^= (poly &x, const ll &y) {return x = x ^ y;}
friend poly operator * (poly x, int y) {for (int &v : x.A) v = 1ll * v * y % mod; return x;}
friend poly operator / (poly x, int y) {y = qpow(y, mod - 2); for (int &v : x.A) v = 1ll * v * y % mod; return x;}
friend poly& operator *= (poly &x, const int &y) {return x = x * y;}
friend poly& operator /= (poly &x, const int &y) {return x = x / y;}
void get_inv() { //mod x^n
int Sn = A.size(); poly _A(1), B(1); B[0] = qpow(A[0], mod - 2);
vector <int> vec; for (int i = Sn; i; i >>= 1) vec.push_back(i); reverse(vec.begin(), vec.end());
for (int n : vec) {
init(n << 1); _A.A.resize(n + 1); for (int i = 0; i <= n; i++) _A[i] = (i >= Sn ? 0 : A[i]);
_A.ntt(G), B.ntt(G); for (int i = 0; i < lim; i++) B[i] = ((2ll * B[i] - 1ll * _A[i] * B[i] % mod * B[i]) % mod + mod) % mod; B.ntt(iG), B.A.resize(n + 1);
}
*this = B, A.resize(Sn);
}
void ln() { //mod x^n
poly IA = *this; IA.get_inv(); dev(); init(A.size() << 1); *this *= IA; ved();
}
void __solve_for_exp(int l, int r, poly &B) {
static poly _A, _B;
if (l == r) {if (!l) B[0] = 1; else B[l] = 1ll * B[l] * qpow(l, mod - 2) % mod; return;}
int mid = (l + r) >> 1; __solve_for_exp(l, mid, B);
_A.A.clear(); for (int i = 0; i <= r - l; i++) _A.A.push_back(A[i]);
_B.A.clear(); for (int i = l; i <= mid; i++) _B.A.push_back(B[i]);
_B *= _A; for (int i = mid + 1; i <= r; i++) (B[i] += _B[i - l - 1]) %= mod;
__solve_for_exp(mid + 1, r, B);
}
void exp() { //mod x^n
int n = A.size(); poly B(n); dev(); __solve_for_exp(0, n - 1, B); *this = B;
}
};
inline poly rev(poly A) {A.rev(); return A;}
inline poly dev(poly A) {A.dev(); return A;}
inline poly ved(poly A) {A.ved(); return A;}
inline poly get_inv(poly A) {A.get_inv(); return A;}
inline poly ln(poly A) {A.ln(); return A;}
inline poly exp(poly A) {A.exp(); return A;}
} using Shit::poly;
浙公网安备 33010602011771号