|
| 1 | +#ifndef CP_ALGO_MATH_XOR_CONVOLUTION_HPP |
| 2 | +#define CP_ALGO_MATH_XOR_CONVOLUTION_HPP |
| 3 | +#include "../number_theory/modint.hpp" |
| 4 | +#include "../util/bit.hpp" |
| 5 | +#include "../util/checkpoint.hpp" |
| 6 | +#include <cassert> |
| 7 | +#include <algorithm> |
| 8 | +#include <vector> |
| 9 | + |
| 10 | +namespace cp_algo::math { |
| 11 | + // Recursive FWHT (XOR) transform for size N (power of two) |
| 12 | + template<auto N> |
| 13 | + void xor_transform(auto &&a) { |
| 14 | + if constexpr (N == 1) { |
| 15 | + return; |
| 16 | + } else { |
| 17 | + constexpr auto half = N / 2; |
| 18 | + xor_transform<half>(&a[0]); |
| 19 | + xor_transform<half>(&a[half]); |
| 20 | + for (uint32_t i = 0; i < half; i++) { |
| 21 | + auto x = a[i] + a[i + half]; |
| 22 | + auto y = a[i] - a[i + half]; |
| 23 | + a[i] = x; |
| 24 | + a[i + half] = y; |
| 25 | + } |
| 26 | + } |
| 27 | + } |
| 28 | + |
| 29 | + // FWHT wrapper that deduces N at compile time via with_bit_floor |
| 30 | + inline void xor_transform(auto &&a, auto n) { |
| 31 | + with_bit_floor(n, [&]<auto NN>() { |
| 32 | + assert(NN == n); |
| 33 | + xor_transform<NN>(a); |
| 34 | + }); |
| 35 | + } |
| 36 | + |
| 37 | + inline void xor_transform(auto &&a) { |
| 38 | + xor_transform(a, std::size(a)); |
| 39 | + } |
| 40 | + |
| 41 | + // In-place XOR convolution on sequences of equal length (power of two) |
| 42 | + void xor_convolution_inplace(auto &a, auto &b) { |
| 43 | + auto N = static_cast<uint32_t>(std::size(a)); |
| 44 | + xor_transform(a); |
| 45 | + xor_transform(b); |
| 46 | + checkpoint("transform"); |
| 47 | + for (uint32_t i = 0; i < N; i++) { |
| 48 | + a[i] *= b[i]; |
| 49 | + } |
| 50 | + checkpoint("dot"); |
| 51 | + xor_transform(a); |
| 52 | + checkpoint("transform"); |
| 53 | + using base = std::decay_t<decltype(a[0])>; |
| 54 | + base ni = base(N).inv(); |
| 55 | + for (auto &it : a) { |
| 56 | + it *= ni; |
| 57 | + } |
| 58 | + checkpoint("mul_inv"); |
| 59 | + } |
| 60 | + |
| 61 | + // Returns XOR convolution of a and b; pads to next power of two |
| 62 | + auto xor_convolution(auto a, auto b) { |
| 63 | + auto n = std::bit_ceil(std::max(std::size(a), std::size(b))); |
| 64 | + a.resize(n); |
| 65 | + b.resize(n); |
| 66 | + xor_convolution_inplace(a, b); |
| 67 | + return a; |
| 68 | + } |
| 69 | +} |
| 70 | +#endif // CP_ALGO_MATH_XOR_CONVOLUTION_HPP |
0 commit comments