多项式复合的 Kinoshita-Li 算法
x383494
·
·
题解
很遗憾，您的《多项式复合的 Kinoshita-Li 算法》不符合推荐标准．原因是：句子末尾应添加句号，且全文使用的句号应一致．
本次修改将打回理由中的句号也修改为实心句号．
以下是第一次修改．
很遗憾，您的《多项式复合的 Kinoshita-Li 算法》不符合推荐标准．原因是：文章不应当由 AI 生成，如有再犯直接撤销使用专栏权限．
文字说明与代码部分均为手写，本次修改加入了一些对具体实现的解释．
总不能是 AFO 之后工程代码写多了导致码风像 AI 吧．
以下是原文．
2024 新出的科技，由 木ノ下恭範（Yasunori Kinoshita）和 李白天（Baitian Li）提出，因此也叫木下-李合成．当时看到之后大为震撼，但是一直没有学．
orz noshi91 & EI．
一些常见记号
例如，$1 + x + x^2 \in \mathbb{R}[x]$，$\dfrac{1}{1-x} = 1 + x + x^2 + \cdots \in \mathbb{R}[[x]]$，$x^{-1} \in \mathbb{R}((x))$．
设 $P(x) \in R[[x]]$，记 $[x^n] P(x)$ 为取 $P(x)$ 中 $x^n$ 项的系数．
多项式 $P \in R[x]$ 的度数 $\deg P$ 定义为 $\max \{ n \in \mathbb{N} : [x^n] P \ne 0 \}$，即最高次项的次数．
## Bostan-Mori 算法
用于分式远端求值．
具体地，给定多项式 $P(x), Q(x) \in \mathbb{R}[x]$，想要求出 $[x^k]\dfrac{P(x)}{Q(x)}$，其中 $\max \{ \deg P, \deg Q\} \le d \ll k$．
注意到 $\dfrac{P(x)}{Q(x)} = \dfrac{P(x) Q(-x)}{Q(x) Q(-x)}$，而分母是偶函数．将分子的奇次项和偶次项拆开，写成 $P(x) Q(-x) = U_0(x^2) + x U_1(x^2)$，根据 $k$ 的奇偶性只需保留一边，于是有子问题 $[x^k]\dfrac{P(x)}{Q(x)} = [x^{\lfloor k/2 \rfloor}] \dfrac{U_{k \bmod 2}(x)}{Q(\sqrt x) Q(- \sqrt x)}$．由于 $Q(x) Q(-x)$ 仅有偶次项，子问题的分母仍然是次数不超过 $d$ 的多项式，分子类似．于是我们用两次多项式乘法的代价将 $k$ 的规模减少了一半．边界情况为 $[x^0] \dfrac{P(x)}{Q(x)} = \dfrac{P(0)}{Q(0)}$．
时间复杂度 $O(d \log d \log k)$．
## Kinoshita-Li 算法
提供一种不需要转置原理的思路．
给定 $f(x), g(x) \in \mathbb{R}[x], \max \{ \deg f, \deg g\} < n$，想要求出 $f(g(x)) \bmod x^n$．令 $Q(x) = f(x^{-1}) \in \mathbb{R}((x))$，即求 $\sum_{i=0}^{\deg f} f_i g^i(x) = \left( \sum_{i=0}^{\deg f} [y^0] f_i y^{-i}\left( \sum_{i = 0} y^i g^i(x) \right) \right) = [y^0]\dfrac{Q(y)}{1 - y g(x)}$．将分母记为 $G(x, y) := 1 - y g(x)$，分子乘上 $y^{\deg f}$ 变成 $P(y) = y^{\deg f} f(y^{-1}) \in \mathbb{R}[x]$，式子变成 $[y^{\deg f}] \dfrac{P(y)}{G(x, y)} \bmod x^n$．
模仿 Bostan-Mori 的做法，即求 $[y^{\deg f}] \dfrac{P(y)}{G(x, y) G(-x, y)} G(-x, y) \bmod x^n$，将问题归为 $\!\bmod x^{\lceil n/2 \rceil}$ 的情况，最后再乘上 $G(-x, y)$．实现时只需取 $y$ 的部分次方，这样每次 $x$ 规模减半，$y$ 规模加倍，总的规模仍是常数，$t = O(\log n)$ 次之后问题转化为 $\! \bmod x$ 的边界情况，即原式代入 $x = 0$．
边界情况有几种解决方法：
- 对初始的 $f$ 做多项式平移变为 $f(x + g(0))$，$f(x) \circ g(x) = f(x+c) \circ (x-c) \circ g(x) = f(x+c) \circ (g(x)-c)$．初始的 $G(x, y) = 1 - yg(x)$，而当 $g(0) = 0$ 时它就为 $1$，因此最后直接返回 $P(y)$ 即可．
- 直接对 $G_t(0, y)$ 求逆，再乘 $P(y)$．
- $G_t(0, y) = (1 - y g(0))^{2^t}$，可以二项式展开，再与 $P(y)$ 做多项式乘法．
现在还剩下一个问题：如何做二元多项式乘法？
设有两多项式 $f(x, y), g(x, y)$，$x$ 的最高次数之和不超过 $p$，那么可以作代换 $x = z^p, y = z$（或者相反），这保证了 $x$ 的进位不会影响到 $y$．然后就是一元的情况．
## 实现
一份 [常数比较大的实现](https://www.luogu.com.cn/record/259450923)，~~现已加入多项式全家桶豪华套餐~~．
网上常见的写法是将 $y$ 代换成 $x$ 的幂，这个实现尝试将 $x$ 代换成 $y$ 的幂．
具体地，维护递归深度 `dep`，每层 $\deg_y G(x, y) \le 2^{dep}$，因此作代换 $x = y^{2^{dep+1}}$ 之后计算 $G(x, y) G(-x, y)$．虽然 $y$ 相乘之后的次数能到 $2^{dep+1}$，但是注意到 $x$ 在相乘之后只有偶数次项，因此仍然可以分离．从下一层递归的视角来看也有 $x = y^{2^{dep+2}}$，相当于同时将 $x$ 的次数减半了．
将 $x$ 代换成 $y$ 的幂之后便于对 $x^n$ 取模，但是需要特判一下 $y$ 相乘之后进位的情况．
```cpp
#include <algorithm>
#include <cassert>
#include <cctype>
#include <cmath>
#include <functional>
#include <iostream>
#include <vector>
#define UP(i, s, e) for(auto i = s; i < e; ++i)
#define DOWN(i, e, s) for(auto i = e; i-- > s;)
using std::cin;
using std::cout;
using u32 = unsigned int;
using u64 = unsigned long long;
constexpr u32 MOD = 998244353, G = 3, TIM = 23;
template<typename T, typename U>
constexpr T qpow(T x, U tim) {
T ans = 1;
while(tim) {
if(tim & 1) ans *= x;
x *= x;
tim >>= 1;
}
return ans;
}
struct m32 {
u32 val;
m32() {}
constexpr m32(u32 x)
: val(x) {}
constexpr bool operator==(m32 b) const { return val == b.val; }
constexpr bool operator!=(m32 b) const { return val != b.val; }
constexpr m32 operator+=(m32 b) {
val += b.val;
val %= MOD;
return *this;
}
constexpr m32 operator*=(m32 b) {
val = (u64)val * b.val % MOD;
return *this;
}
constexpr m32 operator-=(m32 b) {
val += MOD - b.val;
val %= MOD;
return *this;
}
#define ADDOP(op) \
constexpr m32 operator op(m32 b) { \
m32 tmp = *this; \
return tmp op## = b; \
}
ADDOP(+);
ADDOP(-);
ADDOP(*);
};
using CC = m32;
constexpr u32 bitceil(u32 x) {
return x <= 2 ? 2u : (2u << std::__lg(x - 1));
}
constexpr CC getinv(CC x) {
return qpow(x, MOD - 2);
}
CC factorial(u32 x, bool inv = false) {
static std::vector<CC> fac(1, 1), invfac(1, 1);
u32 siz = fac.size();
if(x < fac.size()) {
return inv ? invfac[x] : fac[x];
}
fac.resize(x + 1);
invfac.resize(x + 1);
UP(i, siz, x + 1) {
fac[i] = fac[i - 1] * CC(i);
}
invfac[x] = getinv(fac[x]);
DOWN(i, x + 1, siz) {
invfac[i - 1] = invfac[i] * CC(i);
}
return inv ? invfac[x] : fac[x];
}
CC binom(u32 x, u32 y) {
assert(x >= y);
return factorial(x) * factorial(y, true) * factorial(x - y, true);
}
namespace Poly { // }{{{
constexpr CC ROOT = qpow(CC(G), (MOD - 1) >> TIM), INVROOT = getinv(ROOT);
template<typename T>
void change(T *arr, int len) {
static std::vector<u32> rev;
if((int)rev.size() != len) {
rev.resize(len);
rev[0] = 0;
UP(i, 1, len) {
rev[i] = rev[i >> 1] >> 1;
if(i & 1) rev[i] |= len >> 1;
}
}
UP(i, 0, len)
if((u32)i < rev[i]) std::swap(arr[i], arr[rev[i]]);
}
void ntt(CC *arr, u32 len, bool idft) {
change(arr, len);
for(u32 l = 2; l <= len; l <<= 1) {
CC wn = qpow(idft ? INVROOT : ROOT, (1 << TIM) / l);
for(u32 st = 0; st < len; st += l) {
CC w = 1;
UP(j, st, st + l / 2) {
CC u = arr[j], v = arr[j + l / 2];
v *= w;
arr[j] = u + v;
arr[j + l / 2] = u - v;
w *= wn;
}
}
}
if(idft) {
CC invlen = getinv(len);
UP(i, 0u, len)
arr[i] *= invlen;
}
}
void polymul(CC *a, CC *b, u32 len) {
ntt(a, len, false);
ntt(b, len, false);
UP(i, 0u, len)
a[i] *= b[i];
ntt(a, len, true);
}
void integrate(CC *arr, u32 len) {
DOWN(i, len, 1u)
arr[i] = arr[i - 1] * getinv(CC(i));
arr[0] = 0;
}
void derivative(CC *arr, u32 len) {
UP(i, 0u, len - 1)
arr[i] = arr[i + 1] * CC(i + 1);
arr[len - 1] = 0;
}
// return f(x+c) in O(n log^2 n)
void polyshift(CC *arr, u32 len, CC c, std::vector<CC> &ans) {
if(len == 1) {
ans.clear();
ans.push_back(arr[0]);
return;
}
std::vector<CC> ans1, tmp;
u32 siz = bitceil(len);
ans.reserve(siz);
polyshift(arr, len / 2, c, ans1);
polyshift(arr + len / 2, len - len / 2, c, ans);
ans.resize(siz, 0);
tmp.resize(siz, 0);
CC pow_c = CC(1);
DOWN(i, len / 2 + 1, 0) {
tmp[i] = pow_c * binom(len / 2, i);
pow_c *= c;
}
polymul(ans.data(), tmp.data(), siz);
ans.resize(len);
UP(i, 0u, len / 2)
ans[i] += ans1[i];
}
// G_x_y: x = y^{2 * 2^dep}, deg_y G_x_y = 2^dep
// calculate [y^l ... y^{r-1}]P(y)/G(x, y) mod x^n
// return: x = y^{4 * 2^dep}, mod x^n
std::vector<CC> kinoshita_li(int l, int r, std::vector<CC> &P_y, std::vector<CC> G_x_y, CC g_0, int dep, int n) {
std::vector<CC> ans;
if(n == 1) {
P_y.resize(r, 0);
u32 len = bitceil(2*r);
P_y.resize(len, 0);
ans.reserve(len);
factorial((1<<dep) + r - 2); // for cache
CC g_i = 1;
UP(i, 0, r){
ans.push_back(binom((1<<dep)+i-1, i) * g_i);
g_i *= g_0;
}
ans.resize(len, 0);
polymul(ans.data(), P_y.data(), len);
ans.erase(ans.begin() + r, ans.end());
ans.erase(ans.begin(), ans.begin() + l);
return ans;
}
std::vector<CC> G_pos(G_x_y), G_neg(G_x_y); // G(x, y), G(-x, y), where x = y^{2 * 2^dep}
u32 len = bitceil(n << (dep + 2));
G_pos.resize(len, 0);
G_neg.resize(len, 0);
UP(i, 0, n) {
if(i % 2 == 0) continue;
UP(j, 0, (1 << dep) + 1) {
G_neg[j + (i << (dep + 1))] = CC(0) - G_neg[j + (i << (dep + 1))];
}
}
polymul(G_pos.data(), G_neg.data(), len); // neg is NTTed, pos is passed on
// G_pos: x^2 = y^{2 * 2^{dep+1}} and only x^2 remains
G_pos.resize((n << (dep + 1)) + 1); // mod (y * x^n) = y^(n * 2^(dep+1) + 1) to protect x^{n-1} y^{2^dep} (where x = y^{2^dep})
if(n % 2 == 0) {
// in this case x^{n-1} y^{2^dep} is 0, and we need to erase x^n y^0 term
G_pos[n << (dep + 1)] = 0;
}
int l_nxt = std::max(0, (int)l - (1 << dep));
std::vector<CC> G_rec = kinoshita_li(l_nxt, r, P_y, G_pos, g_0, dep + 1, (n + 1) / 2);
l -= l_nxt, r -= l_nxt;
G_neg.clear();
len = bitceil(n << (dep + 3));
G_neg.resize(len, 0);
UP(i, 0, n) {
UP(j, 0, (1 << dep) + 1) {
CC val = G_x_y[(i << (dep + 1)) + j];
G_neg[(i << (dep + 2)) + j] = (i % 2 == 0 ? val : CC(0) - val);
}
}
G_rec.resize(len, 0);
polymul(G_rec.data(), G_neg.data(), len);
G_rec.resize(n << (dep + 2)); // mod x^n
UP(i, 0, n) {
CC *cur = G_rec.data() + (i << (dep + 2));
UP(j, 0, l) {
cur[j] = 0;
}
UP(j, l, r) {
std::swap(cur[j - l], cur[j]);
}
UP(j, r, 1 + (3 << dep)) {
cur[j] = 0;
}
}
return G_rec;
}
// f(g) mod x^maxlen
void polycompose(CC *f_, CC *g_, u32 len_f, u32 len_g, std::vector<CC> &ans, u32 n) {
std::vector<CC> f, g;
f.resize(n, 0);
g.resize(2 * n, 0);
UP(i, 0u, len_f) f[i] = f_[i];
std::reverse(f.begin(), f.end()); // y^n f(1/y)
g[0] = 1;
UP(i, 0u, len_g) {
g[2 * i + 1] = CC(0) - g_[i];
}
std::vector<CC> tmp = kinoshita_li(n - 1, n, f, g, g_[0], 0, n);
ans.resize(n);
UP(i, 0u, n) {
ans[i] = tmp[i * 4];
}
}
} // {}}}
int main() {
int in, im;
cin >> in >> im;
in++;
im++;
std::vector<CC> f, g;
f.reserve(in);
g.reserve(im);
UP(i, 0, in) {
int x;
cin >> x;
f.push_back(x);
}
UP(i, 0, im) {
int x;
cin >> x;
g.push_back(x);
}
std::vector<CC> ans;
Poly::polycompose(f.data(), g.data(), in, im, ans, in);
UP(i, 0, in) {
cout << ans[i].val << ' ';
}
cout << '\n';
}
```