haar_lib/math/convolution/
subset_conv.rs1use crate::math::convolution::{mobius::*, zeta::*};
6use std::ops::{Add, Mul, Sub};
7
8#[allow(clippy::needless_range_loop, clippy::manual_memcpy)]
13pub fn subset_convolution<T>(f: Vec<T>, g: Vec<T>) -> Vec<T>
14where
15 T: Copy + Default + Add<Output = T> + Sub<Output = T> + Mul<Output = T>,
16{
17 assert_eq!(f.len(), g.len());
18
19 let n = f.len();
20 assert!(n.is_power_of_two());
21
22 let k = n.trailing_zeros() as usize;
23
24 let mut f2 = vec![vec![T::default(); n]; k + 1];
25 let mut g2 = vec![vec![T::default(); n]; k + 1];
26
27 for j in 0..n {
28 let i = j.count_ones() as usize;
29 f2[i][j] = f[j];
30 g2[i][j] = g[j];
31 }
32
33 for i in 0..=k {
34 fast_zeta_subset(&mut f2[i]);
35 fast_zeta_subset(&mut g2[i]);
36 }
37
38 let mut h = vec![vec![T::default(); n]; k + 1];
39
40 for i in 0..=k {
41 for s in 0..=i {
42 for j in 0..n {
43 h[i][j] = h[i][j] + f2[s][j] * g2[i - s][j];
44 }
45 }
46 }
47
48 for i in 0..=k {
49 fast_mobius_subset(&mut h[i]);
50 }
51
52 (0..n)
53 .map(|j| {
54 let i = j.count_ones() as usize;
55 h[i][j]
56 })
57 .collect()
58}