一些奇怪的记号
设 A , B A,B A,B 为两个长度为 n n n 的多项式,在下文中我们做如下约定:
- 当对一个数进行集合操作时,我们认为这个数代表以它二进制表示中为 1 1 1 的位置的下标为元素的集合
- [ x i ] A [x^i]A [xi]A 表示 A A A 的 i i i 次项系数
- A ± B A \pm B A±B 表示 A A A 与 B B B 按位相加减
- A ⋅ B A \cdot B A⋅B 表示 A A A 与 B B B 按位相乘
-
A
∗
∘
B
A *_\circ B
A∗∘B 表示
A
A
A 与
B
B
B 按
∘
\circ
∘ 运算卷积,即若多项式
C
=
A
∗
∘
B
C=A*_\circ B
C=A∗∘B,则
C i = ∑ j ∘ k = i A j B k C_i=\sum_{j\circ k=i}A_jB_k Ci=j∘k=i∑AjBk -
(
A
,
B
)
(A,B)
(A,B) 表示将两个多项式"拼"起来,即若多项式
C
=
(
A
,
B
)
C=(A,B)
C=(A,B),则
C i = { A i i < n B i − n i ⩾ n C_i= \left\{ \begin{array}{lr} A_i &i< n\\ B_{i - n} &i\geqslant n \end{array} \right. Ci={AiBi−ni<ni⩾n
F W T \rm FWT FWT
考虑这样一个式子:
C
k
=
∑
i
∘
j
=
k
A
i
B
j
C_k=\sum_{i \circ j=k}A_iB_j
Ck=i∘j=k∑AiBj
其中 ∘ \circ ∘ 表示一种运算。如果令 ∘ \circ ∘ 表示加法,那么这就是我们熟悉的卷积,用FFT就可以快速计算了。但当 ∘ \circ ∘ 表示位运算的时候,我们要怎么做呢?
回想加法卷积时,我们先是搞了一个变换
D
F
T
{\rm DFT}
DFT 使它满足
D
F
T
(
A
∗
+
B
)
=
D
F
T
(
A
)
⋅
D
F
T
(
B
)
{\rm DFT}(A *_+ B)={\rm DFT(A)\cdot DFT(B)}
DFT(A∗+B)=DFT(A)⋅DFT(B),并且还有能快速计算的逆变换
I
D
F
T
{\rm IDFT}
IDFT。
那么我们能不能找一些新的变换使它们能做位运算呢?
当然可以,这就是 F W T \rm FWT FWT,即 F a s t W a l s h \rm Fast \;Walsh FastWalsh − - − H a d a m a r d T r a n s f o r m \rm Hadamard \;Transform HadamardTransform.
o r \rm or or 卷积
构造
让我们先来看一下
∘
\circ
∘ 是
o
r
\rm or
or 运算的情况:
C
k
=
∑
i
o
r
j
=
k
A
i
B
j
C_k=\sum_{i\,{\rm or}\,j=k}A_iB_j
Ck=iorj=k∑AiBj
如果我们将下标看成它们二进制表示下所代表的集合,容易发现
i
o
r
j
=
k
⇒
i
⊆
k
,
j
⊆
k
i\,{\rm or}\,j=k\Rightarrow i \subseteq k,j\subseteq k
iorj=k⇒i⊆k,j⊆k,而我们希望变换后的两边求和指标是等价的且形式相同,不难发现我们把左边的等于改成包含就能满足这一点,即
i
o
r
j
⊆
k
⇔
i
⊆
k
,
j
⊆
k
i\,{\rm or}\,j\subseteq k\Leftrightarrow i \subseteq k,j\subseteq k
iorj⊆k⇔i⊆k,j⊆k.那么
o
r
\rm or
or 运算时的变换已经非常明显了,我们构造变换
T
{\rm T}
T:
T
(
A
)
i
=
∑
j
⊆
i
A
j
{\rm T}(A)_i=\sum_{j\subseteq i}A_j
T(A)i=j⊆i∑Aj
于是有
T
(
A
∗
o
r
B
)
=
T
(
A
)
⋅
T
(
B
)
{\rm T}(A*_{\rm or}B)={\rm T}(A)\cdot {\rm T}(B)
T(A∗orB)=T(A)⋅T(B)
证明非常简单,其实就是上面的式子
T
(
A
∗
o
r
B
)
i
=
∑
j
⊆
i
∑
j
=
a
o
r
b
A
a
B
b
=
∑
(
a
o
r
b
)
⊆
i
A
a
B
b
=
(
∑
a
⊆
i
A
a
)
(
∑
b
⊆
i
B
b
)
=
T
(
A
)
i
T
(
B
)
i
{\rm T}(A*_{\rm or}B)_i=\sum_{j\subseteq i}\sum_{j=a\,{\rm or}\,b}A_aB_b=\sum_{(a\,{\rm or}\,b)\subseteq i}A_aB_b=\left(\sum_{a\subseteq i}A_a\right)\left(\sum_{b\subseteq i}B_b\right)={\rm T}(A)_i{\rm T}(B)_i
T(A∗orB)i=j⊆i∑j=aorb∑AaBb=(aorb)⊆i∑AaBb=(a⊆i∑Aa)(b⊆i∑Bb)=T(A)iT(B)i
光有
T
{\rm T}
T 不行,我们最后还要将式子变回去,逆变换
T
−
1
{\rm T}^{-1}
T−1 怎么办?
显然,刚才是将所有子集求和,现在只要容斥就行了,即
T
−
1
(
A
)
i
=
∑
j
⊆
i
(
−
1
)
∣
i
∖
j
∣
A
j
{\rm T}^{-1}(A)_i=\sum_{j\subseteq i}(-1)^{|i\setminus j|}A_j
T−1(A)i=j⊆i∑(−1)∣i∖j∣Aj
于是我们就完成了第一步:构造变换
计算
完成构造只是第一步,我们还需要能够快速计算这个变换,否则一切都是白搭。
众所周知位运算具有位独立性,于是我们可以考虑按位分治。
设我们当前要变换的数列为
A
A
A,由于要按位分治,我们先把
A
A
A 的长度补成
2
2
2 的幂,设为
2
k
2^k
2k。再设
A
0
A_0
A0 为
A
A
A 下标最高位为
0
0
0 的部分,
A
1
A_1
A1 为
A
A
A 下标最高位为
1
1
1 的部分(即
A
0
,
A
1
A_0,A_1
A0,A1 分别为
A
A
A 的左半部分与右半部分),考虑如何由
T
(
A
0
)
{\rm T}(A_0)
T(A0) 和
T
(
A
1
)
{\rm T}(A_1)
T(A1) 推出
T
(
A
)
{\rm T}(A)
T(A):
- 对于 i < 2 k − 1 i< 2^{k - 1} i<2k−1,显然右边无法产生贡献(第 k k k 位为 1 1 1),于是有 T ( A ) i {\rm T}(A)_i T(A)i 就是 T ( A 0 ) {\rm T}(A_0) T(A0) 对应位的值
- 而对于 i ⩾ 2 k − 1 i\geqslant 2^{k - 1} i⩾2k−1, T ( A 0 ) , T ( A 1 ) {\rm T}(A_0),{\rm T}(A_1) T(A0),T(A1) 的对应位显然都会第 i i i 位产生贡献,于是 T ( A ) i {\rm T}(A)_i T(A)i 就是两个式子的对应位相加
也就是
T
(
A
)
=
(
T
(
A
0
)
,
T
(
A
0
)
+
T
(
A
1
)
)
{\rm T}(A)=({\rm T}(A_0),{\rm T}(A_0)+{\rm T}(A_1))
T(A)=(T(A0),T(A0)+T(A1))
边界显然,当
k
=
0
k=0
k=0 时,
T
(
A
)
=
A
{\rm T}(A)=A
T(A)=A.
那么我们就可以分治计算变换
T
{\rm T}
T 了,对于
T
−
1
{\rm T}^{-1}
T−1 的计算同理,式子如下:
T
−
1
(
A
)
=
(
T
−
1
(
A
0
)
,
T
−
1
(
A
1
)
−
T
−
1
(
A
0
)
)
{\rm T}^{-1}(A)=({\rm T}^{-1}(A_0),{\rm T}^{-1}(A_1)-{\rm T}^{-1}(A_0))
T−1(A)=(T−1(A0),T−1(A1)−T−1(A0))
代码
/* n := 2^n,下同 */
void FWTor(int *f,int n,int opt){
for(int i = 1; i < n; i <<= 1){
for(int j = 0; j < n; j += i << 1){
for(int k = j; k < j + i; k ++){
if(opt == 1) f[k + i] = add(f[k + i],f[k]);
else f[k + i] = dec(f[k + i],f[k]);
}
}
}
}
a n d \rm and and 卷积
构造
当 ∘ \circ ∘ 是 a n d \rm and and 运算的时候呢?
我们考虑在
∘
\circ
∘ 是
o
r
\rm or
or 时为什么能那么构造,这是因为
i
o
r
j
=
k
⇒
i
⊆
k
,
j
⊆
k
i\,{\rm or}\,j=k\Rightarrow i \subseteq k,j\subseteq k
iorj=k⇒i⊆k,j⊆k.
而我们不难发现
a
n
d
\rm and
and 也具有一个很相似的性质:
i
a
n
d
j
=
k
⇒
i
⊇
k
,
j
⊇
k
i\,{\rm and}\,j=k\Rightarrow i \supseteq k,j\supseteq k
iandj=k⇒i⊇k,j⊇k.
那么我们可以直接仿照上面构造类似的
T
,
T
−
1
{\rm T},{\rm T}^{-1}
T,T−1 与计算方法,这里就不证明了
T ( A ) i = ∑ j ⊇ i A j {\rm T}(A)_i=\sum_{j\supseteq i}A_j T(A)i=j⊇i∑Aj
T − 1 ( A ) i = ∑ j ⊇ i ( − 1 ) ∣ j ∖ i ∣ A j {\rm T}^{-1}(A)_i=\sum_{j\supseteq i}(-1)^{|j\setminus i|}A_j T−1(A)i=j⊇i∑(−1)∣j∖i∣Aj
计算
计算也与
o
r
\rm or
or 类似
T
(
A
)
=
(
T
(
A
0
)
+
T
(
A
1
)
,
T
(
A
1
)
)
{\rm T}(A)=({\rm T}(A_0) + {\rm T}(A_1),{\rm T}(A_1))
T(A)=(T(A0)+T(A1),T(A1))
T − 1 ( A ) = ( T − 1 ( A 0 ) − T − 1 ( A 1 ) , T − 1 ( A 0 ) ) {\rm T}^{-1}(A)=({\rm T}^{-1}(A_0)-{\rm T}^{-1}(A_1),{\rm T}^{-1}(A_0)) T−1(A)=(T−1(A0)−T−1(A1),T−1(A0))
代码
void FWTand(int *f,int n,int opt){
for(int i = 1; i < n; i <<= 1){
for(int j = 0; j < n; j += i << 1){
for(int k = j; k < j + i; k ++){
if(opt == 1) f[k] = add(f[k],f[k + i]);
else f[k] = dec(f[k],f[k + i]);
}
}
}
}
x o r \rm xor xor 卷积
构造
x o r \rm xor xor 可就没有像上面两种运算一样优秀的性质了,我们必须构造一个新的变换 T {\rm T} T.
但是这个变换实在是太过诡异,我也不知道是如何构造出来的…
令
d
(
i
)
d(i)
d(i) 表示
i
i
i 的二进制表示下
1
1
1 的个数,有如下引理:
d
(
k
a
n
d
i
)
x
o
r
d
(
k
a
n
d
j
)
=
d
(
k
a
n
d
(
i
x
o
r
j
)
)
d(k\, {\rm and} \,i) \,{\rm xor}\,d(k \,{\rm and}\,j)=d(k\,{\rm and}\,(i\,{\rm xor}\,j))
d(kandi)xord(kandj)=d(kand(ixorj))
变换
T
{\rm T}
T 为
T
(
A
)
i
=
∑
d
(
i
a
n
d
j
)
=
0
A
j
−
∑
d
(
i
a
n
d
j
)
=
1
A
j
{\rm T}(A)_i=\sum_{d(i \,{\rm and}\, j)=0}A_j-\sum_{d(i \,{\rm and}\, j)=1}A_j
T(A)i=d(iandj)=0∑Aj−d(iandj)=1∑Aj
下面证明
T
(
A
∗
x
o
r
B
)
=
T
(
A
)
⋅
T
(
B
)
{\rm T}(A*_{\rm xor}B)={\rm T}(A)\cdot {\rm T}(B)
T(A∗xorB)=T(A)⋅T(B):
T
(
A
)
i
⋅
T
(
B
)
i
=
(
∑
d
(
i
a
n
d
j
)
=
0
A
j
−
∑
d
(
i
a
n
d
j
)
=
1
A
j
)
(
∑
d
(
i
a
n
d
j
)
=
0
B
j
−
∑
d
(
i
a
n
d
j
)
=
1
B
j
)
=
(
∑
d
(
i
a
n
d
j
)
=
0
A
j
)
(
∑
d
(
i
a
n
d
k
)
=
0
B
k
)
+
(
∑
d
(
i
a
n
d
j
)
=
1
A
j
)
(
∑
d
(
i
a
n
d
k
)
=
1
B
k
)
−
(
∑
d
(
i
a
n
d
j
)
=
0
A
j
)
(
∑
d
(
i
a
n
d
k
)
=
1
B
k
)
−
(
∑
d
(
i
a
n
d
j
)
=
1
A
j
)
(
∑
d
(
i
a
n
d
k
)
=
0
B
k
)
=
(
∑
d
(
i
a
n
d
j
)
x
o
r
d
(
i
a
n
d
k
)
=
0
A
j
B
k
)
−
(
∑
d
(
i
a
n
d
j
)
x
o
r
d
(
i
a
n
d
k
)
=
1
A
j
B
k
)
=
(
∑
d
(
i
a
n
d
(
j
x
o
r
k
)
)
=
0
A
j
B
k
)
−
(
∑
d
(
i
a
n
d
(
j
x
o
r
k
)
)
=
1
A
j
B
k
)
=
T
(
A
∗
x
o
r
B
)
i
\begin{aligned} {\rm T}(A)_i\cdot {\rm T}(B)_i&=\left(\sum_{d(i \,{\rm and}\, j)=0}A_j-\sum_{d(i \,{\rm and}\, j)=1}A_j\right)\left(\sum_{d(i \,{\rm and}\, j)=0}B_j-\sum_{d(i \,{\rm and}\, j)=1}B_j\right)\\ &=\left(\sum_{d(i \,{\rm and}\, j)=0}A_j\right)\left(\sum_{d(i \,{\rm and}\, k)=0}B_k\right)+\left(\sum_{d(i \,{\rm and}\, j)=1}A_j\right)\left(\sum_{d(i \,{\rm and}\,k)=1}B_k\right)\\ &\;\;\;\;\;\;\;-\left(\sum_{d(i \,{\rm and}\, j)=0}A_j\right)\left(\sum_{d(i \,{\rm and}\, k)=1}B_k\right)-\left(\sum_{d(i \,{\rm and}\, j)=1}A_j\right)\left(\sum_{d(i \,{\rm and}\, k)=0}B_k\right)\\ &=\left(\sum_{d(i\, {\rm and} \,j) \,{\rm xor}\,d(i \,{\rm and}\,k)=0}A_jB_k\right)-\left(\sum_{d(i\, {\rm and} \,j) \,{\rm xor}\,d(i \,{\rm and}\,k)=1}A_jB_k\right)\\ &=\left(\sum_{d(i\,{\rm and}\,(j\,{\rm xor}\,k))=0}A_jB_k\right)-\left(\sum_{d(i\,{\rm and}\,(j\,{\rm xor}\,k))=1}A_jB_k\right)\\ &={\rm T}(A*_{\rm xor}B)_i \end{aligned}
T(A)i⋅T(B)i=⎝⎛d(iandj)=0∑Aj−d(iandj)=1∑Aj⎠⎞⎝⎛d(iandj)=0∑Bj−d(iandj)=1∑Bj⎠⎞=⎝⎛d(iandj)=0∑Aj⎠⎞⎝⎛d(iandk)=0∑Bk⎠⎞+⎝⎛d(iandj)=1∑Aj⎠⎞⎝⎛d(iandk)=1∑Bk⎠⎞−⎝⎛d(iandj)=0∑Aj⎠⎞⎝⎛d(iandk)=1∑Bk⎠⎞−⎝⎛d(iandj)=1∑Aj⎠⎞⎝⎛d(iandk)=0∑Bk⎠⎞=⎝⎛d(iandj)xord(iandk)=0∑AjBk⎠⎞−⎝⎛d(iandj)xord(iandk)=1∑AjBk⎠⎞=⎝⎛d(iand(jxork))=0∑AjBk⎠⎞−⎝⎛d(iand(jxork))=1∑AjBk⎠⎞=T(A∗xorB)i
逆变换为
T
−
1
(
A
)
=
1
2
n
T
(
A
)
{\rm T}^{-1}(A)=\frac{1}{2^n}{\rm T}(A)
T−1(A)=2n1T(A)
计算
依然是原来的套路,按位分治,仍然设出
A
0
,
A
1
A_0,A_1
A0,A1。
但这次显然两边都互相有贡献,于是我们略微修改一下,递归求解时不考虑最高位的值。
先递归下去,由于不考虑最高位,边界为
T
(
A
)
=
A
{\rm T_{}}(A)=A
T(A)=A,重点是如何转移。
根据位独立性,不考虑最高位时左右两区间内对应位置的求和指标的对应位置显然也是一样的,故只有最高位会影响贡献。
考虑最高位,不难发现,由于
a
n
d
\rm and
and 运算只有在两边都是
1
1
1的情况下结果为
1
1
1,所以只有
T
(
A
1
)
{\rm T}(A_1)
T(A1) 对
T
(
A
)
{\rm T}(A)
T(A) 右半部分贡献时会使奇偶性翻转,其他奇偶性均不变,故只有这一项是负的,其他都是正的,即
T
(
A
)
=
(
T
(
A
0
)
+
T
(
A
1
)
,
T
(
A
0
)
−
T
(
A
1
)
)
{\rm T}(A)=({\rm T}(A_0)+{\rm T}(A_1),{\rm T}(A_0)-{\rm T}(A_1))
T(A)=(T(A0)+T(A1),T(A0)−T(A1))
而对于逆变换 T − 1 {\rm T}^{-1} T−1,根据其表达式,可以先做一次正变换,再除以 2 n 2^n 2n。
于是我们就莫名其妙地做完了
x
o
r
\rm xor
xor 卷积
代码
void FWTxor(int *f,int n,int opt){
for(int i = 1; i < n; i <<= 1){
for(int j = 0; j < n; j += i << 1){
for(int k = j; k < j + i; k ++){
int x = f[k],y = f[k + i];
f[k] = add(x,y),f[k + i] = dec(x,y);
}
}
}
if(opt == -1){
int t = qpow(n,mod - 2);
for(int i = 0; i < n; i ++) f[i] = 1ll * f[i] * t % mod;
}
}
F W T \rm FWT FWT的应用 ——— 子集卷积
构造
我们考虑这个式子:
C
i
=
∑
j
∪
k
=
i
,
j
∩
k
=
∅
A
j
B
k
C_i=\sum_{j\cup k=i,j\cap k=\varnothing}A_jB_k
Ci=j∪k=i,j∩k=∅∑AjBk
发现
j
∩
k
=
∅
j\cap k=\varnothing
j∩k=∅ 比较难满足,而
j
∪
k
=
i
j\cup k=i
j∪k=i 就是或卷积,我们将它改写一下
C
i
=
∑
j
o
r
k
=
i
,
∣
j
∣
+
∣
k
∣
=
∣
i
∣
A
j
B
k
C_i=\sum_{j{\rm or}k=i,|j|+|k|=|i|}A_jB_k
Ci=jork=i,∣j∣+∣k∣=∣i∣∑AjBk
于是我们可以加一维表示集合大小,即设
[
x
i
]
A
k
′
=
A
i
[
∣
i
∣
=
k
]
[x^i]A'_{k}=A_i[|i|=k]
[xi]Ak′=Ai[∣i∣=k]
进而有
C
i
=
∑
j
o
r
k
=
i
,
∣
j
∣
+
∣
k
∣
=
∣
i
∣
A
j
B
k
=
[
x
i
]
I
F
W
T
o
r
(
∑
j
o
r
k
⊆
i
,
∣
j
∣
+
∣
k
∣
=
∣
i
∣
A
j
B
k
)
=
[
x
i
]
I
F
W
T
o
r
(
∑
p
+
q
=
∣
i
∣
(
∑
j
⊆
i
,
∣
j
∣
=
p
A
j
)
(
∑
k
⊆
i
,
∣
k
∣
=
q
A
j
)
)
=
[
x
i
]
I
F
W
T
o
r
(
∑
p
+
q
=
∣
i
∣
F
W
T
o
r
(
A
p
′
)
⋅
F
W
T
o
r
(
B
q
′
)
)
\begin{aligned} C_i&=\sum_{j\,{\rm or}\,k=i,|j|+|k|=|i|}A_jB_k\\ &=[x^i]{\rm IFWT_{or}}\left(\sum_{j\,{\rm or}\,k\subseteq i,|j|+|k|=|i|}A_jB_k\right)\\ &=[x^i]{\rm IFWT_{or}}\left(\sum_{p+q=|i|}\left(\sum_{j\subseteq i,|j|=p}A_j\right)\left(\sum_{k\subseteq i,|k|=q}A_j\right)\right)\\ &=[x^i]{\rm IFWT_{or}}\left(\sum_{p+q=|i|}{\rm FWT_{or}}\left(A'_p\right)\cdot{\rm FWT_{or}}\left(B'_q\right)\right) \end{aligned}
Ci=jork=i,∣j∣+∣k∣=∣i∣∑AjBk=[xi]IFWTor⎝⎛jork⊆i,∣j∣+∣k∣=∣i∣∑AjBk⎠⎞=[xi]IFWTor⎝⎛p+q=∣i∣∑⎝⎛j⊆i,∣j∣=p∑Aj⎠⎞⎝⎛k⊆i,∣k∣=q∑Aj⎠⎞⎠⎞=[xi]IFWTor⎝⎛p+q=∣i∣∑FWTor(Ap′)⋅FWTor(Bq′)⎠⎞
一共要做 3 n 3n 3n 遍 F W T {\rm FWT} FWT, 2 n 2^n 2n 遍卷积,总时间复杂度为 O ( 4 n 2 2 n ) \Omicron(4n^22^n) O(4n22n).
要注意虽然上面的和式是卷积形式,但 n n n 的规模很小(否则无法承受如此高的复杂度),直接暴力做更快
代码
下面是洛谷模板题的代码
#include <iostream>
#include <cstdio>
using namespace std;
const int maxn = 21,maxm = 1 << 20,mod = 1e9 + 9;
int n,m,cnt[maxm],a[maxn][maxm],b[maxn][maxm],c[maxn][maxm];
int read(){
int x = 0;
char c = getchar();
while(c < '0' || c > '9') c = getchar();
while(c >= '0' && c <= '9') x = x * 10 + (c ^ 48),c = getchar();
return x;
}
inline int add(int x,int y){
if(x + y < mod) return x + y;
else return x + y - mod;
}
inline int dec(int x,int y){
if(x - y >= 0) return x - y;
else return x - y + mod;
}
void FWT(int *f,int n,int opt){
for(int i = 1; i < n; i <<= 1){
for(int j = 0; j < n; j += i << 1){
for(int k = j; k < j + i; k ++){
if(opt == 1) f[k + i] = add(f[k + i],f[k]);
else f[k + i] = dec(f[k + i],f[k]);
}
}
}
}
int main(){
n = read(),m = 1 << n;
for(int i = 0; i < m; i ++) cnt[i] = cnt[i >> 1] + (i & 1);
for(int i = 0; i < m; i ++) a[cnt[i]][i] = read();
for(int i = 0; i < m; i ++) b[cnt[i]][i] = read();
for(int i = 0; i <= n; i ++) FWT(a[i],m,1);
for(int i = 0; i <= n; i ++) FWT(b[i],m,1);
for(int i = 0; i <= n; i ++){
for(int j = 0; i + j <= n; j ++){
for(int k = 0; k < m; k ++){
c[i + j][k] = add(c[i + j][k],1ll * a[i][k] * b[j][k] % mod);
}
}
}
for(int i = 0; i <= n; i ++) FWT(c[i],m,-1);
for(int i = 0; i < m; i ++) printf("%d ",c[cnt[i]][i]);
return 0;
}
例题
先补充一个 F W T \rm FWT FWT 的重要性质,题目中会用到。
任意多项式 A 1 , A 1 , ⋯ A k A_1,A_1,\cdots A_k A1,A1,⋯Ak 均满足:
F W T ( A 1 + A 1 + ⋯ + A k ) = F W T ( A 1 ) + F W T ( A 2 ) + ⋯ + F W T ( A k ) {\rm FWT}(A_1+A_1+\cdots+A_k)={\rm FWT}(A_1)+{\rm FWT}(A_2)+\cdots +{\rm FWT}(A_k) FWT(A1+A1+⋯+Ak)=FWT(A1)+FWT(A2)+⋯+FWT(Ak)其中左式括号内表示多项式按位加,即 “和的 F W T \rm FWT FWT 等于 F W T \rm FWT FWT 的和”
黎明前的巧克力
分析
容易想到两个人最终的异或和相同 ⇔ \Leftrightarrow ⇔ 两人的选数集合异或和为 0 0 0
那么每一个原集合的大小为
k
k
k 的异或和为
0
0
0 的子集都会对答案贡献
2
k
2^k
2k(任意划分成两个集合),于是我们定义一个大小为
k
k
k 的集合的价值为
2
k
2^k
2k.
容易想到一个DP:设
f
i
,
j
f_{i,j}
fi,j 表示由前
i
i
i 个数的子集中异或和为
j
j
j 的子集价值和,转移方程为:
f
i
,
j
=
f
i
−
1
,
j
+
2
f
i
−
1
,
j
x
o
r
a
i
f_{i,j}=f_{i - 1,j} +2f_{i - 1,j \,{\rm xor} \,a_i}
fi,j=fi−1,j+2fi−1,jxorai
边界: f 0 , 0 = 1 f_{0,0}=1 f0,0=1
考虑优化这个DP,发现这个式子本质上就是
f
i
−
1
f_{i - 1}
fi−1 与一个第
0
0
0 位是
1
1
1,第
a
i
a_i
ai 位是
2
2
2 且只有这两个位有值的多项式做
x
o
r
{\rm xor}
xor 卷积。我们不妨设第
i
i
i 个多项式为
A
i
A_i
Ai,那么显然有
F
W
T
(
f
n
)
=
F
W
T
(
A
1
)
⋅
F
W
T
(
A
2
)
⋯
F
W
T
(
A
n
)
{\rm FWT}(f_n)={\rm FWT}(A_1)\cdot {\rm FWT}(A_2)\cdots {\rm FWT}(A_n)
FWT(fn)=FWT(A1)⋅FWT(A2)⋯FWT(An)
那么我们考虑能不能直接求出每一位的所有
F
W
T
(
A
i
)
{\rm FWT}(A_i)
FWT(Ai) 的积:
一方面,因为
A
i
A_i
Ai 只有
0
0
0 与
a
i
a_i
ai 处有值,根据
x
o
r
\rm xor
xor 卷积的
F
W
T
\rm FWT
FWT 定义式,我们可以得到
F
W
T
(
A
i
)
{\rm FWT}(A_i)
FWT(Ai) 的每个位置只能是
−
1
-1
−1 或
3
3
3;另一方面,根据
F
W
T
{\rm FWT}
FWT 的和等于和的
F
W
T
{\rm FWT}
FWT 这一性质,我们只要将
A
i
A_i
Ai 全部相加并做一遍
F
W
T
{\rm FWT}
FWT,就可以知道每一位的
F
W
T
{\rm FWT}
FWT 的总和,又因为值只能是
−
1
-1
−1 或
3
3
3,我们就能直接解出来
−
1
-1
−1 和
3
3
3 分别有几个了!
那么之后就很好办了,我们只要先对 f 0 f_0 f0 做一遍 F W T {\rm FWT} FWT,再按解出来的值对每一位乘上贡献,最后再逆变换回去就能求出 f n f_n fn 了。
最后答案就是 f n , 0 − 1 f_{n,0} - 1 fn,0−1,因为最后还要再减去两个都是空集的情况
代码
#include <iostream>
#include <cstdio>
using namespace std;
const int maxn = 1e6 + 50,maxm = 1.1e6,mod = 998244353,inv2 = 499122177;
int n,m = 1,x,mx,pw[maxn],a[maxm],f[maxm];
int read(){
int x = 0;
char c = getchar();
while(c < '0' || c > '9') c = getchar();
while(c >= '0' && c <= '9') x = x * 10 + (c ^ 48),c = getchar();
return x;
}
inline int add(int x,int y){
if(x + y < mod) return x + y;
else return x + y - mod;
}
inline int dec(int x,int y){
if(x - y >= 0) return x - y;
else return x - y + mod;
}
void FWT(int *f,int n,int opt){
for(int i = 1; i < n; i <<= 1){
for(int j = 0; j < n; j += i << 1){
for(int k = j; k < j + i; k ++){
int t1 = f[k],t2 = f[k + i];
f[k] = add(t1,t2),f[k + i] = dec(t1,t2);
if(opt == -1) f[k] = 1ll * f[k] * inv2 % mod,f[k + i] = 1ll * f[k + i] * inv2 % mod;
}
}
}
}
void init(int n){
pw[0] = 1;
for(int i = 1; i <= n; i ++) pw[i] = 3ll * pw[i - 1] % mod;
for(int i = 1; i <= n; i ++) x = read(),mx = max(mx,x),a[0] ++,a[x] += 2;
while(m <= mx) m <<= 1;
}
int main(){
n = read();
init(n);
FWT(a,m,1);
for(int i = 0; i < m; i ++){
int x = 1ll * dec(3 * n,a[i]) * inv2 % mod * inv2 % mod;
f[i] = x & 1 ? mod - pw[n - x] : pw[n - x];
//由于边界是f[0][0]=1,根据FWT定义式,我们可以得出实际上对f[0]做FWT后每一位都是1,我们就没有必要再做一遍了
}
FWT(f,m,-1);
printf("%d\n",dec(f[0],1));
return 0;
}
[HAOI2015]按位或
上一篇Min-Max容斥
里的题,我就不要脸的搬过来了
为什么上一篇会有
F
W
T
\sout {\rm FWT}
FWT
分析
老套路,我们令 S S S 表示每一位第一次变成一的时间集合,那么我们要求的就是
E ( max ( S ) ) = ∑ ∅ ≠ T ⊆ S ( − 1 ) ∣ T ∣ − 1 E ( min ( T ) ) {\rm E}(\max(S))=\sum_{\varnothing\neq T \subseteq S}(-1)^{|T|-1}{\rm E}(\min(T)) E(max(S))=∅=T⊆S∑(−1)∣T∣−1E(min(T))这题 n n n很小,我们可以直接枚举 T T T,于是问题就变成了如何求 E ( min ( T ) ) {\rm E}(\min(T)) E(min(T))
我们设 f ( S ) f(S) f(S) 表示单次只取到 S S S 中的位置的概率,那么有
E ( min ( T ) ) = ∑ k = 1 + ∞ k f ( S ∖ T ) k − 1 ( 1 − f ( S ∖ T ) ) {\rm E}(\min(T))=\sum_{k = 1}^{+\infty}kf(S \setminus T)^{k - 1}(1-f(S\setminus T)) E(min(T))=k=1∑+∞kf(S∖T)k−1(1−f(S∖T))经过一通爆算,得到(好像也不是很难算)
E ( min ( T ) ) = 1 1 − f ( S ∖ T ) {\rm E}(\min(T))=\frac1{1-f(S\setminus T)} E(min(T))=1−f(S∖T)1(上面两行貌似叫离散随机变量的几何分布)
于是现在问题就转化为了如何求 f ( S ) f(S) f(S),根据定义,显然有
f ( S ) = ∑ T ⊆ S p T f(S)=\sum_{T\subseteq S}p_T f(S)=T⊆S∑pT这不就是 F W T \rm FWT FWT 吗?我们直接 O ( n 2 n ) \Omicron(n2^n) O(n2n) 做就行了
注意不可能时输出
INF
.
代码
#include <iostream>
#include <cstdio>
using namespace std;
const int maxn = 20,maxm = 1 << maxn;
const double eps = 1e-8;
int n,m,cnt[maxm];
double ans,p[maxm];
int main(){
scanf("%d",&n);
m = (1 << n) - 1;
for(int i = 0; i <= m; i ++) scanf("%lf",&p[i]);
for(int i = 1; i <= m; i <<= 1){
for(int j = 0; j <= m; j += i << 1){
for(int k = j; k < j + i; k ++){
p[k + i] += p[k];
}
}
}
for(int i = 1; i <= m; i ++) cnt[i] = cnt[i >> 1] + (i & 1);
for(int i = 1; i <= m; i ++){
double t = 1 - p[i ^ m];
if(t < eps){
printf("INF\n");
return 0;
}
if(cnt[i] & 1) ans += 1 / t;
else ans -= 1 / t;
}
printf("%f\n",ans);
return 0;
}
Graph And Numbers
分析
我们考虑一个简单容斥:
答案
=
=
= 总方案数
−
-
− 没有
0
0
0边
−
-
− 没有
1
1
1边
−
-
− 没有
2
2
2边
+
+
+ 只有
0
0
0边
+
+
+ 只有
1
1
1边
+
+
+ 只有
2
2
2边
−
-
− 没有边
于是我们把可以这些项分开考虑:
- 没有 1 1 1边:也就是每个联通块里只有 1 1 1种数字,记 c n t {cnt} cnt 为联通块个数,则数量为 2 c n t 2^{cnt} 2cnt
- 只有 0 / 2 0/2 0/2边:容易发现这两项的数量是一样的,在大小大于 1 1 1的连通块中只能全选 0 / 1 0/1 0/1,而大小为 1 1 1的连通块中可以任填,故记 c n t 1 cnt1 cnt1 为大小为 1 1 1的连通块个数,则数量为 2 c n t 1 2^{cnt1} 2cnt1
- 只有 1 1 1边:因为只有 1 1 1边,所以相邻点的数字一定是不同的,所以我们将每个连通块黑白染色,如果存在无法染色的情况,则答案为 0 0 0;否则答案为 2 c n t 2^{cnt} 2cnt
- 没有边:显然,只有在 m = 0 m =0 m=0 时数量为 1 1 1,其他时候不存在,但要记得特判
你会发现上面还少了两种情况:没有 0 0 0边 / / / 没有 2 2 2边。这两种情况比较复杂,但显然它们的数量也是相等的,我们只需要计算其中一种即可,不妨算没有 0 0 0边的。
因为 n ⩽ 40 n \leqslant 40 n⩽40,我们可以考虑折半搜索,即搜索出前一半点的状态并在过程中得到当前状态下另一半点应满足的限制条件。显然这些限制条件一定是形如某些位不能为 0 0 0,故可以将这些限制压成二进制。那么不难发现对每一个最终搜出来的限制,它的所有超集所对应的方案都是满足的(只有 0 0 0变成 1 1 1,显然不劣),即它会对它的所有超集做贡献。
那么我们就有了一个大概的想法:先将前一半点的状态都搜索出来,并在每个最终限制位置上 + 1 +1 +1,然后做一遍子集卷积,搜索后一半的状态并统计所有合法状态的答案。
但这样其实是有问题的,因为我们第一次搜索的时候并没有保证后一半互相之间是合法的,且我们做的是子集卷积,所以一些不合法的后一半状态会贡献到它的合法的超集,这些贡献就会被我们统计。
于是我们可以转化一下:我们搜索到最终限制时将下标取反,那么子集卷积就可以改为超集卷积,于是这样我们最终统计的答案就都是合法的了(因为是超集卷积,所有不合法的后一半状态只会贡献到它的子集,而它的子集显然也是不合法的,但我们统计的都是合法的后一半状态,所以这些不合法的不会对我们的答案做贡献)。
代码
#include <iostream>
#include <cstdio>
using namespace std;
const int maxn = 41,maxm = 1 << 20;
int n,m,u,v,cnt1,cnt2,flag = 1,T,lim,vis[maxn];
long long ans,e[maxn],f[maxm];
int read(){
int x = 0;
char c = getchar();
while(c < '0' || c > '9') c = getchar();
while(c >= '0' && c <= '9') x = x * 10 + (c ^ 48),c = getchar();
return x;
}
int dfs(int u,int p){
vis[u] = p;
int sum = 1;
for(int i = 0; i <= n; i ++){
if(e[u] >> i & 1){
if(!vis[i]) sum += dfs(i,-p);
else if(vis[i] == vis[u]) flag = 0;
}
}
return sum;
}
void dfs1(int k,long long S){
if(k == T - 1){
f[(~S) & (lim - 1)] ++;
return;
}
dfs1(k - 1,S);
if(!(S >> k & 1)) dfs1(k - 1,S | e[k]);
}
void dfs2(int k,long long S,long long used){
if(k == T){
ans -= 2 * f[used];
return;
}
dfs2(k + 1,S,used);
if(!(S >> k & 1)) dfs2(k + 1,S | e[k],used | (1 << k));
}
void FWT(long long *f,int n){
for(int i = 1; i < n; i <<= 1){
for(int j = 0; j < n; j += i << 1){
for(int k = j; k < j + i; k ++){
f[k] += f[k + i];
}
}
}
}
int main(){
n = read() - 1,m = read();
if(m == 0){// 记得特判 !!!
cout << 0 << endl;
return 0;
}
for(int i = 1; i <= m; i ++){
u = read() - 1,v = read() - 1;
e[u] |= 1ll << v,e[v] |= 1ll << u;
}
for(int i = 0; i <= n; i ++){
if(!vis[i]){
if(dfs(i,1) == 1) cnt2 ++;
cnt1 ++;
}
}
ans = (1ll << (n + 1)) - (1ll << cnt1) + (2ll << cnt2) + flag * (1ll << cnt1);
T = (n + 1) / 2,lim = 1 << T;
dfs1(n,0);
FWT(f,lim);
dfs2(0,0,0);
printf("%lld\n",ans);
return 0;
}