浅谈FWT的一些进阶模板

· · 算法·理论

学过 FWT 的大佬，也建议看一看本蒟蒻写的 FWT 介绍，后文可能会用其中的一些符号和定义。

卷积形式

我们有两个序列 a 和 b，各有 n 个元素。

求序列 c，满足 c_k=\sum\limits_{i\oplus j=k}a_ib_j。

万能卷积

为了方便，设 n 为 2 的整数次幂。

延续 FFT 的思路，我们把多项式转成点值表达。

我们做一次线性变换，并设出变换系数。

即 A_i=\sum\limits_{j=0}^{n-1}c(i,j)a_j。

那么若 C_i=A_iB_i。

可以推出 c(i,j)c(i,k)=c(i,j\oplus k)。

这是 c 数组的性质一。

能发现此性质只限制了同一个 i 的取值，不同 i 之间的取值关系没有考虑。

为了方便我们定义删最高位操作。

$j$ 的删位性质： $c(i,j)=c(i,j-\text{highbit}(j))c(i,1)$。 # 分治 还是考虑 $A_i=\sum\limits_{j=0}^{n-1}c(i,j)a_j$。 把这 $n$ 个数拆成两半。 $A_i=\sum\limits_{j=0}^{\frac{n}{2}-1}c(i,j)a_j+\sum\limits_{j=\frac{n}{2}}^{n-1}c(i,j)a_j$。 此时考虑将左右两边的 $i$ 和 $j$ 做一次删位操作。 得到： $i$ 最高位为 $0$：$A_i=c(0,0)\sum\limits_{j=0}^{\frac{n}{2}-1}c(i',j')a_j+c(0,1)\sum\limits_{j=\frac{n}{2}}^{n-1}c(i',j')a_j$。 $i$ 最高位为 $1$：$A_i=c(1,0)\sum\limits_{j=0}^{\frac{n}{2}-1}c(i',j')a_j+c(1,1)\sum\limits_{j=\frac{n}{2}}^{n-1}c(i',j')a_j$。 这里可以把上述变换理解为矩阵乘法。 定义 $A=\begin{bmatrix}A_{[0..\frac{n}{2})} \\ A_{[\frac{n}{2}..n)}\end{bmatrix}$。 定义 $B=\begin{bmatrix}c(0,0)&c(0,1)\\c(1,0)&c(1,1)\end{bmatrix}

那么新的 A'=BA。

逆变换

考虑把上述矩阵求逆。

即 B^{-1}A'=A。

那么一个合法 c 就多了一个条件——存在逆矩阵。

合法矩阵

如何构造？

考场上如果忘记矩阵，那么可以考虑解方程。

这个方程存在多个解，于是可以考虑把他们作为矩阵的不同的行。

多尝试几次一般就能成功。

按位或矩阵

1 & 0 \\ 1 & 1 \end{bmatrix}

其逆矩阵：

\begin{bmatrix} 1 & 0 \\ -1 & 1 \end{bmatrix}

按位与矩阵

\begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix}

其逆矩阵：

\begin{bmatrix} 1 & -1 \\ 0 & 1 \end{bmatrix}

按位异或矩阵

\begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix}

其逆矩阵：

\begin{bmatrix} 0.5 & 0.5 \\ 0.5 & -0.5 \end{bmatrix}

模版一：普通 FWT

P4717 【模板】快速莫比乌斯 / 沃尔什变换 (FMT / FWT)

使用 for 而不是递归。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=998244353;
const int N=1e6+10;
int n,len;
ll a[N],b[N],f[N],g[N];
void OR(ll *f,int len,ll v)
{
    for(int h=2;h<=len;h<<=1)
    {
        int mid=h/2;
        for(int i=0;i<len;i+=h)
        {
            for(int j=i;j<i+mid;j++)
            {
                f[j+mid]=(f[j+mid]+f[j]*v)%mod;
            }
        }
    }
}
void AND(ll *f,int len,ll v)
{
    for(int h=2;h<=len;h<<=1)
    {
        int mid=h/2;
        for(int i=0;i<len;i+=h)
        {
            for(int j=i;j<i+mid;j++)
            {
                f[j]=(f[j]+f[j+mid]*v)%mod;
            }
        }
    }
}
void XOR(ll *f,int len,ll v)
{
    for(int h=2;h<=len;h<<=1)
    {
        int mid=h/2;
        for(int i=0;i<len;i+=h)
        {
            for(int j=i;j<i+mid;j++)
            {
                ll a=f[j],b=f[j+mid];
                f[j]=(a+b)*v%mod;f[j+mid]=(a-b)*v%mod;
            }
        }
    }
}
int main()
{
    scanf("%d",&n);len=(1<<n);
    for(int i=0;i<len;i++) scanf("%lld",&a[i]);
    for(int i=0;i<len;i++) scanf("%lld",&b[i]);
    for(int i=0;i<len;i++) f[i]=a[i],g[i]=b[i];
    OR(f,len,1);OR(g,len,1);
    for(int i=0;i<len;i++) f[i]=f[i]*g[i]%mod;
    OR(f,len,-1);
    for(int i=0;i<len;i++) printf("%lld ",(f[i]%mod+mod)%mod);printf("\n");
    for(int i=0;i<len;i++) f[i]=a[i],g[i]=b[i];
    AND(f,len,1);AND(g,len,1);
    for(int i=0;i<len;i++) f[i]=f[i]*g[i]%mod;
    AND(f,len,-1);
    for(int i=0;i<len;i++) printf("%lld ",(f[i]%mod+mod)%mod);printf("\n");
    for(int i=0;i<len;i++) f[i]=a[i],g[i]=b[i];
    XOR(f,len,1);XOR(g,len,1);
    for(int i=0;i<len;i++) f[i]=f[i]*g[i]%mod;
    XOR(f,len,(mod+1)/2);
    for(int i=0;i<len;i++) printf("%lld ",(f[i]%mod+mod)%mod);
    return 0;
}

模板二：子集卷积

问题描述：

有两个长度为 n 的序列 A,B，求 C_k=\sum\limits_{i\&j=0\land i| j=k}A_iB_j。

如果没有 i\& j=k 就是一个普通的 FWT。

能发现若 i\& j=k，那么 |i|+|j|=|i\cup j|。

这提示着我们把集合大小加入状态。

即定义 f_{i,j}=\begin{cases}A_j\,(|j|=i)\\0\end{cases},g_{i,j}=\begin{cases}B_j\,(|j|=i)\\0\end{cases}。

有 h_i=\sum_{k=0}^i f_k * g_{i-k}。

那么答案 C_i=h_{|i|,i}。

为了减少时间复杂度，我们一开始把所有序列 FWT，那么中间的所有卷积都变成直积了。

时间复杂度：O(n^22^n)。

参考代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=1e9+9;
const int N=(1<<20)+10;
void FWT(ll *a,int len,ll v)
{
    for(int h=2;h<=len;h<<=1)
    {
        int mid=h/2;
        for(int i=0;i<len;i+=h)
        {
            for(int j=i;j<i+mid;j++)
            {
                ll x=a[j],y=a[j+mid];
                a[j]=(x+y)*v%mod;a[j+mid]=(x-y+mod)*v%mod;
            }
        }
    }
}
ll f[22][N],g[22][N],h[22][N];
int n,len=0,c[N];
int main()
{
    scanf("%d",&n);len=(1<<n);
    for(int i=0;i<len;i++) c[i]=c[i>>1]+(i&1);
    for(int i=0;i<len;i++)
    {
        int x;scanf("%d",&x);
        f[c[i]][i]=x;
    }
    for(int i=0;i<len;i++)
    {
        int x;scanf("%d",&x);
        g[c[i]][i]=x;
    }
    for(int i=0;i<=n;i++) FWT(f[i],len,1),FWT(g[i],len,1);
    for(int i=0;i<=n;i++)
    {
        for(int j=0;j<=i;j++)
        {
            for(int k=0;k<len;k++) h[i][k]=(h[i][k]+f[j][k]*g[i-j][k])%mod;
        }
    }
    for(int i=0;i<=n;i++) FWT(h[i],len,(mod+1)/2);
    for(int i=0;i<len;i++) printf("%lld ",h[c[i]][i]);
    return 0;
}

核的引入

定义，对于一个序列 F，若存在一个 x，使得若定义 C_i=F_{i\oplus x}，那么 \forall i,j，满足 \text{highbit}(i)=\text{highbit}(j)，有 C_i=C_j。

:::info[引理一] 一个序列若存在核，则他一定最多有 \log_2 n 种数值。 ::: :::success[证明] 根据核的定义即证。 ::: :::info[引理二] 若 A,B 的核为 0，那么 C=A\oplus B（异或卷积）的核也为 0。 ::: :::success[证明] 根据异或卷积的定义有，C_k=\sum\limits_{i\oplus j=k}A_iB_j，证明 C 的核等于 0，等价于要证明 C 取值与 K 无关并与 \beta=\text{highbit}(k) 有关。

我们来分类讨论 i,j 的取值范围。

  2. i<2^{\beta}$ 或 $j<2^{\beta}$，不妨设 $i<2^{\beta}\leq j<2^{\beta+1}$，此时任意 $i$ 所生成的 $j$ 的值都等于 $B_k$，因为 $\text{highbit}(j)=\beta$，不难发现贡献为：$B_k\sum\limits_{i<2^{\beta}} A_i
因此证毕。
:::info[引理三]
若 A 的核为 x，B 的核为 y，那么 C=A\oplus B（异或卷积）的核为 x\oplus y。
:::
:::success[证明]
构造 A'_{i}=A_{i\oplus x},B'_{j}=B_{j\oplus y}，那么 A'\oplus B'=C'，由引理一可知 C' 的核为 0，又因 C_i=C'_{i\oplus x\oplus y}，所以 C 的核为 x\oplus y。
因此 Q.E.D.
:::
# 模板三：有核序列的异或卷积。

根据上述证明，我们可以得到如何求出两个核为 0 的序列的异或卷集（分类讨论 i,j 并分别计算贡献），根据引理一可知我们只需要存储 \log_2n 的数值即可。当然也要存储核。

若核不为 0 呢，那就在算答案的时候异或上核即可（类比引理三的证明），这样求卷集时就可以无视核的存在，只用把核异或起来。

原题连接： #3073. 「2019 集训队互测 Day 2」序列

参考代码：

struct st
{
    ll val[25],vis[25],sp;
    void init(ll x,ll v)
    {
        for(ll i=19;i>=0;i--)
        {
            if(x>>i&1)
            {
                vis[i]=1;ll tmp=x-(x&((1ll<<i)-1));
                val[i]=(tmp*(tmp^v)%mod);
            }
            else vis[i]=0;
        }
    }
    friend st operator*(st a,st b)
    {
        st c;
        ll sum1=a.sp,sum2=b.sp,s=0;
        for(ll i=0;i<=19;i++) sum1=(sum1+a.val[i]*(1ll<<i))%mod,sum2=(sum2+b.val[i]*(1ll<<i))%mod;
        for(ll i=19;i>=0;i--)
        {
            c.vis[i]=a.vis[i]^b.vis[i];
            sum1=(sum1-a.val[i]*(1ll<<i))%mod;sum2=(sum2-b.val[i]*(1ll<<i))%mod;
            c.val[i]=(s+sum1*b.val[i]+sum2*a.val[i])%mod;
            s+=a.val[i]*b.val[i]%mod*(1ll<<i)%mod;s%=mod;
        }
        c.sp=a.sp*b.sp+s;
        return c;
    }
    ll query(ll d)
    {
        for(int i=19;i>=0;i--) if(((d>>i&1)^vis[i])==1) return (val[i]%mod+mod)%mod;
        return (sp%mod+mod)%mod;
    }
}f[N];
int main()
{
    scanf("%d%d%d",&n,&m,&q);
    for(int i=1;i<=n;i++)
    {
        ll x,v;scanf("%lld%lld",&x,&v);
        f[i].init(x,v);
    }
    for(int i=2;i<=n;i++) f[i]=f[i-1]*f[i];
    while(q--)
    {
        int c,d;scanf("%d%d",&c,&d);
        printf("%lld\n",f[c].query(d));
    }
    return 0;
}

模板四：稀疏序列的异或卷积（系数可重集相同）

以一道题目为例。

:::info[题目描述]

有 n 个非负整数三元组 (a_i, b_i, c_i) 和四个非负整数 x,y,z,k。 你利用这 n 个三元组填充了 n 个数组，其中第 i 个数组中有 x 个 a_i，y 个 b_i，z 个 c_i（所以第 i 个数组长度为 (x+y+z)）。 对于 i=0,1,\dots,2^k−1，回答以下询问：

  • 从每个数组中选择恰好一个数，使得这些数的 \mathrm{xor} 和为 i，方案数是多少？
你只需要输出方案数对 998,244,353 取模后得到的结果。
定义 f_i 表示所有序列 FWT 后直积的结果。

能发现稀疏序列的异或卷积有一个前提：所有序列的 x,y,z 相等。

这提示着我们有些序列之间可以划分到一个类中。

先分析 f_i。

显然：f_i=\prod\limits_{i=1}^n (-1)^{u_i}x+(-1)^{v_i}y+(-1)^{w_i}z。

为了简化计算，我们让每个序列的值都异或上 a_i，这样 f_i=\prod\limits_{i=1}^n x+(-1)^{v_i}y+(-1)^{w_i}z 对于两个 (-1) 的分布，一共有 4 种情况。

下面考虑列方程。 最平凡的方程是：$c_1+c_2+c_3+c_4=n$。 接着我们分析 $y$ 前面的系数和，$z$ 前面的系数和，$y,z$ 前面的系数之积的和。这可以列出三个方程，设这三个和分别是 $p_1,p_2,p_3$。 考虑如何计算 $p_1,p_2,p_3$。 1. 计算 $p_1$：$y$ 前面的系数等价于变换系数 $c_{i,b_i}$，为了一次统计所有的 $c_{i,b_i}$，我们可以定义 $f$ 序列，$f_i$ 表示 $b$ 等于 $i$ 的序列个数。那么 $f$ 做 FWT 后的 $f_i$ 就是 $p_1$。 2. $p_2$ 同理。 3. 计算 $p_3$：能发现我们要计算 $c(i,b_i)c(i,c_i)$ 的和，根据变换系数的性质可得，$c(i,b_i)c(i,c_i)=c(i,b_i\oplus c_i)$，那么我们可以定义 $h$ 序列，$h_i$ 表示 $b\oplus c$ 等于 $i$ 的序列个数。那么 $f$ 做 FWT 后的 $h_i$ 就是 $p_3$。 因此我们可以列下方程： $\begin{cases}n=c_1+c_2+c_3+c_4\\p1=c_1+c_2-c_3-c_4\\p2=c_1-c_2+c_3-c_4\\p3=c_1-c_2-c_3+c_4\end{cases}

解得：

\begin{cases}c_1=(n+p1+p2+p3)/4\\c_2=(n+p1-2c_1)/2\\c_3=(n+p2-2c_1)/2\\c_4=(n+p3-2c_1)/2\end{cases}

原题链接：CF1119H Triple

参考代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e6+10;
const ll mod=998244353;
int n,k,len;
ll x,y,z;
void FWT(ll *a,int len,ll v)
{
    for(int h=2;h<=len;h<<=1)
    {
        int mid=h/2;
        for(int i=0;i<len;i+=h)
        {
            for(int j=i;j<i+mid;j++)
            {
                ll x=a[j],y=a[j+mid];
                a[j]=(x+y)*v%mod;a[j+mid]=(x-y)*v%mod;
            }
        }
    }
}
ll fb[N],fc[N],fbc[N],ta=0,f[N];
ll ku(ll a,ll b)
{
    ll ans=1;
    while(b)
    {
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;b>>=1;
    }
    return ans;
}
int main()
{
    scanf("%d%d",&n,&k);len=(1<<k);
    scanf("%lld%lld%lld",&x,&y,&z);
    for(int i=1;i<=n;i++)
    {
        ll a,b,c;scanf("%lld%lld%lld",&a,&b,&c);
        ta^=a;b^=a;c^=a;
        fb[b]++;fc[c]++;fbc[b^c]++;
    }
    FWT(fb,len,1);FWT(fc,len,1);FWT(fbc,len,1);
    for(int i=0;i<len;i++)
    {
        ll p1=fb[i],p2=fc[i],p3=fbc[i];
        ll c1=(n+p1+p2+p3)/4,c2=(n+p1-2*c1)/2,c3=(n+p2-2*c1)/2,c4=(n+p3-2*c1)/2;
        f[i]=ku(x+y+z,c1)*ku(x+y-z,c2)%mod*ku(x-y+z,c3)%mod*ku(x-y-z,c4)%mod;
    }   
    FWT(f,len,(mod+1)/2);
    for(int i=0;i<len;i++) printf("%lld ",(f[i^ta]%mod+mod)%mod);
    return 0;
}

能发现 x,y,z 可以不止有三个。

假设每个序列有 k 个值，数量分别为 x_1,x_2,\dots,x_k，每个序列这 k 个值分别为 a_{i,1},a_{i,2},\dots,a_{i,k}

f_i=\prod_{i=1}^n\sum\limits_{j=1}^k(-1)^{u_{i,j}}a_{i,j}

定义 c_s 表示 s 集合中的值前面的系数为 1，其余的为 -1 的序列的个数。

定义 p_s 表示 s 集合中的值前面的系数之积的和。

我们要先计算出 p_s（计算方法和三个值的差不多），然后列方程。

最后使用高斯消元即可。 时间复杂度：$O((q+n\log_2 n)8^k)

代码就咕掉了（小编太懒了...）。

优化

我们发现可以这个方程刚好是 c_s 的沃尔什变换的结果，所以我们可以直接逆变换回去。

即把高斯消元替换为一次 IFWT。

时间复杂度：O(2^kn\log_2 n+q2^k)

参考代码（用的是EA的代码（依旧小编太懒了...））：

#include <bits/stdc++.h>

using ll = long long;
const int mod = 998244353;
const int inv2 = (mod + 1) / 2;
const int maxn = (1 << 20) + 12;

int n,m,k,w[maxn],z[maxn][12],F[maxn],G[maxn],W[maxn],sum[maxn],lg[maxn];

int qpow(int a,ll b) {
    if (b == 0) return 1;
    ll d = qpow(a,b>>1); d = d * d % mod;
    if (b & 1) d = d * a % mod;
    return d;
}

void DFT(int *a,int lim,int flag){
    for (int i = 1; i < lim; i <<= 1)
    for (int j = 0; j < lim; j += (i << 1))
    for (int k = 0; k < i; ++ k) {
        ll A0 = a[j+k], A1 = a[j+k+i];
        a[j+k] = (A0 + A1) % mod; a[j+k+i] = (A0 - A1 + mod) % mod;
        if (flag == -1) { 
            a[j+k] = (ll) a[j+k] * inv2 % mod;
            a[j+k+i] = (ll) a[j+k+i] * inv2 % mod;
        }
    }
}

int main() {
    scanf("%d%d%d",&n,&m,&k);
    int H[(1<<k)+10][(1<<m)+10];
    std::memset(H,0,sizeof(H));
    for (int i = 0; i < k; ++ i) scanf("%d",&w[i]);
    lg[0] = -1;
    for (int i = 0; i < (1 << k); ++ i) {
        if (i) lg[i] = lg[i/2] + 1;
        for (int j = 0; j < k; ++ j) {
            if (i & (1 << j)) W[i] = (W[i] - w[j] + mod) % mod;
            else W[i] = (W[i] + w[j]) % mod;
        }
    }
    for (int i = 1; i <= n; ++ i) {
        for (int j = 0; j < k; ++ j) 
            scanf("%d",&z[i][j]); 
        H[0][0] += 1;
        for (int mask = 1; mask < (1 << k); ++ mask) {
            int lowbit = mask & -mask;
            int b = lg[lowbit];
            sum[mask] = sum[mask - lowbit] ^ z[i][b];
            H[mask][sum[mask]] += 1;
        }
    }
    for (int mask = 0; mask < (1 << k); ++ mask)
        DFT(H[mask],(1<<m),1);
    for (int i = 0; i < (1 << m); ++ i) {
        for (int j = 0; j < (1 << k); ++ j) 
            G[j] = H[j][i];
        DFT(G,(1<<k),-1);
        F[i] = 1;
        for (int j = 0; j < (1 << k); ++ j) 
            F[i] = (ll) F[i] * qpow(W[j],G[j]) % mod;
    } DFT(F,(1<<m),-1);
    for (int i = 0; i < (1 << m); ++ i)
        printf("%d ",F[i]);
    return 0;
}

模板五：稀疏序列的异或卷积（只有两个非 0 位，）

问题描述：求 m 个序列的异或卷积，每个序列都恰好只有两个非 0 位。

为了方便描述，我们使用形式幂级数：

\prod_{i=1}^k u_ix^{a_i}+v_ix^{b_i}

我们一开始显然可以提出一个 x^{a_i} 来简化形式。

可以变为：

\prod_{i=1}^k x^{a_i}(u_i+v_ix^{a_i\oplus b_i})

单独处理前面的 x^a 部分，后面的部分按照 a \oplus b = c 的 c 分组，每组内乘积仍为 u + v x^c 的形式，且易于计算。

现在问题转化为计算 \prod\limits_{S\subseteq U_n} (c_S + d_S x^S)。

考虑采用分治乘法，分治的每一段区间形如 [a2^b, (a+1)2^b)，可以证明区间内的乘积一定形（归纳证明）：

\sum\limits_{S\subseteq U_b} r_S x^S + x^{a2^b} \sum\limits_{S\subseteq U_b} q_S x^S

于是将 \sum\limits_{S\subseteq U_b} r_S x^S 和 \sum\limits_{S\subseteq U_b} q_S x^S 视作两个定义在 U_b 上的集合幂级数进行维护，每次合并时仅需做 O(1) 次卷积，即可将总复杂度优化至 O(n^2 2^n)。（做法参考归纳证明）

更进一步，我们可以仅维护这两个集合幂级数在 U_b 上做沃尔什变换后的结果，经过简单讨论可验证，所有操作均可在 O(2^b) 的时间复杂度内完成，总复杂度可进一步优化至 O(n 2^n + m)，其中 m 为稀疏多项式的个数。

参考代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=(1<<22)+10;
const ll mod=998244353;
int n,q,len=0;
ll f[N],g[N];
void FWT(int len)
{
    for(int h=2;h<=len;h<<=1)
    {
        int mid=h/2;
        for(int i=0;i<len;i+=h)
        {
            for(int j=i;j<i+mid;j++)
            {
                ll a=f[j],b=f[j+mid],c=g[j],d=g[j+mid];
                f[j]=(a*b+c*d)%mod;f[j+mid]=(a*b-c*d)%mod;
                g[j]=(b*c+a*d)%mod;g[j+mid]=(b*c-a*d)%mod;
            }
        }
    }
}
void IFWT(int len,ll v)
{
    for(int h=2;h<=len;h<<=1)
    {
        int mid=h/2;
        for(int i=0;i<len;i+=h)
        {
            for(int j=i;j<i+mid;j++)
            {
                ll a=f[j],b=f[j+mid];
                f[j]=(a+b)*v%mod;f[j+mid]=(a-b+mod)*v%mod;
            }
        }
    }
}
int main()
{
    scanf("%d%d",&n,&q);len=(1<<n);
    for(int i=0;i<len;i++) f[i]=1,g[i]=0;
    for(int i=1;i<=q;i++)
    {
        ll a,b,c;scanf("%lld%lld%lld",&a,&b,&c);
        ll x=f[c],y=g[c];
        f[c]=(x*a+y*b)%mod;
        g[c]=(y*a+x*b)%mod;
    }
    FWT(len);
    for(int i=0;i<len;i++) f[i]=(f[i]+g[i])%mod;
    IFWT(len,(mod+1)/2);
    ll ans=0;
    for(ll i=0,sum=1;i<len;i++,sum=sum*2ll%mod) ans^=((f[i]*sum%mod+mod)%mod);
    printf("%lld",ans);
    return 0;
}

模板六：k 维异或卷积

众所周知异或是不进位的加法，那如果我们也就可以把模 2 变成模 k。

为了方便，我们先只探讨 k=10 的情况。

先分析贡献系数：c(i,j)c(i,k)=c(i,(j+k)\bmod k)

那么我们的 c(i,j) 可以变成一个有取模性质的值。

可以构造为单位根。

因为维度提升了，所以我们的分治每次就要分成 10 层，所以我们的矩阵就要变成 10\times 10。

即我们需要 c 的 10 组解，且需要存在逆矩阵。

我们可以考虑范德蒙德矩阵。

下面记 \omega=\omega_{10}

V=\begin{pmatrix} \omega^0 & \omega^0 & \cdots & \omega^0 \\ \omega^0 & \omega^1 & \cdots & \omega^9 \\ \vdots & \vdots & \ddots & \vdots \\ \omega^0 & \omega^9 & \cdots & \omega^{9\times9} \end{pmatrix}

即 V_{i,j}=\omega^{ij}

其逆矩阵为：

V^{-1}=\frac{1}{10}\begin{pmatrix} \omega^0 & \omega^0 & \cdots & \omega^0 \\ \omega^0 & \omega^{-1} & \cdots & \omega^{-9} \\ \vdots & \vdots & \ddots & \vdots \\ \omega^0 & \omega^{-9} & \cdots & \omega^{-9\times9} \end{pmatrix}

即 (V^{-1})_{i,j}=\frac{\omega^{-ij}}{10}

这里我们也能看出范德蒙德矩阵的内在性质，他的每一行都是 c_ic_j=c_{(i+j)\bmod k} 的一组解。

但在一些题目中，精度会发生问题，我们需要一个比较通用的一个做法。

我们考虑用一个代数结构来刻画单位根，我们可以让 x 表示 \omega，并且要满足 \omega 的一些性质。

性质一：x^k=1

那么我们不如直接多项式模 x^k-1

这样任何的 x^k 都会被消成 x^k-1

但这真的对吗？

为何不能直接模 x^k-1

我们先思考什么结构能让新的范德蒙德矩阵存在逆矩阵。

我们让范德蒙德矩阵乘他的期望逆矩阵。

对于 i 行乘 i 列，我们期望的结果是得到 1。

乘的结果为： \frac{1}{10}\sum\limits_{j=0}^{9}x^{ij}x^{-ij}=1

这满足条件。

那么 i 行乘 j 列呢？（i\not =j）

乘的结果为 \frac{1}{10}\sum\limits_{p=0}^9 x^{(i-j)p}

这个结果无法整除 x^k-1 。

因为要想整除 x^k-1 ，那么必须整除 \forall i,(x-\omega^i) 。

这等价于代入 \omega^i 的结果为 0。

但能发现这可能做不到。

比如 i-j=4，那么结果为 x^0+x^2+x^4+x^6+x^8

代入 \omega^5 后，能发现结果为 5

所以我们该怎么半？

结论：模分圆多项式是对的

为何可以模分圆多项式

先要证 x^k 会变为 1。

x_k\bmod \Phi_{10}(x)=(x_k-1)+1\bmod \Phi_{10}(x)

根据分圆多项式的性质：x^n-1=\prod\limits_{d|n}\Phi_d(x) 可得：(x_k-1)\bmod \Phi_{10}(x)=0

所以 x_k\bmod \Phi_{10}(x)=1

再证范德蒙德矩阵乘他的期望逆矩阵得到单位矩阵。

和上文一样，对于 i 行乘 i 列，我们的结果显然为 1。

那么 i 行乘 j 列呢？（i\not =j）

乘的结果为 \frac{1}{10}\sum\limits_{p=0}^9 x^{(i-j)p}

然后这相当于每次走 i-j 步，在环上走 n 次。

这次我们只需要代入 \forall \gcd(h,10)=1,\omega^h

能发现对于这样的 h 来说，hp\bmod k 两两不同。

考虑反证法，若存在 i,j,i\not= j 使得 hi\equiv hj\pmod k。

因为 \gcd(h,k)=1，所以 h 存在逆元。

所以 i\equiv hh^{-1}j\pmod k

因此 i\equiv j\pmod k，这与假设矛盾，因此得证。

回到矩阵乘法结果：\frac{1}{10}\sum\limits_{p=0}^9 x^{(i-j)p}=\frac{1}{10}\sum\limits_{p=0}^9 \omega^{(i-j)hp}

根据上文结论： p(i-j)\bmod k 两两不同。

则原式 =\frac{1}{10}\sum\limits_{q=0}^9 \omega^{qh}

下面来分类讨论：

  1. 若 \gcd(h,k)=1，那么 qh\bmod k 两两不同，所以结果等于 0。
  2. 否则，定义 g=\gcd(h,k)，原式 =\frac{1}{10}\sum\limits_{q=0}^9 \omega^{qgp'}=\frac{1}{10}\sum\limits_{q=0}^9 (\omega_k^g)^{qh'}=\frac{1}{10}\sum\limits_{q=0}^9 (\omega_{\frac{k}{g}})^{q\frac{h}{g}}，因为 \gcd(\frac{h}{g},k)=1，所以 q\frac{h}{g}\bmod \frac{k}{g} 两两不同，所以结果等于 0。

因此我们就证明了结果等于 0。

自此我们就证明了结果可以模分圆多项式。

最终实现

但是模分圆多项式复杂度太大了，因为 \Phi_{10}(x)|x^{10}-1，所以我们过程中可以模 x^{10}-1，然后最后再模分圆多项式。

原题链接：CF1103E Radix sum

参考代码：

#include<bits/stdc++.h>
using namespace std;
typedef unsigned long long ll;
const int N=1e5+10;
int n;
ll a[N],c[11][11];
struct poly
{
    ll a[10];
    void clear(){for(int i=0;i<10;i++) a[i]=0;}
    poly operator*(poly b)
    {
        poly c;
        for(int i=0;i<10;i++) c.a[i]=0;
        for(int i=0;i<10;i++)
            for(int j=0;j<10;j++)
                c.a[(i+j)%10]+=a[i]*b.a[j];
        return c;
    } 
    poly operator*(ll b)
    {
        poly c;
        for(int i=0;i<10;i++) c.a[i]=0;
        for(int i=0;i<10;i++) c.a[(i+b)%10]+=a[i];
        return c;
    } 
    poly operator+(poly b)
    {
        poly c;
        for(int i=0;i<10;i++) c.a[i]=a[i]+b.a[i];
        return c;
    }
    poly operator-(poly b)
    {
        poly c;
        for(int i=0;i<10;i++) c.a[i]=a[i]-b.a[i];
        return c;
    } 
}f[N],tmp[N];
ll len=100000;
ll ku(ll a,ll b)
{
    ll ans=1;
    while(b)
    {
        if(b&1) ans=ans*a;
        a=a*a;b>>=1;
    }
    return ans;
}
ll b[N];
void divide(ll *a)
{
    int n=10,m=5;
    ll inv=b[m-1];
    for(int i=n-1;i>=m-1;i--)
    {
        ll d=a[i]*inv;
        for(int j=i;j>=i-m+1;j--) a[j]=a[j]-d*b[m-(i-j)-1];
    }

}
void FWT(poly *f,int op)
{
    for(int h=10;h<=100000;h*=10)
    {
        int mid=h/10;
        for(int i=0;i<100000;i+=h)
        {
            for(int j=i;j<i+mid;j++)
            {
                for(int k=0;k<10;k++) tmp[k]=f[j+k*mid],f[j+k*mid].clear();
                for(int k=0;k<10;k++)
                {
                    for(int p=0;p<10;p++)
                    {
                        f[j+k*mid]=(f[j+k*mid]+(tmp[p]*((op*c[k][p]+10)%10)));
                    }
                }
            }
        }
    }
}
poly ksm(poly a,int b)
{
    poly ans;ans.clear();ans.a[0]=1;
    while(b)
    {
        if(b&1) ans=ans*a;
        a=a*a;b>>=1;
    }
    return ans;
}
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++) scanf("%llu",&a[i]),f[a[i]].a[0]++;
    for(int i=0;i<10;i++)
    {
        for(int j=0;j<10;j++)
        {
            c[i][j]=(i*j%10+10)%10;
        }
    }
    FWT(f,1);
    for(int i=0;i<100000;i++) f[i]=ksm(f[i],n);
    FWT(f,-1);
    b[0]=1;b[1]=-1;b[2]=1;b[3]=-1;b[4]=1;
    for(int i=0;i<n;i++)
    {
        divide(f[i].a);
        ll mod=(1ll<<58);
        ll x=(f[i].a[0]*6723469279985657373)<<1>>6;
        printf("%llu\n",x%mod);
    }
    return 0;
}

致谢和参考资料

感谢 CJX，CXM老师，ZZH与我激情讨论。

参考资料

  1. EternalAlexander：P6097 【模板】子集卷积
  2. 武林的2026年集训队论文：集合幂级数中的稀疏多项式乘法
  3. Fucious_Yin的FWT 小记