m1une's library

This documentation is automatically generated by online-judge-tools/verification-helper

View on GitHub

:heavy_check_mark: Bitwise Convolution
(math/bitwise_convolution.hpp)

Overview

Bitwise convolutions indexed by masks:

The header also exposes the underlying subset, superset, and Walsh-Hadamard transforms. The subset and superset transforms are defined in zeta_mobius_transform.hpp, which this header includes.

#include "math/bitwise_convolution.hpp"

All names are in m1une::math.

Convolutions

For input arrays a and b, the functions compute:

\[c_k = \sum_{i \mathbin{\mathrm{op}} j = k} a_i b_j.\]
template <typename T>
std::vector<T> bitwise_or_convolution(
    std::vector<T> first,
    std::vector<T> second
);

template <typename T>
std::vector<T> bitwise_and_convolution(
    std::vector<T> first,
    std::vector<T> second
);

template <typename T>
std::vector<T> bitwise_xor_convolution(
    std::vector<T> first,
    std::vector<T> second
);

The input vectors are passed by value, so they may be moved into the functions. The originals passed by the caller are otherwise unchanged. The return type is std::vector<T>.

Function Operation Complexity
bitwise_or_convolution(a, b) op is bitwise OR. $O(N \log N)$
bitwise_and_convolution(a, b) op is bitwise AND. $O(N \log N)$
bitwise_xor_convolution(a, b) op is bitwise XOR. $O(N \log N)$

If either input is empty, the result is empty. Otherwise both arrays are zero-padded to the smallest power of two at least max(a.size(), b.size()). The returned vector has that length. Each function uses $O(N)$ additional memory.

The element type must support default construction, addition, subtraction, and multiplication. XOR convolution additionally requires division by the transform size and construction from that size. Integers work because the inverse Walsh-Hadamard transform is exactly divisible; modular types require the power-of-two transform size to be invertible.

Transforms

All transforms operate in place and require a nonempty power-of-two length.

template <typename T>
void walsh_hadamard_transform(
    std::vector<T>& values,
    bool inverse = false
);

The subset and superset transform signatures and their requirements on T are documented in zeta_mobius_transform.hpp. The Walsh-Hadamard transform modifies values directly and returns void.

Function Description Complexity
subset_zeta_transform(values) Replaces each mask with the sum over its submasks. $O(N \log N)$
subset_mobius_transform(values) Inverts the subset zeta transform. $O(N \log N)$
superset_zeta_transform(values) Replaces each mask with the sum over its supermasks. $O(N \log N)$
superset_mobius_transform(values) Inverts the superset zeta transform. $O(N \log N)$
walsh_hadamard_transform(values, inverse) Applies the XOR transform or its inverse. $O(N \log N)$

For a mask $S$, subset zeta computes $F(S) = \sum_{T \subseteq S} f(T)$, while superset zeta computes $F(S) = \sum_{T \supseteq S} f(T)$. Their Möbius transforms recover $f$.

The Walsh-Hadamard transform computes $F(S) = \sum_T (-1)^{\operatorname{popcount}(S \mathbin{\&} T)} f(T)$. Passing true as the second argument applies the same butterflies and divides every element by $N$, producing the inverse transform.

Example

#include "math/bitwise_convolution.hpp"
#include "math/modint.hpp"

#include <vector>

int main() {
    using mint = m1une::math::modint998244353;

    std::vector<mint> a{1, 2, 3, 4};
    std::vector<mint> b{5, 6, 7, 8};
    std::vector<mint> result =
        m1une::math::bitwise_xor_convolution(a, b);
}

Depends on

Required by

Verified with

Code

#ifndef M1UNE_MATH_BITWISE_CONVOLUTION_HPP
#define M1UNE_MATH_BITWISE_CONVOLUTION_HPP 1

#include <cassert>
#include <cstddef>
#include <utility>
#include <vector>

#include "zeta_mobius_transform.hpp"

namespace m1une {
namespace math {

namespace bitwise_convolution_detail {

inline std::size_t common_size(
    std::size_t first_size,
    std::size_t second_size
) {
    std::size_t required = first_size > second_size
        ? first_size
        : second_size;
    std::size_t size = 1;
    while (size < required) size <<= 1;
    return size;
}

template <typename T>
std::vector<T> pointwise_product(
    std::vector<T> first,
    const std::vector<T>& second
) {
    assert(first.size() == second.size());
    for (std::size_t index = 0; index < first.size(); ++index) {
        first[index] *= second[index];
    }
    return first;
}

}  // namespace bitwise_convolution_detail

template <typename T>
void walsh_hadamard_transform(
    std::vector<T>& values,
    bool inverse = false
) {
    assert(zeta_mobius_transform_detail::is_power_of_two(values.size()));
    for (std::size_t length = 1; length < values.size(); length <<= 1) {
        for (
            std::size_t block = 0;
            block < values.size();
            block += length << 1
        ) {
            for (std::size_t offset = 0; offset < length; ++offset) {
                T first = values[block + offset];
                T second = values[block + offset + length];
                values[block + offset] = first + second;
                values[block + offset + length] = first - second;
            }
        }
    }
    if (inverse) {
        T size = T(static_cast<long long>(values.size()));
        for (T& value : values) value /= size;
    }
}

template <typename T>
std::vector<T> bitwise_or_convolution(
    std::vector<T> first,
    std::vector<T> second
) {
    if (first.empty() || second.empty()) return {};
    std::size_t size = bitwise_convolution_detail::common_size(
        first.size(),
        second.size()
    );
    first.resize(size);
    second.resize(size);
    subset_zeta_transform(first);
    subset_zeta_transform(second);
    first = bitwise_convolution_detail::pointwise_product(
        std::move(first),
        second
    );
    subset_mobius_transform(first);
    return first;
}

template <typename T>
std::vector<T> bitwise_and_convolution(
    std::vector<T> first,
    std::vector<T> second
) {
    if (first.empty() || second.empty()) return {};
    std::size_t size = bitwise_convolution_detail::common_size(
        first.size(),
        second.size()
    );
    first.resize(size);
    second.resize(size);
    superset_zeta_transform(first);
    superset_zeta_transform(second);
    first = bitwise_convolution_detail::pointwise_product(
        std::move(first),
        second
    );
    superset_mobius_transform(first);
    return first;
}

template <typename T>
std::vector<T> bitwise_xor_convolution(
    std::vector<T> first,
    std::vector<T> second
) {
    if (first.empty() || second.empty()) return {};
    std::size_t size = bitwise_convolution_detail::common_size(
        first.size(),
        second.size()
    );
    first.resize(size);
    second.resize(size);
    walsh_hadamard_transform(first);
    walsh_hadamard_transform(second);
    first = bitwise_convolution_detail::pointwise_product(
        std::move(first),
        second
    );
    walsh_hadamard_transform(first, true);
    return first;
}

}  // namespace math
}  // namespace m1une

#endif  // M1UNE_MATH_BITWISE_CONVOLUTION_HPP
#line 1 "math/bitwise_convolution.hpp"



#include <cassert>
#include <cstddef>
#include <utility>
#include <vector>

#line 1 "math/zeta_mobius_transform.hpp"



#line 7 "math/zeta_mobius_transform.hpp"

namespace m1une {
namespace math {

namespace zeta_mobius_transform_detail {

inline bool is_power_of_two(std::size_t size) noexcept {
    return size != 0 && (size & (size - 1)) == 0;
}

inline std::vector<std::size_t> primes_up_to(std::size_t limit) {
    std::vector<std::size_t> primes;
    std::vector<bool> is_prime(limit + 1, true);
    if (!is_prime.empty()) is_prime[0] = false;
    if (limit >= 1) is_prime[1] = false;
    for (std::size_t value = 2; value <= limit; ++value) {
        if (!is_prime[value]) continue;
        primes.emplace_back(value);
        if (value > limit / value) continue;
        for (
            std::size_t multiple = value * value;
            multiple <= limit;
            multiple += value
        ) {
            is_prime[multiple] = false;
        }
    }
    return primes;
}

}  // namespace zeta_mobius_transform_detail

template <typename T>
void subset_zeta_transform(std::vector<T>& values) {
    assert(zeta_mobius_transform_detail::is_power_of_two(values.size()));
    for (std::size_t bit = 1; bit < values.size(); bit <<= 1) {
        for (
            std::size_t block = 0;
            block < values.size();
            block += bit << 1
        ) {
            for (std::size_t offset = 0; offset < bit; ++offset) {
                values[block + bit + offset] += values[block + offset];
            }
        }
    }
}

template <typename T>
void subset_mobius_transform(std::vector<T>& values) {
    assert(zeta_mobius_transform_detail::is_power_of_two(values.size()));
    for (std::size_t bit = 1; bit < values.size(); bit <<= 1) {
        for (
            std::size_t block = 0;
            block < values.size();
            block += bit << 1
        ) {
            for (std::size_t offset = 0; offset < bit; ++offset) {
                values[block + bit + offset] -= values[block + offset];
            }
        }
    }
}

template <typename T>
void superset_zeta_transform(std::vector<T>& values) {
    assert(zeta_mobius_transform_detail::is_power_of_two(values.size()));
    for (std::size_t bit = 1; bit < values.size(); bit <<= 1) {
        for (
            std::size_t block = 0;
            block < values.size();
            block += bit << 1
        ) {
            for (std::size_t offset = 0; offset < bit; ++offset) {
                values[block + offset] += values[block + bit + offset];
            }
        }
    }
}

template <typename T>
void superset_mobius_transform(std::vector<T>& values) {
    assert(zeta_mobius_transform_detail::is_power_of_two(values.size()));
    for (std::size_t bit = 1; bit < values.size(); bit <<= 1) {
        for (
            std::size_t block = 0;
            block < values.size();
            block += bit << 1
        ) {
            for (std::size_t offset = 0; offset < bit; ++offset) {
                values[block + offset] -= values[block + bit + offset];
            }
        }
    }
}

template <typename T>
void divisor_zeta_transform(std::vector<T>& values) {
    if (values.size() <= 2) return;
    const std::size_t limit = values.size() - 1;
    const std::vector<std::size_t> primes =
        zeta_mobius_transform_detail::primes_up_to(limit);
    for (std::size_t prime : primes) {
        for (std::size_t value = 1; value <= limit / prime; ++value) {
            values[value * prime] += values[value];
        }
    }
}

template <typename T>
void divisor_mobius_transform(std::vector<T>& values) {
    if (values.size() <= 2) return;
    const std::size_t limit = values.size() - 1;
    const std::vector<std::size_t> primes =
        zeta_mobius_transform_detail::primes_up_to(limit);
    for (std::size_t prime : primes) {
        for (
            std::size_t value = limit / prime;
            value >= 1;
            --value
        ) {
            values[value * prime] -= values[value];
        }
    }
}

template <typename T>
void multiple_zeta_transform(std::vector<T>& values) {
    if (values.size() <= 2) return;
    const std::size_t limit = values.size() - 1;
    const std::vector<std::size_t> primes =
        zeta_mobius_transform_detail::primes_up_to(limit);
    for (std::size_t prime : primes) {
        for (
            std::size_t value = limit / prime;
            value >= 1;
            --value
        ) {
            values[value] += values[value * prime];
        }
    }
}

template <typename T>
void multiple_mobius_transform(std::vector<T>& values) {
    if (values.size() <= 2) return;
    const std::size_t limit = values.size() - 1;
    const std::vector<std::size_t> primes =
        zeta_mobius_transform_detail::primes_up_to(limit);
    for (std::size_t prime : primes) {
        for (std::size_t value = 1; value <= limit / prime; ++value) {
            values[value] -= values[value * prime];
        }
    }
}

}  // namespace math
}  // namespace m1une


#line 10 "math/bitwise_convolution.hpp"

namespace m1une {
namespace math {

namespace bitwise_convolution_detail {

inline std::size_t common_size(
    std::size_t first_size,
    std::size_t second_size
) {
    std::size_t required = first_size > second_size
        ? first_size
        : second_size;
    std::size_t size = 1;
    while (size < required) size <<= 1;
    return size;
}

template <typename T>
std::vector<T> pointwise_product(
    std::vector<T> first,
    const std::vector<T>& second
) {
    assert(first.size() == second.size());
    for (std::size_t index = 0; index < first.size(); ++index) {
        first[index] *= second[index];
    }
    return first;
}

}  // namespace bitwise_convolution_detail

template <typename T>
void walsh_hadamard_transform(
    std::vector<T>& values,
    bool inverse = false
) {
    assert(zeta_mobius_transform_detail::is_power_of_two(values.size()));
    for (std::size_t length = 1; length < values.size(); length <<= 1) {
        for (
            std::size_t block = 0;
            block < values.size();
            block += length << 1
        ) {
            for (std::size_t offset = 0; offset < length; ++offset) {
                T first = values[block + offset];
                T second = values[block + offset + length];
                values[block + offset] = first + second;
                values[block + offset + length] = first - second;
            }
        }
    }
    if (inverse) {
        T size = T(static_cast<long long>(values.size()));
        for (T& value : values) value /= size;
    }
}

template <typename T>
std::vector<T> bitwise_or_convolution(
    std::vector<T> first,
    std::vector<T> second
) {
    if (first.empty() || second.empty()) return {};
    std::size_t size = bitwise_convolution_detail::common_size(
        first.size(),
        second.size()
    );
    first.resize(size);
    second.resize(size);
    subset_zeta_transform(first);
    subset_zeta_transform(second);
    first = bitwise_convolution_detail::pointwise_product(
        std::move(first),
        second
    );
    subset_mobius_transform(first);
    return first;
}

template <typename T>
std::vector<T> bitwise_and_convolution(
    std::vector<T> first,
    std::vector<T> second
) {
    if (first.empty() || second.empty()) return {};
    std::size_t size = bitwise_convolution_detail::common_size(
        first.size(),
        second.size()
    );
    first.resize(size);
    second.resize(size);
    superset_zeta_transform(first);
    superset_zeta_transform(second);
    first = bitwise_convolution_detail::pointwise_product(
        std::move(first),
        second
    );
    superset_mobius_transform(first);
    return first;
}

template <typename T>
std::vector<T> bitwise_xor_convolution(
    std::vector<T> first,
    std::vector<T> second
) {
    if (first.empty() || second.empty()) return {};
    std::size_t size = bitwise_convolution_detail::common_size(
        first.size(),
        second.size()
    );
    first.resize(size);
    second.resize(size);
    walsh_hadamard_transform(first);
    walsh_hadamard_transform(second);
    first = bitwise_convolution_detail::pointwise_product(
        std::move(first),
        second
    );
    walsh_hadamard_transform(first, true);
    return first;
}

}  // namespace math
}  // namespace m1une
Back to top page