k进制快速沃尔什变换（k进制FWT）

引入

我们已经知道了多项式版的二进制 FWT，但是我们似乎不是很好将其扩展到 $k$ 进制下，于是我们来考虑另一种方式的 FWT。

原理

二进制

规定三个多项式的长度为 $2^n$，如果不够则往后面补 $0$。

还是先考虑二进制，我们假设存在矩阵使得 $TA\ast TB=TC$，考虑矩阵长什么样。

我们可以得到如下式子

$$FWT(A{)}_x=\sum _{i=0}^{2^n-1}{w}_{x,i}\ast A_i$$

因为 $FWT(A)\ast FWT(B)=FWT(C)$，所以 $FWT(A{)}_x\ast FWT(B{)}_x=FWT(C{)}_x$，可写为

$$\sum _{i=0}^{2^n-1}{w}_{x,i}\ast A_i\ast \sum _{j=0}^{2^n-1}{w}_{x,j}\ast B_j=\sum _{k=0}^{2^n-1}{w}_{x,k}\ast C_k$$

即

$$\sum _{i=0}^{2^n-1}\sum _{j=0}^{2^n-1}{w}_{x,i}\ast w_{x,j}\ast A_i\ast B_j=\sum _{k=0}^{2^n-1}{w}_{x,k}\ast C_k$$

因为我们要使得

$$\sum _{i\oplus j=k}{A}_i\ast B_j=C_k$$

即要使得以下式子成立

$$\sum _{i\oplus j=k}{w}_{x,i}\ast w_{x,j}\ast A_i\ast B_j=w_{x,k}\ast C_k$$

所以我们可以得出

$$w_{x,i}\ast w_{x,j}=w_{x,k}(i\oplus j=k)$$

此时的矩阵就是范德蒙德矩阵。（就是我们再莫比乌斯反演里面学的哪个单位根矩阵）

我们考虑如何快速求出 $FWT(A{)}_x$。

我们定义 $A_0$ 表示前半段，$A_1$ 表示后半段，$x_1$ 表示 $x$ 的最高位（如果最高位为 $1$ 则 $x_1$ 为 $1$，否则为 $0$），$x_0$ 表示 $x$ 除了最高位其他的位，$i_0$ 也表示 $i$ 除了最高位其它的位。

则我们可以推出

$$\begin{array}{rl}FWT(A{)}_x& =\sum _{i=0}^{2^n-1}{w}_{x,i}\ast A_i\\ & =\sum _{i=0}^{2^{n-1}-1}{w}_{x,i}\ast A_i+\sum _{i=2^n}^{2^n-1}{w}_{x,i}\ast A_i\\ & =\sum _{i=0}^{2^{n-1}-1}{w}_{x_1,0}\ast w_{x_0,i_0}\ast A_i+\sum _{i=2^n}^{2^n-1}{w}_{x_1,1}\ast w_{x_0,i_0}\ast A_i\\ & =w_{x_1,0}\ast FWT(A_0{)}_{x_0}+w_{x_1,1}\ast FWT(A_1{)}_{x_0}\end{array}$$

于是就跟多项式下的 FWT 一样可以递归求了。

k进制

于是我们进入正题，来看看 $k$ 进制下的 FWT。（其实跟二进制下很像，只不过有一些细节上的区别）

同样规定多项式长度为 $k^n$。

还是要构造矩阵使得 $TA\ast TB=TC$。

同样有

$$\begin{array}{rl}FWT(A{)}_x& =\sum _{i=0}^{k^n-1}{w}_{x,i}\ast A_i\\ & =\sum _{i=0}^{k-1}{w}_{x_1,i}\ast FWT(A_i{)}_{x_0}\end{array}$$

上面的定义和二进制类似，证明也类似，就直接写出来算了。（绝对不是因为懒）

此时一样递归求就可以了。

矩阵就是下面的矩阵，其中 $w_k$ 代表 $k$ 下的单位根。（跟 FFT 里的矩阵一样）

$$\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& {\omega }_n^1& {\omega }_n^2& \cdots & {\omega }_n^{n-1}\\ 1& {\omega }_n^2& {\omega }_n^4& \cdots & {\omega }_n^{2(n-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& \cdots & {\omega }_n^{(n-1)(n-1)}\end{array}\right]$$

可以看看代码加深理解！！！

代码

我这份代码是这道题的，三进制下的。

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
int n, k, len = 1;
LL ans;
complex<double> a[1000005];
const complex<double> w = { -0.5, 0.5 * sqrt(3) }, w2 = { -0.5, -0.5 * sqrt(3) };
int in() {
    char ch = getchar();
    int s = 0;
    while (ch < '0' || ch > '9') ch = getchar();
    while (ch <= '9' && ch >= '0') s = s * 3 + ch - '1', ch = getchar();
    return s;
}
void FWT(complex<double> *f, int flag) {
    for (int mid = 1; mid < len; mid = mid * 3) {
        for (int i = 0; i < len; i = i + mid * 3) {
            for (int j = i; j < i + mid; j++) {
                complex<double> t0 = f[j], t1 = f[j + mid], t2 = f[j + mid * 2];
                if (flag == 1) {
                    f[j] = t0 + t1 + t2;
                    f[j + mid] = t0 + t1 * w + t2 * w2;
                    f[j + mid * 2] = t0 + t1 * w2 + t2 * w;
                } else {
                    f[j] = t0 + t1 + t2;
                    f[j + mid] = t0 + t1 * w2 + t2 * w;
                    f[j + mid * 2] = t0 + t1 * w + t2 * w2;
                    double t = f[j].real();
                    f[j].real(t / 3);

                    t = f[j + mid].real();
                    f[j + mid].real(t / 3);

                    t = f[j + mid * 2].real();
                    f[j + mid * 2].real(t / 3);

                    t = f[j].imag();
                    f[j].imag(t / 3);

                    t = f[j + mid].imag();
                    f[j + mid].imag(t / 3);

                    t = f[j + mid * 2].imag();
                    f[j + mid * 2].imag(t / 3);
                }
            }
        }
    }
}
int main() {
    scanf("%d%d", &n, &k);
    for (int t = 0; t < k; t++) len = len * 3;
    for (int i = 0; i < n; i++) a[in()].real(1);
    FWT(a, 1);
    for (int i = 0; i < len; i++) a[i] = a[i] * a[i] * a[i];
    FWT(a, -1);
    ans = a[0].real() + 0.5;
    printf("%lld\n", (ans - n) / 6);
    return 0;
}