集合幂级数“全”家桶（转载自洛谷）
集合幂级数全家桶（？）
1. 集合幂级数
对于一个全集 \(U\)，给 \(U\) 的每个子集 \(S\) 赋予一个权值 \(f_S\)。我们也可以通过类似“多项式”的方法描述：
由此可以进行一些运算。
2. 基本操作
- 子集卷积
- 集合幂级数求逆
- 集合幂级数 \(\ln\)
- 集合幂级数 \(\exp\)
3. FWT, FMT
3.1 基本操作
集合幂级数的基本运算是：
其中 \(\odot\) 是一种常见的集合运算符，包括 \(\cup\)（并），\(\cap\)（交），\(\Delta\)（差），\(\sqcup\)（不交并）。
前面三种是基础的，最后一种基于前面的算法。
显然，我们可以 \(O(4^n)\) 来枚举 \(R, T\) 来得到 \(H\)，但这太慢了。这就需要引入算法：快速沃尔什变换（FWT）
3.2 恰当的中间值
在 FFT 中，我们先对 \(F, G\) 进行变换 \(\hat{F}, \hat{G}\)，然后逐位相乘得到 \(\hat{H}\)，再逆变换得到 \(H\)。在集合幂级数的运算中，我们同样有类似的中间值。此时，对于不同的运算 \(\cup, \cap, \Delta\)，我们有：
其中 \(T \circ S\) 表示 \(|T \cap S| \bmod 2\)。我们一一验证一下它的正确性。对于 \(\cup\)，
对于 \(\cap\)，我们每一个 \(S\) 取补集就和 \(\cup\) 一样了。对于 \(\Delta\)，证明稍微复杂：
3.3 高维前缀和
由于在这些运算中，位之间都是相对独立的。我们可以逐位处理。
对于 \(\cup\)，我们在处理一个位 \(i\) 时，我们执行 \(f_S \overset{+}{\to} f_{S \cup \{i\}}\)，这样执行完后，对于 \(S\) 的每个子集 \(T\)，其恰好有一种贡献路径达到 \(S\)。至于怎么处理，可以类似 FFT 中的分治处理：
void fwt_or(ll *f, int len) // len 是集合个数
{
for (int h = 2; h <= len; h <<= 1) // 枚举处理的位 * 2
{
for (int j = 0; j < len; j += h)
{
for (int i = 0; i < (h >> 1); i++)
{
(f[j + i + (h >> 1)] += f[j + i]) %= MOD;
}
}
}
}
我们也将处理 \(\cup\) 操作的 FWT 称为快速莫比乌斯变换（FMT）。
将其反过来处理（交换第 \(9\) 行），就可以得到 \(\cap\) 的 FWT 结果了。
而对于 \(\Delta\) 操作，我们只要只需进行 \((f'_S, f'_{S + \{i\}}) \gets (f'_{S} + f'_{S + \{i\}}, f'_{S} - f'_{S + \{i\}})\) 即可。这是因为，对没有 \(i\) 的 \(S\)，添加一个 \(i\) 之后，两集合的交大小奇偶性不会变化，而对于有 \(i\) 的 \(S\) 会变，所以就可以得到这个式子了。
3.4 逆变换
这其实相当的简单。例如 \(\cup\) 操作，我们在 FWT 时进行了 \(f_S \overset{+}{\to} f_{S \cup \{i\}}\)，IFWT 时进行 \(f_{S \cup \{i\}} \gets f_{S \cup \{i\}} - f_S\) 即可。
其他操作也类似。在 \(\Delta\) 操作的 IFWT 中，我们还要除以 \(2\)（类似于三角函数的积化和差）。可以发现，这些过程和 FWT 过程非常类似，我们可以放在一个函数内写。最终可以得到：
void fwt_or(ll *f, int len, bool on)
{
ll type = on ? 1 : -1;
for (int h = 2; h <= len; h <<= 1)
{
for (int j = 0; j < len; j += h)
{
for (int i = 0; i < (h >> 1); i++)
{
(f[j + i + (h >> 1)] += MOD + f[j + i] * type) %= MOD;
}
}
}
}
void fwt_and(ll *f, int len, bool on)
{
ll type = on ? 1 : -1;
for (int h = 2; h <= len; h <<= 1)
{
for (int j = 0; j < len; j += h)
{
for (int i = 0; i < (h >> 1); i++)
{
(f[j + i] += MOD + f[j + i + (h >> 1)] * type) %= MOD;
}
}
}
}
void fwt_xor(ll *f, int len, bool on)
{
ll type = on ? 1 : I2;
for (int h = 2; h <= len; h <<= 1)
{
for (int j = 0; j < len; j += h)
{
for (int i = 0; i < (h >> 1); i++)
{
ll u = f[j + i], t = f[j + i + (h >> 1)];
f[j + i] = (u + t) * type % MOD;
f[j + i + (h >> 1)] = (u - t + MOD) * type % MOD;
}
}
}
}
4. 子集卷积
我们还需要应对不交并 \(\sqcup\) 操作。
4.1 多项式套幂级数？
对于集合 \(S, T\)，我们观察它们的大小 \(|S|, |T|\)。当进行不交并时，若有 \(|S \cup T| = |S| + |T|\)，我们才会真正乘一起。而当我们真正进行这个过程时，相当于将大小为 \(|S|\) 的集合与大小为 \(|T|\) 的集合权值乘一起，贡献到 \(|S| + |T|\) 的集合去了。能够发现，这和多项式里的 \(x^{|S|} \times x^{|T|} = x^{|S| + |T|}\) 很像。
因此，我们可以定义一个新的（不知道用什么词）式来表达我们初始的 \(f_S\)：
当我们给 \(F(x, y)\) 与 \(G(x, y)\) 相乘时，所有的 \(f_S\) 与 \(g_T\) 会被乘进 \([x^{S \cup T}y^{|S| + |T|}](F \times G)(x, y)\)。而我们最终只取 \([x^Sy^{|S|}](F \times G)(x, y)\)，也就保证了当 \(|S| + |T| = |S \cup T|\) 时，乘法的结果才会被我们取到。所以我们得到了正确的不交并卷积，也就是子集卷积。
4.2 多套幂的操作
我们可以将 \(F\) 存进一个二维数组 ll f[MAXL][MAXN]，其中第一维存 \(y\)，第二维存 \(x\)。
考虑到乘法的实质，我们要做的是先对 \(x\) 这一维进行 FWT，然后对 \(y\) 这一维进行 FFT，然后再分别位置相乘。这个可以达到理论复杂度 \(O(2^nn\log n)\)。
然而在实践中，由于 \(y\) 的次数较小，我们不会对它进行 FFT 优化。也就是用 \(O(n^2)\) 的多项式乘法，达到 \(O(2^nn^2)\)。
/**
* calculate h(x) = f(x) * g(x), where x^S * x^T = { 0 if S ∩ T != ∅
* { x^{S ∪ T} else
* please offer f, g, h in ll[log len][len], where [x^S]f(x) is stored in f[|S|][S].
*/
template <typename ArrayType>
void subset_conv(__restrict__ ArrayType f, __restrict__ ArrayType g, int loglen, __restrict__ ArrayType h)
{
int len = 1 << loglen;
for (int i = 0; i <= loglen; i++)
{
fwt_or(f[i], len, true);
fwt_or(g[i], len, true);
}
for (int i = 0; i <= loglen; i++) // 这一段实际上就是多项式乘法
{
for (int j = 0; j < len; j++)
{
h[i][j] = 0;
for (int k = 0; k <= i; k++)
{
(h[i][j] += f[k][j] * g[i - k][j] % MOD) %= MOD;
}
}
}
for (int i = 0; i <= loglen; i++)
{
fwt_or(f[i], len, false);
fwt_or(g[i], len, false);
fwt_or(h[i], len, false);
}
}
5. 初等函数
5.1 多项式的初等函数搬进 \(\sqcup\)
在实际情况中，我们遇到的基本都是 \(\sqcup\) 操作。既然我们定义了 \(\sqcup\) 上的乘法（也就是子集卷积），我们同样可以衍生出别的函数，例如 \(\exp, \ln\) 和求逆。
在这些函数中，可以发现它们在 \(y\) 这一维度上都和多项式一样。这可能很违反直觉，也可能难以理解，这是因为它和我们常见的操作差异巨大，且比较难找到合适的组合意义。
这类操作都可以封装在一个模板中：
/**
* calculate f(x) = oper(g(x)), where x^S * x^T = { 0 if S ∩ T != ∅
* { x^{S ∪ T} else
* please offer f, g in ll[log len][len], where [x^S]g(x) is stored in g[|S|][S].
*
* @param oper The operation on sets.
* Its type should be void (*) (ll *f, ll *g, int loglen) where f is the destination
*/
template <typename ArrayType, typename Oper> // 理论上 ArrayType 应当是类似 (ll **) 或 (ll *[])
void subset_oper(__restrict__ ArrayType f, __restrict__ ArrayType g, int loglen, const Oper oper)
{
int len = 1 << loglen;
for (int i = 0; i <= loglen; i++)
{
fwt_or(g[i], len, true);
}
static ll F[MAXL], G[MAXL];
for (int j = 0; j < len; j++)
{
for (int i = 0; i <= loglen; i++)
{
G[i] = g[i][j];
}
oper(F, G, loglen);
for (int i = 0; i <= loglen; i++)
{
f[i][j] = F[i];
}
}
ll inl = pow(len, MOD - 2, MOD);
for (int i = 0; i <= loglen; i++)
{
fwt_or(g[i], len, false);
fwt_or(f[i], len, false);
}
}
5.2 \(\exp\) 的平方解
\(\exp\) 算是最好理解的操作了。\(F(x) = \exp(G(x))\) 中，\(f_S\) 相当于将 \(S\) 划分成任意多个子集 \(S_1 \sqcup S_1 \sqcup \dots \sqcup S_k = S\)，然后对每个 \(S_i\) 取 \(g_{S_i}\) 然后乘起来得到结果。这个过程在 \(|S|\) 意义下就是普通的多项式 \(\exp\)。我们要得到一个好写的 \(O(n^2)\) 求 \(\exp\) 的方法（下文称计算 \(\exp (f)\)）。
事实上，我们已经在求多项式 \(\exp\) 时得到了这么一个式子（见多项式全家桶初等函数段）：
这个式子显然是可以 \(O(n^2)\) 转移的。
5.3 \(\ln\) 操作
首先肯定可以通过 \(\int \frac{f'(x)}{f(x)}\) 得到，但过于复杂。考虑到 \(\ln\) 是 \(\exp\) 的反函数，即：
5.4 求逆
这相当于 \(\sum f_if^{-1}_{n - i} = [n = 0]\)，直接移项就可以得到。我们可以写出以上三者的代码：
void exp_bf(ll *f, ll *g, int loglen) // 建议不要和我自己写的 POLY::poly_exp 混用，因为 POLY::poly_exp 只支持 2^n 长度
{
const ll *inv = get_inv();
f[0] = 1;
for (int i = 1; i <= loglen; i++)
{
f[i] = 0;
for (int j = 1; j <= i; j++)
{
(f[i] += f[i - j] * g[j] % MOD * j) %= MOD;
}
(f[i] *= inv[i]) %= MOD;
}
}
void ln_bf(ll *f, ll *g, int loglen)
{
const ll *inv = get_inv();
f[0] = 0;
for (int i = 1; i <= loglen; i++)
{
f[i] = 0;
for (int j = 1; j < i; j++)
{
(f[i] += f[j] * g[i - j] % MOD * j) %= MOD;
}
f[i] = (g[i] - inv[i] * f[i] % MOD + MOD) % MOD;
}
}
void inv_bf(ll *f, ll *g, int loglen)
{
f[0] = pow(g[0], MOD - 2, MOD);
for (int i = 1; i <= loglen; i++)
{
f[i] = 0;
for (int j = 1; j <= i; j++)
{
(f[i] += g[j] * f[i - j]) %= MOD;
}
f[i] = (MOD - f[0] * f[i] % MOD) % MOD;
}
}
6. 完整代码
#define MAXL 22
#define MAXN 1100005
#define MOD 998244353ll
#define I2 499122177ll
#define popcnt __builtin_popcount
ll pow(ll b, ll p, ll m)
{
ll r = 1;
while (p)
{
if (p & 1)
{
r = r * b % m;
}
b = b * b % m;
p >>= 1;
}
return r;
}
namespace SSPOLY // subset poly
{
void fwt_or(ll *f, int len, bool on)
{
ll type = on ? 1 : -1;
for (int h = 2; h <= len; h <<= 1)
{
for (int j = 0; j < len; j += h)
{
for (int i = 0; i < (h >> 1); i++)
{
(f[j + i + (h >> 1)] += MOD + f[j + i] * type) %= MOD;
}
}
}
}
void fwt_and(ll *f, int len, bool on)
{
ll type = on ? 1 : -1;
for (int h = 2; h <= len; h <<= 1)
{
for (int j = 0; j < len; j += h)
{
for (int i = 0; i < (h >> 1); i++)
{
(f[j + i] += MOD + f[j + i + (h >> 1)] * type) %= MOD;
}
}
}
}
void fwt_xor(ll *f, int len, bool on)
{
ll type = on ? 1 : I2;
for (int h = 2; h <= len; h <<= 1)
{
for (int j = 0; j < len; j += h)
{
for (int i = 0; i < (h >> 1); i++)
{
ll u = f[j + i], t = f[j + i + (h >> 1)];
f[j + i] = (u + t) * type % MOD;
f[j + i + (h >> 1)] = (u - t + MOD) * type % MOD;
}
}
}
}
/**
* @param fwt The FWT function.
* Its type should be void (*) (ll *f, int len, bool on).
* Provided functions are fwt_or, fwt_and, fwt_xor.
*/
template <typename FWT_Type>
void conv(ll __restrict__ *f, ll __restrict__ *g, int len, ll __restrict__ *h, const FWT_Type fwt)
{
fwt(f, len, true);
fwt(g, len, true);
for (int i = 0; i < len; i++)
{
h[i] = f[i] * g[i] % MOD;
}
fwt(f, len, false);
fwt(g, len, false);
fwt(h, len, false);
}
/**
* calculate h(x) = f(x) * g(x), where x^S * x^T = x^{S ∪ T}
*/
void conv_or(ll __restrict__ *f, ll __restrict__ *g, int len, __restrict__ ll *h)
{
conv(f, g, len, h, fwt_or);
}
/**
* calculate h(x) = f(x) * g(x), where x^S * x^T = x^{S ∩ T}
*/
void conv_and(ll __restrict__ *f, ll __restrict__ *g, int len, __restrict__ ll *h)
{
conv(f, g, len, h, fwt_and);
}
/**
* calculate h(x) = f(x) * g(x), where x^S * x^T = x^{S Δ T}
*/
void conv_xor(ll __restrict__ *f, ll __restrict__ *g, int len, __restrict__ ll *h)
{
conv(f, g, len, h, fwt_xor);
}
/**
* calculate h(x) = f(x) * g(x), where x^S * x^T = { 0 if S ∩ T != ∅
* { x^{S ∪ T} else
* please offer f, g, h in ll[log len][len], where [x^S]f(x) is stored in f[|S|][S].
*/
template <typename ArrayType>
void subset_conv(__restrict__ ArrayType f, __restrict__ ArrayType g, int loglen, __restrict__ ArrayType h)
{
int len = 1 << loglen;
for (int i = 0; i <= loglen; i++)
{
fwt_or(f[i], len, true);
fwt_or(g[i], len, true);
}
for (int i = 0; i <= loglen; i++)
{
for (int j = 0; j < len; j++)
{
h[i][j] = 0;
for (int k = 0; k <= i; k++)
{
(h[i][j] += f[k][j] * g[i - k][j] % MOD) %= MOD;
}
}
}
for (int i = 0; i <= loglen; i++)
{
fwt_or(f[i], len, false);
fwt_or(g[i], len, false);
fwt_or(h[i], len, false);
}
}
/**
* only inverse element under log len.
*/
const ll *get_inv()
{
static ll inv[MAXL];
if (!inv[1])
{
for (int i = 1; i < MAXL; i++)
{
inv[i] = pow(i, MOD - 2, MOD);
}
}
return inv;
}
void exp_bf(ll *f, ll *g, int loglen)
{
const ll *inv = get_inv();
f[0] = 1;
for (int i = 1; i <= loglen; i++)
{
f[i] = 0;
for (int j = 1; j <= i; j++)
{
(f[i] += f[i - j] * g[j] % MOD * j) %= MOD;
}
(f[i] *= inv[i]) %= MOD;
}
}
void ln_bf(ll *f, ll *g, int loglen)
{
const ll *inv = get_inv();
f[0] = 0;
for (int i = 1; i <= loglen; i++)
{
f[i] = 0;
for (int j = 1; j < i; j++)
{
(f[i] += f[j] * g[i - j] % MOD * j) %= MOD;
}
f[i] = (g[i] - inv[i] * f[i] % MOD + MOD) % MOD;
}
}
void inv_bf(ll *f, ll *g, int loglen)
{
f[0] = pow(g[0], MOD - 2, MOD);
for (int i = 1; i <= loglen; i++)
{
f[i] = 0;
for (int j = 1; j <= i; j++)
{
(f[i] += g[j] * f[i - j]) %= MOD;
}
f[i] = (MOD - f[0] * f[i] % MOD) % MOD;
}
}
/**
* calculate f(x) = oper(g(x)), where x^S * x^T = { 0 if S ∩ T != ∅
* { x^{S ∪ T} else
* please offer f, g in ll[log len][len], where [x^S]g(x) is stored in g[|S|][S].
*
* @param oper The operation on sets.
* Its type should be void (*) (ll *f, ll *g, int loglen) where f is the destination
* Provided functions are exp_bf, ln_bf, inv_bf
*/
template <typename ArrayType, typename Oper>
void subset_oper(__restrict__ ArrayType f, __restrict__ ArrayType g, int loglen, const Oper oper)
{
int len = 1 << loglen;
for (int i = 0; i <= loglen; i++)
{
fwt_or(g[i], len, true);
}
static ll F[MAXL], G[MAXL];
for (int j = 0; j < len; j++)
{
for (int i = 0; i <= loglen; i++)
{
G[i] = g[i][j];
}
oper(F, G, loglen);
for (int i = 0; i <= loglen; i++)
{
f[i][j] = F[i];
}
}
ll inl = pow(len, MOD - 2, MOD);
for (int i = 0; i <= loglen; i++)
{
fwt_or(g[i], len, false);
fwt_or(f[i], len, false);
}
}
/**
* calculate f(x) = exp(g(x))
* see subset_oper for more detail
*/
template <typename ArrayType>
void subset_exp(__restrict__ ArrayType f, __restrict__ ArrayType g, int loglen)
{
subset_oper(f, g, loglen, exp_bf);
}
/**
* calculate f(x) = ln(g(x))
* see subset_oper for more detail
*/
template <typename ArrayType>
void subset_ln(__restrict__ ArrayType f, __restrict__ ArrayType g, int loglen)
{
subset_oper(f, g, loglen, ln_bf);
}
/**
* calculate f(x) = 1 / g(x)
* see subset_oper for more detail
*/
template <typename ArrayType>
void subset_inv(__restrict__ ArrayType f, __restrict__ ArrayType g, int loglen)
{
subset_oper(f, g, loglen, inv_bf);
}
}
7. 后记
我也不知道我为什么要学这个，也不知道我能不能记住，实战时更想不到。反正 CCF 不会考，但是 ATC 还是会有的。
实际上多项式全家桶和集合幂级数全家桶还没写完。
已知没写的：多项式复合、复合逆、多项式复合集合幂级数、……
浙公网安备 33010602011771号