competitive/math/
subset_convolve.rs1use super::{BitwiseorConvolve, ConvolveSteps, Invertible, Ring};
2use std::marker::PhantomData;
3
4pub struct SubsetConvolve<M> {
5 _marker: PhantomData<fn() -> M>,
6}
7
8impl<R> ConvolveSteps for SubsetConvolve<R>
9where
10 R: Ring<T: PartialEq, Additive: Invertible>,
11{
12 type T = Vec<R::T>;
13 type F = Vec<Vec<R::T>>;
14
15 fn length(t: &Self::T) -> usize {
16 t.len()
17 }
18
19 fn transform(t: Self::T, len: usize) -> Self::F {
20 let k = len.trailing_zeros() as usize;
21 let mut f = vec![vec![R::zero(); len]; k + 1];
22 for (i, t) in t.iter().enumerate() {
23 let f = &mut f[i.count_ones() as usize][i];
24 *f = R::add(f, t);
25 }
26 for f in f.iter_mut() {
27 BitwiseorConvolve::<R::Additive>::zeta_transform(f);
28 }
29 f
30 }
31
32 fn inverse_transform(mut f: Self::F, len: usize) -> Self::T {
33 for f in f.iter_mut() {
34 BitwiseorConvolve::<R::Additive>::mobius_transform(f);
35 }
36 let mut t = vec![R::zero(); len];
37 for (i, t) in t.iter_mut().enumerate() {
38 *t = R::add(t, &f[i.count_ones() as usize][i]);
39 }
40 t
41 }
42
43 fn multiply(f: &mut Self::F, g: &Self::F) {
44 for i in (0..f.len()).rev() {
45 let (lf, rf) = f.split_at_mut(i);
46 for (x, y) in rf[0].iter_mut().zip(&g[0]) {
47 *x = R::mul(x, y);
48 }
49 for (x, y) in lf.iter().rev().zip(g.iter().skip(1)) {
50 for ((x, y), z) in x.iter().zip(y).zip(&mut rf[0]) {
51 *z = R::add(z, &R::mul(x, y));
52 }
53 }
54 }
55 }
56
57 fn convolve(a: Self::T, b: Self::T) -> Self::T {
58 assert_eq!(a.len(), b.len());
59 let len = a.len();
60 let same = a == b;
61 let mut a = Self::transform(a, len);
62 let b = if same {
63 a.clone()
64 } else {
65 Self::transform(b, len)
66 };
67 Self::multiply(&mut a, &b);
68 Self::inverse_transform(a, len)
69 }
70}
71
72#[cfg(test)]
73mod tests {
74 use super::*;
75 use crate::{algebra::AddMulOperation, rand, tools::Xorshift};
76
77 const A: i64 = 100_000;
78
79 #[test]
80 fn test_subset_convolve() {
81 let mut rng = Xorshift::default();
82
83 for k in 0..12 {
84 let n = 1 << k;
85 rand!(rng, f: [-A..A; n], g: [-A..A; n]);
86 let mut h = vec![0i64; n];
87 for i in 0..n {
88 for j in 0..n {
89 if i & j == 0 {
90 h[i | j] += f[i] * g[j];
91 }
92 }
93 }
94 let i = SubsetConvolve::<AddMulOperation<i64>>::convolve(f, g);
95 assert_eq!(h, i);
96 }
97 }
98}