数列幂级数与快速数列变换(Fast Array Transform)

前言：方括号表示艾弗森括号，括号中的命题为真则值为 \(1\)，否则是 \(0\)。（例如 \([1\leq 2]=1, [2<2]=0\)）

数列幂级数

考虑有 \(n\) 个集合 \(S_1, S_2, ......, S_n\)，设 \(U=S_1\times S_2\times ......\times S_n\)（乘法是笛卡尔积）。

\(U\) 中的每一个元素都是一个 \(n\) 元组，或者也可以理解为一个长度为 \(n\) 的数列 \(\{a_1, a_2, ..., a_n\}\)，其中 \(a_i\in S_i\)。

定义 \(f(x)=\sum_{a\in U}f_ax^a\) 是数列幂级数。

定义数列幂级数的加法 \(f(x)+g(x)=\sum_{a\in U}(f_a+g_a)x^a\)。

为了定义数列幂级数的乘法，首先必须有数列加函数 \(\mu_1(x, y), \mu_2(x, y), ..., \mu_n(x, y)\)，定义数列加法 \((a+b)_i=\mu_i(a_i, b_i)\)。

由此定义数列幂级数乘法 \(f(x)g(x)=\sum_{a\in U}(\sum_{b+c=a}f_bg_c)x^a\)。

不难发现，集合幂级数可以看作特殊的数列幂级数，满足 \(S_i=\{0, 1\}\) 表示这个 \(i\) 是否在集合中。

FAT（快速数列变换）

效仿 FWT 的做法，考虑映射 \(\text{FAT}\) 和其逆映射 \(\text{FAT}^{-1}\)，满足

\[\text{FAT}(f(x))\cdot \text{FAT}(g(x))=\text{FAT}(f(x)\times g(x)) \]

即将乘法变为逐点相乘。

为此需要构造特征函数 \(\phi(a, b)\) 满足同态性 \(\phi(a, c)\phi(b, c)=\phi(a+b, c)\)，如果能找到这样的函数，则可以定义 $$\text{FAT}(F(x))=\sum_{a\in U}(\sum_{b\in U}\phi(b, a)f_b)x^a$$
显然这样的 \(\text{FAT}\) 满足条件。

为了构造特征函数，考虑首先构造 \(n\) 个子特征函数 \(\phi_i(x, y)\)，满足单位同态性:$$\phi_i(x, z)\phi_i(y, z)=\phi(\mu_i(x, y), z)$$

设 \(\phi(a, b)=\prod_{i=1}^n\phi_i(a_i, b_i)\)。容易证明如果子特征函数都满足单位同态性，则 \(\phi(a, b)\) 满足同态性。

计算方法

考虑逐位计算，设\(f_0\) 是初始的数列幂级数, \(f_i\) 是计算完前 \(i\) 位后的幂级数，满足 $$[x^a]f_i(x)=\underset{b\in U, b_{i+1}=a_{i+1}, b_{i+2}=a_{i+2}......, b_n=a_n}{\sum}(\prod_{j=1}^i\phi_j(b_j, a_j))[x^b]f_0(x)$$

显然 \(f_n=\text{FAT}(f_0)\)。

考虑利用 \(f_i\) 计算 \(f_{i+1}\) 的方法，注意到

\[[x^a]f_{i+1}(x)=\underset{b\in U, b_{i+2}=a_{i+2}, b_{i+3}=a_{i+3}......, b_n=a_n}{\sum}(\prod_{j=1}^{i+1}\phi_j(b_j, a_j))[x^b]f_0(x) \]

\[=\sum_{p\in S_{i+1}}\phi_{i+1}(p, a_{i+1})\underset{b\in U, b_{i+1}=p, b_{i+2}=a_{i+2}, b_{i+3}=a_{i+3}......, b_n=a_n}{\sum}(\prod_{j=1}^{i}\phi_j(b_j, a_j))[x^b]f_0(x) \]

（这一步枚举 \(b\) 的第 \(i+1\) 项 \(p\)，并将 \(\phi_{i+1}\) 提前）

\[=\sum_{p\in S_{i+1}}\phi_{i+1}(p, a_{i+1})([x^{a'}]f_{i}(x)) \]

（\(a'\) 是将数列 \(a\) 的第 \(i+1\) 项替换为 \(p\) 后得到的数列）

从 \(f_i\) 计算 \(f_{i+1}\) 的复杂度是 \(O(|U||S_{i+1}|)\)，所以总复杂度 \(O(|U|\sum_{i=1}^n|S_i|)=O((\prod_{i=1}^n|S_i|)(\sum_{i=1}^n|S_i|))\)

IFAT（快速数列逆变换）

IFAT 是 FAT 的逆映射，为了计算 \(\text{IFAT}(f(x))\)，考虑反转上述过程，从 \(f_n(x)\) 逐步求到 \(f_0(x)\)，考虑如何从 \(f_{i+1}(x)\) 求出 \(f_{i}(x)\)。

根据公式 \([x^a]f_{i+1}(x)=\sum_{p\in S_{i+1}}\phi_{i+1}(p, a_{i+1})([x^{a'}]f_{i}(x))\)，这实际上是一个向量乘以矩阵的形式，因此只需求出该矩阵的逆矩阵，然后执行矩阵乘法即可。

单次推导的复杂度 \(O(|U||S_{i+1}|+|S_{i+1}|^3)\)，总复杂度

\[O((\prod_{i=1}^n|S_i|)(\sum_{i=1}^n|S_i|)+\sum_{i=1}^n|S_i|^3) \]

FAT 的一些常见特化

FWT

考虑 \(S_i=\{0, 1\}\)，\(\mu_i(x, y)=\max(x, y), \phi_i(x, y)=[x\leq y]\)，此时的数列幂级数乘法相当于集合并卷积。
\(\mu_i(x, y)=\min(x, y), \phi_i(x, y)=[x\geq y]\)，此时是集合交卷积。
\(\mu_i(x, y)=(x+y)\mod 2, \phi_i(x, y)=(-1)^{\min(x, y)}\)，此时是集合对称差（异或）卷积。
并/交卷积的逆变换很简单，对称差卷积的逆变换如何求？当然可以套用 IFAT，但更好的办法是把 \(2\times 2\) 的矩阵写下来，手动求逆。

上述三种算法的复杂度均为

\[O(2^n\times 2n+n\times 2^3)=O(n2^n) \]

注：容易扩展到 k-FWT 的与/或形式。

二进制混合卷积

给出长为 \(\log_2 V\) 的序列，每一项都是 |，&，^ 中的一个。根据这个运算，定义两个小于 \(V\) 的自然数的混合积 \(a\oplus b\) ，表示对于二进制下第 \(i\) 位使用序列第 \(i\) 项对应的运算。
定义两个形式幂级数的卷积为 \([x^k](f(x)g(x))=\sum_{a\oplus b=k}f_ag_b\)。现在给定 \(f(x), g(x)\)，求 \(f(x)g(x)\)。（显然，如果序列中所有的运算都相同，那么混合卷积退化为普通的集合并/交/异或卷积）
做法很简单，令 \(S_i=\{0, 1\}\)，\(\mu_i\) 和 \(\phi_i\) 根据第 \(i\) 位的运算从三种运算中选择，然后直接套用 FAT 的模板即可。

代码

const ll MOD=998244353, N=(1LL<<13)+5, MOD2=499122177;//499122177是0.5模998244353
ll n, m;
string s;
ll a[N], b[N];
void FAT(ll a[], ll cnt){
    rep(i, 0, cnt-1){
        if(s[i]=='|'){
            rep(x, 0, (1LL<<cnt)-1){
                if((x>>i)&1)(a[x]+=a[x^(1<<i)])%=MOD;
            }
        }else if(s[i]=='&'){
            rep(x, 0, (1LL<<cnt)-1){
                if(((x>>i)&1)==0)(a[x]+=a[x^(1<<i)])%=MOD;
            }
        }else{
            rep(x, 0, (1LL<<cnt)-1){
                (b[x]=a[x]*(((x>>i)&1)?-1:1)+a[x^(1<<i)])%=MOD;
            }
            rep(x, 0, (1LL<<cnt)-1)a[x]=b[x];
        }
    }
}
void IFAT(ll a[], ll cnt){
    rrep(i, cnt-1, 0){
        if(s[i]=='|'){
            rep(x, 0, (1LL<<cnt)-1){
                if((x>>i)&1)(a[x]-=a[x^(1<<i)])%=MOD;
            }
        }else if(s[i]=='&'){
            rep(x, 0, (1LL<<cnt)-1){
                if(((x>>i)&1)==0)(a[x]-=a[x^(1<<i)])%=MOD;
            }
        }else{
            rep(x, 0, (1LL<<cnt)-1){
                if((x>>i)&1)b[x]=-a[x]+a[x^(1<<i)], (b[x]*=MOD2)%=MOD;
                else b[x]=a[x]+a[x^(1<<i)], (b[x]*=MOD2)%=MOD;
            }
            rep(x, 0, (1LL<<cnt)-1)a[x]=b[x];
        }
    }
}

例题 [THUPC 2019] 找树
考虑对于 \(v\in[0, 2^w-1]\)，算出权值为 \(v\) 的生成树个数。这个问题显然严格强于原问题。
解决这个问题只需要令权值为 \(v\) 的边的新权值为 \(x^v\)，然后用新权值直接跑矩阵树定理即可，算出的行列式就是一个形式幂级数，\(x^v\) 的系数就是权值为 \(v\) 的生成树个数。
复杂度 \(O(2^wn^3)\)。

子集卷积

直接套用 FAT 可以得到如下的 \(O(n3^n)\) 做法。

\[S_i=\{0, 1, \text{DNE}\} \]

\[\mu_i(x, y)= \begin{cases} x+y, x, y\ne \text{DNE}, x+y \leq 1\\ \text{DNE}, x+y=2\\ \text{DNE}, x=\text{DNE}\text{或}y=\text{DNE} \end{cases}\]

此时计算的数列幂级数乘法相当于子集卷积。这里 \(0\) 和 \(1\) 的含义不变, \(\text{DNE}\) 表示不存在。这样如果两个集合有交（存在某一位，在两边都是 \(1\)），那么就会产生一个包含了 \(\text{DNE}\) 的数列，这样的数列自然不会被我们统计到。

所以更好的特化方法是在 FWT 这样的特化基础上加上占位多项式，从而达到 \(O(n\log n 2^n)\)。

快速 GCD/LCM 变换

例题:有 \(n\) 个集合 \(A_1, A_2, ......, A_n\subset \mathbb{Z}^+\)，和 \(n\) 个数 \(a_1, a_2, ..., a_n\)。\(q\) 次询问，每次询问从 \(A_i\) 中选择 \(a_i\) 个（可以相同）数，总计 \(\sum_{i=1}^na_i\) 个数的 \(\text{LCM}\) 为一个给定的值 \(k\) 的方案数。
\(n\leq 40, \forall x\in A_i, x\leq 10^6, q\leq 10^6, k\leq 10^6\)。

解法：看到计数题首先想幂级数。如果定义 \(f_i(x)=\sum_{a\in A_i}x^a\)，并且 \(x^a\times x^b=x^{\text{lcm}(a, b)}\)，那么只需要计算 \(\prod_{i=1}^nf_i(x)\) 即可。
这显然是数列幂级数的特化，只需 \(n=1, S_1=\{1, 2, ..., 10^6\}, \mu_1(x, y)=\text{lcm}(x, y)\) 即可。
容易发现 \(\phi_1(x, y)=[x|y]\) 能够保证同态性，所以 \(\text{FAT}\) 相当于狄利克雷前缀和，而 \(\text{IFAT}\) 就是狄利克雷差分。

注：可以扩展到 \(n>1\) 的情况，此时是高维狄利克雷前缀和/差分。