首页 > 其他分享 >快速沃尔什变换与子集卷积

快速沃尔什变换与子集卷积

时间:2024-02-08 10:24:40浏览次数:35  
标签:nonumber 卷积 mathscr sum int 子集 沃尔什 operatorname MOD

前置知识:FFT(快速傅里叶变换)。

快速沃尔什变换

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

快速沃尔什变换(Fast Walsh–Hadamard transform)解决二进制运算下的卷积。

给定序列 \(f,g\),求以下三个序列 \(A,B,C\):

\[A_i=\sum_{j \operatorname{or} k = i} f_j \times g_k \nonumber \]

\[A_i=\sum_{j \operatorname{and} k = i} f_j \times g_k \nonumber \]

\[A_i=\sum_{j \operatorname{xor} k = i} f_j \times g_k \nonumber \]

其中 \(\operatorname{or}\)、\(\operatorname{and}\)、\(\operatorname{xor}\) 为二进制的或、与、异或运算。

显然可以 \(\mathcal{O}(3^n)\) 做平凡的子集枚举,但这不是我们想要的,我们的目标是 \(\mathcal{O}(n \log n)\) 求,其中 \(n\) 是序列长度。

FWT运用了和FFT相同的思想,将序列变换后做运算再逆变换回去来加速卷积。

或卷积

考虑构造一个运算 \(\mathscr{F}\) 使得若

\[A_i=\sum_{j \operatorname{or} k = i} f_j \times g_k \nonumber \]

\[(\mathscr{F}A)_i=(\mathscr{F}f)_i \times (\mathscr{F}g)_i \nonumber \]

并且给定序列 \(f\),我们可以快速地计算出 \(\mathscr{F}A\) 以及 \(\mathscr{F}^{-1}A\),其中 \(\mathscr{F}^{-1}\) 表示 \(\mathscr{F}\) 的逆变换,即 \(\mathscr{F}^{-1} (\mathscr{F}A)=A\)。

构造算子 \(\mathscr{F}\) 满足:

\[(\mathscr{F}A)_i=\sum_{j \operatorname{or} i = i}A_j \nonumber \]

也就是 \((\mathscr{F}A)_i\) 是所有 \(A_j\) 的和,其中下标 \(j\) 是 \(i\) 表示的集合的子集。

容易验证 \((\mathscr{F}A)_i=(\mathscr{F}f)_i \times (\mathscr{F}g)_i\):

\[\begin{aligned} (\mathscr{F}f)_i \times (\mathscr{F}g)_i=\sum_{j \operatorname{or} i = i}\sum_{k \operatorname{or} i = i} f_jg_k=\sum_{(j \operatorname{or} k) \operatorname{or}i = i} f_jg_k=(\mathscr{F}A)_i \end{aligned} \nonumber \]

于是我们只需要快速对一个序列进行变换及逆变换。

对一个序列 \(A\) 进行变换是容易的,先将序列补齐为 \(2^k\) 项,考虑分治:设 \(A_0\) 表示 \(A\) 的前 \(2^{k-1}\) 项(即下标首位为 \(0\) 的位置)构成的序列,\(A_1\) 表示 \(A\) 的后 \(2^{k-1}\) 项即下标首位为 \(1\) 的位置构成的序列,根据 \(\mathscr{F}A_0\) 和 \(\mathscr{F}A_1\) 求出 \(\mathscr{F}A\)。

对于 \(\mathscr{F}A\) 的前 \(2^{k-1}\) 项,它只能由前 \(2^{k-1}\) 项贡献过来,因为后 \(2^{k-1}\) 项首位都是 \(1\),不可能或出 \(0\) 来;对于后 \(2^{k-1}\) 项,整个序列都可以贡献,因此得到

\[\mathscr{F}A=\operatorname{merge}\Big( \mathscr{F}A_0,\mathscr{F}A_0+\mathscr{F}A_1 \Big) \nonumber \]

其中 \(\operatorname{merge}(P,Q)\) 就是将两个序列像字符串那样拼起来,加法就是对两个序列对应位置做加法得到的新序列。

逆变换也就呼之欲出:

\[\mathscr{F}^{-1}A=\operatorname{merge}\Big( \mathscr{F}^{-1}A_0,\mathscr{F}A_1-\mathscr{F}^{-1}A_0 \Big) \nonumber \]

时间复杂度 \(\mathcal{O}(n \log n)\),代码如下:

inline void OR_FWT(int A[], int n, int type) {
// type = -1 表示逆变换,否则type = 1
// n 是一个 2^k
    for(int h = 2, H = 1; h <= n; h <<= 1, H <<= 1) {
        for(int i = 0; i < n; i += h) 
            for(int j = 0; j < H; j ++) 
                A[i + j + H] = ((A[i + j + H] + type * A[i + j]) % MOD + MOD) % MOD;
    }
}

与卷积

构造和或卷积类似的算子 \(\mathscr{F}\) 使得

\[(\mathscr{F}A)_i=\sum_{j \operatorname{and} i = i} A_j \nonumber \]

容易验证若 \(A_i=\sum\limits_{j \operatorname{and} k = i} f_jg_k\),那么 \((\mathscr{F}A)_i=(\mathscr{F}f)_i \times (\mathscr{F}g)_i\)。

同样的,将 \(A\) 补齐到 \(2^k\) 项,设 \(A_0\) 表示前 \(2^{k-1}\) 项构成的序列,\(A_1\) 表示后 \(2^{k-1}\) 项构成的序列,不难得到

\[\begin{aligned} &\mathscr{F}A=\operatorname{merge}\Big( \mathscr{F}A_0+\mathscr{F}A_1,\mathscr{F}A_1 \Big)\newline &\mathscr{F}^{-1}A=\operatorname{merge}\Big( \mathscr{F}^{-1}A_0-\mathscr{F}^{-1}A_1,\mathscr{F}^{-1}A_1 \Big) \end{aligned} \nonumber \]

时间复杂度同样是 \(\mathcal{O}(n \log n)\),代码如下:

inline void AND_FWT(int A[], int n, int type) {
// type = -1 表示逆变换,否则type = 1
// n 是一个 2^k
    for(int h = 2, H = 1; h <= n; h <<= 1, H <<= 1) {
        for(int i = 0; i < n; i += h)
            for(int j = 0; j < H; j ++)
                A[i + j] = ((A[i + j] + type * A[i + j + H]) % MOD + MOD) % MOD;
    }
}

异或卷积

异或的情况就比较不同了,我们不能照搬之前与和或的形式来做变换(\(j \operatorname{xor} i = i \Rightarrow j = 0\),这样的变换肯定是不可行的)。

接下来给出异或的正确变换,为了方便,对于数 \(x\),记 \(|x|\) 表示 \(x\) 的二进制表示中 \(1\) 的个数,即 \(|x|=\operatorname{popcount}(x)\)。

基于奇偶性,由于 \(\big((a \operatorname{and} b) \bmod 2\big) \operatorname{xor} \big((b \operatorname{and} c)]\bmod 2 \big)=\big( (a \operatorname{and} c) \bmod 2 \big)\),得到变换 \(\mathscr{F}\):

\[(\mathscr{F}A)_i=\sum_{|j \operatorname{and} i| \bmod 2 = 0}A_j-\sum_{|j \operatorname{and} i|\bmod 2=1}A_j \nonumber \]

为了方便,记 \(S_i\) 表示满足 \(|j \operatorname{and} i| \bmod 2 = 0\) 的 \(j\) 构成的集合,\(T_i\) 表示满足 \(|j \operatorname{and} i| \bmod 2 = 1\) 的 \(j\) 构成的集合,那么上式可以就写成

\[(\mathscr{F}A)_i=\sum_{j \in S_i}A_j - \sum_{j \in T_i}A_j \nonumber \]

下证若 \(A_i=\sum\limits_{j \operatorname{xor} k=i}f_jg_k\),那么 \((\mathscr{F}A)_i=(\mathscr{F}f)_i \times (\mathscr{F}g)_i\)。

证明:

\[\begin{aligned} (\mathscr{F}f)_i \times (\mathscr{F}g)_i&=\left(\sum_{j \in S_i}f_j-\sum_{j \in T_i}f_j\right) \times \left(\sum_{j \in S_i}g_j-\sum_{j \in T_i}g_j\right)\newline &=\sum_{j,k \in S_i}f_jg_k + \sum_{j,k \in T_i}f_jg_k-\sum_{j \in S_i, k \in T_i}f_jg_k - \sum_{j \in T_i, k \in S_i}f_jg_k\newline &=\sum_{(j \operatorname{xor} k) \in S_i}f_jg_k-\sum_{(j \operatorname{xor} k) \in T_i}f_jg_k\newline &=(\mathscr{F}A)_i \end{aligned} \nonumber \]

然后还是分治求 \(\mathscr{F}A\) 和 \(\mathscr{F}^{-1}A\),将 \(A\) 补齐为 \(2^k\) 项,设 \(A_0\) 为前 \(2^{k-1}\) 项构成的序列,\(A_1\) 为后 \(2^{k-1}\) 项构成的序列,那么

\[\begin{aligned} &\mathscr{F}A=\operatorname{merge}(\mathscr{F}A_0+\mathscr{F}A_1,\mathscr{F}A_0-\mathscr{F}A_1)\newline &\mathscr{F}^{-1}A=\operatorname{merge}(\frac{\mathscr{F}^{-1}A_0+\mathscr{F}^{-1}A_1}{2},\frac{\mathscr{F}^{-1}A_0-\mathscr{F}^{-1}A_1}{2}) \end{aligned} \nonumber \]

时间复杂度 \(\mathcal{O}(n \log n)\),代码如下:

inline void XOR_FWT(int A[], int n, int type) {
// type = 2的逆元 表示逆变换,否则type = 1
// n 是一个 2^k
    for(int h = 2, H = 1, x, y; h <= n; h <<= 1, H <<= 1) {
        for(int i = 0; i < n; i += h)
            for(int j = 0; j < H; j ++) {
                x = A[i + j], y = A[i + j + H]; 
                A[i + j] = 1ll * (x + y) % MOD * type % MOD, A[i + j + H] = 1ll * ((x - y) % MOD + MOD) % MOD * type % MOD;
            }
    }
}

三个FWT合并起来的代码如下:

code for Luogu P4717
#include <bits/stdc++.h>

using namespace std;

const int maxn = (1 << 18) + 10, MOD = 998244353, inv2 = 499122177;
int n, A[maxn], B[maxn], C[maxn];

inline void print(int A[], int n) {
    for(int i = 0; i < n; i ++) cout << A[i] << ' ';
    cout << '\n';
}

inline void OR_FWT(int A[], int n, int type) {
// type = -1 表示逆变换,否则type = 1
// n 是一个 2^k
    for(int h = 2, H = 1; h <= n; h <<= 1, H <<= 1) {
        for(int i = 0; i < n; i += h) 
            for(int j = 0; j < H; j ++) 
                A[i + j + H] = ((A[i + j + H] + type * A[i + j]) % MOD + MOD) % MOD;
    }
}
inline void AND_FWT(int A[], int n, int type) {
// type = -1 表示逆变换,否则type = 1
// n 是一个 2^k
    for(int h = 2, H = 1; h <= n; h <<= 1, H <<= 1) {
        for(int i = 0; i < n; i += h)
            for(int j = 0; j < H; j ++)
                A[i + j] = ((A[i + j] + type * A[i + j + H]) % MOD + MOD) % MOD;
    }
}
inline void XOR_FWT(int A[], int n, int type) {
// type = 2的逆元 表示逆变换,否则type = 1
// n 是一个 2^k
    for(int h = 2, H = 1, x, y; h <= n; h <<= 1, H <<= 1) {
        for(int i = 0; i < n; i += h)
            for(int j = 0; j < H; j ++) {
                x = A[i + j], y = A[i + j + H]; 
                A[i + j] = 1ll * (x + y) % MOD * type % MOD, A[i + j + H] = 1ll * ((x - y) % MOD + MOD) % MOD * type % MOD;
            }
    }
}

int main() {
    ios::sync_with_stdio(false), cin.tie(0);

    cin >> n; n = (1 << n);
    for(int i = 0; i < n; i ++) cin >> A[i];
    for(int i = 0; i < n; i ++) cin >> B[i];

    OR_FWT(A, n, 1), OR_FWT(B, n, 1);
    for(int i = 0; i < n; i ++) C[i] = 1ll * A[i] * B[i] % MOD;
    OR_FWT(C, n, -1); print(C, n); OR_FWT(A, n, -1), OR_FWT(B, n, -1);

    AND_FWT(A, n, 1), AND_FWT(B, n, 1);
    for(int i = 0; i < n; i ++) C[i] = 1ll * A[i] * B[i] % MOD;
    AND_FWT(C, n, -1); print(C, n); AND_FWT(A, n, -1), AND_FWT(B, n, -1);

    XOR_FWT(A, n, 1), XOR_FWT(B, n, 1);
    for(int i = 0; i < n; i ++) C[i] = 1ll * A[i] * B[i] % MOD;
    XOR_FWT(C, n, inv2); print(C, n); XOR_FWT(A, n, inv2), XOR_FWT(B, n, inv2);
    
    return 0;
}

集合幂级数

就是将集合放到指数位置上,例如 \(f(x)=x^{\emptyset}+2x^{\{1\}}+4x^{\{2\}}+3x^{\{1,2\}}\)。

设全集为 \(U\),那么集合幂级数可以表示为

\[f(x)=\sum_{S \subseteq U} f_Sx^S \nonumber \]

定义集合幂级数的加、减法为

\[f(x) \pm g(x)=\sum_{S \subseteq U} (f_S \pm g_S)x^S \nonumber \]

子集卷积

Luogu P6097 【模板】子集卷积

基本形式

给定两个长度为 \(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 \operatorname{and} j=0}\newline{i \operatorname{or} j=k} } } a_i b_j \nonumber \]

其中 \(\operatorname{and}\) 和 \(\operatorname{or}\) 为位运算。

尝试将它转变成我们熟悉的与、或卷积的形式:注意到若 \(i \operatorname{or} j=k\) 并且 \(|i|+|j|=|k|\)(其中 \(|x|\) 表示 \(\operatorname{popcount}(x)\),即 \(x\) 的二进制中 \(1\) 的个数),那么就自然满足 \(i \operatorname{and} j=k\),于是设

\[A_{i,j}=\begin{cases} a_j, &|j|=i\newline 0, &|j|\ne i \end{cases} \nonumber \]

同理定义 \(B_{i,j},C_{i,j}\),上式就可以改写为

\[C_{l,k}=\sum_{i=0}^{l}\left( \sum_{x \operatorname{or} y = k} A_{i,x} \times B_{l-i,y} \right) \nonumber \]

后面括号里面是或卷积的形式,可以 \(\mathcal{O}(n2^n)\) 算,对于每个 \(i \in [0,|k|]\) 都算一次,时间复杂度 \(\mathcal{O}(n^2 2^n)\)。

集合幂级数下的解释

在集合幂级数的意义下,相当与加一个变量 \(z^{|S|}\),两个集合幂级数的乘积为

\[\left( \sum_{S \subseteq U} f_Sx^Sz^{|S|} \right) \times \left( \sum_{S \subseteq U} g_Sx^Sz^{|S|} \right) \nonumber \]

最终提取出 \(x^{T}z^{|T|}\) 的系数就是答案。

code for Luogu P6097
#include <bits/stdc++.h>

using namespace std;

const int N = (1 << 20) + 10, MOD = 1e9 + 9;
inline int Plus(int a, int b) {return a + b >= MOD ? a + b - MOD : a + b; }
inline int Minus(int a, int b) {return a - b < 0 ? a - b + MOD : a - b; }
int n, a[N], b[N], c[N];
int A[21][N], B[21][N], C[21][N];

inline void FWT(int A[], int n, int type) {
    for(int h = 2; h <= n; h <<= 1) 
        for(int i = 0; i < n; i += h)
            for(int j = i; j < i + (h >> 1); j ++) {
                if(type == 1) A[j + (h >> 1)] = Plus(A[j + (h >> 1)], A[j]);
                else A[j + (h >> 1)] = Minus(A[j + (h >> 1)], A[j]);
            }
}

int main() {
    ios::sync_with_stdio(false), cin.tie(0);

    cin >> n;
    for(int i = 0; i < (1 << n); i ++)
        cin >> a[i], A[__builtin_popcount(i)][i] = a[i];
    for(int i = 0; i < (1 << n); i ++)
        cin >> b[i], B[__builtin_popcount(i)][i] = b[i];
    for(int i = 0; i <= n; i ++)
        FWT(A[i], 1 << n, 1), FWT(B[i], 1 << n, 1);
    
    for(int i = 0; i <= n; i ++) {
        for(int k = 0; k <= i; k ++)
            for(int j = 0; j < (1 << n); j ++)
                C[i][j] = Plus(C[i][j], 1ll * A[k][j] * B[i - k][j] % MOD);
        FWT(C[i], 1 << n, -1);
    }
    for(int i = 0; i < (1 << n); i ++)
        c[i] = C[__builtin_popcount(i)][i];
    for(int i = 0; i < (1 << n); i ++)
        cout << c[i] << ' ';
    cout << '\n';
    
    return 0;
}

注意看这个式子

\[C_{l,k}=\sum_{i=0}^{l}\left( \sum_{x \operatorname{or} y = k} A_{i,x} \times B_{l-i,y} \right) \nonumber \]

改写成

\[C_i=\sum_{i_1+i_2=i}A_{i_1} *_{\operatorname{or}} B_{i_2} \nonumber \]

一个卷积?是不是也可以做求逆、\(\exp\)、\(\ln\)!

定义集合幂级数 \(\exp f\) 为满足 \(\exp f =\sum_{i=0}^{+\infty}\frac{f^i}{i!}\) 的集合幂级数,\(\ln f\) 为满足 \(\exp (\ln f) = f\) 的集合幂级数。

实际上是当然可以的,对每个 \(A_i,B_i\) 都FWT一下,现在对每个 \(x\),\(C_{i, x}\) 是互相独立的,因此可以对每个 \(x \in [0,2^n-1]\) 求出这一维的答案,然后最终合并,问题转化为长度为 \(n\) 的序列的求逆、\(\ln\)、\(\exp\),当然可以上多项式全家桶做到 \(\mathcal{O}(2^n n \log n)\),但是这是完全没必要的,\(n\) 很小,也就 \(20\) 左右,多项式全家桶的大常数是得不偿失的,可以直接 \(\mathcal{O}(n^2)\) 暴力做这些操作,时间复杂度是 \(\mathcal{O}(n^2 2^n)\)。

再细说一下的话,求逆用递推 \(g_n=\frac{1-\sum_{j=0}^{i-1}g_jf_{i-j}}{f_0}\),\(\exp\) 用 \(g'=gf'\),写成系数形式就是 \(g_n=\frac{1}{n}\sum_{i=1}^{n}if_ig_{n-i}\)。\(\ln\) 可以利用它是 \(\exp\) 的逆变换得到 \(g_n=f_n-\frac{1}{n}\sum_{i=1}^{n-1}ig_if_{n-i}\)。

求逆
for(int i = 0; i <= n; i ++)
    FWT(f[i], 1 << n, 1);
for(int mask = 0; mask < (1 << n); mask ++) {
    g[0][mask] = 1;
    for(int i = 1; i <= n; i ++)
        for(int j = 1; j <= i; j ++)
            g[i][mask] = Minus(g[i][mask], 1ll * f[j][mask] * g[i - j][mask] % MOD);
}
for(int i = 0; i <= n; i ++)
    FWT(g[i], 1 << n, -1);
$\exp$
inline void exp(int f[], int g[], int n) {
    for(int i = 0; i <= n; i ++) g[i] = 0;
    g[0] = 1;
    for(int i = 1; i <= n; i ++) {
        for(int j = 1; j <= i; j ++)
            g[i] = Plus(g[i], 1ll * f[j] * j % MOD * g[i - j] % MOD);
        g[i] = 1ll * g[i] * inv[i] % MOD;   // inv[i] 是 i 的逆元
    }
}

// 主函数中
for(int i = 0; i <= n; i ++)
    FWT(f[i], 1 << n, 1);
for(int mask = 0; mask < (1 << n); mask ++) {
    static int F[N + 1], G[N + 1];
    for(int i = 0; i <= n; i ++) 
        F[i] = f[i][mask];
    exp(F, G, n);
    for(int i = 0; i <= n; i ++)
        g[i][mask] = G[i];
}
for(int i = 0; i <= n; i ++)
    FWT(g[i], 1 << n, -1);
$\ln$
inline void ln(int f[], int g[], int n) {
    for(int i = 0; i <= n; i ++) g[i] = 0;
    for(int i = 1; i <= n; i ++) {
        for(int j = 1; j < i; j ++)
            g[i] = Plus(g[i], 1ll * g[j] * j % MOD * f[i - j] % MOD);
        g[i] = Minus(f[i], 1ll * g[i] * inv[i] % MOD); // inv[i] 是 i 的逆元
    }
}

// 主函数中
for(int i = 0; i <= n; i ++)
    FWT(f[i], 1 << n, 1);
for(int mask = 0; mask < (1 << n); mask ++) {
    static int F[N + 1], G[N + 1];
    for(int i = 0; i <= n; i ++)
        F[i] = f[i][mask];
    ln(F, G, n);
    for(int i = 0; i <= n; i ++)
        g[i][mask] = G[i];
}
for(int i = 0; i <= n; i ++)
    FWT(g[i], 1 << n, -1);

例题

「WC2018」 州区划分

题目描述

小 S 现在拥有 \(n\) 座城市,第 \(i\) 座城市的人口为 \(w_i\),城市与城市之间可能有双向道路相连。

现在小 S 要将这 \(n\) 座城市划分成若干个州,每个州由至少一个城市组成,每个城市在恰好一个州内。

假设小 S 将这些城市划分成了 \(k\) 个州,设 \(V_i\) 是第 \(i\) 个州包含的所有城市组成的集合。定义一条道路是一个州的内部道路,当且仅当这条道路的两个端点城市都在这个州内。如果一个州内部存在一条起点终点相同,不经过任何不属于这个州的城市,且经过这个州的所有内部道路都恰好一次并且经过这个州的所有城市至少一次的路径(路径长度可以为 \(0\)),则称这个州是不合法的。

定义第 \(i\) 个州的满意度为:第 \(i\) 个州的人口在前 \(i\) 个州的人口中所占比例的 \(p\) 次幂,即:

\[\left(\dfrac{\sum _ {x \in V _ i} w _ x}{\sum _ {j = 1} ^ i \sum _ {x \in V _ j} w _ x}\right) ^ p \]

定义一个划分的满意度为所有州的满意度的乘积。

求所有合法的划分方案的满意度之和。

答案对 \(998244353\) 取模。
两个划分 \(\{V_1, V _ 2, \cdots, V_k\}\) 和 \(\{C_1, C _ 2, \cdots, C_s\}\) 是不同的,当且仅当 \(k \neq s\),或存在某个 \(1 \leq i \leq k\),使得 \(V_i \neq C_i\)。

保证对于所有数据有:\(0 \leq n \leq 21\), \(0 \leq m \leq \dfrac{n\times (n-1)}{2}\) , \(0 \leq p \leq 2\), \(1 \leq w_i \leq 100\)。

题解

注意到一个州不合法当且仅当它连通且存在欧拉回路,状压dp,设 \(f_S\) 表示集合 \(S\) 的答案,\(g_S\) 表示集合 \(S\) 内所有点的 \(w\) 和,就有转移

\[f_S=\sum_{T \subset S}f_{S-T}\left(\frac{g_T}{g_S}\right)^p=\frac{1}{g_S^p}\sum_{T \subset S} f_{S-T}g_T^p \nonumber \]

这是一个子集卷积的形式,按照 \(|S|\) 升序,依次算出 \(f\) 即可。

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

code
#include <bits/stdc++.h>

using namespace std;

const int N = 21, M = (1 << N), MOD = 998244353;
inline int Plus(int a, int b) {return a + b >= MOD ? a + b - MOD : a + b; }
inline int Minus(int a, int b) {return a - b < 0 ? a - b + MOD : a - b; }
inline int ksm(long long a, int b) {
    long long r = 1;
    for(; b; b >>= 1, a = a * a % MOD)
        if(b & 1) r = r * a % MOD;
    return r;
}

int n, m, p, w[N];
vector<int> G[N];
int f[N + 1][M], g[N + 1][M], ig[M];

bool vis[N];
void dfs(int u, int mask, int c) {
    vis[u] = c;
    for(auto v : G[u]) 
        if((mask >> v & 1) && (vis[v] != c))
            dfs(v, mask, c);
}
inline int check(int mask) {
    static int deg[N];
    bool flag = false;
    for(int i = 0; i < n; i ++) if(mask >> i & 1) {
        for(auto j : G[i]) if(mask >> j & 1) deg[i] ++;
        if(deg[i] & 1) flag = true;
    }
    for(int i = 0; i < n; i ++) if(mask >> i & 1) {
        for(auto j : G[i]) if(mask >> j & 1) deg[i] --;
    }
    for(int i = 0; i < n; i ++) 
        if(mask >> i & 1) {dfs(i, mask, 1); break; }
    int siz = 0;
    for(int i = 0; i < n; i ++) 
        if(mask >> i & 1) siz += vis[i];
    if(siz != __builtin_popcount(mask)) flag = true;
    for(int i = 0; i < n; i ++)
        if(mask >> i & 1) {dfs(i, mask, 0); break; }
    if(!flag) return 0;
    int sum = 0;
    for(int i = 0; i < n; i ++)
        if(mask >> i & 1) sum = Plus(sum, w[i]);
    return sum;
}

inline void FWT(int A[], int n, int type) {
    for(int h = 2; h <= n; h <<= 1) 
        for(int i = 0; i < n; i += h)
            for(int j = i; j < i + (h >> 1); j ++) {
                if(type == 1) A[j + (h >> 1)] = Plus(A[j + (h >> 1)], A[j]);
                else A[j + (h >> 1)] = Minus(A[j + (h >> 1)], A[j]);
            }
}

int main() {
    ios::sync_with_stdio(false), cin.tie(0);

    cin >> n >> m >> p;
    for(int i = 0; i < m; i ++) {
        int a, b; cin >> a >> b;
        a --, b --; G[a].emplace_back(b), G[b].emplace_back(a);
    }
    for(int i = 0; i < n; i ++) cin >> w[i];

    for(int mask = 1; mask < (1 << n); mask ++) {
        int x = check(mask);
        if(x != 0) {
            if(p == 0) g[__builtin_popcount(mask)][mask] = 1;
            else if(p == 1) g[__builtin_popcount(mask)][mask] = x;
            else g[__builtin_popcount(mask)][mask] = 1ll * x * x % MOD;
        }
        int sum = 0;
        for(int i = 0; i < n; i ++)
            if(mask >> i & 1) sum = Plus(sum, w[i]);
        ig[mask] = ksm(ksm(sum, MOD - 2), p);
    }
    for(int i = 0; i <= n; i ++)
        FWT(g[i], 1 << n, 1);

    f[0][0] = 1; FWT(f[0], 1 << n, 1);
    for(int i = 1; i <= n; i ++) {
        for(int j = 0; j < i; j ++)
            for(int k = 0; k < (1 << n); k ++)
                f[i][k] = Plus(f[i][k], 1ll * f[j][k] * g[i - j][k] % MOD);
        FWT(f[i], 1 << n, -1);
        for(int k = 0; k < (1 << n); k ++)
            f[i][k] = __builtin_popcount(k) == i ? 1ll * f[i][k] * ig[k] % MOD : 0;
        FWT(f[i], 1 << n, 1);
    }
    FWT(f[n], 1 << n, -1);
    cout << f[n][(1 << n) - 1] << '\n';

    return 0;
}

「CEOI2019」 Amusement Park

题目描述

有一个含 \(n\) 个点,\(m\) 条边的有向图,图无重边,无自环,两点之间不成环。

现在我们想改变一些边的方向,使得该有向图无环。

您需要求出,每一种改变方向后使得该有向图无环的方案的需改变边的数量之和 \(\bmod\ 998244353\) 之后的答案。

对于 \(100\%\) 的数据,保证 \(1\le n\le 18\),\(0\le m\le \frac{n\times (n-1)}{2}\),\(1\le a_i,b_i\le n\),\(a_i\not=b_i\),对于 \(i\not=j\),均有 \(a_i\not=a_j\) 或者 \(b_i\not=b_j\),无序数对 \(\{a_i,b_i\}\) 互不相同。

题解

主要到一个DAG的所有边反向后仍然是DAG,因此答案就是方案数乘 \(\frac{m}{2}\)。

设 \(f_S\) 表示集合 \(S\) 中点的方案数,转移时钦定一些度数为 \(0\) 的点后容斥出结果,具体的有

\[f_S=\sum_{\substack{ {T \subseteq S}\newline\newline{T\text{是独立集} } } }(-1)^{|T|-1}f_{S-T} \nonumber \]

于是就可以枚举子集做到 \(\mathcal{O}(3^n)\) 了,应该可以通过,但接下来我们将它优化到 \(\mathcal{O}(n^22^n)\)。

设 \(g_S=(-1)^{|S|-1}[S\text{是独立集}]\),那么 \(f=f \times g + 1\),这里的乘法是子集卷积。解得 \(f=\frac{1}{1-g}\),做一个子集卷积求逆即可,时间复杂度 \(\mathcal{O}(n^22^n)\)。

code
#include <bits/stdc++.h>

using namespace std;

const int N = 18, M = (1 << N), MOD = 998244353;
inline int Plus(int a, int b) {return a + b >= MOD ? a + b - MOD : a + b; }
inline int Minus(int a, int b) {return a - b < 0 ? a - b + MOD : a - b; }
inline int ksm(long long a, int b) {
    long long r = 1;
    for(; b; b >>= 1, a = a * a % MOD)
        if(b & 1) r = r * a % MOD;
    return r;
}
int n, m;
int flag[M];

inline void FWT(int A[], int n, int type) {
    for(int h = 2; h <= n; h <<= 1)
        for(int i = 0; i < n; i += h)
            for(int j = i; j < i + (h >> 1); j ++) {
                if(type == 1) A[j + (h >> 1)] = Plus(A[j + (h >> 1)], A[j]);
                else A[j + (h >> 1)] = Minus(A[j + (h >> 1)], A[j]);
            }
}

int f[N + 1][M], g[N + 1][M];

int main() {
    ios::sync_with_stdio(false), cin.tie(0);

    cin >> n >> m;
    for(int i = 0; i < m; i ++) {
        int a, b; cin >> a >> b; a--, b --;
        flag[(1 << a) | (1 << b)] = 1;
    }
    FWT(flag, 1 << n, 1);
    for(int mask = 0; mask < (1 << n); mask ++) {
        if(flag[mask]) continue;
        if(__builtin_popcount(mask) & 1) f[__builtin_popcount(mask)][mask] = 1;
        else f[__builtin_popcount(mask)][mask] = MOD - 1;
    }

    for(int i = 0; i <= n; i ++)
        FWT(f[i], 1 << n, 1);
    for(int mask = 0; mask < (1 << n); mask ++) {
        for(int i = 1; i <= n; i ++)
            f[i][mask] = Minus(0, f[i][mask]);
        g[0][mask] = f[0][mask] = 1;
        for(int i = 1; i <= n; i ++)
            for(int j = 1; j <= i; j ++)
                g[i][mask] = Minus(g[i][mask], 1ll * f[j][mask] * g[i - j][mask] % MOD);
    }
    FWT(g[n], 1 << n, -1);
    cout << 1ll * g[n][(1 << n) - 1] * m % MOD * ksm(2, MOD - 2) % MOD << '\n';

    return 0;
}

ABC236Ex Distinct Multiples

题目描述

给定两个正整数 \(N,M\) 和一个正整数序列 \(D\),询问满足条件的序列 \(A\) 的个数:

  • \(1\leq A_i\leq M(1\leq i\leq N)\)
  • \(A_i\neq A_j(1\leq i<j\leq N)\)
  • \(D_i|A_i\)

满足 \(2\leq N\leq 16\),\(1\leq M\leq 10^{18}\),\(1\leq D_i\leq M\)。

题解

AtCoder ABC Ex乱写

loj154. 集合划分计数

题目描述

给定一个集合 \(S = \{ {x_1}, {x_2}, \dots, {x_n} \}\) 和一个 \(S\) 上的集合族 \(\mathcal F = \{ {S_0}, {S_1}, \dots, {S_{m-1}} \}\)。

一个划分 \(\mathcal P\) 是 \(\mathcal F\) 的一个子族,满足 \(\mathcal P\) 中所有集合的并为 \(S\) ,任意两个集合不相交。

求大小不大于 \(k\) 的划分的数量 \(\text{mod } 998244353\)。

两个划分 \({\mathcal P_1}\), \({\mathcal P_2}\) 不同,当且仅当存在 \(i\) 使 \(S_i \in {\mathcal P_1} \land S_i \notin {\mathcal P_2}\) 或 \(S_i \notin {\mathcal P_1} \land S_i \in {\mathcal P_2}\)。\(S_i\) 和 \(S_j\) 不同当且仅当 \(i \neq j\)。

题解

这就是一个不完整的子集卷积 \(\exp\) 。记 \(f(x)\) 为所给集族 \(\mathcal{F}\),所求即为

\[\sum_{i=0}^{k}\frac{f^i(x)}{i!} \nonumber \]

这里的乘法是用子集卷积定义的,所以不用担心转移时集合之间有交的问题。

怎么算这个东西呢?记上面式子为 \(g\),求导得到

\[g'=\sum_{i=0}^{k}\frac{if^{i-1}f'}{i!}=f'\sum_{i=0}^{k-1}\frac{f^i}{i!} \nonumber \]

\[g'=f'(g-\frac{f^k}{k!}) \nonumber \]

用先 \(\ln\) 再乘 \(k\) 再 \(\exp\) 的方法容易求得 \(h=\frac{f^k}{k!}\),之后系数形式就可以写作

\[g_n=\frac{1}{n}\sum_{i=1}^{n}if_i\left(g_{n-i}-h_{n-i}\right) \nonumber \]

然后就可以直接求了。

code
#include <bits/stdc++.h>

using namespace std;

const int N = 21, MOD = 998244353;
inline int Plus(int a, int b) {return a + b >= MOD ? a + b - MOD : a + b; }
inline int Minus(int a, int b) {return a - b < 0 ? a - b + MOD : a - b; }
inline int ksm(long long a, int b) {
    long long r = 1;
    for(; b; b >>= 1, a = a * a % MOD)
        if(b & 1) r = r * a % MOD;
    return r;
}
int n, m, k, f[N + 1][1 << N], g[N + 1][1 << N];
int inv[N + 10], fac[N + 1], ifac[N + 1];

inline void FWT(int A[], int n, int type) {
    for(int h = 2; h <= n; h <<= 1)
        for(int i = 0; i < n; i += h)
            for(int j = i; j < i + (h >> 1); j ++) {
                if(type == 1) A[j + (h >> 1)] = Plus(A[j + (h >> 1)], A[j]);
                else A[j + (h >> 1)] = Minus(A[j + (h >> 1)], A[j]);
            }
}

inline void exp(int f[], int g[], int n) {
    for(int i = 0; i <= n; i ++) g[i] = 0;
    g[0] = 1;
    for(int i = 1; i <= n; i ++) {
        for(int j = 1; j <= i; j ++)
            g[i] = Plus(g[i], 1ll * f[j] * j % MOD * g[i - j] % MOD);
        g[i] = 1ll * g[i] * inv[i] % MOD;   // inv[i] 是 i 的逆元
    }
}
inline void ln(int f[], int g[], int n) {
    for(int i = 0; i <= n; i ++) g[i] = 0;
    for(int i = 1; i <= n; i ++) {
        for(int j = 1; j < i; j ++)
            g[i] = Plus(g[i], 1ll * g[j] * j % MOD * f[i - j] % MOD);
        g[i] = Minus(f[i], 1ll * g[i] * inv[i] % MOD); // inv[i] 是 i 的逆元
    }
}

inline void power(int f[], int n, int k) {
    static int g[N];
    for(int i = 0; i <= n; i ++) g[i] = 0;
    int pos = 0;
    while(pos <= n && f[pos] == 0) pos ++;
    if(pos > n || pos * k > n) {
        for(int i = 0; i <= n; i ++) f[i] = 0;
        return;
    }
    int d = f[pos], id = ksm(d, MOD - 2);
    for(int i = 0; i <= n; i ++)
        if(i + pos <= n) f[i] = 1ll * f[i + pos] * id % MOD;
        else f[i] = 0;
    ln(f, g, n - pos);
    for(int i = 0; i <= n - pos; i ++)
        f[i] = 1ll * g[i] * k % MOD;
    exp(f, g, n - pos);
    for(int i = n, mul = ksm(d, k); i >= 0; i --)
        if(i - pos * k >= 0) f[i] = 1ll * g[i - pos * k] * mul % MOD;
        else f[i] = 0;
}

int main() {
    ios::sync_with_stdio(false), cin.tie(0);

    cin >> n >> m >> k;
    inv[1] = 1;
    for(int i = 2; i <= n; i ++)
        inv[i] = Minus(0, 1ll * (MOD / i) * inv[MOD % i] % MOD);
    fac[0] = 1; for(int i = 1; i <= n; i ++) fac[i] = 1ll * fac[i - 1] * i % MOD;
    ifac[n] = ksm(fac[n], MOD - 2); for(int i = n; i >= 1; i --) ifac[i - 1] = 1ll * ifac[i] * i % MOD;
    while(m --) {
        int x; cin >> x;
        f[__builtin_popcount(x)][x] ++;
    }

    const int Lim = (1 << n);
    for(int i = 0; i <= n; i ++)
        FWT(f[i], 1 << n, 1);
    for(int mask = 0; mask < Lim; mask ++) {
        static int F[N + 1], G[N + 1], H[N + 1];
        for(int i = 0; i <= n; i ++)
            H[i] = F[i] = f[i][mask], G[i] = 0;
        power(H, n, k);
        for(int i = 0; i <= n; i ++)
            H[i] = 1ll * H[i] * ifac[k] % MOD;
        for(int i = 0, mul = 1; i <= k; i ++) {
            G[0] = Plus(G[0], 1ll * mul * ifac[i] % MOD);
            mul = 1ll * mul * F[0] % MOD;
        }
        for(int i = 1; i <= n; i ++) {
            for(int j = 1; j <= i; j ++)
                G[i] = Plus(G[i], 1ll * j * F[j] % MOD * Minus(G[i - j], H[i - j]) % MOD);
            G[i] = 1ll * G[i] * inv[i] % MOD;
        }
        for(int i = 0; i <= n; i ++) g[i][mask] = G[i];
    }
    FWT(g[n], Lim, - 1);
    cout << g[n][Lim - 1] << '\n';

    return 0;
}

标签:nonumber,卷积,mathscr,sum,int,子集,沃尔什,operatorname,MOD
From: https://www.cnblogs.com/313la/p/18010434

相关文章

  • 代码随想录 day42 背包问题 分割等和子集
    01背包leetcode没有原题这里是解法importjava.util.Arrays;publicclassBagProblem{publicstaticvoidmain(String[]args){int[]weight={1,3,4};int[]value={15,20,30};intbagSize=4;testWeightBagProblem(weight,value,bagSize);}/***初始化dp......
  • 机器学习中一维卷积的作用是什么
    一维卷积在机器学习中的应用特别适合处理时间序列数据或者是一维信号数据。其作用主要体现在以下几个方面:特征提取:一维卷积通过在数据上滑动一个较小的窗口(卷积核),并计算窗口内数据的加权和(可能还包括偏置项),从而在局部区域内提取特征。这种操作有助于识别一维数据中的局部模式和特征......
  • 深度学习网络的感受野与卷积核
    https://www.bilibili.com/read/cv27451493/?jump_opus=1https://zhuanlan.zhihu.com/p/484653541?utm_id=0一般认为,网络越深,卷积核越大,感受野也就越大。同时,也会丢失一定的小尺度捕捉能力。在《Residualnetworksbehavelikeensemblesofrelativelyshallownetworks》中,说......
  • 卷积神经网络理解(5)
    1、参数计算         在一次卷积过程中,卷积核进行共享,即每个通道采用一个卷积核即可。 在Pytorch的nn模块中,封装了nn.Conv2d()类作为二维卷积的实现。参数如下图所示    in_channels:输入张量的channels数out_channels:输出张量的channels数kernel_s......
  • 如何优雅的处理特殊的子集 dp 问题
    sosdp&高维前缀和求\[g_i=\sum_{j\&i>0}f_j(i\leq2^n-1)\]我们将\(i,j\)进行二进制拆分,拆成\(n\)个维度。类似于:\[g_{a_1,a_2,a_3,a_4,a_5...a_n}=\sum_{a_k\leqb_k}f_{b_1,b_2,b_3,b_4,b_5...b_n}(a_i,b_i\subseteq\{0,1\}......
  • 积性函数与狄利克雷卷积
    积性函数【定义】若对于一个数论函数\(f\),有:对\((a,b)=1\),有\(f(a\timesb)=f(a)\timesf(b)\),称\(f\)是一个积性函数。特别地,若对于任意数\(a,b\),有\(f(a\timesb)=f(a)\timesf(b)\),称\(f\)是一个完全积性函数。积性函数举例:欧拉函数,因数个数,因数和……完......
  • 卷积神经网络理解(4)
    1、CNN中常见的名词padding:padding(填充)参数的作用是决定在进行卷积或池化操作时,是否对输入的图像矩阵边缘补0stride:滑动卷积核时的步长stride(例如每次滑动一个或两个)kernal:卷积核,通常为3x3或者5x5filter:卷积核的数量(神经元的数量)。这个地方怎么理解呢,一个3x3的卷积核有9个参......
  • 卷积神经网络理解(3)
    1、定义LeNet是深度学习领域的一个经典卷积神经网络模型,由YannLeCun等人于1998年提出,被广泛应用于手写数字识别和其他图像识别任务。LeNet的网络结构相对简单,包含两个卷积层和三个全连接层,是卷积神经网络的基础。LeNet对于现代的图像识别任务来说可能过于简单,但其对于深度学习......
  • 卷积神经网络理解(二)
    1、卷积神经网络的特点卷积神经网络相对于普通神经网络在于以下四个特点:局部感知域:CNN的神经元只与输入数据的一小部分区域相连接,这使得CNN对数据的局部结构具有强大的敏感性,可以自动学习到图像的特征。参数共享:在CNN中,同一个卷积核(filter)在整个输入图像上滑动,共享权重和偏置......
  • Pdfium.Net.Free 一个免费的Pdfium的 .net包装器--创建字符子集
    项目地址:Pdfium.Net:https://github.com/1000374/Pdfium.NetPdfiumViewer:https://github.com/1000374/PdfiumViewerPdfium.Net.Free一个免费的Pdfium的.net包装器--加载字体 接上篇,怎么创建字符子集呢?获取字符集内的字形符号需要引用wpf下PresentationCore.dll,根据比对传入......