FWT 学习笔记

标签：begin return int end 笔记学习 pmatrix FWT ModInt

快速沃尔什变换

模板题

给定长度为 \(2^n\) 两个序列 \(A,B\)，设

\[C_i=\sum_{j\oplus k = i}A_j \times B_k \]

分别当 \(\oplus\) 是 or, and, xor 时求出 \(C\)。

FWT，中文名称：快速沃尔什变换。因为已经有 FFT 和 NTT 基础，所以直接考虑构造 FWT 的变换。

不失一般性，先考虑 \(n = 1\)，即数组长度为 \(2\) 的情况（运算为 and）：

\[\begin{cases} C_0 = A_0 \cdot B_0 + A_1 \cdot B_0 + A_0 \cdot B_1 \\ C_1 = A_1 \cdot B_1 \end{cases} \]

那么，如果将 \(A\) 视作一个长度为 \(2^n\) 的向量，并构造下列变换：

对于 \(|A| = 2\)：

\[f\begin{pmatrix}A_0 \\ A_1\end{pmatrix} = f\begin{pmatrix}A_0 + A_1 \\ A_1\end{pmatrix} \\ \]

对于 \(|A| > 2\)：

\[f\begin{pmatrix}A_0 \\ A_1 \\ \vdots \\ A_{2^n - 1}\end{pmatrix} = f\begin{pmatrix}f\begin{pmatrix}A_0 \\ A_1 \\ \vdots \\ A_{2^{n - 1} - 1}\end{pmatrix} \\ f\begin{pmatrix}A_{2^{n - 1}} \\ A_{2^{n - 1} + 1} \\ \vdots \\ A_{2^n - 1}\end{pmatrix}\end{pmatrix} \]

作为一个变换，它是否满足 \(f(A) \times f(B) = f(C)\)，让我们验证一下。容易发现，根据数学归纳法，如果对于 \(n = 1\) 成立，那么对于 \(n > 1\) 都成立。下面是 \(n = 1\) 的式子：

\[\begin{aligned} f(A) &= \begin{pmatrix}A_0 + A_1 \\ A_1\end{pmatrix} \\ f(B) &= \begin{pmatrix}B_0 + B_1 \\ B_1\end{pmatrix} \\ f(A \cdot B) &= \begin{pmatrix}A_0B_0 + A_1B_0 + A_0B_1 + A_1B_1 \\ A_1B_1\end{pmatrix}\\ f(C) &= \begin{pmatrix}C_0 + C_1 \\ C_1\end{pmatrix} = \begin{pmatrix}A_0B_0 + A_1B_0 + A_0B_1 + A_1B_1 \\ A_1B_1\end{pmatrix} = f(A \cdot B) \end{aligned} \]

由此可见，该构造满足 FWT 在 and 运算下的变换。

对于 or 和 xor 运算类似构造，只考虑更改 \(n = 1\) 时的变换，那么有：

对于 or 运算：
\[f(A) = \begin{pmatrix}A_0 \\ A_0 + A_1 \end{pmatrix} \]
对于 xor 运算：
\[f(A) = \begin{pmatrix}A_0 + A_1 \\ A_0 - A_1 \end{pmatrix} \]

有了 FWT 自然有 IFWT，其变换构造方式很显然，不再详细赘述。特别地，对于 xor 运算，可以发现其逆变换为：

\[f'(A) = \begin{pmatrix} \frac{A_0 + A_1}{2} \\ \frac{A_0 - A_1}{2} \end{pmatrix} = \frac{f(A)}{2} \]

在代码实现时，不用专门写其 IFWT，只需重新 FWT 一遍后，将 \(C_i \gets \frac{C_i}{2^n}\) 即可。

代码（运算为 xor）

template<class T>
void FWT(T *A, int N) {
    for (int i = 1; i < N; i *= 2) {
        for (int j = 0; j < N; j += 2 * i) {
            for (int k = 0; k < i; ++k) {
                auto a = A[j + k];
                auto b = A[i + j + k];
                A[j + k] = a + b;
                A[i + j + k] = a - b;
            }
        }
    }
}

时间复杂度为 \(O(n2^n)\)。

一些性质：

假设 \(A\) 序列经过 FWT 后变为了序列 \(A'\)，那么有：

若运算为 or，则：
- 正变换：\(A'_i = \sum\limits_{j \subseteq i} A_j\)；
- 逆变换：\(A_i = \sum\limits_{j \subseteq i} (-1)^{|i| - |j|} A'_j\)。
若运算为 and，则：
- 正变换：\(A_i' = \sum_{j \supseteq i} A_i\)；
- 逆变换：\(A_i = \sum_{j \supseteq i} (-1)^{|j| - |i|} A'_j\)。
若运算为 xor，则：
- 正变换：\(A_i' = \sum_j (-1)^{|i \cup j|} A_j\)；
- 逆变换：\(A_i = \frac{1}{2^n} \sum_j (-1)^{|i \cup j|} A'_j\)。

这些都是重要性质，常常用到。

子集卷积

模板题

给定两个长度为 \(2^n\) 的序列 \(a_0,a_1,\cdots,a_{2^n-1}\) 和 \(b_0,b_1,\cdots,b_{2^n-1}\)，你需要求出一个序列 \(c_0,c_1,\cdots,c_{2^n-1}\)，其中 \(c_k\) 满足：

\[c_k=\sum_{\substack{{i \& j=0}\\{i~\mid~ j=k}}} a_i b_j \]

其中\(~\mid~\)表示按位或，\(\&\)表示按位与。

答案对 \(10^9+9\) 取模。

如果只有一个条件 \(i \mid j = k\) 就是 FWT 板题。考虑多加一维得到 \(f(i, S)\)，定义是 \(f(i, S) = \sum\limits_{|T| = i \land T \subseteq S} a_T\)，而 \(g(i, S)\) 是基于 \(b\) 序列的类似的函数，显然这两个函数可以通过在同一层 FWT 一遍得到。又设 \(h(i, S) = \sum\limits_{j = 0}^i f(j, S)g(i - j, S)\)，发现 \(h(i, S)\) 其实就等于 \(\sum\limits_{|T| = i \land T \subseteq S} c_T\)，于是做一遍 IFWT 就得到了 \(c\)。

代码

#include <bits/stdc++.h>

using u32 = unsigned;
using i64 = long long;
using u64 = unsigned long long;

template<typename T>
T power(T a, i64 b) {
    T res = 1;
    for (; b; b /= 2, a = a * a) {
        if (b & 1) {
            res = res * a;
        }
    }
    return res;
}

template<const i64 P>
class ModInt {
public:
    i64 x;

    ModInt() : x{0} {}
    ModInt(int n) : x{(n % getMod() + getMod()) % getMod()} {}
    ModInt(i64 n) : x{(n % getMod() + getMod()) % getMod()} {}

    static i64 Mod;
    constexpr static void setMod(i64 _Mod) {
        Mod = _Mod;
    }
    constexpr static i64 getMod() {
        return P > 0 ? P : Mod;
    }

    constexpr ModInt &operator+=(ModInt k) & {
        x = (x + k.x >= getMod() ? x + k.x - getMod() : x + k.x);
        return *this;
    }
    constexpr ModInt &operator-=(ModInt k) & {
        x = (x - k.x < 0 ? x - k.x + getMod() : x - k.x);
        return *this;
    }
    constexpr ModInt &operator*=(ModInt k) & {
        x = 1LL * x * k.x % getMod();
        return *this;
    }

    friend constexpr ModInt operator+(ModInt lhs, ModInt rhs) {
        return lhs += rhs;
    }
    friend constexpr ModInt operator-(ModInt lhs, ModInt rhs) {
        return lhs -= rhs;
    }
    friend constexpr ModInt operator*(ModInt lhs, ModInt rhs) {
        return lhs *= rhs;
    }
    friend constexpr ModInt operator/(ModInt lhs, ModInt rhs) {
        return lhs /= rhs;
    }

    constexpr ModInt inv() {
        return power(*this, getMod() - 2);
    }
    constexpr ModInt &operator/=(ModInt k) & {
        return (*this) *= k.inv();
    }

    friend constexpr std::istream &operator>>(std::istream &is, ModInt &k) {
        i64 val;
        is >> val;
        k = val;
        return is;
    }
    friend constexpr std::ostream &operator<<(std::ostream &os, ModInt k) {
        return os << k.x;
    }

    friend constexpr bool operator==(ModInt lhs, ModInt rhs) {
        return lhs.x == rhs.x;
    }
    friend constexpr bool operator!=(ModInt lhs, ModInt rhs) {
        return lhs.x != rhs.x;
    }
    constexpr bool operator!() {
        return !x;
    }

    constexpr ModInt &operator++() & {
        return (*this) += 1;
    }
    constexpr ModInt &operator++(int) & {
        ModInt temp = *this;
        (*this) += 1;
        return temp;
    }
    constexpr ModInt &operator--() & {
        return (*this) -= 1;
    }
    constexpr ModInt &operator--(int) & {
        ModInt temp = *this;
        (*this) -= 1;
        return temp;
    }
} ;

template<>
i64 ModInt<0>::Mod = 1E9 + 7;
constexpr int P = 1E9 + 9;
using Z = ModInt<P>;

struct Comb {
    std::vector<Z> _inv, _invfac, _fac;
    int n;

    Comb() {
        n = 0;
        _inv.push_back(0);
        _invfac.push_back(1);
        _fac.push_back(1);
    }
    void Init(int m) {
        if (m <= n) {
            return ;
        }

        _inv.resize(m + 1);
        _fac.resize(m + 1);
        _invfac.resize(m + 1);

        for (int i = n + 1; i <= m; ++i) {
            _fac[i] = _fac[i - 1] * i;
        }
        _invfac[m] = _fac[m].inv();
        for (int i = m; i > n; --i) {
            _invfac[i - 1] = _invfac[i] * i;
            _inv[i] = _invfac[i] * _fac[i - 1];
        }
    }

    Z fac(int m) {
        Init(2 * m);
        return _fac[m];
    }
    Z inv(int m) {
        Init(2 * m);
        return _inv[m];
    }
    Z invfac(int m) {
        Init(2 * m);
        return _invfac[m];
    }
    Z binom(int n, int m) {
        if (m < 0 || m > n) {
            return 0;
        } else {
            return fac(n) * invfac(m) * invfac(n - m);
        }
    }
} comb;

void FWT(auto &A, int N, int coef) {
    for (int i = 1; i < N; i *= 2) {
        for (int j = 0; j < N; j += 2 * i) {
            for (int k = 0; k < i; ++k) {
                A[i + j + k] += coef * A[j + k];
            }
        }
    }
}

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int n;
    std::cin >> n;

    std::vector f(n + 1, std::vector<Z>(1 << n, 0)), g(f), h(g);
    for (int i = 0; i < (1 << n); ++i) {
        std::cin >> f[std::popcount((u32)i)][i];
    }
    for (int i = 0; i < (1 << n); ++i) {
        std::cin >> g[std::popcount((u32)i)][i];
    }

    for (int i = 0; i <= n; ++i) {
        FWT(f[i], (1 << n), 1);
        FWT(g[i], (1 << n), 1);
    }

    for (int i = 0; i <= n; ++i) {
        for (int j = 0; j <= i; ++j) {
            for (int S = 0; S < (1 << n); ++S) {
                h[i][S] += f[j][S] * g[i - j][S];
            }
        }
    }

    for (int i = 0; i <= n; ++i) {
        FWT(h[i], (1 << n), -1);
    }
    for (int i = 0; i < (1 << n); ++i) {
        std::cout << h[std::popcount((u32)i)][i] << " ";
    }

    return 0;
}

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

标签：begin,return,int,end,笔记,学习,pmatrix,FWT,ModInt
From： https://www.cnblogs.com/CTHOOH/p/18367604

快速沃尔什变换

子集卷积

相关文章

赞助商

阅读排行