Bitwise XOR Convolution
(Src/Sequence/BitwiseXORConvolution.hpp)

View this file on GitHub
Last update: 2025-06-02 16:41:21+09:00
Include: #include "Src/Sequence/BitwiseXORConvolution.hpp"

概要

長さ $N = 2^{K}$ の数列 $a, b$ に対して $\displaystyle c_{i} = \sum_{p, q, p\oplus q = i} a_{p}\times b_{q}$ を $O(N\log N)$ で計算する。

内部でWalsh-Hadamard Transformというものを用いている。愚直なTransformは低速で、これを用いても畳み込みの計算量は $\Theta (N^{2})$ から改善しない。しかし、Transformの式をぐっと睨むことにより高速化に成功する(俗に高速アダマール変換と呼ばれる)。

計算結果の値が小さいと見積もることができても、計算途中でオーバーフローする場合があるので、返り値の型には注意すること(例: ABC396-G)

参考

Depends on

標準データ型のエイリアス (Src/Template/TypeAlias.hpp)

Verified with

Code

#pragma once

#include "../Template/TypeAlias.hpp"

#include <bit>
#include <concepts>
#include <cassert>
#include <vector>

namespace zawa {

// note: 返り値の各点の値は真の値より(2^{k/2})倍されている
template <class T>
void FastWalshHadamardTransform(std::vector<T>& A) {
    if (A.empty()) return;
    while (!std::has_single_bit(A.size())) A.push_back(T{0});
    const usize k = std::bit_width(A.size()) - 1, n = A.size();
    for (usize i = 0 ; i < k ; i++) {
        const usize bit = 1 << i;
        for (usize j = 0 ; j < n ; j += bit << 1) {
            for (usize k = 0 ; k < bit ; k++) {
                const T p = A[j + k], q = A[j + k + bit];
                A[j + k] = p + q;
                A[j + k + bit] = p - q;
            }
        }
    }
}

template <class T>
std::vector<T> BitwiseXORConvolution(std::vector<T> A, std::vector<T> B) {
    FastWalshHadamardTransform(A);
    FastWalshHadamardTransform(B);
    if (A.size() > B.size()) std::swap(A, B);
    for (usize i = 0 ; i < A.size() ; i++) A[i] *= B[i];
    FastWalshHadamardTransform(A);
    if (A.empty()) return A;
    assert(std::has_single_bit(A.size()));
    if constexpr (std::integral<T>) {
        const usize d = std::bit_width(A.size()) - 1;
        for (auto& a : A) a >>= d;
    }
    else {
        const T d = [&]() {
            assert(std::has_single_bit(A.size()));
            usize exp = std::bit_width(A.size()) - 1;
            T v = T{1} / T{2}, res = T{1};
            while (exp) {
                if (exp & 1) res = res * v;
                v = v * v;
                exp >>= 1;
            }
            return res;
        }();
        for (T& a : A) a *= d;
    }
    return A;
}

template <class T, class U>
requires (!std::same_as<T, U> and std::convertible_to<U, T>)
std::vector<T> BitwiseXORConvolution(std::vector<U> A, std::vector<U> B) {
    std::vector<T> a(A.size()), b(B.size());
    for (usize i = 0 ; i < A.size() ; i++) a[i] = static_cast<T>(std::move(A[i]));
    for (usize i = 0 ; i < B.size() ; i++) b[i] = static_cast<T>(std::move(B[i]));
    return BitwiseXORConvolution<T>(a, b);
}

} // namespace zawa

#line 2 "Src/Sequence/BitwiseXORConvolution.hpp"

#line 2 "Src/Template/TypeAlias.hpp"

#include <cstdint>
#include <cstddef>

namespace zawa {

using i16 = std::int16_t;
using i32 = std::int32_t;
using i64 = std::int64_t;
using i128 = __int128_t;

using u8 = std::uint8_t;
using u16 = std::uint16_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;

using usize = std::size_t;

} // namespace zawa
#line 4 "Src/Sequence/BitwiseXORConvolution.hpp"

#include <bit>
#include <concepts>
#include <cassert>
#include <vector>

namespace zawa {

// note: 返り値の各点の値は真の値より(2^{k/2})倍されている
template <class T>
void FastWalshHadamardTransform(std::vector<T>& A) {
    if (A.empty()) return;
    while (!std::has_single_bit(A.size())) A.push_back(T{0});
    const usize k = std::bit_width(A.size()) - 1, n = A.size();
    for (usize i = 0 ; i < k ; i++) {
        const usize bit = 1 << i;
        for (usize j = 0 ; j < n ; j += bit << 1) {
            for (usize k = 0 ; k < bit ; k++) {
                const T p = A[j + k], q = A[j + k + bit];
                A[j + k] = p + q;
                A[j + k + bit] = p - q;
            }
        }
    }
}

template <class T>
std::vector<T> BitwiseXORConvolution(std::vector<T> A, std::vector<T> B) {
    FastWalshHadamardTransform(A);
    FastWalshHadamardTransform(B);
    if (A.size() > B.size()) std::swap(A, B);
    for (usize i = 0 ; i < A.size() ; i++) A[i] *= B[i];
    FastWalshHadamardTransform(A);
    if (A.empty()) return A;
    assert(std::has_single_bit(A.size()));
    if constexpr (std::integral<T>) {
        const usize d = std::bit_width(A.size()) - 1;
        for (auto& a : A) a >>= d;
    }
    else {
        const T d = [&]() {
            assert(std::has_single_bit(A.size()));
            usize exp = std::bit_width(A.size()) - 1;
            T v = T{1} / T{2}, res = T{1};
            while (exp) {
                if (exp & 1) res = res * v;
                v = v * v;
                exp >>= 1;
            }
            return res;
        }();
        for (T& a : A) a *= d;
    }
    return A;
}

template <class T, class U>
requires (!std::same_as<T, U> and std::convertible_to<U, T>)
std::vector<T> BitwiseXORConvolution(std::vector<U> A, std::vector<U> B) {
    std::vector<T> a(A.size()), b(B.size());
    for (usize i = 0 ; i < A.size() ; i++) a[i] = static_cast<T>(std::move(A[i]));
    for (usize i = 0 ; i < B.size() ; i++) b[i] = static_cast<T>(std::move(B[i]));
    return BitwiseXORConvolution<T>(a, b);
}

} // namespace zawa

Bitwise XOR Convolution (Src/Sequence/BitwiseXORConvolution.hpp)

概要

参考

Depends on

Verified with

Code

Bitwise XOR Convolution
(Src/Sequence/BitwiseXORConvolution.hpp)